体素邻域索引-使用线性索引检测26个相邻体素访问中的"越界"情况。

Question

嗯，我不知道如何用标题描述我的问题，希望我得到的是正确的。

我有一个矩阵（在下面的示例中为M），它是一个3D图像，由11x11x11体素组成（我只是为了方便而逻辑化，并且大小也只是一个示例）。

在我的代码中，我需要访问一些体素的26个邻居，为此我使用了一些花哨的线性索引，这些索引可以在http://www.mathworks.com/matlabcentral/answers/86900-how-to-find-all-neighbours-of-an-element-in-n-dimensional-matrix中找到。

问题是，如果point处于M的“边界”上，则将尝试访问一些越界值，这将生成错误。

要解决这个问题，一个很好的方法是在M周围创建一个边界，在每个维度上增加+2大小，并以零填充，但是我真的想避免改变M，因为我的代码比示例中的代码复杂得多。

我找不到任何方法来做到这一点，我有点卡在这里。 有什么建议吗？

编辑：@Dan的答案有效，但是我想看看是否存在可能使用此线性索引方法的解决方案。

% Example data
M=round(randn(11,11,11))~=0;

% Fancy way of storing 26 neigh indices for good accesing 
s=size(M);
N=length(s);
[c1{1:N}]=ndgrid(1:3);
c2(1:N)={2};
neigh26=sub2ind(s,c1{:}) - sub2ind(s,c2{:});

point=[5 1 6];

% This will work unless the point is in the boundary (like in this example)
neighbours=M(sub2ind(s,point(1),point(2),point(3))+neigh26) 

- Ander Biguri

---

@Divakar 理想情况是在这些情况下有不同的大小（neigh26）。因此，如果该点在一个边界内，neigh26（忘记变量名）将是一个2x3x3矩阵，在边缘的情况下为2x2x3矩阵，在角落的情况下为2x2x2矩阵。 - Ander Biguri

@AnderBiguri 你是不是指的边缘情况是 2x3x3，而角落情况是 2x2x3？此外，如果你正在处理多个体素，由于大小不同，因此如果你想一次性获取所有体素的邻居而不使用循环，我认为你肯定需要使用单元数组。 - Divakar

@Divakar 更具体地说，在理想情况下，如果 neigh26 仅与 X 平面接触，则为 2x3x3，如果在 Y 平面中，则为 3x2x3，如果在 Z 中，则为 3x3x2（对于 2 平面情况进行推断）。然而，这可能过于复杂。只要我知道哪些点没有被访问，如果 neigh26 始终是 3x3x3，并且在不使用 neigh 的位置上有 0，那么这不是问题。 - Ander Biguri

@AnderBiguri 抱歉，我必须离开我的系统。那么我们用NaN填充，这样我们将拥有M_padded作为13x13x13，并索引到M_padded而不是M。索引后，NaN元素将指示边界point的"无效"元素。这种方法适合你吗？ - Divakar

@AnderBiguri 没关系，我想我能处理 :) 好吧，我有一段使用 bsxfun 的有趣代码，我想我可以和你一起运行。 - Divakar

Answer 1

那个线性索引的东西很重要吗？如果你使用下标索引和min和max，处理边界条件就很容易了：

p = [5, 1, 6];

neighbourhood = M(max(1,p(1)-1)):min(p(1)+1,end),
                  max(1,p(2)-1)):min(p(2)+1,end),
                  max(1,p(3)-1)):min(p(3)+1,end))

%// Get rid of the point it self (i.e. the center)
neighbours = neighbourhood([1:13, 15:end])

如果你想要更广泛的邻域，也可以轻松地推广这种方法：

p = [5, 1, 6];
n = 2;
neighbourhood = M(max(1,p(1)-n)):min(p(1)+n,end),
                  max(1,p(2)-n)):min(p(2)+n,end),
                  max(1,p(3)-n)):min(p(3)+n,end))

%// Get rid of the point it self (i.e. the center)
mid = ceil(numel(neigbourhood)/2);
neighbours = neighbourhood([1:mid-2, mid+1:end])

