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 =

enter image description here

Get the mesh region:

myMeshRegion = ImageMesh[myImage]

enter image description here

Generate a bunch of points inside such region:

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

enter image description here

Get the Voronoi mesh from the points:

myVoronoi = VoronoiMesh[myPoints]

enter image description here

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

enter image description here

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?


P.S. I'm using MMA 12.2


2 Answers 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 = Internal`PartitionRagged[Range[Total[lengths]], lengths];

meshregion = MeshRegion[allpts, Polygon[allpolygons]]

enter image description here

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]

enter image description here

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


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

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

  • $\begingroup$ @kglr Ah, using the cloud is a great idea. Thank you for the edit! $\endgroup$ Commented Mar 10, 2021 at 12:15
  • $\begingroup$ 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. $\endgroup$ Commented Mar 10, 2021 at 13:08
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Mar 10, 2021 at 13:39
  • 1
    $\begingroup$ Maybe MeshConnectivityGraph has a bug. See my edit; I showed a workaround there. $\endgroup$ Commented Mar 10, 2021 at 19:27
  • 1
    $\begingroup$ @TumbiSapichu Thank you for the bounty! $\endgroup$ Commented Mar 19, 2021 at 7:32

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

polygons = Graphics`PolygonUtils`PolygonCombine @ 
       MeshPrimitives[myMeshRegion, 2]] & /@ MeshPrimitives[myVoronoi, 2];

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

enter image description here

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]

enter image description here

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]

enter image description here

  • $\begingroup$ 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). $\endgroup$ Commented Mar 10, 2021 at 19:07
  • $\begingroup$ I'm using MMA 12.2 $\endgroup$ Commented Mar 10, 2021 at 19:37
  • $\begingroup$ @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. $\endgroup$
    – kglr
    Commented Mar 10, 2021 at 19:38
  • $\begingroup$ 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. $\endgroup$ Commented Mar 10, 2021 at 19:43
  • $\begingroup$ Also, I always get an error for Graphics`PolygonUtils`PolygonCombine @ Graphics`PolygonUtils`PolygonIntersection, how were you able to run this? I have MMA 12.2, maybe a compatibility issue? $\endgroup$ Commented Mar 24, 2021 at 23:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.