18
$\begingroup$

In addition to the accepted answer, see also the answer by Chip Hurst. This functionality is built in, but not documented.


Given an arbitrary mesh region, how can I efficiently obtain the graph describing the adjacency structure of mesh cells?

For example, given the following mesh,

SeedRandom[123]
pts = RandomReal[1, {10, 2}];
mesh = VoronoiMesh[pts, MeshCellLabel -> {2 -> "Index"}]

enter image description here

I need this Graph:

enter image description here

It tells me that e.g. cells 4 and 5 are neighbours.

This example was for 2D cells, but the problem can be stated generally for cells of any dimension:

  • two edges (1D cells) are adjacent if they have a common point
  • two faces (2D cells) are adjacent if they have a common edge
  • two 3D cells are adjacent if they have a common face
  • ...

I am looking both for specialized methods to obtain the face-adjacency graph of a 2D or 3D mesh, and general methods to obtain the $k$-dimensional cell adjacency graph of a $d > k$ dimensional mesh.


Here's a naïve method for face-adjacency. It is too slow for practical use due to its $O(n^2)$ complexity.

SimpleGraph[
 RelationGraph[
  Length@Intersection[First@MeshCells[mesh, #1], First@MeshCells[mesh, #2]] >= 2 &,
  MeshCellIndex[mesh, 2]
  ],
 VertexLabels -> "Name"
 ]

This method exploits the fact that if two polygons (faces) are adjacent, then they will have (at least) two points in common. For example, Polygon[{14, 8, 7, 11}] and Polygon[{11, 7, 2, 4, 13}] have points 7 and 11 in common. It generalizes to higher dimensions as well: two 3D cells are adjacent if they have at least 3 points in common.

However, it is quite slow because RelationGraph will test each pair of cells for adjacency.

SeedRandom[123]
pts = RandomReal[1, {500, 2}];
mesh = VoronoiMesh[pts];

RelationGraph[
   Length@Intersection[First@MeshCells[mesh, #1], 
       First@MeshCells[mesh, #2]] >= 2 &,
   MeshCellIndex[mesh, 2]
   ]; // AbsoluteTiming

(* {2.36978, Null} *)

While this method can be sped up by a constant factor with minor tweaks, this won't fix the core problem: quadratic complexity.

cells = MeshCells[mesh, 2][[All, 1]];        
RelationGraph[Length@Intersection[#1, #2] >= 2 &, cells]; // AbsoluteTiming
(* {0.815857, Null} *)

1 second for only 500 cells is still too slow. Can we do significantly better?

$\endgroup$
  • $\begingroup$ You really need only the connectivity graph, not the vertex positions? I think I can do something about it... $\endgroup$ – Henrik Schumacher Nov 22 '17 at 10:00
  • $\begingroup$ If you have the initial points, wouldn't it be easier to use DelaunayMesh and MeshPrimitives? As in pi = PositionIndex[MeshPrimitives[DelaunayMesh[pts], 0][[;; , 1]]][[;; , 1]], Graph[Values[pi], MeshPrimitives[DelaunayMesh[pts], 1][[;; , 1]] /. pi, VertexLabels -> "Name"]. $\endgroup$ – aardvark2012 Nov 22 '17 at 10:02
  • $\begingroup$ @HenrikSchumacher I only need the connectivity graph. Vertex positions are irrelevant. $\endgroup$ – Szabolcs Nov 22 '17 at 10:08
  • $\begingroup$ @aardvark2012 The VoronoiMesh was just an example. I am looking for a method that works for any mesh. I used a Voronoi mesh as an example because, unlike most other meshes one might encounter in Mathematica, its cells are not triangles (some methods may work only for triangles). $\endgroup$ – Szabolcs Nov 22 '17 at 10:09
15
$\begingroup$

I think, I found a general and even faster way, but I haven't tested it for $1$- or $3$-dimensional MeshRegions.

The following function computes the cell-vertex-adjacency matrix A. Two $d$-dimensional cells ($d>0$) are adjacent if they share at least $d$ common points. We can find these pairs by looking for entries $\geq d$ in A.Transpose[A].

ToPack = Developer`ToPackedArray;
ClearAll[getCellCellAdjacencyList];
getCellCellAdjacencyList[R_MeshRegion, d_] := 
 Module[{pts, cells, A, lens, n, m, nn},
  pts = MeshCoordinates[R];
  cells = ToPack[MeshCells[R, d][[All, 1]]];
  lens = Length /@ cells;
  n = Length[pts];
  m = Length[cells];
  nn = Total[lens];
  A = SparseArray @@ {Automatic, {m, n}, 0, {1, {
       ToPack[Join[{0}, Accumulate[lens]]],
       ArrayReshape[Flatten[Sort /@ cells], {nn, 1}]
       },
      ConstantArray[1, nn]}};
  SparseArray[
    UnitStep[UpperTriangularize[A.Transpose[A], 1] -  d]
  ]["NonzeroPositions"]
  ]

A special treatment is necessary for 0-dimensional cells; it's just the edges that we need.

getCellCellAdjacencyList[R_MeshRegion, 0] := ToPack[MeshCells[R, 1][[All, 1]]]

Here are some examples:

SeedRandom[123]
pts = RandomReal[1, {10, 2}];
R = VoronoiMesh[pts];

GraphicsGrid[Table[
  {VoronoiMesh[pts, MeshCellLabel -> {d -> "Index"}],
   Graph[getCellCellAdjacencyList[R, d], VertexLabels -> "Name"]
   }, {d, 0, 2}], ImageSize -> Large]

enter image description here

And some timings for comparison:

SeedRandom[123]
pts = RandomReal[1, {10000, 2}];
R = VoronoiMesh[pts]; // RepeatedTiming
getCellCellAdjacencyList[R, 0]; // RepeatedTiming
getCellCellAdjacencyList[R, 1]; // RepeatedTiming
getCellCellAdjacencyList[R, 2]; // RepeatedTiming

{0.636, Null}

{0.015, Null}

{0.031, Null}

{0.041, Null}

Edit

It's now rather straight-forward to write methods for the various adjacency matrices, lists, and graphs, even for cells of different dimensions (see below).

Edit 2

As Chip Hurst pointed out, the adjacency matrix of a MeshRegion R for distinct dimensions d1, d2 can be found as pattern SparseArray under R["ConnectivityMatrix"[d1,d2]]. (Its "RowPointers" and "ColumnIndices" must have been computed immediately when the MeshRegion was built.)

Many applications of adjacency matrices, in particular in finite elements, need 1 instead of Pattern as nonzero entries. Even computing vertex rings in a graph by using MatrixPowers of the adjacency matrix is considerably faster with (real) numeric matrices. A remedy could be the function SparseArrayFromPatternArray below. As Chip Hurst has pointed out, we can turn a pattern array into a numerical one with Unitize. I updated my old code to utilize this observation, leading to a tremendous performance boost. Somewhat surprisingly, even the old implementation of CellAdjacencyMatrix[R, 1, 2] tends to be faster than R["ConnectivityMatrix"[1,2]], so that I decided to use the new approach only for the case when either d1 or d2 is equal to 0.

CellAdjacencyMatrix[R_MeshRegion, d_, 0] := If[MeshCellCount[R, d] > 0,
   Unitize[R["ConnectivityMatrix"[d, 0]]],
   {}
   ];

CellAdjacencyMatrix[R_MeshRegion, 0, d_] := If[MeshCellCount[R, d] > 0,
   Unitize[R["ConnectivityMatrix"[0, d]]],
   {}
   ];

CellAdjacencyMatrix[R_MeshRegion, 0, 0] := 
  If[MeshCellCount[R, 1] > 0,
   With[{A = CellAdjacencyMatrix[R, 0, 1]},
    With[{B = A.Transpose[A]},
     SparseArray[B - DiagonalMatrix[Diagonal[B]]]
     ]
    ],
   {}
   ];

CellAdjacencyMatrix[R_MeshRegion, d1_, d2_] := 
  If[(MeshCellCount[R, d1] > 0) && (MeshCellCount[R, d2] > 0), 
   With[{B = CellAdjacencyMatrix[R, d1, 0].CellAdjacencyMatrix[R, 0, d2]},
    SparseArray[
     If[d1 == d2,
      UnitStep[B - DiagonalMatrix[Diagonal[B]] - d1],
      UnitStep[B - (Min[d1, d2] + 1)]
      ]
     ]
    ],
   {}
   ];

CellAdjacencyLists[R_MeshRegion, d1_, d2_] := 
  If[(MeshCellCount[R, d1] > 0) && (MeshCellCount[R, d2] > 0),
   Module[{i1, i2, data},
    data = If[d1 == d2,
      UpperTriangularize[CellAdjacencyMatrix[R, d1, d2], 1]["NonzeroPositions"], 
      CellAdjacencyMatrix[R, d1, d2]["NonzeroPositions"]
      ];
    If[Length[data] > 0,
     {i1, i2} = Transpose[data];
     Transpose[
      {
       Transpose[{ConstantArray[d1, {Length[i1]}], i1}],
       Transpose[{ConstantArray[d2, {Length[i2]}], i2}]
       }
      ],
     {}
     ]
    ],
   {}
   ];

CellAdjacencyGraph[R_MeshRegion, d1_, d2_] := Graph[
   Join[MeshCellIndex[R, d1], MeshCellIndex[R, d2]],
   UndirectedEdge @@@ CellAdjacencyLists[R, d1, d2],
   VertexLabels -> "Name"
   ];

Note that CellAdjacencyLists and CellAdjacencyGraph use labels that are compatible with those obtained from MeshCellIndex. Applied to Szabolcs's example MeshRegion, theses graphs look as follows:

GraphicsGrid[
 Table[CellAdjacencyGraph[R, d1, d2], {d1, 0, 2}, {d2, 0, 2}], 
 ImageSize -> Full]

enter image description here

As for comparing the performance of these new implementations to getCellCellAdjacencyList:

{
 getCellCellAdjacencyList[R, 0]; // RepeatedTiming // First,
 getCellCellAdjacencyList[R, 1]; // RepeatedTiming // First,
 getCellCellAdjacencyList[R, 2]; // RepeatedTiming // First
 }
{
 CellAdjacencyLists[R, 0, 0]; // RepeatedTiming // First,
 CellAdjacencyLists[R, 1, 1]; // RepeatedTiming // First,
 CellAdjacencyLists[R, 2, 2]; // RepeatedTiming // First
 }

{0.015, 0.030, 0.037}

{0.0068, 0.011, 0.0066}

$\endgroup$
  • $\begingroup$ Thank you for this! I haven't forgotten, and I will accept an answer in due time. $\endgroup$ – Szabolcs Nov 28 '17 at 17:40
  • $\begingroup$ Would you mind if I incorporated some of this into IGraph/M (with credits of course)? github.com/szhorvat/IGraphM $\endgroup$ – Szabolcs Nov 28 '17 at 17:43
  • $\begingroup$ @Szabolcs No, I would not. To the contrary: I would feel honored! I would also be thankful for any improvements you might find as I use these things for my own work on surface meshes. $\endgroup$ – Henrik Schumacher Nov 28 '17 at 20:18
  • $\begingroup$ Just a note that while Transpose[Transpose[edges][[{2, 1}]] is faster than Reverse[edges, {2}], it fails on {}. I think that given that computing the edge list in that function takes two orders of magnitude longer for big meshes, the slowdown due to Reverse is acceptable. $\endgroup$ – Szabolcs Nov 30 '17 at 12:13
  • $\begingroup$ @Szabolcs Thank you! Good point. I will add a check against empty edge lists. You are right that MeshCells[R, 1] takes pretty long. In an ideal word, the edge list would be cached within the MeshRegion for it is used for quite a number of things. Actually, that's what I do in my own data type for simplicial surfaces. If you subtract the time for getting the edges, then Reverse slows down building the sparse matrix by a factor of two. This is why I prefer a check against emptiness of the edge list. As long as MeshRegions don't cache the edge list, this is a matter of taste, though. $\endgroup$ – Henrik Schumacher Nov 30 '17 at 12:40
12
$\begingroup$

I need three compiled helper functions:

getEdgesFromPolygons = Compile[{{f, _Integer, 1}},
   Table[
    {
     Min[Compile`GetElement[f, i], Compile`GetElement[f, Mod[i + 1, Length[f], 1]]], 
     Max[Compile`GetElement[f, i], Compile`GetElement[f, Mod[i + 1, Length[f], 1]]]
     },
    {i, 1, Length[f]}
    ],
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   CompilationTarget -> "C"
   ];
takeSortedThread = Compile[{{data, _Integer, 1}, {ran, _Integer, 1}},
   Sort[Part[data, ran[[1]] ;; ran[[2]]]],
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   CompilationTarget -> "C"
   ];
extractIntegerFromSparseMatrix = Compile[
   {{vals, _Integer, 1}, {rp, _Integer, 1}, {ci, _Integer, 
     1}, {background, _Integer},
    {i, _Integer}, {j, _Integer}},
   Block[{k},
    k = rp[[i]] + 1;
    While[k < rp[[i + 1]] + 1 && ci[[k]] != j, ++k];
    If[k == rp[[i + 1]] + 1, background, vals[[k]]]
    ],
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   CompilationTarget -> "C"
   ];

The following functions takes a MeshRegion and finds all pairs of neighboring two-dimensional MeshCells. First, it generates a list of all edges (with sorted indices) and creates a lookup table for edges in form of a SparseArray. With the lookup table, we can find the indices of all edges bounding a given polygon, so that we can build the SparseArray edgepolygonadjacencymatrix, whose "AdjacencyLists" is what we are looking for. The method should have linear complexity.

ToPack = Developer`ToPackedArray;
getPolygonPolygonAdjacencyList[R_MeshRegion] := 
 Module[{pts, polygons, edgesfrompolygons, edges, edgelookupcontainer,
    polyranges, polygonsneighedges, edgepolygonadjacencymatrix, acc},
  pts = MeshCoordinates[R];
  polygons = ToPack[MeshCells[R, 2][[All, 1]]];
  edgesfrompolygons =  ToPack[Flatten[getEdgesFromPolygons[polygons], 1]];
  edges = DeleteDuplicates[edgesfrompolygons];
  edgelookupcontainer = 
   SparseArray[
    Rule[Join[edges, Transpose[Transpose[edges][[{2, 1}]]]], 
     Join[Range[1, Length[edges]], Range[1, Length[edges]]]], {Length[
      pts], Length[pts]}];
  acc = Join[{0}, Accumulate[ToPack[Length /@ polygons]]];
  polyranges = Transpose[{Most[acc] + 1, Rest[acc]}];
  polygonsneighedges = takeSortedThread[extractIntegerFromSparseMatrix[
      edgelookupcontainer["NonzeroValues"],
      edgelookupcontainer["RowPointers"],
      Flatten@edgelookupcontainer["ColumnIndices"],
      edgelookupcontainer["Background"],
      edgesfrompolygons[[All, 1]],
      edgesfrompolygons[[All, 2]]],
     polyranges];
  edgepolygonadjacencymatrix = Transpose@With[{
      n = Length[edges], m = Length[polygons],
      data = ToPack[Flatten[polygonsneighedges]]
      },
     SparseArray @@ {Automatic, {m, n}, 
       0, {1, {acc, Transpose[{data}]}, ConstantArray[1, Length[data]]}}
     ];
  Select[(edgepolygonadjacencymatrix["AdjacencyLists"]), Length[#] == 2 &]
  ]

Testing with OP's example:

SeedRandom[123]
pts = RandomReal[1, {10, 2}];
R = VoronoiMesh[pts, MeshCellLabel -> {2 -> "Index"}]
Graph[
 UndirectedEdge @@@ getPolygonPolygonAdjacencyList[R],
 VertexLabels -> "Name"
 ]

enter image description here

enter image description here

Speed test

SeedRandom[123]
pts = RandomReal[1, {10000, 2}];
R = VoronoiMesh[pts, 
    MeshCellLabel -> {2 -> "Index"}]; // RepeatedTiming
getPolygonPolygonAdjacencyList[R]; // RepeatedTiming

{0.625, Null}

{0.086, Null}

Edit

Slight improvement by merging Sort into takeThread (takeThread replaced by takeSortedThread).

Slight improvement by replacing Extract with extractIntegerFromSparseMatrix.

$\endgroup$
10
$\begingroup$

Using a few undocumented properties of MeshRegion[] objects, we have the following:

BlockRandom[SeedRandom[123];
            vm = VoronoiMesh[RandomReal[1, {10, 2}]]];

Show[vm,
     Graph[Range[vm["FaceCount"]], 
           Union[Sort /@ Flatten[MapIndexed[Thread[UndirectedEdge[#2[[1]], #1]] &, 
                                            vm["FaceFaceConnectivity"]]]],
           PlotTheme -> "ClassicDiagram", 
           VertexCoordinates -> Map[Mean, vm["FaceCoordinates"]]]]

mesh and its face-face connectivity graph

I'm not sure why the labels from this version are not consistent with the labels from MeshCellLabel, tho.

$\endgroup$
  • $\begingroup$ vm["EdgeEdgeConnectivity"] and vm["VertexVertexConnectivity"] might also be of interest. $\endgroup$ – J. M. is away Nov 22 '17 at 14:12
  • $\begingroup$ Ah, "FaceFaceConnectivity" works for MeshRegions in the plane! I tried these properties several times for two-dimensional MeshRegions in $\mathbb{R}^3$ and it did not work. Seemingly, there is still a lot to be done internally... $\endgroup$ – Henrik Schumacher Nov 22 '17 at 14:19
8
$\begingroup$

Here's another way.

Data from OP:

SeedRandom[123]
pts = RandomReal[1, {10, 2}];
mesh = VoronoiMesh[pts];

Get the adjacency matrix:

conn = mesh["ConnectivityMatrix"[2, 1]];
adj = conn.Transpose[conn];

Find the cell centroids for visualization purposes:

centers = PropertyValue[{mesh, 2}, MeshCellCentroid];

g = AdjacencyGraph[adj, PlotTheme -> "Scientific", VertexCoordinates -> centers];

Show[mesh, g]

enter image description here

Using the same profiling code as Henrik, we have

SeedRandom[123]
pts = RandomReal[1, {10000, 2}];
R = VoronoiMesh[pts]; // RepeatedTiming

getCellCellAdjacencyList[R, 2]; // RepeatedTiming

RepeatedTiming[
  conn = R["ConnectivityMatrix"[2, 1]];
  conn . Transpose[conn];
]

{0.632, Null}

{0.042, Null}

{0.012, Null}

$\endgroup$
  • 2
    $\begingroup$ So this was already built in ... it's really a pity that all this stuff is undocumented!! $\endgroup$ – Szabolcs Feb 23 '18 at 19:23
  • $\begingroup$ I imagine it's partly because this is not a good design for documented functionality. "ConnectivityMatrix"[2, 1] is quite strange looking. $\endgroup$ – Chip Hurst Feb 23 '18 at 19:26
  • 1
    $\begingroup$ Just write it as mesh@"ConnectivityMatrix"[2, 1] and it looks fine. It's exactly the same notation that J/Link uses for method calls. The @ replaces the . that most object-oriented languages would use. $\endgroup$ – Szabolcs Feb 23 '18 at 19:29
  • 1
    $\begingroup$ Oh my, that's so embarassing. Why aren't these things documented? @Szabolcs Of course I don't mind if you make this the accepted answer. $\endgroup$ – Henrik Schumacher Feb 24 '18 at 12:47
  • 1
    $\begingroup$ @ChipHurst Do you know of an efficient way to turn R["ConnectivityMatrix"[2, 1]] into a numerical array with 1 instead of Pattern as nonzero values? That would be very useful for many applications... $\endgroup$ – Henrik Schumacher Mar 8 '18 at 11:53
1
$\begingroup$

Note VoronoiMesh is the dual of the DelaunayMesh, so

Show[mesh, 
 AdjacencyGraph[DelaunayMesh[pts]["AdjacencyMatrix"], 
  VertexCoordinates -> MeshCoordinates[DelaunayMesh[pts]], 
  EdgeStyle -> Red]]

But I don't know how to get the right vertex label..

$\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.