当附近scipy.spatial.Delaunay
当我比较scipy(0.9.0)和matplotlib(1.0.1)Delaunay三角测量例程时,我注意到一个无法解释的行为。我的观点是存储在numpy.array([[easting, northing], [easting, northing], [easting, northing]])
中的UTM坐标。 Scipy的边缘缺少我的一些观点,而matplotlib都在那里。有没有修复,或者我做错了什么?当附近scipy.spatial.Delaunay
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
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章,标题为复制与impresicion。在本章中,处理舍入误差的算法描述:
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]
该移动计算到浮点数是相当密集的区域。
感谢您的回答,一个很好的解决方案。 – Benjamin