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

Question

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

Answer 1

2016 年 7 月更新：

感谢许多志愿者（特别是 Nikolai Nowaczyk 和我），现在有了更强大/正确的代码来处理 Python 中球体表面的 Voronoi 图。 从 scipy 版本 0.18 开始，它以 scipy.spatial.SphericalVoronoi 形式正式提供。 官方文档中有一个使用和绘图的工作示例。

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

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

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

• Nikolai PyData 伦敦 2016 演讲
• < a href="https://www.youtube.com/watch?v=gxNa9BD5CnQ" rel="noreferrer">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 的网站 (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

Answer 2

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

或者，如果您精通 Fortran（天啊！），那么这个网站。

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

Answer 3

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

Answer 4

简而言之，尝试 cssgrid 来自 NCAR 图形。 我在codereview.stackexchange.com上为类似问题写了一个更长的答案。

Answer 5

INRIA 有一篇关于球体上点的 Delaunay 三角测量 (DT) 的论文：CAROLI，Manuel，等人。 对球体上或球体附近的点进行稳健且高效的 Delaunay 三角剖分。 2009 年。 他们讨论了 CGAL 中的实现。

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

引用论文：

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

为了计算凸包，论文建议：

1. Hull，一个凸包程序。
2. Qhull。
3. 三维凸包。 FORTRAN 格式。三维凸包。
4. FORTRAN 语言中的 STRIPACK。

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

根据这篇文章< /a> 作者：Monique Teillaud（上述论文的作者之一） 在我看来，2012 年 11 月，实施尚未准备就绪。

Answer 6

这个问题已经有一段时间了，但我发现了两篇实现Fortune 算法< /a> （效率 O(N lg N)，内存 O(N)）在球体表面。 也许未来的观看者会发现这些信息很有用。

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

我现在正在自己研究它们，所以我无法很好地解释。 基本思想是，只要正确计算点的边界抛物线，《财富》算法就适用于球体表面。 由于球体的表面是包裹的，因此您还可以使用圆形列表来包含海滩线，而不必担心处理矩形空间边缘的单元格。 这样，您可以从球体的北极向南扫过，然后再次返回，跳至向海滩线引入新点（向海滩线添加抛物线）或引入单元顶点（删除单元顶点）的站点。海滩线的抛物线）。

这两篇论文都希望读者能够对线性代数有较高的了解，以理解这些概念，但在开始解释算法本身时，它们都让我迷失了方向。 不幸的是，两者都没有提供源代码。

Answer 7

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

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

Answer 8

这里有一个很好的 Voronoi 图示例程序（包括 Delphi 5/的源代码） 6).

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

另请注意，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 正在开发“spherical kernel”包，它可以精确地计算这类事情。 不幸的是，尚未发布，但也许会在下一个版本中发布，因为他们已经在 3 月份的 Google 技术演讲中提到过

Answer 10

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

Answer 11

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

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