Язык R: преобразование изображений с обычными точками во фрейм данных

Question

У меня много ультрафиолетовых изображений (* .png), каждое из которых имеет 96 "лунок", которые отображаются в виде кругов. Эти колодцы расположены в формате 8 на 12. Пожалуйста, смотрите пример ниже.

На каждом изображении некоторые из 96 кругов подсвечиваются (реагируя на УФ), а некоторые нет. Я хотел бы дать номер каждому кружку, а затем идентифицировать их как «подсвеченный» или «не подсвеченный» (с заранее определенным отсечением).

Как вы думаете, какой самый простой способ достичь этого в R? Я играл с пакетом imager без особого успеха.

Еще одно замечание: не все мои изображения имеют одинаковое увеличение (т. Е. Круги не всегда имеют одинаковый размер в разных файлах, но внутри файла они примерно одинакового размера).

Answer 1

Это может быть немного поздно, но этот вопрос решает проблему, с которой я часто сталкиваюсь.Возможно, более общее решение может представлять интерес для других.Вот предлагаемое решение, которое использует еще несколько функций пакета EBImage и взаимодействует с пользователем.Его можно применять к многолуночным блюдам любой геометрии.

# EBImage needs to be available
  if (!require(EBImage)) {
    source("https://bioconductor.org/biocLite.R")
    biocLite("EBImage")
    library(EBImage)
  }

# Read the image and examine it
  fn <- file.choose() # select saved image
  img <- readImage(fn)
  img <- channel(img, "gray") # gray scale for simplicity
  plot(img)

# Define the geometry of the plate
  nx <- 12 # number of columns
  ny <- 8 # number of rows
  nwells <- seq_len(nx*ny) # index of each well in plate

# Use locator() to interactively define the upper left "corner"
# of the top left and bottom right well with a red mark to confirm
  p <- locator(2, type = "p", col = 2)

# Calculate distance between wells and create coordinates in xx, yy
  dx <- abs(diff(p$x))/(nx - 1)
  dy <- abs(diff(p$y))/(ny - 1)
  ix <- (nwells - 1) %% nx  # index for rows
  iy <- (nwells - 1) %/% nx # index for columns
  xx <- ix*dx + p$x[1] # the first point must be the upper left well
  yy <- iy*dy + p$y[1]

# Confirm the calculated coordinates with red dots
  points(xx, yy, pch = 16, col = 2)

Теперь будет создана маска объекта для использованияс семейством функций computeFeatures.Маска создается пошагово со следующим кодом.

# First a pure black image
  mask <- Image(0, dim = dim(img)) # black image same size as image

# Add white rectangles allowing 5% space between rectangles
  for (x in xx)
    for (y in yy) 
      mask[x:(x + 0.95*dx), y:(y + 0.95*dy)] <- 1

# Convert image to an 'object mask' 
  mask <- bwlabel(mask)

# One could draw circular masks with drawCircle() but this is very,
# very inefficient and doesn't improve the quality of the data
#
# mask <- Image(0, dim = dim(img)) # black image same size as image
# for (x in xx)
#   for (y in yy) 
#     mask <- drawCircle(mask, x + 0.5*dx, y + 0.5*dy,
#       radius = 0.95*dx/2, col = 1, fill = TRUE)

# Show mask on top of original image
  plot(paintObjects(mask, img, col = "white"))

# Use the computeFeatures family of functions to measure the mean intensity
  M <- computeFeatures.basic(mask, img) # a matrix of values is returned
  options(digits = 4)
  head(M)
>   b.mean   b.sd   b.mad  b.q001  b.q005  b.q05 b.q095 b.q099
> 1 0.2800 0.2114 0.18605 0.04706 0.05882 0.2118 0.7451 0.8959
> 2 0.1778 0.0981 0.09303 0.05490 0.06275 0.1569 0.3882 0.4745
> 3 0.1817 0.1028 0.10465 0.05098 0.05882 0.1569 0.4039 0.4824
> 4 0.1854 0.1029 0.10465 0.05490 0.06275 0.1608 0.3961 0.4880
> 5 0.3300 0.2425 0.23256 0.05490 0.06667 0.2510 0.7882 0.9451
> 6 0.1967 0.1076 0.11047 0.05490 0.06275 0.1765 0.4157 0.5216

# Extract the mean intensity and examine with a simple barplot
  val <- M[,"b.mean"]
  barplot(val)

# Combination of median and mad provides a fair estimate of the cutoff
  bgnd <- median(val) + mad(val)
  abline(h = bgnd, col = 2)

