scipy.spatial.Delaunay中被遗漏的邻近点

Question

我发现了一个奇怪的情况，比较了scipy（0.9.0）和matplotlib（1.0.1）中的Delaunay三角剖分功能。我的点是UTM坐标，存储在numpy.array([[easting, northing], [easting, northing], [easting, northing]])里。scipy的边缘缺少了一些我的点，而matplotlib的则都包含了。请问有什么解决办法，还是我哪里做错了？

import scipy
import numpy
from scipy.spatial import Delaunay
import matplotlib.delaunay


def delaunay_edges(points):
    d = scipy.spatial.Delaunay(points)
    s = d.vertices
    return numpy.vstack((s[:,:2], s[:,1:], s[:,::-2]))


def delaunay_edges_matplotlib(points):
        cens, edges, tri, neig = matplotlib.delaunay.delaunay(points[:,0], points[:,1])
        return edges


points = numpy.array([[500000.25, 6220000.25],[500000.5, 6220000.5],[500001.0, 6220001.0],[500002.0, 6220003.0],[500003.0, 6220005.0]])

edges1 = delaunay_edges(points)
edges2 = delaunay_edges_matplotlib(points)

numpy.unique(edges1).shape # Some points missing, presumably nearby ones
numpy.unique(edges2).shape # Includes all points

Answer 1

这个 scipy.spatial.Delaunay 的行为可能和浮点运算的不精确性有关。

你可能知道，scipy.spatial.Delaunay 是通过 C 语言的 qhull 库来计算 Delaunay 三角剖分的。而 qhull 实际上是实现了一个叫做 Quickhull 算法，这个算法的详细信息可以在 这篇论文 (1) 中找到。你可能也知道，计算机使用的浮点运算是按照 IEEE 754 标准进行的（你可以在 维基百科 上了解更多）。根据这个标准，每个有限的数字可以用三个整数来简单描述：s= 符号（0 或 1），c= 有效数字（或称为“系数”），q= 指数。用来表示这些整数的位数根据数据类型的不同而变化。因此，很明显，浮点数在数轴上的分布密度并不是恒定的——数字越大，它们的分布就越稀疏。你甚至可以在谷歌计算器上看到这一点——你可以用 3333333333333334 减去 3333333333333333，结果竟然是 0。这是因为 3333333333333333 和 3333333333333334 都被四舍五入成了同一个浮点数。

现在，了解了舍入误差后，我们可能想看看这篇论文的第 4 章，标题是 应对不精确性。在这一章中，描述了一种处理舍入误差的算法：

Quickhull partitions a point and determines its horizon facets by computing 
whether the point is above or below a hyperplane. We have assumed that 
computations return consistent results ... With floating-point arithmetic, we 
cannot prevent errors from occurring, but we can repair the damage after 
processing a point. We use brute force: if adjacent facets are nonconvex, one of 
the facets is merged into a neighbor. Quickhull merges the facet that minimizes 
the maximum distance of a vertex to the neighbor.

所以可能发生的情况是——Quickhull 无法区分两个相近的点，因此它将两个面合并，从而成功地消除了其中一个。为了消除这些错误，你可以尝试移动你的坐标原点：

co = points[0]
points = points - co

edges1 = delaunay_edges(points)
edges2 = delaunay_edges_matplotlib(points)

print numpy.unique(edges1) 
>>> [0 1 2 3 4]
print numpy.unique(edges2)
>>> [0 1 2 3 4]

这样可以将计算移动到浮点数相对密集的区域。