# How to find the MeshConnectivityGraph[] from a list of independent cell mesh regions?

I'm generating a bunch of points inside a BoundaryMeshRegion. Then I generate the Voronoi mesh of the points, and take the intersection of the Voronoi cells into the region. How could I recover the adjacency of the resulting cells? The problem is that at the end I end up with a list of independent mesh cell regions, so I don't know how could I get the adjacency of the cells from the list.

What I'm doing below.

Take the following binary image:

myImage =


Get the mesh region:

myMeshRegion = ImageMesh[myImage]


Generate a bunch of points inside such region:

myPoints = RandomPoint[myMeshRegion, 100];
Show[{myMeshRegion, Graphics[Point[myPoints]]}]


Get the Voronoi mesh from the points:

myVoronoi = VoronoiMesh[myPoints]


Get the intersection of myMeshRegion and the cells in myVoronoi, using the method in this question.

myIntersectionCells = DeleteCases[
RegionIntersection[DiscretizeGraphics@#, myMeshRegion] & /@
MeshPrimitives[myVoronoi, 2], _RegionIntersection];
Show[myIntersectionCells]


By inspection, we can see that myIntersectionCells is a list of BoundaryMeshRegion elements, representing every "cell" in this mesh. So, from myIntersectionCells, can I get the MeshConnectivityGraph[]? Related to this question, how could I make myIntersectionCells a single mesh, instead of a list of mesh cells?

Thanks

P.S. I'm using MMA 12.2

This is surprisingly cumbersome. Even if one first converts all the Boundary MeshRegion s, then RegionUnion (which was my first guess) creates just a MeshRegion with only the boundary of the geometry stored as a single polygon.

The following should work however. First we extract all the polygons in myIntersectionCells. Then we concatenate all the point lists (in the given order!) and generate a new list allpolygons of the polygon index lists, and hand these over to MeshRegion. The latter is (hopefully) clever enough to deal with duplicate vertices.

allprimitives = (Join @@ (MeshPrimitives[#, 2] & /@ myIntersectionCells))[[All, 1]];
allpts = Join @@ allprimitives;
lengths = Length /@ allprimitives;
allpolygons = InternalPartitionRagged[Range[Total[lengths]], lengths];

meshregion = MeshRegion[allpts, Polygon[allpolygons]]


Afterwards, you should be able to apply

MeshConnectivityGraph[ meshregion, 2]


I have only version 12.0, so I have to simulate MeshConnectivityGraph with with the help of the IGraphM package:

Needs["IGraphM "]


I seem to obtain always connected graphs. So maybe MeshConnectivityGraph is a bit buggy. It would not be the first region-related function whose first version is buggy...

Update

One can generate myIntersectionCells a bit faster (5 times as fast on my machine) by using the following piece of code:

myIntersectionCells =
With[{R = MeshPrimitives[myMeshRegion, 2][[1]]},
Map[
BoundaryMeshRegion[
GraphicsPolygonUtilsPolygonIntersection[{R, #}]] &,
MeshPrimitives[myVoronoi, 2]
]
];


It is still orders of magnitude slower than it could be.

• @kglr Ah, using the cloud is a great idea. Thank you for the edit! – Henrik Schumacher Mar 10 at 12:15
• Very cool, thanks! I've tried it a few times, and it seems that in many cases it generates isolated components. Is this because this particular shape could have cells that span out-and-into the shape? I've tried increasing the point number, which reduces this effect, but ideally it would be cool having a single connected mesh. Even if that would mean adding "additional" points to the seeded. – TumbiSapichu Mar 10 at 13:08
• Ah, I see. I had made an error in the definition of allprimitives because I had falsely assumed that each mesh in myIntersectionCells would consists of only a single polygon. I corrected it. Please have a look whether it works better for you now. – Henrik Schumacher Mar 10 at 13:39
• Maybe MeshConnectivityGraph has a bug. See my edit; I showed a workaround there. – Henrik Schumacher Mar 10 at 19:27
• @TumbiSapichu Thank you for the bounty! – Henrik Schumacher Mar 19 at 7:32

An alternative approach to get intersection cells and face connectivity graph:

polygons = GraphicsPolygonUtilsPolygonCombine @
GraphicsPolygonUtilsPolygonIntersection[#,
MeshPrimitives[myMeshRegion, 2]] & /@ MeshPrimitives[myVoronoi, 2];

Show[myMeshRegion,
Graphics[{EdgeForm[White], RandomColor[], #} & /@ polygons],
ImageSize -> 700]


adjm = (1 - IdentityMatrix[Length@polygons])
Outer[Boole @ IntersectingQ[Join @@ #[[1]],Join @@ #2[[1]]] &, polygons, polygons];

VertexCoordinates -> (RegionCentroid /@ polygons),
EdgeStyle -> Directive[Thick, Black],
Prolog -> ({Opacity[.7], EdgeForm[White], RandomColor[], #} & /@ polygons),
ImageSize -> 700]


Update: In comparison, using meshregion from Henrik's answer with MeshConnectivityGraph we miss some connections:

MeshConnectivityGraph[meshregion, 2,
VertexCoordinates -> (RegionCentroid /@ MeshPrimitives[ meshregion, 2]),
EdgeStyle -> Directive[Thick, Black],
Prolog -> ({Opacity[.7], EdgeForm[White], RandomColor[], #} & /@
MeshPrimitives[ meshregion, 2]), ImageSize -> 700]


• What does Outer[] do? This method is not working well for me, the last part, when creating the adjacency matrix, usually crashes the kernel (I'm using 1000 points). The couple of times it has not crashed, it generates a massive amount of edges between the nodes (half a million from the 1000 nodes, when the other method -as I would expect- generates about 3K edges total). – TumbiSapichu Mar 10 at 19:07
• I'm using MMA 12.2 – TumbiSapichu Mar 10 at 19:37
• @TumbiSapichu, Outer checks if two polygons (possibly with multiple cells) share a coordinate for every pair polygons to create a binary adjacency matrix. I tried this method mainly to see if we can capture the connections missed by meshregion+MeshConnectivity approach in the original example. Outer is essentially a brute force method to compute the face connectivity matrix; and as is it does twice as many comparisons as necessary. Definitely not recommended for large inputs. – kglr Mar 10 at 19:38
• Oh, I see. Maybe that's why the kernel is crashing. But it's strange that when it does not crash, there are many more extra edges. – TumbiSapichu Mar 10 at 19:43
• Also, I always get an error for GraphicsPolygonUtilsPolygonCombine @ GraphicsPolygonUtilsPolygonIntersection, how were you able to run this? I have MMA 12.2, maybe a compatibility issue? – TumbiSapichu Mar 24 at 23:45