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

Question

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

Answer 1

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 = Internal`PartitionRagged[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 `"]
IGMeshCellAdjacencyGraph[meshregion, 2]

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[
      Graphics`PolygonUtils`PolygonIntersection[{R, #}]] &,
    MeshPrimitives[myVoronoi, 2]
    ]
   ];

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

Answer 2

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

polygons = Graphics`PolygonUtils`PolygonCombine @ 
    Graphics`PolygonUtils`PolygonIntersection[#, 
       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];
AdjacencyGraph[adjm, 
 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]