BUG: trying to fix #1605, griddata spitting out nans #212

merged 6 commits into from May 27, 2012


None yet

2 participants

SciPy member
pv commented May 12, 2012

Now comes somewhat tricky stuff: ticket http://projects.scipy.org/scipy/ticket/1605 reports a bug in interpolate.griddata, which turns out to be due to numerical precision issues in "is-this-point-inside-this-simplex" checks. These turn up when the triangulation contains nearly-zero-volume simplices, which is actually not so rare in rectangular grids.

Unfortunately, this is sort of a release blocker, and should be solved for 0.11.0.

This patch set does the following:

  • Detect which simplices are "bad", i.e., for which we cannot reliably detect whether points are inside or not

  • This leaves in nearly-zero-volume "holes" inside the triangulation (which is always convex in reality), in which it is numerically unclear whether a point belongs to a simplex or not

  • To resolve this, we add sqrt(eps) sized margins for each simplex, so that the holes are blocked. Or, are they?

The problem here is that the broadening of the simplices by sqrt(eps) is obviously a fudge factor, and it's probably not guaranteed that even this is large enough in all cases. It also introduces errors of the same relative order of magnitude, sqrt(eps), for interpolation --- although this is likely not so bad.

I think this problem cannot be avoided just by telling Qhull to use its QJ joggling mode, as badly conditioned simplices probably still occur e.g. in manually perturbed rectangular grids.

I'll try to think up a better solution, but in case nothing comes up, this PR should be a reasonable workaround for practical purposes.

pv added some commits Feb 22, 2012
@pv pv BUG: spatial/qhull: discard degenerate barycentric transforms (fixes #…

For some point sets, the Delaunay triangulation produces nearly
degenerate simplices. These have singular barycentric transforms.

Previously, we discarded exactly singular transforms.  However,
computing the transform may succeed, but the result may still be invalid
and cause problems later on e.g. in _find_simplex_directed.

This commit adds a better check for validity of the transform: we
explicitly check the condition number, and drop those for which it is
large compared to machine precision.
@pv pv TST: spatial: move test data under tests/data 77e5b57
@pv pv TST: spatial/qhull: add tests for finding simplices in nearly degener…
…ate cases

Related to #1605
@pv pv GEN: spatial: regenerate qhull.c 45a8c75
SciPy member

The sqrt(eps) errors are introduced for all results, or only in case bad simplices are present? The latter should be OK as a workaround.

SciPy member

An alternative could be to check the output for nans where there shouldn't be any (just a thought, not sure if this is easy), and then fix them up.

Something like this nan inpainting function could be useful: http://www.mathworks.com/matlabcentral/fileexchange/4551

SciPy member
pv commented May 12, 2012

The sqrt(eps) are introduced for all simplices, but only sqrt(eps) close to the boundaries. However, the is-inside check for bad simplices could be replaced by a is-inside-neighbor check, which could be more robust.

In-painting cannot be done, since the output is also unstructured set of data points.

pv added some commits May 12, 2012
@pv pv BUG: spatial/qhull: use larger tolerance in find_simplex only for deg…
…enerate simplices

Use a larger tolerance in is-inside-check for simplices next to
degenerate simplices. This makes the code behave as previously in
well-formed triangulations, but makes it cope better with triangulations
with (nearly) degenerate simplices.
@pv pv GEN: regenerate interpnd.c and qhull.c 25ecea2
SciPy member
pv commented May 12, 2012

Ok, now it uses a larger epsilon only for faces neighboring degenerate simplices. This should be more controlled, and seems to work as well.

SciPy member

Sounds good.

All tests pass on OS X / python 2.6.

@rgommers rgommers merged commit d77fae0 into scipy:master May 27, 2012
SciPy member

OK, merged this. Time for 0.11 I think.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment