How to obtain the cell-adjacency graph of a mesh?
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 MatrixPower
s of the adjacency matrix is considerably faster with (real) numeric matrices. A remedy could be the function As Chip Hurst has pointed out, we can turn a pattern array into a numerical one with SparseArrayFromPatternArray
below.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}
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"
]
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
.
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}