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

I'm looking for a simple (if exists) algorithm to find the Voronoi diagram for a set of points on the surface of a sphere. Source code would be great. I'm a Delphi man (yes, I know...), but I eat C-code too.

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

11 个解决方案

#1

Here's a paper on spherical Voronoi diagrams.

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

Or if you grok Fortran (bleah!) there's this site.

或者,如果你了解Fortran(bleah!),那就是这个网站。

#2

Update in July 2016:

Thanks to a number of volunteers (especially Nikolai Nowaczyk and I), there is now far more robust / correct code for handling Voronoi diagrams on the surface of a sphere in Python. This is officially available as scipy.spatial.SphericalVoronoi from version 0.18 of scipy onwards. There's a working example of usage and plotting in the official docs.

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

The algorithm follows quadratic time complexity. While loglinear is the theoretical optimum for Voronoi diagrams on the surfaces of spheres, this is currently the best we've been able to implement. If you'd like to find out more and help with the development effort there are some open issues related to improving the way Python handles spherical Voronoi diagrams and the related data structures:

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

• Effort to improve the plotting of spherical polygons in matplotlib

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

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

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

For further background on the theory / development / challenges related to this Python code and related computational geometry efforts you can also check out some talks from Nikolai and I:

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

• Nikolai PyData London 2016 talk

Nikolai PyData伦敦2016年演讲

• Tyler PyData London 2015 talk

Tyler PyData伦敦2015谈话

• Tyler PyCon 2016 Computational Geometry tutorial

Tyler PyCon 2016计算几何教程

---

Original Answer:

I've actually recently written some open source Python code for Voronoi diagrams on the surface of a sphere: https://github.com/tylerjereddy/py_sphere_Voronoi

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

The usage, algorithm, and limitations are documented on readthedocs (http://py-sphere-voronoi.readthedocs.org/en/latest/voronoi_utility.html). There are some detailed examples there but I'll place one or two below as well. The module also handles the calculation of the Voronoi region surface areas, albeit with some numerical weaknesses in the current development version.

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

I haven't seen many well-documented open source implementations for spherical Voronoi diagrams, but there has been a bit of buzz about the JavaScript implementation on Jason Davies' website (http://www.jasondavies.com/maps/voronoi/). I don't think his code is open though. I also saw a blog post about using Python to deal with part of the problem (http://jellymatter.com/2014/01/29/voronoi-tessellation-on-the-surface-of-a-sphere-python-code/). Many of the primary literature sources cited in the above posts seemed very challenging to implement (I tried some of them) but maybe some people will find my implementation useful or even suggest ways to improve it.

我还没有看到很多关于球形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 /)。以上帖子中引用的许多主要文献资料似乎都很难实现(我尝试了其中一些),但也许有些人会发现我的实现很有用,甚至建议改进它的方法。

Examples:

1) Produce a Voronoi diagram for a pseudo-random set of points on the unit sphere:

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) Calculate the surface areas of the Voronoi region polygons and verify that the reconstituted surface area is sensible:

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

#3

Notice that Delaunay triangulation on a sphere is just the convex hull. Thus you can compute the 3D convex hull (e.g. using CGAL) and take the dual.

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

#4

In short, try cssgrid from NCAR Graphics. I wrote a longer answer for a similar question at codereview.stackexchange.com.

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

#5

There is a paper from INRIA about the Delaunay Triangulation (DT) of points lying on a sphere: CAROLI, Manuel, et al. Robust and Efficient Delaunay triangulations of points on or close to a sphere. 2009. where they talk about an implementation in CGAL.

有一篇来自INRIA的论文,关于位于球体上的点的Delaunay三角剖分(DT):CAROLI,Manuel,et al。鲁棒和有效的Delaunay三角剖分在球体上或靠近球体。他们谈论CGAL的实施。

The paper refers to various available implementation of DT algorithms.

本文提到了DT算法的各种可用实现。

Quoting from the paper:

引自论文:

An easy and standard answer consists in computing the 3D convex hull of the points, which is notoriously equivalent.

一个简单而标准的答案在于计算点的三维凸包,这是众所周知的等价物。

for computing the convex hull the paper suggests:

为了计算凸包,本文建议:

1. Hull, a program for convex hulls.

赫尔,凸壳的程序。

2. Qhull.
3. Three-dimensional convex hulls. in FORTRAN.Three-dimensional convex hulls.

三维凸壳。在FORTRAN.Three维凸壳。

4. STRIPACK in FORTRAN.

FORTRAN中的STRIPACK。

The DT C++ class of CGAL has the method dual to get the Voronoi Diagram.

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

According to this post by Monique Teillaud (one of the author of the above mentioned paper) it seems to me that in November 2012 the implementation was not still ready.

根据Monique Teillaud(上述论文的作者之一)的这篇文章,在我看来,2012年11月的实施还没有准备好。

#6

It's been a while since the question has been answered, but I've found two papers that implement Fortune's algorithm (efficiency O(N lg N), memory O(N)) over the surface of the sphere. Maybe a future viewer will find this information useful.

问题已经回答了一段时间,但我发现两篇文章在球体表面实现了Fortune的算法(效率O(N lg N),记忆O(N))。也许未来的观众会发现这些信息很有用。

