當附近scipy.spatial.Delaunay

Question

當我比較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 

Answer 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章，標題爲複製與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] 

該移動計算到浮點數是相當密集的區域。