Mesh cell count for Voronoi mesh too low

Question

I was trying to solve this question out of interest and thought perhaps creating a Voronoi mesh, cropping it to a circle, and colouring the mesh cells might work. However, if I ask VoronoiMesh to create cells for too many points MeshCellCount[mesh, 2] (or equivalently Length@MeshCells[mesh]) it returns a number that is smaller than the number of points provided initially.

I've tried using different functions to generate the points around which the cells should be built, used both exact and real numbers, and checked out the documentation for VoronoiMesh and MeshRegion, but I'm still not sure what's causing this. Are my points simply too close together for VoronoiMesh to uniquely determine a cell for some of them?

The simplest code that reproduces this is:

MeshCellCount[
  VoronoiMesh[
    Flatten[Quiet[Thread[CirclePoints[Range[100], 360]]], 1]
  ], 
2]

which should return 36,000 since it is 100 radial points and 360 azimuthal points, but instead returns 35,985. For this code it seems to start when there's around 32,000 elements. If the radial points inside Range are set to 87, I get the expected result. If the radial points are set to 88 (with the same 360 azimuthal points) I get an unexpected result. For all smaller numbers it seems to work as expected.

For some reason, if I use the following code to determine the number of cells, this discrepancy shows up at even smaller numbers of cells.

generate[i_] := 
 Table[
   {r Sin[θ], r Cos[θ]}, 
   {θ, 0, 359 π/180, π/180}, 
   {r, 1/2, (i - 1) + 1/2}
 ]
66*360 - MeshCellCount[VoronoiMesh[Flatten[generate[66], 1]], 2]

The result of this code is 2 where I would expect it to be zero for all values passed to generate.

Does anyone know what I'm doing wrong or if there is a workaround? Or am I simply asking too much of VoronoiMesh?

Answer 1

Let's look at the centroids of the faces to see if we can figure where the issue lies.

pts = Flatten[Quiet[Thread[CirclePoints[Range[100], 360]]], 1];

vor = VoronoiMesh[pts];

centroids = PropertyValue[{vor, 2}, MeshCellCentroid];
norms = Sort[Norm /@ centroids];

KeySelect[Counts[Round[norms]], LessThan[100]]

<|1 -> 345, 2 -> 360, 3 -> 360, 4 -> 360, 5 -> 360, 6 -> 360, ...|>

So it looks like there are 15 missing inner most faces. Let's take a look:

MeshRegion[
  MeshCoordinates[vor], 
  Pick[MeshCells[vor, 2], RegionMember[Disk[{0, 0}, 1], centroids]]
]

I don't know what went wrong, nor do I know how to fix the builtin behavior. But we can find a workaround by adapting the answer here:

Block[{Print}, <<IGraphM`];

Voronoi2D[pts_] :=
  Block[{minmax, padding, vpts, dm, prims, vnodes, conn, adj, vlines, mr1d, g, faces, lens},
    minmax = MinMax /@ Transpose[pts];
    padding = Abs[Subtract @@@ minmax];
    vpts = Join[pts, (minmax[[All, 1]] - padding)IdentityMatrix[2], (minmax[[All, 2]] + padding)IdentityMatrix[2]];

    dm = DelaunayMesh[vpts];
    prims = MeshPrimitives[dm, 2, "Multicells" -> True][[1, 1]];
    vnodes = circumCenter2D[prims];


    conn = dm["ConnectivityMatrix"[2, 1]];
    adj = conn.Transpose[conn];
    vlines = UpperTriangularize[adj, 1]["NonzeroPositions"];

    mr1d = Quiet @ MeshRegion[vnodes, Line[vlines]];

    g = IGMeshGraph[mr1d];
    faces = IGFaces[g];

    (* delete outer face *)
    lens = Length /@ faces;
    faces = Pick[faces, UnitStep[lens - Max[lens]], 0];

    MeshRegion[MeshCoordinates[mr1d], Polygon[faces]]
  ]

(* speed up from calling Circumsphere... but some rounding error could be introduced *)
With[{
  a = Det[{{x1, y1, 1}, {x2, y2, 1}, {x3, y3, 1}}],
  bx = Det[{{x1^2+y1^2, y1, 1}, {x2^2+y2^2, y2, 1}, {x3^2+y3^2, y3, 1}}],
  by = Det[{{x1^2+y1^2, x1, 1}, {x2^2+y2^2, x2, 1}, {x3^2+y3^2, x3, 1}}],
  ε = 2^22 * $MachineEpsilon
  },
  circumCenter2D = Compile[{{pts, _Real, 2}},
    Block[{x1, y1, x2, y2, x3, y3},
      x1 = pts[[1, 1]];
      y1 = pts[[1, 2]];
      x2 = pts[[2, 1]];
      y2 = pts[[2, 2]];
      x3 = pts[[3, 1]];
      y3 = pts[[3, 2]];

      Round[.5Divide[{bx, -by}, a], ε]
    ],
    CompilationTarget -> "C",
    Parallelization -> True,
    RuntimeOptions -> "Speed",
    RuntimeAttributes -> {Listable}
  ]
];

Your example:

Voronoi2D[pts] // MeshCellCount

{36691, 72752, 36000}