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?


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

enter image description here

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.


enter image description here

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.


3 Answers 3


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
 Cases[pieces2, _BoundaryMeshRegion],
    Position[pieces2, Except[_BoundaryMeshRegion], 1, Heads -> False]]


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

enter image description here

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

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]]]]]] &, 
  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}}, 
       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


enter image description here

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


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

enter image description here

can lead to

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

enter image description here

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

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

enter image description here

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


(vs 33.0193 for pieces = ...)


enter image description here

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

enter image description here


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), 
    ] /. Dispatch[Thread[pts -> Range[Length[pts]]]];

    MeshRegion[pts, cells]

A test:

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

enter image description here


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.