Bdtjtk

Why are some Voronoi regions empty?

I have a point cloud of the surface of a human that I'm trying to find the Voronoi diagram of. So I found scipy.spatial.Voronoi(), including some pretty good documentation. But sometimes vor.regions includes empty regions, and the documentation doesn't quite explain why they're there. At best it explains the relation to its Delaunay calculation "Note here that due to similar numerical precision issues as in Delaunay triangulation above, there may be fewer Voronoi regions than input points." But if so, it's only because qhull calculates Voronoi using a similar (or the exact same) algorithm as it uses for Delaunay Triangulation that this happens. Mathematically, the only reason a Voronoi region should be empty is if its points are duplicated

Specific goals:

I'm looking to replicate the results from this paper that use the Voronoi diagrams to calculate normal vectors at each input point. So I could do a workaround where I use a similar normal to a nearby point's properly calculated normal, or I could ignore any point that gives a zero Voronoi region. But I'd really like to know how to solve these empty regions so I can use as much of the input data as possible.

Background and other research

This answer suggests Voronoi(points, qhull_options='Qbb Qc Qx') and in the docs scipy suggests qhull_options="QJ Pp"

input data: http://columbia.edu/~nxb2101/minimal_stickfig_.npy

Code:

import numpy as np; np.set_err(all='raise')

import scipy.spatial

pt_cloud=np.load('minimal_stickfig_.npy')

vor=scipy.spatial.Voronoi(pt_cloud)



for idx in range(len(vor.regions)):

  region=vor.regions[idx]

  vertices = vor.vertices[region]

  hull=scipy.spatial.ConvexHull(vertices)

  volume=hull.volume

  triangle_mesh_hull=vertices[hull.simplices] # (n,3,3)

  # triang mesh calculation taken from

  #   https://stackoverflow.com/questions/26434726/return-surface-triangle-of-3d-scipy-spatial-delaunay/26516915



  inner_pt = np.mean(vertices[:2],axis=0).reshape((1,3))

  CoM=np.zeros((3,))

  pif("inner_pt is {0}".format(inner_pt))

  for triangle in triangle_mesh_hull:

    pif("triangle is: n{0}".format(triangle))

    tetra=np.concatenate((inner_pt,triangle),axis=0)

    CoM_tetra=np.mean(tetra,axis=0)

    vol_tetra=volume_tetra(tetra)

    CoM+=(CoM_tetra*vol_tetra)

  CoM /= volume

  # do more with CoM, and volume

Any kind of help is appreciated; if you're not sure why qhull does this or how to get around it, please let me know what you did to make a high-quality mesh from a point cloud