[SciPy-User] ND interpolation with Qhull

[Pauli Virtanen]

Mon, 06 Jun 2011 20:40:10 +1000, Wolfgang Kerzendorf wrote:
[clip]
> I understand that one can use qhull to construct the convex hull in a 
> d-dimensional space. If I want the delauney triangulation of n points
> in d dimensions I just need to project these points on a paraboloid
> in d+1 dimensional space build the convex hull and reproject this
> onto d-dimensions.

Qhull actually does that for you -- you can ask it to directly provide
a Delaunay triangulation. But yes, that's how it works it out internally.

[clip]
> I know that in the barycentric coordinate system I have three
> coefficents and if the sum of two of them is less than 1 to reproduce my
> point, then I found my triangle.

You also need to require that the two coordinates are in the range [0, 1].

> But this requires me to go through all triangles. I'm sure there 
> is a faster way (which is probably used by scipy).

You can read the code to find out how it works:

https://github.com/scipy/scipy/blob/master/scipy/spatial/qhull.pyx#L771

Basically, you do a directed walk in the triangle neigbourhood graph.
There are two things you can do: first, you walk on the convex hull
in (d+1)-dim to the correct direction until you see the target point
over the horizon, and then you walk towards the target point in d-dim.

Or, you just do the directed walk in d-dim. Qhull itself uses the
"walk towards the horizon" approach, but in practice this doesn't
seem to be much better than the directed walk.

If the walk ends up in a degenerate triangle in the triangulation,
these approaches either fail or enter an infinite cycle, so you need
to fall back to brute force search through all the triangles.
Getting degenerate triangles in the triangulation seems difficult to
avoid in practice (people like to use these routines for data on
rectangular grids), so you just have to live with them.

Also, when you do interpolation, you start the walk from the point at which
you were last at, because the next point to interpolate is usually close
to the last one.

> Once I have the triangle I just determine the two coefficients (in two
> dimensions) and add the vectors up to get the interpolation?

You can use the three barycentric coordinates [c3 = 1 - c1 - c2] to
compute the weights you want. Barycentric interpolation is simply

	v = c1 v1 + c2 v3 + c3 v3

But if you want *smooth* spline interpolation rather than linear,
things get hairy, as you need to ensure that C1 continuity is satisfied
at the triangle boundaries. In 2D these conditions are not too difficult
to satisfy, but things start to get substantially more hairy in higher
dimensions. First, there are more matching conditions, and satisfying
them is more difficult. Second, you need higher-degree splines, and so
have more free parameters -- so in 3D you need to estimate not only
gradients but also the Hessians from the set of data points, etc etc.
As far as I know, a general solution for N dimensions is not known
so far. Instead, you have a number of cooking recipes in 3D and 4D.

To my understanding, in higher dimensions, natural neighbor interpolation
becomes easier to implement than spline interpolation, and IIRC, if done
correctly, it does provide global C1 smoothness. However, for natural
neighbor the computational complexity goes up fast with dimensions,
as you in the end need to do local modifications to the triangulation.
In principle, one could reuse Qhull here, but this is not done in Scipy
yet.

-- 
Pauli Virtanen