从R中的二进制矩阵中提取子矩阵

Question

说二进制矩阵m：

      # [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
 # [1,]    0    0    0    0    0    0    0    0    0
 # [2,]    0    0    0    0    0    0    0    0    0
 # [3,]    0    0    0    1    1    1    1    0    0
 # [4,]    0    0    0    1    1    1    1    0    0
 # [5,]    0    0    0    1    1    1    1    0    0
 # [6,]    0    0    0    0    0    0    0    0    0
 # [7,]    0    1    1    0    0    0    0    1    1
 # [8,]    0    1    1    0    1    1    0    1    1
 # [9,]    0    0    0    0    1    1    0    1    1
# [10,]    0    0    0    0    1    1    0    0    0

m <- structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 
1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 
0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 
1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 
0, 0, 0, 0, 0, 0, 1, 1, 1, 0), .Dim = c(10L, 9L))

我们如何提取那些1-valued 子矩阵？例如

m[7:9,8:9]

#     [,1] [,2]
#[1,]    1    1
#[2,]    1    1
#[3,]    1    1

关键是我想通过算法提取它们，而不是像 m[7:9,8:9] 那样明确索引它们。

• 输入是二进制矩阵
• 作为输出的子矩阵列表（因此，四个昏暗矩阵的列表3*4、2*2、3*2 和3*2）
• 子矩阵是1-valued 矩形
• 子矩阵的边界用零保护。

Answer 1

我会将其视为一个空间问题，您有一个栅格并想要检测连接单元的区域。

library(raster)
r <- raster(m)

library(igraph)
rc <- clump(r)

plot(rc, col = rainbow(rc@data@max))

m1 <- as.matrix(rc)

lapply(seq_len(rc@data@max), function(x) {
  inds <- which(m1 == x, arr.ind = TRUE)
  nrow <- diff(range(inds[, "row"])) + 1
  ncol <- diff(range(inds[, "col"])) + 1
  matrix(1, ncol = ncol, nrow = nrow)
})
#[[1]]
#     [,1] [,2] [,3] [,4]
#[1,]    1    1    1    1
#[2,]    1    1    1    1
#[3,]    1    1    1    1
#
#[[2]]
#     [,1] [,2]
#[1,]    1    1
#[2,]    1    1
#
#[[3]]
#     [,1] [,2]
#[1,]    1    1
#[2,]    1    1
#[3,]    1    1
#
#[[4]]
#     [,1] [,2]
#[1,]    1    1
#[2,]    1    1
#[3,]    1    1

Answer 2

在光栅包中使用focal 和适当的权重矩阵w。它。将w 与m 卷积，得到一个与m 相同维度的矩阵，每个左上角的值big 和其他位置的其他值，因此将其与big 进行比较得到一个逻辑矩阵，该矩阵在上部为真矩形的左角。使用which，我们得到rc，每个矩形有一行，两列代表该矩形左上角的i和j坐标。 Map 调用遍历左上角坐标，在每个坐标上调用 genmap。 genmap 使用 rle（在 rl 函数中定义）来查找每个坐标方向上的行的长度，并返回具有这些维度的矩阵。

library(raster)

big <- 100
r <- raster(m)
w <- matrix(0, 3, 3); w[1:2, 1:2] <- 1; w[2, 2] <- big
rc <- which(as.matrix(focal(r, w, pad = TRUE, padValue = 0)) == big, arr = TRUE)

rl <- function(x) rle(x)$lengths[1]
genmat <- function(i, j) matrix(1, rl(m[i:nrow(m), j]), rl(m[i, j:ncol(m)]))
Map(genmat, rc[, 1], rc[, 2])

给予：

[[1]]
     [,1] [,2]
[1,]    1    1
[2,]    1    1

[[2]]
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    1    1    1
[3,]    1    1    1    1

[[3]]
     [,1] [,2]
[1,]    1    1
[2,]    1    1
[3,]    1    1

[[4]]
     [,1] [,2]
[1,]    1    1
[2,]    1    1
[3,]    1    1

更新简化代码。

Answer 3

一个相当冗长的答案，但您可以像我在SO answer 中所做的那样通过图像标记来做到这一点。这将很好地扩展到 1 的非矩形斑点。

