[SciPy-User] Scipy's Delaunay triangulation

[Pauli Virtanen]

Hi,

21.02.2012 15:26, barbara.padova kirjoitti:
> I am noticing an unexplained behaviour when use scipy's (0.9.0)
> Delaunay triangulation routine. My points are UTM coordinates stored
> in numpy.array([[easting, northing], [easting, northing], [easting,
> northing]]).
> If I use the true coordinates Scipy's edges are missing some of my
> points.
[clip]

Delaunay triangulations are not well-defined for all input data sets,
and are often even less well-defined numerically.

The behavior of Scipy's triangulation routine is inherited from the
computational geometry library it uses, Qhull (http://qhull.org/). In
some cases, Qhull can exclude points that are ambiguous due to numerical
precision from the triangulation, and I think this is what you observe.

Scipy runs qhull with options "d Qz Qbb Qt" plus the defaults; you can
look up the precise meanings here: http://qhull.org/html/qh-quick.htm

One could argue that Scipy should raise an error in ambiguous cases, and
that we should not pass "Qbb" and "Qz" by default, but rather let the
user customize how they want to deal with the triangulation failure.
These "numerically bad" cases are not uncommon, though, and I'm not 100%
sure it is possible to tell Qhull to always include all points (see [1]).

> There is a solution to use my original data with scipy's Delaunay
> triangulation routine?

It's a limitation of the algorithm.

One option is to deduce which points are missing, and live with the fact
(e.g. considering them "merged" them with the nearest points in the
triangulation in your own code, or just ignoring them).

Qhull also has an option "QJ" which adds numerical noise to the least
significant digits of the input data, so that results become numerically
better defined. Unfortunately, at the moment the Scipy interface does
not allow passing this option down to Qhull.

You should however be able to achieve a similar effect by adding noise
to the input data also yourself, e.g.

	numpy.random.seed(1234)
	joggled = points * (
		1 + 1e-8*(numpy.random.rand(points.shape) - 0.5))
	tri = Delaunay(joggled)

it's not exactly the same as the "QJ" option, but should behave in the
same way.

-- 
Pauli Virtanen


[1] Try feeding the following to qdelaunay:

3
10
  -0.5 -0.5 -0.5
  -0.5 -0.5  0.5
  -0.5  0.5 -0.5
  -0.5  0.5  0.5
   0.5 -0.5 -0.5
   0.5 -0.5  0.5
   0.5  0.5 -0.5
   0.5  0.5  0.5000000000000003
   0.5 -0.5  0.500000000000001
   0.5  0.5  0.5