分割环，即 R 中的非完整对象（在 EBIimage 或其他中）

Question

我依靠边缘检测（而不是颜色检测）从血细胞中提取特征。原图如下：

我正在使用 R EBImage 包来运行 sobel + 低通滤波器来达到这样的效果：

library(EBImage)
library(data.table)
img <- readImage("6hr-007-DIC.tif")

#plot(img)


#print(img, short = T)

# 1. define filter for edge detection
hfilt <- matrix(c(1, 2, 1, 0, 0, 0, -1, -2, -1), nrow = 3) # sobel

# rotate horizontal filter to obtain vertical filter
vfilt <- t(hfilt)

# get horizontal and vertical edges
imgH <- filter2(img, hfilt, boundary="replicate")
imgV <- filter2(img, vfilt, boundary="replicate")

# combine edge pixel data to get overall edge data
hdata <- imageData(imgH)
vdata <- imageData(imgV)
edata <- sqrt(hdata^2 + vdata^2)

# transform edge data to image
imgE <- Image(edata)
#print(display(combine(img, imgH, imgV, imgE), method = "raster", all = T))

display(imgE, method = "raster", all = T)


# 2. Enhance edges with low pass filter

hfilt <- matrix(c(1, 1, 1, 1, 1, 1, 1, 1, 1), nrow = 3) # low pass

# rotate horizontal filter to obtain vertical filter
vfilt <- t(hfilt)

# get horizontal and vertical edges
imgH <- filter2(imgE, hfilt, boundary="replicate")
imgV <- filter2(imgE, vfilt, boundary="replicate")

# combine edge pixel data to get overall edge data
hdata <- imageData(imgH)
vdata <- imageData(imgV)
edata <- sqrt(hdata^2 + vdata^2)


# transform edge data to image
imgE <- Image(edata)
plot(imgE)

我想知道是否有任何方法可以填充大环（血细胞）中的孔，使其成为固体，有点像：

（显然这不是同一张图片，但想象一下最后一张图片只是从边缘开始的。）

然后我想使用 EBImage 包中的 computeFeatures() 之类的方法（据我所知，它仅适用于实体）

Answer 1

编辑更多的代码来提取具有“连接”到边框的对象内部。附加代码包括定义分段单元的凸包和创建填充蒙版。

简短的回答是fillHull 和floodFill 可能有助于填充具有明确边界的单元格。

下面较长（已编辑）的答案建议使用floodFill 的方法可能有用。您在从低对比度 DIC 图像中提取信息方面做得很好，但更多的图像处理可能会有所帮助，例如对嘈杂的 DIC 图像进行“平场校正”。这个原理在Wikipedia page 中有描述，但是一个简单的实现就可以创造奇迹。此处建议的编码解决方案需要用户交互来选择单元格。这不是一个强大的方法。尽管如此，也许更多的图像处理结合代码来定位细胞可能会起作用。最后，细胞内部被分割出来，可以通过computeFeatures进行分析。

代码从阈值图像开始（修剪边缘并转换为二进制）。

# Set up plots for 96 dpi images
  library(EBImage)
  dm <- dim(img2)/96
  dev.new(width = dm[1], height = dm[2])
# Low pass filter with gblur and make binary
  xb <- gblur(img2, 3)
  xt <- thresh(xb, offset = 0.0001)
  plot(xt) # thresh.jpg
# dev.print(jpeg, "thresh.jpg", width = dm[1], unit = "in", res = 96) 

# Keep only "large" objects
  xm <- bwlabel(xt)
  FS <- computeFeatures.shape(xm)
  sel <- which(FS[,"s.area"] < 800)
  xe <- rmObjects(xm, sel)

# Make binary again and plot
  xe <- thresh(xe)
  plot(xe) # trimmed.jpg
#  dev.print(jpeg, "trimmed.jpg", width = dm[1], unit = "in", res = 96)

# Choose cells with intact interiors
# This is done by hand here but with more pre-processing, it may be
# possible to have the image suitable for more automated analysis...
  pp <- locator(type = "p", pch = 3, col = 2) # marked.jpg
#  dev.print(jpeg, "marked.jpg", width = dm[1], unit = "in", res = 96)

# Fill interior of each cell with a unique integer
  myCol <- seq_along(pp$x) + 1
  xf1 <- floodFill(xe, do.call(rbind, pp), col = myCol)

# Discard original objects from threshold (value = 1) and see
  cells1 <- rmObjects(xf1, 1)
  plot(colorLabels(cells1))
# dev.print(jpeg, "cells1.jpg", width = dm[1], unit = "in", res = 96)