find.contiguous <- function(img, x, bg) {
  ## we need to deal with a single (row,col) matrix index
  ## versus a collection of them in a two column matrix separately.
  if (length(x) > 2) {
    lbl <- img[x][1]
    img[x] <- bg
    xc <- x[,1]
    yc <- x[,2]
  } else {
    lbl <- img[x[1],x[2]]
    img[x[1],x[2]] <- bg
    xc <- x[1]
    yc <- x[2]
  }    
  ## find all neighbors of x
  xmin <- ifelse((xc-1) < 1, 1, (xc-1))
  xmax <- ifelse((xc+1) > nrow(img), nrow(img), (xc+1))
  ymin <- ifelse((yc-1) < 1, 1, (yc-1))
  ymax <- ifelse((yc+1) > ncol(img), ncol(img), (yc+1))
  ## find all neighbors of x
  x <- rbind(cbind(xmin, ymin),
             cbind(xc  , ymin),
             cbind(xmax, ymin),
             cbind(xmin, yc),
             cbind(xmax, yc),
             cbind(xmin, ymax),
             cbind(xc  , ymax),
             cbind(xmax, ymax))
  ## that have the same label as the original x
  x <- x[img[x] == lbl,]
  ## if there is none, we stop and return the updated image
  if (length(x)==0) return(img);
  ## otherwise, we call this function recursively
  find.contiguous(img,x,bg)
}

find.contiguous 是一个递归函数，对于它收到的每个调用：

1. 图片img的工作副本。
2. 像素（矩阵）索引x (row,col) 的集合，属于图像img 中的对象。
3. 背景值bg

find.contiguous 然后继续：

1. 将img 中x 处的所有像素设置为bg 颜色。这标志着我们已经访问了像素。
2. 查找x 的所有相邻像素，其标签（值）与x 中的标签（值）相同。这增加了同一对象的区域。请注意，由于x 不一定是单个像素，所以x 会几何地增长，因此，事实上，这个函数并没有懈怠。
3. 如果没有更多的邻居属于同一个对象，我们返回更新后的图像；否则，我们进行递归调用。

从对应于对象的单个像素开始，对find.contiguous 的调用将扩大区域以包含对象的所有像素，并返回更新后的图像，其中对象被背景替换。然后可以循环重复此过程，直到图像中不再有对象，从而能够提取所有 1 的子矩阵。

使用您的数据：

m <- structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 
                 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 
                 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 
                 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 
                 0, 0, 0, 0, 0, 0, 1, 1, 1, 0), .Dim = c(10L, 9L))
## make a copy to img which will be converted to all-zeros in the process 
## as matrices of 1's are extracted by the process
img <- m
## get all pixel coordinates that are objects
x <- which(img==1, arr.ind=TRUE)
## loop until there are no more pixels that are objects
##the output is in the list out
count <- 0
out <- list()
while (length(x) > 0) {
  ## choose a single (e.g., first) pixel location. This belongs to the current
  ## object that we will grow and remove from the image using find.contiguous
  if (length(x) > 2) {
    x1 <- x[1,]
  }
  ## make the call to remove the object from img
  img <- find.contiguous(img, x1, 0)
  ## find the remaining pixel locations belonging to objects
  xnew <- which(img==1, arr.ind=TRUE)
  count <- count + 1
  ## extract the indices for the 1's found by diffing new with x 
  out.ind <- x[!(x[,1] %in% xnew[,1] & x[,2] %in% xnew[,2]),]
  ## set it as a matrix in the output
  out[[count]] <- matrix(m[out.ind],nrow=length(unique(out.ind[,1])),ncol=length(unique(out.ind[,2])))
  x <- xnew
}

您的输出是列表out：

print(out)
##[[1]]
##     [,1] [,2]
##[1,]    1    1
##[2,]    1    1
##
##[[2]]
##     [,1] [,2] [,3] [,4]
##[1,]    1    1    1    1
##[2,]    1    1    1    1
##[3,]    1    1    1    1
##
##[[3]]
##     [,1] [,2]
##[1,]    1    1
##[2,]    1    1
##[3,]    1    1
##
##[[4]]
##     [,1] [,2]
##[1,]    1    1
##[2,]    1    1
##[3,]    1    1

请注意，您可以轻松地从out.ind 输出提取的 1 的位置：