使用 R 加速循环和条件

Question

我想在 R 中加速这段代码。

输入是一个包含整数的 3x3x3 数组，并基于邻居，如果它们为零，则将它们替换为相应的数字。

输出是带有新值的数组“mask_roi”。

###### Start here

list_neig = array(0, dim = c(3,3,3))

mask_roi = array(sample(c(0,1,2),27,replace=T), dim = c(3,3,3))

values_mask = array(1:27, dim = c(3,3,3))

values_mask_melted = melt(values_mask, varnames=c("x","y","z"))

### Tranform the 3D Matrix in a data.table wit 4 columns position and value
image_melted  <- melt(mask_roi, varnames=c("x","y","z"))  # 4 columns: x, y, z, value

image_melted$box = rownames(image_melted)

image_melted_non_zeros<-image_melted[!(image_melted$value==0),]

box_neigbors = vector("list", nrow(image_melted))

for (i in 1:(nrow(image_melted_non_zeros))){
  cat(i,"\n")
  x = image_melted_non_zeros[i,1]
  y = image_melted_non_zeros[i,2]
  z = image_melted_non_zeros[i,3] 

  box_neigbors[[image_melted_non_zeros[i,5]]] <- list(nearestNeighbors(values_mask, elem = c(x,y,z), dist = 1,dim = c(3,3,3)))

}

我已经完成了“box_neighbors”向量，只是将它包含在此处以说明如何获取它，我们需要从这里到结束更快。这个想法是，检查所有不同于零的体素并检查他的所有邻居。如果他的邻居为零，他将具有相同的值，如果不为零，则保持原来的值。

for (i in 1:(nrow(image_melted_non_zeros))){
  cat(i,"\n")  
  x = image_melted_non_zeros[i,1]
  y = image_melted_non_zeros[i,2]
  z = image_melted_non_zeros[i,3] 

  number_of_nei = length(box_neigbors[[image_melted_non_zeros[i,5]]][[1]] )
  value_vozel = mask_roi[x,y,z]  # it will give this new value 

  for (j in 1:number_of_nei){
    nei_number = box_neigbors[[image_melted_non_zeros[i,5]]][[1]][j]

    xx = image_melted[nei_number,1]
    yy = image_melted[nei_number,2]
    zz = image_melted[nei_number,3]     

    value_nei = mask_roi[xx,yy,zz]

    if(value_nei == 0){      
      mask_roi[xx,yy,zz] = value_vozel
    }  
  }  
}

我需要为 256x256x256 阵列而不是 3x3x3 执行此操作。

非常感谢！

nearestNeighbors <- function(ary, elem, dist, dims){
  usedims <- mapply(function(el, d) {
    seq(max(1, el - dist), min(d, el + dist))
  }, elem, dims, SIMPLIFY=FALSE)
  df <- as.matrix(do.call('expand.grid', usedims))
  ndist <- sqrt(apply(df, 1, function(x) sum((x - elem)^2)))
  ret <- df[which(ndist > 0 & ndist <= dist),,drop = FALSE]

  return(ary[ret])

}

Answer 1

我整理了一个使用 K-d 树的实现。在配备 16GB RAM 和 2.3 GHz i.7 处理器的 MacBookPro 上运行时，它可以在大约 13 秒内处理 256x256x256 阵列。您没有给出任何具体的基准，但我认为 13s 足以发布答案。我在下面概述了我的步骤。如果我误解了问题的一部分，请告诉我。

设置：

我们有一个边长为 n 的盒子，里面装满了点。 方框中的一个点由坐标 i,j,k 确定，它可以 范围从 1 到 n。总共，该框包含 n^3 个唯一点。 每个点都有一个相关的整数值 0、1 或 2。

问题：

带有 n = 256 的框。 对于每个具有 0 值的点 P，找到其最近的 k 非零值邻居并使用该邻居的值更新 P。 更新后方框中的每个点都应该是非零的。

解决方案：

我们的盒子有 16,777,216 (256^3) 个点，所以蛮力方法不可用。 幸运的是，这正是 K-d 树的用途 https://en.wikipedia.org/wiki/K-d_tree。 有一些 R 库专注于度量数据结构。 我在此示例中使用 FNN，因为我认为它具有更强大的 API 比替代品https://cran.r-project.org/web/packages/FNN/index.html.

守则：

框表示为具有列名（i、j、k、值）的矩阵。 每行代表框中的一个点。

set.seed(256)
library(FNN)
len = 256
values = c(0, 1, 2)
createBox = function(n, vals) {
    index = 1:len^3
    value = sample(vals, length(index), replace = T)
    box = as.matrix(cbind(index, index, index, value))
    dimnames(box) = list(NULL, c("i", "j", "k", "value"))
    box
}
box= createBox(len, values)

knnx.index 函数接受框矩阵和查询矩阵（框矩阵的子集） 作为参数并返回查​​询中每个点的最近邻索引。

updateZeroValuedPoints = function(box, kval) {
  zeroPointIndx = which(box[ , "value"] == 0)
  nonZeroPoints = box[-1 * zeroPointIndx, ]
  zeroPoints = box[zeroPointIndx, ]
  nnIdx = knnx.index(nonZeroPoints, zeroPoints, k = kval, algorithm = "kd_tree")
  zeroPoints[, "value"] = nonZeroPoints[nnIdx[ , ncol(nnIdx)], "value"]
  zeroPoints
}

一旦你有了邻居索引，就可以直接交换更新值，不需要 for 循环。

system.time(updateZeroValuedPoints(box, 1))
# > system.time(updateZeroValuedPoints(box, 1))
# user  system elapsed
# 13.517   1.162  14.676

希望这很有用，并且接近您的性能预期。