我需要引入算法来连接顶点之间的整数点并填充一个凸多边形。此处的代码实现了 Bresenham 算法，并使用了仅适用于凸（简单）多边形的简单多边形填充例程。

#
# Bresenham's balanced integer line drawing algorithm
#
bresenham <- function(x, y = NULL, close = TRUE)
{
# accept any coordinate structure
  v <- xy.coords(x = x, y = y, recycle = TRUE, setLab = FALSE)
  if (!all(is.finite(v$x), is.finite(v$y)))
    stop("finite coordinates required")

  v[1:2] <- lapply(v[1:2], round) # Bresenham's algorithm IS for integers
  nx <- length(v$x)
  if (nx == 1) return(list(x = v$x, y = v$y)) # just one point
  if (nx > 2 && close == TRUE) { # close polygon by replicating 1st point
    v$x <- c(v$x, v$x[1])
    v$y <- c(v$y, v$y[1])
    nx <- nx + 1
  }
# collect result in 'ans, staring with 1st point
  ans <- lapply(v[1:2], "[", 1)

# process all vertices in pairs
  for (i in seq.int(nx - 1)) {
    x <- v$x[i] # coordinates updated in x, y
    y <- v$y[i]
    x.end <- v$x[i + 1]
    y.end <- v$y[i + 1]

    dx <- abs(x.end - x); dy <- -abs(y.end - y)
    sx <- ifelse(x < x.end, 1, -1)
    sy <- ifelse(y < y.end, 1, -1)
    err <- dx + dy

  # process one segment
    while(!(isTRUE(all.equal(x, x.end)) && isTRUE(all.equal(y, y.end)))) {
      e2 <- 2 * err
      if (e2 >= dy) { # increment x
        err <- err + dy
        x <- x + sx
      }
      if (e2 <= dx) { # increment y
        err <- err + dx
        y <- y + sy
      }
      ans$x <- c(ans$x, x)
      ans$y <- c(ans$y, y)
    }
  }
# remove duplicated points (typically 1st and last)
  dups <- duplicated(do.call(cbind, ans), MARGIN = 1) 
  return(lapply(ans, "[", !dups))
}

还有一个简单的例程来查找简单多边形的内部点。

#
# Return x,y integer coordinates of the interior of a CONVEX polygon
#
cPolyFill <- function(x, y = NULL) 
{
  p <- xy.coords(x, y = y, recycle = TRUE, setLab = FALSE)
  p[1:2] <- lapply(p[1:2], round)
  nx <- length(p$x)
  if (any(!is.finite(p$x), !is.finite(p$y)))
    stop("finite coordinates are needed")

  yc <- seq.int(min(p$y), max(p$y))
  xlist <- lapply(yc, function(y) sort(seq.int(min(p$x[p$y == y]), max(p$x[p$y == y]))))
  ylist <- Map(rep, yc, lengths(xlist))
  ans <- cbind(x = unlist(xlist), y = unlist(ylist))
  return(ans)
}

现在这些可以与ocontour() 和chull() 一起使用来创建和填充每个分段单元格周围的凸包。这会“修复”那些被入侵的单元格。

# Create convex hull mask
  oc <- ocontour(cells1) # for all points along perimeter
  oc <- lapply(oc, function(v) v + 1) # off-by-one flaw in ocontour
  sel <- lapply(oc, chull) # find points that define convex hull
  xh <- Map(function(v, i) rbind(v[i,]), oc, sel) # new vertices for convex hull
  oc2 <- lapply(xh, bresenham) # perimeter points along convex hull

# Collect interior coordinates and fill
  coords <- lapply(oc2, cPolyFill)
  cells2 <- Image(0, dim = dim(cells1))
  for(i in seq_along(coords))
    cells2[coords[[i]]] <- i # blank image for mask
  xf2 <- xe
  for (i in seq_along(coords))
    xf2[coords[[i]]] <- i # early binary mask

# Compare before and after
  img <- combine(colorLabels(xf1), colorLabels(cells1),
    colorLabels(xf2), colorLabels(cells2))
  plot(img, all = T, nx = 2)
  labs <- c("xf1", "cells1", "xf2", "cells2")
  ix <- c(0, 1, 0, 1)
  iy <- c(0, 0, 1, 1)
  text(dm[1]*96*(ix + 0.05), 96*dm[2]*(iy + 0.05), labels = labs,
    col = "white", adj = c(0.05,1))
 # dev.print(jpeg, "final.jpg", width = dm[1], unit = "in", res = 96)