Crop a Voronoi diagram and get a proper MeshRegion

Question

How can I crop a VoronoiMesh (or any similar MeshRegion) to a non-rectangular region and receive the output as a MeshRegion that properly encodes cell-adjacency information?

Example:

voronoi = VoronoiMesh@RandomPoint[Disk[], 100];

I would like to crop this to the following disk, or an approximation of it as a polygon (i.e. BoundaryDiscretizeRegion@Disk[])

Show[voronoi, Graphics@Circle[]]

The output should look like this, but be a correct MeshRegion:

pieces = BoundaryDiscretizeRegion@RegionIntersection[Disk[], #] & /@ MeshPrimitives[voronoi, 2]; // AbsoluteTiming
(* {29.715, Null} *)

This is much too slow at 30 seconds and only 100 cells. The solution should work for thousands of cells.

Show[pieces]

This result is not a MeshRegion. The pieces could be assembled into one, but doing so is non-trivial as it is important that no point is duplicated, i.e. when vertices of neighbouring polygons coincide, they should be represented with a single point. One way to do this is the technique based on Nearest shown by Henrik Schumacher.

---

A more general formulation of the question would be:

Given a MeshRegion in 2D or 3D and a BoundaryMeshRegion of the same dimensionality, how can we crop all cells of the MeshRegion to be contained within the BoundaryMeshRegion?

Note that the naïve

RegionIntersection[voronoi, BoundaryDiscretizeRegion@Disk[]]

will not preserve the internal cells of the MeshRegion.

Answer 1

First Try (doesn't work due to a bug)

That's very frustrating. I thought that it might be a good idea to discretize the disk first and compute the intersections afterwards. That improves the timing by a factor of 20 but it gets some of the intersections wrong.

disk = BoundaryDiscretizeRegion[Disk[], MaxCellMeasure -> 0.01];
cells = MeshPrimitives[voronoi, 2];
pieces2 = BoundaryDiscretizeRegion@RegionIntersection[disk, #] & /@ cells; // AbsoluteTiming // First
Show[
 Cases[pieces2, _BoundaryMeshRegion],
 Graphics[{Red,
   Extract[cells, 
    Position[pieces2, Except[_BoundaryMeshRegion], 1, Heads -> False]]
   }]
 ]

1.45424

This is the result by Mmm 11.3. The intersections of the red cells with the circle could not be computed (probably a bug?):

Second Try

Not as simple anymore but quite efficient and maybe sufficiently general.

Let's say R is a MeshRegion and B is BoundaryMeshRegion. We want to crop R to be contained in B.

The strategy: Cells of R are divided into three groups: internalcells which are entirely contained in B, externalcells which do not intersect with B, and the bndcells which intersect the boundary. The internalcells can just be copied, the externalcells get discarded and the bndcells are first discretized by BoundaryDiscretizeRegion and then intersected with B. In the end, everythng is merged by the method that you fortunately mentioned. (I completely forgot about that post). Note that also IGMeshCellAdjacencyMatrix comes in handy here, in counting the number of vertices of each cell that lie within B. Note however that this is guaranteed to work correctly only if each cell of R is also convex (a mild condition) and if B is also convex (a more restrictive condition).

Needs["IGraphM`"];
cropRegion[R_MeshRegion, B_BoundaryMeshRegion] := 
 Module[{f, pts, cells, booleans, A, intersec, internalcells, 
   externalcells, bndcells, bndpieces, ϵ, plist, newpts, nf, 
   cellprimitives, celldata},
  f = RegionMember[B];
  pts = MeshCoordinates[R];
  cellprimitives = MeshPrimitives[R, 2];
  booleans = Map[Boole@*f, pts];
  
  A = IGMeshCellAdjacencyMatrix[R, 2, 0];
  intersec = (A.booleans);
  
  internalcells = Flatten[Position[Subtract[intersec, Total[A, {2}]], 0]];
  externalcells = Flatten[Position[intersec, 0]];
  bndcells = Complement[Range[Length[A]], internalcells, externalcells];
  
  bndpieces = ParallelMap[
     RegionIntersection[B, With[{p = #[[1]]}, BoundaryMesh@MeshRegion[p, Head[#][Range[Length[p]]]]]] &, 
     cellprimitives[[bndcells]]
     ];
  
  celldata = Join[
    Join @@ (MeshPrimitives[#, 2] & /@ bndpieces)[[All, All, 1]],
    cellprimitives[[internalcells, 1]]
    ];
  ϵ = 0.5 Min[
     PropertyValue[{R, 1}, MeshCellMeasure],
     PropertyValue[{#, 1}, MeshCellMeasure] & /@ bndpieces
     ]; 
  pts = Flatten[celldata, 1]; 
  plist = Sort[DeleteDuplicates[
     Compile[{{idx, _Integer, 1}}, 
       Min[idx], 
       RuntimeAttributes -> {Listable}, 
       Parallelization -> True
       ][Nearest[pts -> Automatic, pts, {∞, ϵ}]]
     ]]; 
  newpts = pts[[plist]]; nf = Nearest[newpts -> Automatic]; 
  cells = Internal`PartitionRagged[Flatten[nf[pts]], Length /@ celldata];
  MeshRegion[newpts, Polygon[cells]]
  ];

Usage example and timing

voronoi = VoronoiMesh@RandomPoint[Disk[], 10000];
r = 1;
disk = BoundaryDiscretizeRegion[Disk[{0, 0}, r], MaxCellMeasure -> 0.01];
S = cropRegion[voronoi, disk]; // RepeatedTiming // First
S

0.416

The original Voronoi mesh with 100 cells gets processed in about 0.04 seconds. That's about 700 times as fast as what we have started from.

There is still some optimization potential since we do not really need to search the interior vertices of the interior region for duplicates...

Warning

If B is not convex, then some less pleasant things can happen. For example, using this star fish

B = BoundaryMeshRegion[
  Map[t \[Function] (2 + Cos[5 t])/3 {Cos[t], Sin[t]}, Most@Subdivide[0., 2. Pi, 2000]],
  Line[Partition[Range[2000], 2, 1, 1]]
  ]

can lead to

SeedRandom[17];
cropRegion[VoronoiMesh@RandomPoint[Disk[], 100], B]

But this becomes more and more negligible, the smaller the cells in R are:

cropRegion[VoronoiMesh@RandomPoint[Disk[], 10000], B]

Answer 2

pieces2 = BoundaryDiscretizeRegion /@ 
  (Graphics`PolygonUtils`PolygonIntersection[#, Polygon[CirclePoints[100]]] &/@ 
    MeshPrimitives[voronoi, 2]); // AbsoluteTiming // First

1.09499

(vs 33.0193 for pieces = ...)

Show[pieces2]

With 1000 points it takes 11.3679 seconds to get pieces2, and Show[pieces2] gives:

Answer 3

Here's another implementation based on the idea of grabbing cells fully contained in the disk and intersecting the ones with partial overlap.

maskedMeshRegion[mr_?MeshRegionQ, mask_?ConstantRegionQ, δ_:0.01] :=
  Block[{prims, inside, frontier, bmr, fpieces, pts, cells},
    prims = MeshPrimitives[mr, 2];

    inside = Select[prims, RegionWithin[mask, #]&];

    frontier = Select[Complement[prims, inside], !RegionDisjoint[mask, #]&];
    bmr = BoundaryDiscretizeRegion[mask, MaxCellMeasure -> δ];
    fpieces = RegionIntersection[bmr, 
      BoundaryMeshRegion[#, Line[Append[Range[Length[#]], 1]]]]& /@ frontier[[All, 1]];

    pts = Union[Join @@ MeshCoordinates /@ fpieces, Join @@ inside[[All, 1]]];
    cells = Join[
      Join @@ (MeshPrimitives[#, 2]& /@ fpieces), 
      inside
    ] /. Dispatch[Thread[pts -> Range[Length[pts]]]];

    MeshRegion[pts, cells]
]

A test:

voronoi = VoronoiMesh@RandomPoint[Disk[], 10000];
maskedMeshRegion[voronoi, Disk[]] // AbsoluteTiming