• "Sweeping the Sphere" by Dinis and Mamede, published in the 2010 International Symposium on Voronoi Diagrams in Science and Engineering. Can be purchased at http://dx.doi.org/10.1109/ISVD.2010.32

Dinis和Mamede的“席卷整个领域”,发表于2010年Voronoi图表科学与工程国际研讨会。可以在http://dx.doi.org/10.1109/ISVD.2010.32购买

• "A Plane Sweep Algorithm for the Voronoi Tessellation of the Sphere" by Zheng et al. I'm not sure it was published because of the first one, but it is dated 13 December 2011. It is available for free at http://www.e-lc.org/tmp/Xiaoyu__Zheng_2011_12_05_14_35_11.pdf

郑等人的“球面Voronoi镶嵌的平面扫描算法”。我不确定它是否因为第一个而被公布,但是它的日期是2011年12月13日。它可以在http://www.e-lc.org/tmp/Xiaoyu__Zheng_2011_12_05_14_35_11.pdf免费获得。

I'm working through them myself at the moment, so I can't explain it well. The basic idea is that Fortune's algorithm works on the surface of the sphere so long as you calculate the points' bounding parabolas correctly. Because the surface of the sphere wraps, you can also use a circular list to contain the beach line and not worry about handling cells at the edge of rectangular space. With that, you can sweep from the north pole of the sphere to the south and back up again, skipping to sites that introduce new points to the beach line (adding a parabola to the beach line) or the introduction of cell vertices (removing a parabola from the beach line).

我现在正在自己解决这些问题,所以我无法解释清楚。基本思想是,只要你正确地计算点的边界抛物线,Fortune的算法就可以在球体的表面上工作。由于球体表面包裹,您还可以使用圆形列表来包含海滩线,而不必担心处理矩形空间边缘的单元格。通过这种方式,您可以从球体的北极向南扫掠并再次向后移动,跳过到向海滩线引入新点的位置(向海滩线添加抛物线)或引入单元顶点(移除来自海滩线的抛物线)。

Both papers expect a high level of comfort with linear algebra to understand the concepts, and they both keep losing me at the point they start explaining the algorithm itself. Neither provide source code, unfortunately.

这两篇论文都期望线性代数具有高水平的舒适度来理解这些概念,并且他们在开始解释算法本身的时候都会让我失去理智。不幸的是,它们都没有提供源代码。

#7

I think the Voronoi plane for each point can be constructed using non-Euclidian geometry. What was normally a line on a 2d plane, is now a 'great circle' on the sphere (see Wikipedia:elliptic geometry). It is easy to find which points are on the wrong side of any great circle between two points, by simply rotating the sphere such that the dividing great circle is the equator, and then it's all points on the other hemisphere than the point you're constructing the Voronoi plane for.

我认为每个点的Voronoi平面可以使用非欧几里德几何构造。通常是2d平面上的线,现在是球体上的“大圆”(参见*:椭圆几何)。通过简单地旋转球体使得分割的大圆是赤道,很容易找到哪个点在两个点之间的任何大圆的错误一侧,然后它是另一个半球上的所有点而不是你的点。构建Voronoi平面。

This is not the entire answer, but this is where I'd start..

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

#8

There's a nice Voronoi diagram example program here (including source code for Delphi 5/6).

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

I think "points on the surface of a sphere" means that you first have to remap them to 2D-coordinates, create the Voronoi diagram and then remap them to sphere surface coordinates. Are the two formulas from Wikipedia UV mapping article working here?

我认为“球体表面上的点”意味着您首先必须将它们重新映射到2D坐标,创建Voronoi图,然后将它们重新映射到球体表面坐标。*UV映射文章中的两个公式是否在这里工作?

Also notice that the Voronoi diagram will have the wrong topology (it is inside a rectangle and does not "wrap around"), here it could help to copy all the points from (0,0)-(x, y) to the neighbour regions above (0, -y * 2)-(x, 0), below (0, y)-(x, y * 2), left (-x, 0)-(0, y) and right (x, 0)-(x*2, y). I hope you know what I mean, feel free to ask :)

另请注意,Voronoi图将具有错误的拓扑结构(它位于矩形内部并且不会“环绕”),这里它可以帮助将所有点从(0,0) - (x,y)复制到邻居区域(0,-y * 2) - (x,0),低于(0,y) - (x,y * 2),左( - x,0) - (0,y)和右(x, 0) - (x * 2,y)。我希望你知道我的意思,随便问:)

#9

CGAL is working on the "spherical kernel" package, which would allow to compute exactly these kind of things. Unfortunately, is not released yet, but maybe it will be in their next release, since they already mentioned it in a google tech talk in march

CGAL正在开发“球形内核”软件包,它可以准确地计算出这些东西。不幸的是,尚未发布,但也许它将在他们的下一个版本中,因为他们已经在3月份的Google技术谈话中提到了它

#10

Quoting from this reference: http://www.qhull.org/html/qdelaun.htm

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

To compute the Delaunay triangulation of points on a sphere, compute their convex hull. If the sphere is the unit sphere at the origin, the facet normals are the Voronoi vertices of the input.

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

#11

If your points are within one hemisphere, you could do a gnomonic projection from spherical to planar coordinates, and then triangulate, since great-circles become straight lines of shortest distance.

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