在 scipy.spatial.Delaunay 附近时忽略的点答案

【问题标题】：Points left out when nearby in scipy.spatial.Delaunay在 scipy.spatial.Delaunay 附近时忽略的点
【发布时间】：2011-12-25 15:54:16
【问题描述】：

在比较 scipy 的 (0.9.0) 和 matplotlib 的 (1.0.1) Delaunay 三角测量例程时，我注意到一个无法解释的行为。我的点是存储在numpy.array([[easting, northing], [easting, northing], [easting, northing]]) 中的 UTM 坐标。 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

【问题讨论】：

标签： python numpy matplotlib scipy delaunay

【解决方案1】：

scipy.spatial.Delaunay 的这种行为可能与浮点运算的印象有关。

您可能知道，scipy.spatial.Delaunay 使用 C 语言 qhull 库来计算 Delaunay 三角剖分。 Qhull 又是Quickhull algorithm 的实现，作者在this 论文(1) 中有详细描述。您也可能知道计算机中使用的浮点运算是使用 IEEE 754 标准执行的（例如，您可以在 Wikipedia 中了解它）。根据标准，每个有限数最简单地用三个整数描述：s= 一个符号（零或一），c= 一个有效数（或“系数”），q= 一个指数。用于表示这些整数的位数因数据类型而异。所以，很明显，浮点数在数值轴上的分布密度不是恒定的——数字越大，分布越松散。即使使用 Google 计算器也可以看到 - 您可以从 3333333333333334 和 get 0 中减去 3333333333333333。这是因为 3333333333333333 和 3333333333333334 都被舍入到相同的浮点数。

现在，知道了舍入误差，我们可能想阅读论文 (1) 的第 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 无法区分 2 个附近的点，因此它合并了两个方面，从而成功地消除了其中一个。为了消除这些错误，您可以尝试移动坐标原点：

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]

这会将计算转移到浮点数相当密集的区域。

【讨论】：