在球体上计算 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

Answer 2

这是一篇关于 spherical Voronoi diagrams 的论文。

或者，如果您熟悉 Fortran（哎呀！），那就是 this site。

原链接（失效）：https://people.sc.fsu.edu/~jburkardt/f_src/sxyz_voronoi/sxyz_voronoi.html

Answer 3

请注意，球体上的 Delaunay 三角剖分只是凸包。 因此，您可以计算 3D 凸包（例如，使用 CGAL） 并采取对偶。

Answer 4

简而言之，从NCAR Graphics 尝试cssgrid。我写了a longer answer for a similar question at codereview.stackexchange.com。

Answer 5

INRIA 有一篇关于球体上点的 Delaunay 三角剖分 (DT) 的论文：CAROLI, Manuel, et al. Robust and Efficient Delaunay triangulations of points on or close to a sphere. 2009.，他们在其中讨论了 CGAL 中的实现。

本文参考了 DT 算法的各种可用实现。

引自论文：

一个简单而标准的答案是计算 3D 凸包 的点，这是众所周知的等价。

论文建议计算凸包：

1. Hull, a program for convex hulls.
2. Qhull。
3. Three-dimensional convex hulls. 在 FORTRAN 中。三维凸包。
4. STRIPACK 在 FORTRAN 中。

CGAL的DT C++类有方法dual获取Voronoi图。

根据 Monique Teillaud（上述论文的作者之一）的 this post，在我看来，2012 年 11 月的实施还没有准备好。

Answer 6

这个问题已经有一段时间没有回答了，但是我发现了两篇论文在球体表面上实现了Fortune's algorithm（效率 O(N lg N)，内存 O(N)）。也许未来的观众会发现这些信息很有用。

• Dinis 和 Mamede 的“Sweeping the Sphere”，发表在 2010 年科学与工程中的 Voronoi 图国际研讨会上。可以在http://dx.doi.org/10.1109/ISVD.2010.32购买
• Zheng 等人的“用于球体 Voronoi 镶嵌的平面扫描算法”。由于第一个，我不确定它是否已发布，但它的日期为 2011 年 12 月 13 日。它可在http://www.e-lc.org/tmp/Xiaoyu__Zheng_2011_12_05_14_35_11.pdf免费获得

目前我自己正在处理它们，所以我无法很好地解释它。基本思想是，只要您正确计算点的边界抛物线，Fortune 的算法就可以在球体表面上工作。由于球体表面包裹，您还可以使用圆形列表来包含海滩线，而不必担心处理矩形空间边缘的单元格。有了这个，您可以从球体的北极向南扫过并再次返回，跳到向海滩线引入新点的站点（向海滩线添加抛物线）或引入单元顶点（删除一个海滩线的抛物线）。

两篇论文都希望通过线性代数来理解这些概念，并且在他们开始解释算法本身时都让我感到很困惑。不幸的是，两者都没有提供源代码。

Answer 7

我认为每个点的 Voronoi 平面都可以使用非欧几里得几何构造。通常是二维平面上的一条线，现在是球体上的一个“大圆圈”（参见维基百科：elliptic geometry）。很容易找到两个点之间的任何大圆的错误一侧的点，只需旋转球体，使分割大圆是赤道，然后它是另一个半球上的所有点而不是你所在的点构建 Voronoi 平面。

这不是完整的答案，但这是我要开始的地方..

Answer 8

有一个很好的 Voronoi 图示例程序 here（包括 Delphi 5/6 的源代码）。

我认为“球体表面上的点”意味着您首先必须将它们重新映射到 2D 坐标，创建 Voronoi 图，然后将它们重新映射到球体表面坐标。 Wikipedia UV mapping article 的两个公式在这里有效吗？

还要注意，Voronoi 图会有错误的拓扑结构（它在一个矩形内并且没有“环绕”），在这里它可以帮助复制 (0,0)-(x, y) 中的所有点到 (0, -y * 2)-(x, 0) 上方、(0, y)-(x, y * 2) 下方、左侧 (-x, 0)-(0, y) 和右侧的相邻区域(x, 0)-(x*2, y)。我希望你知道我的意思，随时问:)

Answer 9

CGAL 正在开发“球形内核”包，它可以精确计算这类东西。不幸的是，还没有发布，但可能会在他们的下一个版本中，因为他们已经mentioned it in a google tech talk in march

Answer 10

如果您的点在一个半球内，您可以进行从球面坐标到平面坐标的日光投影，然后进行三角测量，因为大圆变成最短距离的直线。

Answer 11

引用此参考：http://www.qhull.org/html/qdelaun.htm

要计算球体上点的 Delaunay 三角剖分，请计算它们的凸包。如果球体是原点的单位球体，则面法线是输入的 Voronoi 顶点。