21
$\begingroup$

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

Mathematica graphics

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

Mathematica graphics

which nicely shows how the three planes are not connected.

$\endgroup$
1
  • $\begingroup$ chris, can you try with the following modification of arrayGraph: add before the last line edges2= UndirectedEdge @@@ DeleteDuplicates[Sort /@ Flatten[Thread /@ ComponentMeasurements[Transpose/@m, "Neighbors"]]];edges=DeleteDuplicates[Join[edges,edges2]]; (I am using the free online version of V10, and the input matrix is too large given the limitations of the free version.) $\endgroup$
    – kglr
    Oct 17 '14 at 20:31
10
$\begingroup$

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

$\endgroup$
2
  • $\begingroup$ I knew that you would come out with something funky, +1! $\endgroup$
    – Öskå
    Oct 18 '14 at 1:07
  • $\begingroup$ thank you @Oska. Got lucky here, V10 happened to add to ComponentMeasurements just what we needed for this question:) $\endgroup$
    – kglr
    Oct 18 '14 at 1:20
10
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Nice answer too; sorry I could not accept both! $\endgroup$
    – chris
    Oct 18 '14 at 8:34
4
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

Here's an alternative solution with RelationGraph.

Here's an array:

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

Image3D[arr]

enter image description here

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
]

enter image description here

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.