# Collect values in labeled data.frame and score the values
  df <- data.frame(row = factor(LETTERS[iy+1]), column = factor(ix+1),
    val = val, positive = val > bgnd)
  head(df)
>   row column    val positive
> 1   A      1 0.2800     TRUE
> 2   A      2 0.1778    FALSE
> 3   A      3 0.1817    FALSE
> 4   A      4 0.1854    FALSE
> 5   A      5 0.3300     TRUE
> 6   A      6 0.1967    FALSE

# And for visual confirmation...or to check the background value
  plot(img)
  text(xx[df$positive], yy[df$positive], "+", col = 2, cex = 2)

Answer 2

source("https://bioconductor.org/biocLite.R")
biocLite("EBImage")
library(EBImage)

fn = YOURFILE
img <- readImage(fn)

# remove outer frame
border <- 5 # 5px
dims <- dim(img)
img <- img[border:(dims[1]-border), border:(dims[2]-border),1:dims[3]]

# some despectling
img <- medianFilter(img,size=10) # blur image 
# display(img) # check if blurring is fine

# get plate mask.
highpass_rect <- .02 # if images are darker, lower this
rect <- bwlabel(dilate(img>highpass_rect, makeBrush(size=3)))
rect <- rect > 0 # remove background
# display(rect) # check if plate is recognized correctly

highpass_light <- .4 # again, darker images need lower values here for the light objects
# get single light objects
img <- bwlabel(img > highpass_light)
# display(img) # # check if all lights are displayed

# get dimensions of rectangle
rect_mid <- round(dim(rect@.Data[,,1]) / 2)
x_range <- c(min(which(rect@.Data[,rect_mid[2],1])), max(which(rect@.Data[,rect_mid[2],1])))
y_range <- c(min(which(rect@.Data[rect_mid[1],,1])), max(which(rect@.Data[rect_mid[1],,1])))

# now substract border of plate
# measured from the image you provided. 'should' scale to other images as well
# x: left 69px, right: 74px, ROI: 1085
x_range[1] <- x_range[1] + round((diff(x_range)) * (74/(1085+74+69)))
x_range[2] <- x_range[2] - round((diff(x_range)) * (69/(1085+74+69)))

# y: top 17px, bottom: 12px, ROI: 722
y_range[1] <- y_range[1] + round((diff(y_range)) * (17/(722+17+12)))
y_range[2] <- y_range[2] - round((diff(y_range)) * (12/(722+17+12)))

# get pixel ranges for the 12x8 cells
# we will use this as indexes in df
x_cuts <- c(rep(NA,x_range[1]), cut(x_range[1]:x_range[2],12,dig.lab = 0,include.lowest = T,labels=F))
y_cuts <- c(rep(NA,y_range[1]), cut(y_range[1]:y_range[2],8,dig.lab = 0,include.lowest = T,labels=F))

# create 12X8 matrix
df <- matrix(rep(0,12*8),nrow=8,dimnames = list(levels(cut(y_range,8,dig.lab = 0)),
                                                levels(cut(x_range,12,dig.lab = 0,include.lowest = T))))
# now go through lighted objects
for (i in 1:(dim(table(img))-1)) { 
  # img == i is light nr i
  # get position of object
  pos <- which(img@.Data[,,1] == i, arr.ind = T)

  # add up enlighted pixels in df
  for (row in 1:nrow(pos)) {
    df[y_cuts[pos[row,2]], x_cuts[pos[row,1]]] <- df[y_cuts[pos[row,2]], x_cuts[pos[row,1]]] + 1
  }
}

print(df)

>                  [143,235] (235,326] (326,417] (417,508] (508,599] (599,690] (690,781] (781,872] (872,963] (963,1.05e+03] (1.05e+03,1.14e+03] (1.14e+03,1.24e+03]
> (90,1.8e+02]           1116         0         0         0      1974         0         0         0         0              0                   0                   0
> (1.8e+02,2.7e+02]         0      2528         0         0         0      2528         0         0         0              0                1938                 845
> (2.7e+02,3.6e+02]         0      2479      2518      2430         0         0         0         0         0           2477                2409                   0
> (3.6e+02,4.5e+02]      1731      2339         0      2601         0      2775         0         0         0              0                   0                   0
> (4.5e+02,5.4e+02]         0         0         0      2549         0      1033         0         0         0              0                   0                   0
> (5.4e+02,6.3e+02]      2370      2449      2570         0         0         0         0      2555         0              0                   0                   0
> (6.3e+02,7.3e+02]         0         0         0         0         0         0         0         0         0              0                1836                2348
> (7.3e+02,8.2e+02]         0      1447      1760         0         0         0         0         0         0              0                2303                   0