R语言：将带有规则点的图像转换为数据框
R language: Converting images with regular dots to data frame
我有很多紫外线图像 (*.png)，每个图像有 96 个 "wells" 显示为圆圈。这些孔以 8 x 12 的方式排列。请看下面的例子。
在每张图片中，96 个圆圈中的一些点亮了（对紫外线有反应），而另一些则没有。我想给每个圆圈一个数字，然后将它们标识为 "lighted up" 或 "not lighted up"（具有预定义的截止值）。
您认为在 R 中实现此目的的最简单方法是什么？
我一直在玩这个包 imager 但没有成功。
请注意：并非我所有的图像都具有相同的放大倍率（即，圆圈在文件中的大小并不总是相同，但在一个文件中，它们的大小大致相同）。
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
这可能有点晚了，但这个问题解决了我经常遇到的问题。也许其他人可能会对更通用的解决方案感兴趣。这是一个建议的解决方案，它使用 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)
我有很多紫外线图像 (*.png)，每个图像有 96 个 "wells" 显示为圆圈。这些孔以 8 x 12 的方式排列。请看下面的例子。
在每张图片中，96 个圆圈中的一些点亮了（对紫外线有反应），而另一些则没有。我想给每个圆圈一个数字，然后将它们标识为 "lighted up" 或 "not lighted up"（具有预定义的截止值）。
您认为在 R 中实现此目的的最简单方法是什么？
我一直在玩这个包 imager 但没有成功。
请注意：并非我所有的图像都具有相同的放大倍率（即，圆圈在文件中的大小并不总是相同，但在一个文件中，它们的大小大致相同）。
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
这可能有点晚了，但这个问题解决了我经常遇到的问题。也许其他人可能会对更通用的解决方案感兴趣。这是一个建议的解决方案，它使用 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)