Robert Kern wrote:
> On 2009-07-13 13:20, Robert Cimrman wrote:
>> Hi all,
>>
>> I would like to use griddata() to interpolate a function given at
>> specified points of a bunch of other points. While the method works
>> well, it slows down considerably as the number of points to interpolate
>> to increases.
>>
>> The dependence of time/(number of points) is nonlinear (see the
>> attachment) - it seems that while the Delaunay trinagulation itself is
>> fast, I wonder how to speed-up the interpolation. The docstring says,
>> that it is based on "natural neighbor interpolation" - how are the
>> neighbors searched?
>
> Using the Delaunay triangulation. The "natural neighbors" of an interpolation
> point are those points participating in triangles in the Delaunay triangulation
> whose circumcircles include the interpolation point. The triangle that encloses
> the interpolation point is found by a standard walking procedure, then the
> neighboring triangles (natural or otherwise) are explored in a breadth-first
> search around the starting triangle to find the natural neighbors.
I see, thanks for the explanation. The walking procedure is what is
described e.g. in [1], right? (summary; starting from a random triangle,
a line is made connecting that triangle with the interpolation point,
and triangles along that line are probed.)
[1] http://www.geom.uiuc.edu/software/cglist/GeomDir/ptloc96.ps.gz
> Unfortunately, griddata() uses the unstructured-interpolation-points API rather
> than the more efficient grid-interpolation-points API. In the former, each
> interpolation point uses the last-found enclosing triangle as the start of the
> walking search. This works well where adjacent interpolation points are close to
> each other. This is not the case at the ends of the grid rows. The latter API is
> smarter and starts a new row of the grid with the triangle from the triangle
> from the *start* of the previous row rather than the end. I suspect this is
> largely the cause of the poor performance.
Good to know, I will try to pass the points in groups of close points.
>> Does it use the kd-trees like scipy.spatial? I have
>> a very good experience with scipy.spatial performance.
>>
>> Also, is there a way of reusing the triangulation when interpolating
>> several times using the same grid?
>
> One would construct a Triangulation() object with the (x,y) data points, get a
> new NNInterpolator() object using the .nn_interpolator(z) method for each new z
> data set, and then interpolate your grid on the NNInterpolator.
So if the above fails, I can bypass griddata() by using the delaunay
module directly, good.
thank you,
r.
|