[SciPy-User] ND interpolation with Qhull
Wolfgang Kerzendorf
I'm interested in learning how the LinearNDInterpolator actually works.
So I read up on qhull, convex hulls and delauney triangulation:
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.
So now that I have the triangles I just need to find the triangle
containing the point to be interpolated. And that is where I'm a little
bit unclear: How do I find the point?
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. But this requires me to go through all
triangles. I'm sure there is a faster way (which is probably used by scipy).
Once I have the triangle I just determine the two coefficients (in two
dimensions) and add the vectors up to get the interpolation?
Help is much appreciated
Wolfgang
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
Pauli Virtanen
[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
Pauli Virtanen
ke, 2011-06-08 kello 03:47 -0700, denis kirjoitti:
> A trick for smoothing barycentric interpolation is to warp big
> coefficients toward 1,
> so that each vertex Vj attracts nearby X more strongly:
>
> In: X = sum cj Vj, inside hull of Vj
> Zj = value at Vj
> Out: sum( warp(cj) Zj ) / sum( warp(cj) )
> instead of sum( cj Zj )
Yeah, you can smooth things inside the triangle. However, you will still
get only C0 continuity, as there are discontinuities in the gradient at
the triangle boundaries. Of course, whether this matters depends on the
problem.
Pauli