How to obtain the cell-adjacency graph of a mesh?

Question

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

I need this Graph:

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?

Answer 1

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]

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]

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}

Answer 2

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

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.

Answer 3

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]

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}

Answer 4

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

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

Answer 5

In version 12.1, the function MeshConnectivityGraph[] is now built-in. Using the examples in Henrik's answer:

Table[Show[MeshConnectivityGraph[mesh, k, VertexLabels -> "Index"]],
      {k, 0, 2}] // GraphicsRow

Answer 6

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