在 MATLAB 中构造 3D 点阵图
Construct 3D lattice graph in MATLAB
我对这个 question/answer () 的扩展很感兴趣，以将 4-connected 案例扩展到三维。
问题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);
答案 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 扩展到三维 ()：
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)
此方法使用 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，因此我无法测试哪个生成图表的速度更快。
我对这个 question/answer () 的扩展很感兴趣，以将 4-connected 案例扩展到三维。
问题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);
答案 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 扩展到三维 ()：
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)
此方法使用 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，因此我无法测试哪个生成图表的速度更快。