Очки не учитываются, когда рядом в scipy.spatial. Delaunay - PullRequest
Новая
       16

Очки не учитываются, когда рядом в scipy.spatial. Delaunay

3 голосов
/ 10 ноября 2011

Я замечаю необъяснимое поведение при сравнении процедур триангуляции Делоне для scipy (0.9.0) и matplotlib (1.0.1).Мои точки - это UTM координаты, хранящиеся в numpy.array([[easting, northing], [easting, northing], [easting, northing]]).Края Сципи пропускают некоторые из моих очков, в то время как у 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

1 Ответ

6 голосов
/ 21 ноября 2011

Такое поведение scipy.spatial.Delaunay может быть связано с неточностью арифметики с плавающей запятой.

Как вы, возможно, знаете, scipy.spatial.Delaunay использует библиотеку C qhull для расчета триангуляции Делоне. Qhull, в свою очередь, является реализацией алгоритма Quickhull , который подробно описан авторами в этой статье (1). Вы также можете знать, что арифметика с плавающей точкой, используемая в компьютерах, выполняется с использованием стандарта IEEE 754 (об этом вы можете прочитать, например, в Wikipedia ). Согласно стандарту, каждое конечное число проще всего описывается тремя целыми числами: s = знак (ноль или один), c = значение и (или «коэффициент»), q = показатель степени. Количество битов, используемых для представления этих целых чисел, зависит от типа данных. Итак, очевидно, что плотность распределения чисел с плавающей запятой по числовой оси не постоянна - чем больше числа, тем более свободно они распределены. Это можно увидеть даже с помощью калькулятора Google - вы можете вычесть 3333333333333333 из 3333333333333334 и получить 0 . Это происходит потому, что 3333333333333333 и 3333333333333334 оба округлены до одного и того же числа с плавающей запятой.

Теперь, зная об ошибках округления, мы могли бы прочитать главу 4 статьи (1), озаглавленную Копирование с неточностью . В этой главе описан алгоритм работы с ошибками округления: 

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]

Это перемещает вычисления в область, где числа с плавающей запятой достаточно плотные.

...