如何修改R语言DFS算法识别矩阵中最大独立连通簇？

一、DFS算法的无嵌套大簇识别优化 #

我在R语言中生成了如下10×10二进制矩阵：

set.seed(123)
matrix_1 <- matrix(rbinom(100, 1, 0.5), nrow = 10, ncol = 10)

我编写了一个基于深度优先搜索（DFS）的算法来识别矩阵中值为1的连通簇，定义簇为满足8连通性（含对角线）且最小大小为3的连续单元。该算法运行速度快，但会返回所有符合条件的簇（共44个，包含嵌套在大簇内的小簇），而我通过可视化仅观察到3个符合大小要求的大簇。请问如何修改该函数，使其仅识别矩阵中最大的无嵌套连通簇？

解决方案 #

要实现仅识别无嵌套的最大连通簇，核心是避免重复标记已属于大簇的单元格，具体修改思路如下：

• 标记已访问单元格：遍历过程中，一旦某个单元格被归入簇，就标记为已访问，后续不再处理，彻底避免嵌套簇生成；
• 筛选最大簇：遍历完成后，先保留所有符合最小大小的簇，再从中提取尺寸最大的簇（多个同尺寸最大簇全部保留）。

修改后的实现示例（用BFS实现，与DFS逻辑等价，仅遍历顺序不同）：

find_largest_non_nested_clusters <- function(mat, min_size = 3, directions = 8) {
  rows <- nrow(mat)
  cols <- ncol(mat)
  visited <- matrix(FALSE, nrow = rows, ncol = cols)
  clusters <- list()
  
  # 定义8连通邻域的坐标偏移（排除自身）
  offsets <- expand.grid(c(-1,0,1), c(-1,0,1))[-5,]
  
  for (i in 1:rows) {
    for (j in 1:cols) {
      if (mat[i,j] == 1 && !visited[i,j]) {
        current_cluster <- list(c(i,j))
        visited[i,j] <- TRUE
        queue <- list(c(i,j))
        
        # 遍历当前单元格的所有邻域
        while (length(queue) > 0) {
          cell <- queue[[1]]
          queue <- queue[-1]
          
          for (k in 1:nrow(offsets)) {
            ni <- cell[1] + offsets[k,1]
            nj <- cell[2] + offsets[k,2]
            
            # 检查坐标合法性、未访问且值为1
            if (ni >=1 && ni <= rows && nj >=1 && nj <= cols && 
                !visited[ni,nj] && mat[ni,nj] == 1) {
              visited[ni,nj] <- TRUE
              current_cluster <- c(current_cluster, list(c(ni,nj)))
              queue <- c(queue, list(c(ni,nj)))
            }
          }
        }
        
        # 仅保留符合最小大小要求的簇
        if (length(current_cluster) >= min_size) {
          clusters <- c(clusters, list(current_cluster))
        }
      }
    }
  }
  
  # 筛选最大簇
  if (length(clusters) == 0) return(list())
  cluster_sizes <- sapply(clusters, length)
  max_size <- max(cluster_sizes)
  
  return(clusters[cluster_sizes == max_size])
}

# 测试调用
set.seed(123)
matrix_1 <- matrix(rbinom(100, 1, 0.5), nrow = 10, ncol = 10)
largest_clusters <- find_largest_non_nested_clusters(matrix_1)

二、terra包双通道簇识别的性能优化 #

我后续采纳建议使用terra包实现了针对两种通道类型的簇识别函数，处理10000个100×100矩阵耗时约3.9分钟：

library(dplyr)
library(terra)

# M: 整数类型矩阵 

find_clusters_2chan <- function(M) {
  
  # 提取值为1的单元格
  ones <- M == 1
  
  # 转换为栅格对象
  raster_ones <- ones |> rast()
  
  # 识别连通簇（将0视为NA，即断开连通）
  clusters_ones <- patches(raster_ones,
                           directions = 8,
                           zeroAsNA = TRUE)
  
  # 生成簇大小频率表
  ones_freq <- clusters_ones |> freq()
  
  # 筛选大小≥3的簇
  ONES <- ones_freq$count[ones_freq$count >= 3]
  
  #-------------------------------------------------------------------------------
  
  # 提取值为2的单元格
  twos <- M == 2
  
  # 转换为栅格对象
  raster_twos <- twos |> rast()
  
  # 识别连通簇（将0视为NA，即断开连通）
  clusters_twos <- patches(raster_twos,
                           directions = 8,
                           zeroAsNA = TRUE)
  
  # 生成簇大小频率表
  twos_freq <- clusters_twos |> freq()
  
  # 筛选大小≥3的簇
  TWOS <- twos_freq$count[twos_freq$count >= 3]
  
  clusters_list <- list(channel_1 = ONES,
                        channel_2 = TWOS)
  
  return(clusters_list)
  
}

start <- Sys.time()

clusters_big_list <- lapply(list_of_matrices, find_clusters_2chan)

end <- Sys.time()

end - start

# 运行时间 = 3.902859 分钟

性能优化建议 #

针对批量处理场景，可通过以下方式显著降低耗时：

1. 复用栅格模板
   如果所有矩阵维度一致，预先创建栅格模板，避免重复初始化：

# 预先创建100×100的栅格模板
template_rast <- rast(nrow = 100, ncol = 100)

# 修改函数内的转栅格步骤
raster_ones <- template_rast
values(raster_ones) <- as.vector(ones)

2. 并行处理
   利用多核CPU并行处理多个矩阵，替代单线程的lapply：

library(parallel)
# 创建并行集群（保留1个核心给系统）
cl <- makeCluster(detectCores() - 1)
# 导出所需资源到集群
clusterExport(cl, c("find_clusters_2chan", "template_rast"))
clusterEvalQ(cl, {library(terra); library(dplyr)})

# 并行批量处理
clusters_big_list <- parLapply(cl, list_of_matrices, find_clusters_2chan)

# 关闭集群
stopCluster(cl)

3. 简化频率筛选逻辑
   用基础R函数替代dplyr管道，减少额外开销：

# 替代原筛选步骤
ONES <- ones_freq[ones_freq$count >= 3, "count", drop = TRUE]

4. 合并重复逻辑
   将双通道的处理逻辑封装为内部函数，减少代码冗余：

find_clusters_2chan <- function(M, template_rast) {
  process_channel <- function(val) {
    channel <- M == val
    raster_channel <- template_rast
    values(raster_channel) <- as.vector(channel)
    clusters <- patches(raster_channel, directions = 8, zeroAsNA = TRUE)
    freq_tab <- freq(clusters)
    freq_tab[freq_tab$count >= 3, "count", drop = TRUE]
  }
  
  list(
    channel_1 = process_channel(1),
    channel_2 = process_channel(2)
  )
}

内容的提问来源于stack exchange，提问作者Tavaro Evanis