在球体上计算 Voronoi 图的算法？

Question

我正在寻找一种简单的(如果存在的话)算法来查找球体表面上一组点的 Voronoi 图。源代码会很棒。我是 Delphi 人(是的，我知道...)，但我也吃 C 代码。

Answer 1

2016 年 7 月更新:

感谢许多志愿者(尤其是 Nikolai Nowaczyk 和我)，现在有了更强大/更正确的代码来处理 Python 中球体表面的 Voronoi 图。从 scipy 的 0.18 版本开始，它作为 scipy.spatial.SphericalVoronoi 正式可用。官方 docs 中有一个使用和绘图的工作示例.

该算法遵循二次时间复杂度。虽然对数线性是球体表面 Voronoi 图的理论最优值，但这是目前我们能够实现的最佳值。如果您想了解更多信息并帮助开发工作，可以找到一些与改进 Python 处理球形 Voronoi 图和相关数据结构的方式相关的未解决问题:

• Effort to improve the plotting of spherical polygons in matplotlib
• Effort to improve the handling of spherical polygon surface area calculations in scipy

有关与此 Python 代码和相关计算几何工作相关的理论/开发/挑战的更多背景信息，您还可以查看 Nikolai 和我的一些谈话:

• Nikolai PyData London 2016 talk
• Tyler PyData London 2015 talk
• Tyler PyCon 2016 Computational Geometry tutorial

---

原答案:

实际上，我最近为球体表面的 Voronoi 图编写了一些开源 Python 代码:https://github.com/tylerjereddy/py_sphere_Voronoi

用法、算法和限制记录在 readthedocs ( http://py-sphere-voronoi.readthedocs.org/en/latest/voronoi_utility.html ) 上。那里有一些详细的例子，但我也会在下面放一两个例子。该模块还处理 Voronoi 区域表面积的计算，尽管在当前开发版本中存在一些数值缺陷。

我还没有看到很多关于球形 Voronoi 图的有据可查的开源实现，但是在 Jason Davies 的网站 (http://www.jasondavies.com/maps/voronoi/) 上有一些关于 JavaScript 实现的讨论。我不认为他的代码是开放的。我还看到一篇关于使用 Python 处理部分问题的博文 (http://jellymatter.com/2014/01/29/voronoi-tessellation-on-the-surface-of-a-sphere-python-code/)。上述帖子中引用的许多主要文献来源似乎很难实现(我尝试了其中一些)，但也许有些人会发现我的实现很有用，甚至会提出改进方法。

示例:

1) 为单位球体上的一组伪随机点生成 Voronoi 图:

import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import numpy as np
import scipy as sp
import voronoi_utility
#pin down the pseudo random number generator (prng) object to avoid certain pathological generator sets
prng = np.random.RandomState(117) #otherwise, would need to filter the random data to ensure Voronoi diagram is possible
#produce 1000 random points on the unit sphere using the above seed
random_coordinate_array = voronoi_utility.generate_random_array_spherical_generators(1000,1.0,prng)
#produce the Voronoi diagram data
voronoi_instance = voronoi_utility.Voronoi_Sphere_Surface(random_coordinate_array,1.0)
dictionary_voronoi_polygon_vertices = voronoi_instance.voronoi_region_vertices_spherical_surface()
#plot the Voronoi diagram
fig = plt.figure()
fig.set_size_inches(2,2)
ax = fig.add_subplot(111, projection='3d')
for generator_index, voronoi_region in dictionary_voronoi_polygon_vertices.iteritems():
   random_color = colors.rgb2hex(sp.rand(3))
   #fill in the Voronoi region (polygon) that contains the generator:
   polygon = Poly3DCollection([voronoi_region],alpha=1.0)
   polygon.set_color(random_color)
   ax.add_collection3d(polygon)
ax.set_xlim(-1,1);ax.set_ylim(-1,1);ax.set_zlim(-1,1);
ax.set_xticks([-1,1]);ax.set_yticks([-1,1]);ax.set_zticks([-1,1]); 
plt.tick_params(axis='both', which='major', labelsize=6)

2) 计算 Voronoi 区域多边形的表面积并验证重构的表面积是否合理:

import math
dictionary_voronoi_polygon_surface_areas = voronoi_instance.voronoi_region_surface_areas_spherical_surface()
theoretical_surface_area_unit_sphere = 4 * math.pi
reconstituted_surface_area_Voronoi_regions = sum(dictionary_voronoi_polygon_surface_areas.itervalues())
percent_area_recovery = round((reconstituted_surface_area_Voronoi_regions / theoretical_surface_area_unit_sphere) * 100., 5)
print percent_area_recovery
97.87551 #that seems reasonable for now