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

Question

Following this and that question and the corresponding answers, Mathematica can build 2D graph from 2D matrices.

Question

How to make the corresponding 3D graphs from 3D cubes?

Attempt

Modifying ever so slightly this great 2D answer from @kguler

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

Then I get

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

producing the following funky plot

which nicely shows how the three planes are not connected.

Answer 1

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]

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

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.

Answer 2

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]}]]

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

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]

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]

And if diagonal links are not required,

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

It works also on the skeletons

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

Answer 3

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]

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

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

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.

---

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.

Answer 4

Here's an alternative solution with RelationGraph.

Here's an array:

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

Image3D[arr]

Take the coordinates of the marked cells:

pts = Keys@Most@ArrayRules[arr];

Connect points whose distance is no larger than $\sqrt{3}$, then remove self-connection with SimpleGraph.

SimpleGraph[
 RelationGraph[Norm[#1 - #2] <= Sqrt[3] &, pts],
 VertexCoordinates -> pts
]

If you only want connections through a common edge or face (instead of a common corner point), use $\sqrt{2}$ or $1$ as the distance threshold.

This solution is quite slow (much slow than the other one I posted), but it is by far the simplest.