加快循環和其中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確定，座標i可以是 ，範圍從1到n。總共，該框包含n^3個獨特點。 每個點具有相關聯的整數值的問題0，1或2。

：

隨着盒具有n = 256 對於具有0值的每個點P，發現其K最近 NON -ZERO-VALUED鄰居並用該鄰居的值更新P. 更新後，框中的每個點都應該不爲零。

解決方案：

我們框中有16777216（256^3）分那麼暴力方法都出來了。 幸運的是，這正是K-d樹木可用於 https://en.wikipedia.org/wiki/K-d_tree。 有幾個R庫專注於度量數據結構。 我在這個例子中使用FNN，因爲我認爲它比替代品https://cran.r-project.org/web/packages/FNN/index.html擁有更強大的API 。

的代碼：

箱子被表示爲具有列名（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。索引函數接受盒矩陣和查詢矩陣（盒矩陣的子集） 作爲參數，並返回查詢中每個點的最近鄰索引。

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 

希望這是有用的和接近您的性能預期。