我对扩展这个问题/答案(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);
答案1:下面的答案3中的循环效率更高一些。与其建立巨大的边缘和权重矩阵,不如直接在循环中迭代地建立图形。这可以很快完成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 = digraph();
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 % loop through all "horizontal" (column) connections in each z plane
% example:
% 1 - 3
% 2 - 4
for z = 1:zdim
if y ~= 1
nodes1 = imat(:,y,z); % grab all rows (nodes) of column y
nodes2 = imat(:,y-1,z); % grab all rows (nodes) of column y-1
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
if y ~= ydim
nodes1 = imat(:,y,z);
nodes2 = imat(:,y+1,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
end
for x = 1:xdim % loop through all "vertical" (row) connections within each z plane
% example:
% 1 3
% | |
% 2 4
for z = 1:zdim
if x ~= 1
nodes1 = imat(x,:,z); % grab all rows (nodes) of column y
nodes2 = imat(x-1,:,z); % grab all rows (nodes) of column y-1
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
if x ~= ydim
nodes1 = imat(x,:,z); % grab all rows (nodes) of column y
nodes2 = imat(x+1,:,z); % grab all rows (nodes) of column y-1
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
end
for z = 1:zdim % loop through all "deep" connections in each z plane
% example:
% 5 7
% / /
% 1 3
for x = 1:xdim
if z ~= 1
nodes1 = imat(x,:,z); % grab all rows (nodes) of column y
nodes2 = imat(x,:,z-1); % grab all rows (nodes) of column y-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
if z ~= zdim
nodes1 = imat(x,:,z); % grab all rows (nodes) of column y
nodes2 = imat(x,:,z+1); % grab all rows (nodes) of column y-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
end
disp('done.')
figure
spy(adjacency(g))
答案2:将链接的答案的4个连接的情况扩展到第三维(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:一个(非常慢的低效率)for循环版本,它使用权重计算创建一个边列表:
% input matrix of values
X=3;Y=3;Z=2;
A = reshape([1:(X*Y*Z)],X,Y,Z);
% # nodes is equal to numel(A)
% # edges is equal to number of pairs of 4-connected nodes (cityblock)
nodes1 = [];
nodes2 = [];
wts = [];
xdim = size(A,1);
ydim = size(A,2);
zdim = size(A,3);
% create adjacency matrix
k=1;
disp('finding connected nodes for edge list with weights...')
for z1 = 1:zdim
for x1 = 1:xdim
for y1 = 1:ydim
for z2 = 1:zdim
for x2 = 1:xdim
for y2 = 1:ydim
node1 = xdim*(x1-1) + y1 + (xdim*ydim)*(z1-1);
node2 = xdim*(x2-1) + y2 + (xdim*ydim)*(z2-1);
d = sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2);
if d==1
% add first node index:
nodes1(k) = node1;
% add second nodename:
nodes2(k) = node2;
% add edge weight between these nodes:
wts(k) = nanmean([A(x1,y1,z1),A(x2,y2,z2)]);
k=k+1;
end
end
end
end
end
end
end
g = digraph(nodes1, nodes2, wts);
figure
spy(adjacency(g))
disp('done.')
此人实际上可能适用于大箱子,但我没有耐心寻找。