如果你想保持魔方的形状，那么也许可以这样做：

neighbours = neighbourhood;
neighbours(mid) = NaN;

如果您想在代码中多次使用此功能，最好将其重构为一个m文件函数，该函数仅返回索引：

function ind = getNeighbours(M,p,n)
    M = zeros(size(M));
    M(max(1,p(1)-n)):min(p(1)+n,end), max(1,p(2)-n)):min(p(2)+n,end), max(1,p(3)-n)):min(p(3)+n,end)) = 1;
    M(p(1), p(2), p(3)) = 0;
    ind = find(M);
end

- Dan

---

然而，'end' 的使用非常棒，你应该知道 ;) - Ander Biguri

@AnderBiguri 添加了匿名函数，希望没有打错字！ - Dan

@AnderBiguri 哎呀，我把 n 加成了一个参数。我也会再加上一个只有26个邻居的版本。 - Dan

我一直在仔细查看我的代码，尝试将你的答案融入其中。但是，匿名函数对我没用。我真的需要与这些邻居进行相当多的交互，因此使用线性索引方法说M1（p + neigh26）= M2(p + neigh26)+1; （显然只是一个例子）太容易了。下标索引会使事情变得更加复杂（但当然可以工作），但匿名函数在这里行不通。但是我真的很感激你为给我这个答案所做的努力！ - Ander Biguri

@AnderBiguri 我明白了。最后一次尝试，我添加了一个帮助器m文件函数。你可以像这样使用它：M1(getNeighbours(M1,p,1)) = ... - Dan

Answer 2

基本理论：将输入数组扩展到左右、上下、第三维的每一侧，使用NaN填充。这将使我们能够使用统一的3x3x3网格，然后稍后使用那些NaN来检测超出输入数组边界的元素，并因此丢弃它们。

代码：

%// Initializations
sz_ext = size(M)+2; %// Get size of padded/extended input 3D array
M_ext = NaN(sz_ext); %// Initialize extended array
M_ext(2:end-1,2:end-1,2:end-1) = M; %// Insert values from M into it

%// Important stuff here : Calculate linear offset indices within one 3D slice
%// then for neighboring 3D slices too
offset2D = bsxfun(@plus,[-1:1]',[-1:1]*sz_ext(1)); %//'
offset3D = bsxfun(@plus,offset2D,permute([-1:1]*sz_ext(1)*sz_ext(2),[1 3 2]));

%// Get linear indices for all points
points_linear_idx = sub2ind(size(M_ext),point(:,1)+1,point(:,2)+1,point(:,3)+1);
%// Linear indices for all neighboring elements for all points; index into M_ext
neigh26 = M_ext(bsxfun(@plus,offset3D,permute(points_linear_idx,[4 3 2 1])))

如何使用： 因此，第四维中的每个切片表示27个元素（相邻元素加上元素本身），形式为3x3x3数组。因此，neigh26将是一个3x3x3xN数组，其中N是点阵列中点的数量。

示例： 例如，假设在M和Point中有一些随机值 -

M=rand(11,11,11);
point = [
    1 1 4;
    1 7 1]

在使用先前的代码和这些输入运行时，我会得到类似于这样的结果 -

neigh26(:,:,1,1) =
       NaN       NaN       NaN
       NaN    0.5859    0.4917
       NaN    0.6733    0.6688
neigh26(:,:,2,1) =
       NaN       NaN       NaN
       NaN    0.0663    0.5544
       NaN    0.3440    0.3664
neigh26(:,:,3,1) =
       NaN       NaN       NaN
       NaN    0.3555    0.1257
       NaN    0.4424    0.9577
neigh26(:,:,1,2) =
   NaN   NaN   NaN
   NaN   NaN   NaN
   NaN   NaN   NaN
neigh26(:,:,2,2) =
       NaN       NaN       NaN
    0.7708    0.3712    0.2866
    0.7088    0.3743    0.2326
neigh26(:,:,3,2) =
       NaN       NaN       NaN
    0.4938    0.5051    0.9416
    0.1966    0.0213    0.8036

- Divakar