Construct 3D lattice graph in MATLAB

Question

I was interested in an extension of this question/answer (https://stackoverflow.com/a/3283732/2371031) to extend the 4-connected case into a third dimension.

Question 1: given an X x Y x Z dimensional matrix, how does one construct the 6-connection adjacency matrix (or edge list)?

Question 2: given the adjacency matrix (or edge list), how does one compute the edge weights as some function (e.g., mean) of the values of connected nodes?

Question 3: How does one solve Question 1 and 2 for very large matrices without running out of memory? Sparse matrices or edge lists seem like probable solutions, but how to implement?

% 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);

Do you want the 4-connected neighbors of the 3d matrix, or the 6-connected neighbors so that the planes are connected as well? — beaker, Jun 14 '20 at 15:07
What will you use this graph for? An explicit representation is usually not necessary, for example in image processing all graph methods can be implemented with implicit graphs much more efficiently. If you do need an explicit graph representation, it is often worth while to use the edge list representation rather than an adjacency matrix. Adjacency matrices are only efficient with densely connected graphs. For neighbors in a matrix (image?) the graph is always very sparse, making the adjacency matrix inefficient both computationally and memory-wise. — Cris Luengo, Jun 14 '20 at 16:49
Brian, If you're writing your own, all you need is a `neighbors()` function which, as @CrisLuengo said, will be much more memory (and time) efficient. — beaker, Jun 15 '20 at 15:21
In my Answer 1 below, I build the graph using an edge list. But it takes a really long time to loop through all the edges in the large case. — Brian D, Jun 15 '20 at 15:32
Why are you using `digraph` instead of `graph`? Are your edges actually directed? Adding edges to the graph one-by-one is going to be very slow. It would be much better to create the edge list/adjacency matrix and pass that to the graph constructor. I'll clean up some code to generate the edge list and post it. — beaker, Jun 15 '20 at 15:35
ah good point. I was thinking that i'd use digraph to allow the shortest path to go in either direction... but i suppose it doesn't matter as long as the edge weight is the same in either direction. In the case where I want a solution to always travel one direction in one of the dimensions then I could use a digraph to set the weights to Inf in the "wrong" direction, or just remove those connections... probably simpler to just use an undirected graph though. — Brian D, Jun 15 '20 at 15:42
If you want to compute the shortest path I think test best option is using [graydist](https://www.mathworks.com/help/images/ref/graydist.html#bvih8ch). As explained in the examples you need to perform two distance transforms one from the source and the other from the destination and sum the two transforms. Now you can use a simple loop to find the path. — rahnema1, Jun 17 '20 at 01:04
@rahnema1 Sure, Steve Eddins did a [series of articles](https://blogs.mathworks.com/steve/2011/11/01/exploring-shortest-paths-part-1/) on this approach a while back. But by the time you've done the distance transforms and found a single path, I think you're probably just better off implementing Dijkstra using a custom function that reads the neighbors of a node directly from the grid rather than constructing an adjacency matrix or list. — beaker, Jun 17 '20 at 15:09
@beaker Thanks for the link! Probably Steve Eddins's last article from that series would be using `graydist` to compute _weighted shortest path_ on an image. "Reading neighbors of a node directly from the grid" is good idea. However it is just re-implementation of graydist but in MATLAB language that lacks efficient data structures like priority queue. I think graydist is simplest approach. We should wait for @ BrianD to provide detailed answer comparing different methods. — rahnema1, Jun 17 '20 at 16:24
@beaker I wanted to convert the matrix to a graph so i would have access to the graph `shortestpath` method which is extremely fast even on a very large graph (less than a second). However, it takes about a minute to (1) build the edge list and (2) create the graph. Are you suggesting that it would be faster to write a shortest path algorithm that takes in (a) a matrix (instead of a graph), (b) start and end node, and (c) connectivity, and run the algorithm directly on the matrix without bothering to convert to a graph? I have not tried the `greydist` method. — Brian D, Jun 18 '20 at 17:42
@BrianD ***If*** you have a fast implementation of, say, Dijkstra's algorithm and can replace the function that returns the neighbors of a node, finding the neighbors directly from the grid is no more costly than finding the neighbors from an adjacency list (i.e., *O(1)*, in the case of a constant number of neighbors). The advantage is that you do not have to spend the time converting the grid to some other graph representation beforehand. — beaker, Jun 18 '20 at 18:51

score 1 · Accepted Answer · answered Jun 15 '20 at 16:28

This approach uses ndgrid to find each node's neighbor in one direction. For example, in dimension 2 we find the (valid) neighbor to the right of each node.

Here, the blue nodes are the source of the edge, and the red nodes to their right are the target nodes. So the grids for the source and target nodes are offset by 1 from each other.

% 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);

Running this in Octave with your large testcase:

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

takes about 45 seconds for the (undirected) edge list, and an additional 20 seconds or so to generate the sparse adjacency matrix. Since I use Octave, I can't test which one will be faster in generating the graph.

Nice! I like this approach. It definitely works for the large case. That part was ~35 seconds for me and creating the graph from S, T, W was ~40 seconds (I didn't try multiple runs for benchmarking, just a quick n=1 test). Quite promising as I expect I'll be able to greatly reduce the total matrix/edge list search space. — Brian D, Jun 17 '20 at 20:43

Brian D · Answer 2 · 2021-07-06T16:23:27.727

Answer 1: Rather than building up a huge matrix of edges and weights... iteratively build the graph directly in the loop using the edge list, looping over each dimension (rows, cols, planes). This one does 64x64x10 pretty quick (seconds), but for the large example (256x256x1000) runs MUCH longer. Would be nice to work out a faster version.

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))

Answer 2: extends the linked answer's 4-connected case into a third dimension (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)

For the large case Matlab complains with an error:

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

Update

Answer 3:
(updated version of @beaker's answer)

I was curious whether I could get an improved shortest path result by using the 26-connected case (all 9 linearly or diagonally adjacent points in the superior plane, 9 in the inferior plane, and 8 in the coincident plane).

% 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)

Construct 3D lattice graph in MATLAB

2 Answers2