How to build a 3D graph from a 3D matrix?

Here is how it works. If you have a volume in 3d it is essential, that you use connected component labeling in 3d so that components that are connected over layers stick together and get the same label. Lucky for us that MorphologicalComponents can do this. Let's create a test volume

data = With[{init = RandomChoice[{0, 0, 1}, {10, 10}]},
   NestList[
    If[RandomChoice[{True, True, False}], 
      RotateLeft /@ #, #*Transpose[#]] &, init, 9]
   ];
cmp = MorphologicalComponents[data, 
  CornerNeighbors -> False]; 
Colorize[cmp]

Mathematica graphics

The colors indicate that objects are correctly recognized throughout the layers. Note that I used the $N_6$-foreground neighborhood here because this is a bit more intuitive when watching a 3d volume.

What you need to do now is to extract all voxel positions a label that are directly connected within the component. This can be done by simply checking Norm[Subtract[p1,p2]]===1. With Outer I compare all position combinations and create edges for directly neighboring positions. Note that if I would need to implement this fast, this is not the way to go.

With these edges I create a temporary Graph3D for one component. I use this to extract the correctly sorted VertexList which I need later to set the VertexCoordinates correctly. This all is done in the following function

builtGraph[components_, label_] :=
 Module[{g},
  g = Graph3D[
    DeleteDuplicates@
     Flatten@Last@Reap@Block[{labelOnly = Position[components, label]},
         Outer[
          If[Norm[Subtract[##]] === 1, 
            Sow[UndirectedEdge @@ Sort[{##}]]] &, labelOnly, 
          labelOnly, 1]
         ]];
  {EdgeList[g], VertexList[g]}
  ]

After that, I need to apply this function to each component label in the volume and built a final graph from it. Combining the single edge- and vertex-lists is a bit playing with Transpose and Flatten, but nothing really hard. The rest is simple

Graph3D[#1, VertexCoordinates -> #2] & @@ (Flatten[#, 1] & /@ 
   Transpose@Table[builtGraph[cmp, i], {i, Max[Flatten[cmp]]}])

Mathematica graphics

If you look close, you see the stair-structure on the right which is on the left in the image above and the green upside-down "T" of the image is in the upper left corner of the graph. This comes from Position which reverses the position in the matrix, but I'm sure you can fix this yourself.

One final note, since I was extracting edges rather than vertices and I haven't allowed self-referencing edges, you won't find components in the graph which consist only of one voxel.


Using ComponentMeasurements twice, on the original matrix m and on Transpose/@m we can get all Neighbors:

mat = RandomInteger[{0, 1}, {3, 3, 3}];
m = Module[{i = 1}, mat /. 1 :> i++];

v = ComponentMeasurements[m, "Label"][[All, 1]]

{1,2,3,4,5,6,7,8,9,10,11,12,13,14}

vcoords = ComponentMeasurements[m, "Centroid"][[All, -1]]

{{0.5,2.5,2.5},{1.5,2.5,2.5},{0.5,1.5,2.5},{1.5,1.5,2.5},{0.5,0.5,2. 5},
{2.5,0.5,2.5},{2.5,2.5,1.5},{0.5,1.5,1.5},{0.5,0.5,1.5},{1.5,0.5,1. 5},
{1.5,2.5,0.5},{2.5,2.5,0.5},{2.5,1.5,0.5},{1.5,0.5,0.5}}

linksa = List @@@ DeleteDuplicates[Sort /@ Flatten[Thread /@ 
   ComponentMeasurements[m, "Neighbors", CornerNeighbors -> False]]]

{{1,3},{2,4},{3,5},{3,8},{5,9},{7,12},{8,9},{10,14},{12,13}}

linksb = List @@@ DeleteDuplicates[Sort /@ Flatten[Thread /@ 
  ComponentMeasurements[Transpose /@ m, "Neighbors", CornerNeighbors -> False]]]

{{1,2},{3,4},{3,8},{5,9},{7,12},{9,10},{10,14},{11,12}}

alllinks = DeleteDuplicates[Join[linksa, linksb]]

{{1,3},{2,4},{3,5},{3,8},{5,9},{7,12},{8,9},{10,14},{12,13},{1,2},{3,4},{9,10},{11,12}}

Graphics3D[GraphicsComplex[vcoords, {PointSize[Large], Red, 
   Sphere[#, .1] & /@ v, Thick, Blue, Line[alllinks]}]]

enter image description here

Update: a larger input matrix

mat = RandomInteger[{0, 1}, {10, 10, 10}];

(* ... same calculations for v, vcoords and alllinks as above *)

Graphics3D[GraphicsComplex[vcoords, {PointSize[Large],
   Red, Sphere[#, .15] & /@ v, Opacity[.7], Orange, Tube[alllinks]}],
 ImageSize -> 500, Background -> Black, Boxed -> False]

we get

enter image description here

Graphics3D[GraphicsComplex[vcoords,
 {PointSize[Large], Red, Sphere[#, .15] & /@ v, 
  Opacity[.1], White, Cuboid[vcoords[[#]] - .5, vcoords[[#]] + .5] & /@ v,
  Opacity[.9], Orange, Tube[alllinks]}],
 ImageSize -> 500, Background -> Black, Boxed -> False]

enter image description here

Note: Had to use Graphics3D to produce the output instead of Graph3D because somehow

Graph3D[UndirectedEdge@@@alllinks, VertexCoordinates->Thread[v->vcoords]]

gives Null on the free version of Wolfram Programming Cloud.

Update

This can be encapsulated as follows

ClearAll[arrayGraph];
arrayGraph[mat_, opts : OptionsPattern[]] := 
Module[{m = Module[{i = 1}, mat /. 1 :> i++], edges, edges2, vcs, v},
v = ComponentMeasurements[m, "Label"][[All, 1]];
vcs = ComponentMeasurements[m, "Centroid"];
edges = 
UndirectedEdge @@@ 
DeleteDuplicates[
 Sort /@ Flatten[Thread /@ ComponentMeasurements[m, "Neighbors",
     FilterRules[Flatten[{opts}], 
      Options[ComponentMeasurements]]]]];
edges2 = 
UndirectedEdge @@@ 
DeleteDuplicates[
 Sort /@ Flatten[
   Thread /@ ComponentMeasurements[Transpose /@ m, "Neighbors",
     FilterRules[Flatten[{opts}], 
      Options[ComponentMeasurements]]]]]; 
edges = DeleteDuplicates[Join[edges, edges2]];
Graph3D[v, edges, VertexCoordinates -> vcs, 
FilterRules[Flatten[{opts}], Options[Graph3D]]]]

which works as follows:

 mat = RandomInteger[{0, 1}, {3, 3, 3}];
 arrayGraph[mat, VertexSize -> .3, EdgeStyle -> Directive[Thick, Red], 
 Boxed -> True]

Mathematica graphics

And if diagonal links are not required,

 arrayGraph[mat, VertexSize -> .3, EdgeStyle -> Directive[Thick, Red], 
 Boxed -> True, CornerNeighbors-> False]

Mathematica graphics

It works also on the skeletons

 (pl = {(skl = skel[dat]) // Image3D[#, ImageSize -> 300] & // 
 Rasterize,Image[arrayGraph[skl // Normal, CornerNeighbors -> False], 
 ImageSize -> 300]}); ImageMultiply@@pl

Mathematica graphics


IGraph/M makes this really easy using its mesh/graph conversion functions along with ArrayMesh.

Here's an array:

SeedRandom[42]
arr = RandomInteger[1, {8, 8, 8}];

Make a mesh:

mesh = ArrayMesh[arr]

enter image description here

Construct the adjacency graph of 3D cells:

IGMeshCellAdjacencyGraph[ArrayMesh[arr], 3]

We can add vertex coordinates manually, if needed:

IGMeshCellAdjacencyGraph[mesh, 3, 
 VertexCoordinates -> PropertyValue[{mesh, 3}, MeshCellCentroid]]

enter image description here

This graph connects cells with a common adjacent face (von Neumann neighbourhood). If instead you want to connect cells with a common adjacent edge (1-dimensional cell), you can do this:

mat = IGMeshCellAdjacencyMatrix[mesh, 3, 1];

SimpleGraph[
 AdjacencyGraph@Unitize[mat.Transpose[mat]],
 VertexCoordinates -> PropertyValue[{mesh, 3}, MeshCellCentroid]
]

enter image description here

Or perhaps you want to connect cells with a common vertex (0-dimensional cell), i.e. also include the diagonals of the cubes as connections.

enter image description here


You can also use IGBipartiteProjections to get the same graph like this:

SetProperty[
 First@IGBipartiteProjections@IGMeshCellAdjacencyGraph[mesh, 3, 1],
 VertexCoordinates -> PropertyValue[{mesh, 3}, MeshCellCentroid]
]

This last solution is a bit sloppy because it determines the bipartite partitions automatically, and there is more than one valid partitioning. A robust way is to specify partitions manually:

SetProperty[
 First@IGBipartiteProjections[IGMeshCellAdjacencyGraph[mesh, 3, 1], {MeshCellIndex[mesh, 3], MeshCellIndex[mesh, 1]}],
 VertexCoordinates -> PropertyValue[{mesh, 3}, MeshCellCentroid]
]

Some background:

IGMeshCellAdjacencyGraph[mesh, 3, 1] gives the bipartite graph describing the connectivity of dimension-3 and dimension-1 mesh cells. A bipartite graph is a graph with two groups of vertices, $A$ and $B$, where connections may only run between the two groups (but not within a single group). A bipartite projection is a graph whose vertex set is $A$, and two of its vertices $a_1, a_2 \in A$ are connected iff there is some vertex $b \in B$ so that the connections $a_1 \leftrightarrow b$ and $a_2 \leftrightarrow b$ both exist.