在 MATLAB 中构建 3D 点阵图

Question

我对扩展此问题/答案 (https://stackoverflow.com/a/3283732/2371031) 以将 4 连接案例扩展到第三维感兴趣。

问题 1：给定一个 X x Y x Z 维矩阵，如何构造 6 连接邻接矩阵（或边列表）？

问题 2：给定邻接矩阵（或边列表），如何将边权重计算为连接节点值的某个函数（例如均值）？

问题 3：对于非常大的矩阵，如何解决问题 1 和 2 而不会耗尽内存？稀疏矩阵或边缘列表似乎是可能的解决方案，但如何实现？

% simple case
X = 2; Y = 2; Z = 2;
mat = reshape([1:(X*Y*Z)],X,Y,Z);
>> mat

mat(:,:,1) =

     1     3
     2     4


mat(:,:,2) =

     5     7
     6     8


% large case
X = 256; Y = 256; Z = 1000;
mat = reshape([1:(X*Y*Z)],X,Y,Z);

Answer 1

这种方法使用ndgrid 在一个方向上查找每个节点的邻居。例如，在维度 2 中，我们找到每个节点右侧的（有效）邻居。

这里，蓝色节点是边缘的源，它们右边的红色节点是目标节点。因此源节点和目标节点的网格相互偏移 1。

% Given a 3d matrix M, generate a graph representing
% the 6-connected neighbors.
% Node numbers are equal to node indices in M.
% Edge weights are given by mean of adjacent node values.
[m, n, p] = size(M);

% Initialize Source and Target node vectors and corresponding Weight
S = T = W = [];

% List neighboring nodes where S(i) < T(i)
% Neighbors along dimension 1
[X1, X2, X3] = ndgrid(1:m-1, 1:n, 1:p);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(2:m, 1:n, 1:p);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Neighbors along dimension 2
[X1, X2, X3] = ndgrid(1:m, 1:n-1, 1:p);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(1:m, 2:n, 1:p);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Neighbors along dimension 3
[X1, X2, X3] = ndgrid(1:m, 1:n, 1:p-1);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(1:m, 1:n, 2:p);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Calculate the weight for each edge 
W = mean(M([S, T]), 2);

%% Adjacency matrix A is for an undirected graph
%%    therefore only edges where S(i) < T(i) are listed.
%%    If the backward edges are also necessary:
% tempS = S;
% S = [S;T];
% T = [T;tempS];
% W = [W;W];

% [S T W] gives the graph's edge list.
% If the sparse adjacency matrix is required:
A = sparse(S, T, W);

使用您的大型测试用例在 Octave 中运行它：

X = 256; Y = 256; Z = 1000;
M = reshape([1:(X*Y*Z)],X,Y,Z);

（无向）边缘列表大约需要 45 秒，另外需要 20 秒左右来生成稀疏邻接矩阵。由于我使用 Octave，我无法测试哪个会更快地生成图表。

Answer 2

答案 1： 与其建立一个巨大的边和权重矩阵......使用边列表直接在循环中迭代地构建图形，循环每个维度（行、列、平面）。这个 64x64x10 非常快（秒），但对于大型示例（256x256x1000）运行时间要长得多。要是能开发出一个更快的版本就好了。

X=256;Y=256;Z=1000;
% input matrix of 0-1 values 
mat = rand(X,Y,Z);
% matrix of indices for input matrix coordinates
imat = reshape([1:(X*Y*Z)],X,Y,Z);

% initialize empty graph
g = graph();

xdim = size(imat,1);
ydim = size(imat,2);
zdim = size(imat,3);

% create edge list with weights
disp('calculating start/end nodes for edges and associated weights...')

for y = 1:(ydim-1)  % loop through all "horizontal" (column) connections in each z plane
    % example:
    % 1 - 3
    % 2 - 4
    for z = 1:zdim
        nodes1 = imat(:,y,z);
        nodes2 = imat(:,y+1,z); % grab column y+1 nodes of plane z
        wts = ((1-mat(:,y,z)) + (1-mat(:,y+1,z))) ./ 2;  % average of 1-value for each node
        g = addedge(g,nodes1,nodes2,wts); % add to graph
    end
end

for x = 1:(xdim-1)  % loop through all "vertical" (row) connections within each z plane
    % example:
    % 1  3
    % |  |
    % 2  4

    for z = 1:zdim
        nodes1 = imat(x,:,z);
        nodes2 = imat(x+1,:,z); % grab row x+1 nodes of plane z
        wts = ((1-mat(x,:,z)) + (1-mat(x+1,:,z))) ./ 2;  % average of 1-value for each node
        g = addedge(g,nodes1,nodes2,wts); % add to graph
    end
end

for z = 1:(zdim-1)  % loop through all "deep" connections across z planes
    % example:
    %   5  7
    %  /  /
    % 1  3

    for x = 1:xdim
        nodes1 = imat(x,:,z);
        nodes2 = imat(x,:,z+1); % grab row x nodes of plane z+1
        wts = ((1-mat(x,:,z)) + (1-mat(x,:,z+1))) ./ 2;  % average of 1-value for each node
        g = addedge(g,nodes1,nodes2,wts); % add to graph
    end
end
disp('done.')

figure
spy(adjacency(g))

答案 2： 将链接答案的 4-connected case 扩展到第三维 (https://stackoverflow.com/a/3283732/2371031)：

disp("calculating adjacency matrix...")

X=3; Y=3; Z=2;
mat = reshape([1:(X*Y*Z)],X,Y,Z);
[x, y, z] = size(mat);                                     % Get the matrix size
diagVec1 = repmat(repmat([ones(y-1, 1); 0], x, 1), z, 1);  % Make the first diagonal vector 
                                                           % (for y-planar connections)
diagVec1 = diagVec1(1:end-1);                              % Remove the last value

diagVec2 = repmat([ones(y*(x-1), 1); zeros(x,1)], z, 1);   % Make the second diagonal vector 
                                                           % (for x-planar connections)
diagVec2 = diagVec2(1:end-x);                              % drop the last x values (zeros)

diagVec3 = ones(x*y*(z-1), 1);                             % make the third diagonal vector 
                                                           % (for z-planar connections)

adj = diag(diagVec1, 1)+diag(diagVec2, y)+diag(diagVec3, x*y);   % Add the diagonals to a zero matrix
adj = adj+adj.';                                           % Add the matrix to a transposed copy of
                                                           %   itself to make it symmetric

disp("done.")

figure
spy(adj)

对于大案例Matlab报错：

Error using diag
Requested 65536000x65536000 (32000000.0GB) array exceeds maximum array size preference. Creation of
arrays greater than this limit may take a long time and cause MATLAB to become unresponsive.

Error in blah (line X)
adj = diag(diagVec1, 1) + diag(diagVec2, y) + diag(diagVec3, x*y);   % Add the diagonals to a zero matrix

---

更新

答案 3：
（@beaker 答案的更新版本）

我很好奇是否可以通过使用 26 连通情况（所有 9 个线性或对角相邻点在上平面，9 个在下平面，8 个在重合平面）来获得改进的最短路径结果。

% 3x3x3 'cube'
X=3;Y=3;Z=3;

% input matrix of 0-1 values 
M_wts = rand(X,Y,Z);
M_wts = 1-M_wts;

% matrix of indices for input matrix coordinates
M = reshape([1:(X*Y*Z)],X,Y,Z);

% initialize empty graph
g = graph();

%% Given a 3d matrix M, generate a graph representing
% the 6-connected neighbors.
% Node numbers are equal to node indices in M.
% Edge weights are given by mean of adjacent node values.
[x, y, z] = size(M);

% Initialize Source and Target node vectors and corresponding Weight
S = []; 
T = [];

%% 6-connected case: 
% List neighboring nodes where S(i) < T(i)
% Neighbors along dimension 1 (x/lateral)
[X1, X2, X3] = ndgrid(1:x-1, 1:y, 1:z);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(2:x, 1:y, 1:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Neighbors along dimension 2 (y/vertical)
[X1, X2, X3] = ndgrid(1:x, 1:y-1, 1:z);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(1:x, 2:y, 1:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Neighbors along dimension 3 (z/planar)
[X1, X2, X3] = ndgrid(1:x, 1:y, 1:z-1);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(1:x, 1:y, 2:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

%% 
% Calculate the weight for each edge 
W = mean(M_wts([S, T]), 2);

g = graph(S, T, W);
gplot = plot(g,'EdgeLabel',g.Edges.Weight);
layout(gplot,'force3')
view(3)

%% diagonal connections
% Neighbors within dimension 3 (diagonal xy1)
[X1, X2, X3] = ndgrid(1:x-1, 1:y-1, 1:z);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(2:x, 2:y, 1:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Neighbors within dimension 3 (diagonal xy2)
[X1, X2, X3] = ndgrid(1:x-1, 2:y, 1:z);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(2:x, 1:y-1, 1:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];


%% across-plane diagonal connections
% Neighbors along dimension 3 (across-plane diagonal xz1)
[X1, X2, X3] = ndgrid(1:x, 1:y-1, 1:z-1);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(1:x, 2:y, 2:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Neighbors along dimension 3 (across-plane diagonal xz2)
[X1, X2, X3] = ndgrid(1:x, 2:y, 1:z-1);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(1:x, 1:y-1, 2:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Neighbors along dimension 3 (across-plane diagonal yz1)
[X1, X2, X3] = ndgrid(1:x-1, 1:y, 1:z-1);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(2:x, 1:y, 2:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Neighbors along dimension 3 (across-plane diagonal yz2)
[X1, X2, X3] = ndgrid(2:x, 1:y, 1:z-1);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(1:x-1, 1:y, 2:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Neighbors along dimension 3 (across-plane diagonal xyz1)
[X1, X2, X3] = ndgrid(1:x-1, 1:y-1, 1:z-1);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(2:x, 2:y, 2:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Neighbors along dimension 3 (across-plane diagonal xyz2)
[X1, X2, X3] = ndgrid(2:x, 1:y-1, 1:z-1);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(1:x-1, 2:y, 2:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Neighbors along dimension 3 (across-plane diagonal xyz3)
[X1, X2, X3] = ndgrid(1:x-1, 2:y, 1:z-1);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(2:x, 1:y-1, 2:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];

% Neighbors along dimension 3 (across-plane diagonal xyz4)
[X1, X2, X3] = ndgrid(2:x, 2:y, 1:z-1);
S = [S; sub2ind(size(M), X1(:), X2(:), X3(:))];

[X1, X2, X3] = ndgrid(1:x-1, 1:y-1, 2:z);
T = [T; sub2ind(size(M), X1(:), X2(:), X3(:))];


%% 3D plot it
% Calculate the weight for each edge 
W = mean(M_wts([S, T]), 2);

g = graph(S, T, W);
gplot = plot(g,'EdgeLabel',g.Edges.Weight);
layout(gplot,'force3')
view(3)