算法计算球体上的Voronoi图？

我正在寻找一个简单的(如果存在的)算法来找到球体表面上一组点的Voronoi图.源代码会很棒.我是德尔福人(是的,我知道......),但我也吃C代码.

2016年7月更新:

感谢许多志愿者(尤其是Nikolai Nowaczyk和我),现在有更强大/正确的代码用于在Python中处理球体表面上的Voronoi图.这是scipy.spatial.SphericalVoronoi从0.18scipy 版本开始正式提供的.在官方文档中有一个使用和绘图的工作示例.

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

努力改进matplotlib中球形多边形的绘制

努力改进scipy中球面多边形表面积计算的处理

有关与此Python代码和相关计算几何工作相关的理论/开发/挑战的进一步背景,您还可以查看来自Nikolai和我的一些演讲:

Nikolai PyData伦敦2016年演讲

Tyler PyData伦敦2015谈话

Tyler PyCon 2016计算几何教程

原答案:

我实际上最近在球体表面上为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的网站上有一些关于JavaScript实现的讨论(http://www.jasondavies.com/maps/voronoi/) .我不认为他的代码是开放的.我还看到一篇关于使用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

这是一篇关于球形Voronoi图的论文.

或者如果你想要Fortran(bleah!)就有这个网站.

请注意,球体上的Delaunay三角剖分仅仅是凸包.因此,您可以计算3D凸包(例如使用CGAL)并采用双重.