Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to obtain the gradient of a scipy interpolant directly?

I have a largish 3D numpy array of scalar values (OK call it a "volume" if you must). I want to interpolate a smooth scalar field over this at a succession of irregular, not all known up-front non-integral xyz coordinates.

Now Scipy's support for this is just excellent: I filter the volume with

filtered_volume = scipy.ndimage.interpolation.spline_filter(volume)

and invoke

scipy.ndimage.interpolation.map_coordinates(
    filtered_volume,
    [[z],[y],[x]],
    prefilter=False)

for (x,y,z) of interest to obtain apparently nicely behaved (smooth etc) interpolated values.

So far so good. However, my application also needs the local derivatives of the interpolated field. Currently I obtain these by central-differencing: I also sample the volume at 6 additional points (this can at least be done with just one call to map_coordinates) and calculate e.g the x derivative from (i(x+h,y,z)-i(x-h,y,z))/(2*h). (Yes I know I could reduce the number of additional taps to 3 and do "one sided" differences, but the asymmetry would annoy me.)

My instinct is that there ought to be a more direct way of obtaining the gradient but I don't know enough spline math (yet) to figure it out, or understand what's going on in the guts of the Scipy implementation: scipy/scipy/ndimage/src/ni_interpolation.c.

Is there a better way of obtaining my gradients "more directly" than central differencing ? Preferably one which allows them to be obtained using the existing functionality rather than hacking on Scipy's innards.

like image 560
timday Avatar asked Aug 01 '11 21:08

timday


People also ask

How does Scipy interpolation work?

The interp1d class in the scipy. interpolate is a convenient method to create a function based on fixed data points, which can be evaluated anywhere within the domain defined by the given data using linear interpolation. By using the above data, let us create a interpolate function and draw a new interpolated graph.

What does Scipy interpolate interp1d return?

Interpolate a 1-D function. x and y are arrays of values used to approximate some function f: y = f(x). This class returns a function whose call method uses interpolation to find the value of new points.

What is interp1d in Python?

The interp1d() function of scipy. interpolate package is used to interpolate a 1-D function. It takes arrays of values such as x and y to approximate some function y = f(x) and then uses interpolation to find the value of new points.

Is interp1d linear interpolation?

vq = interp1( x , v , xq ) returns interpolated values of a 1-D function at specific query points using linear interpolation.


1 Answers

Aha: according to the classic paper on splines cited in the numpy code, splines of order n and their derivatives are related by

   n          n-1           n-1
 dB (x)/dx = B   (x+1/2) - B   (x-1/2)

So using SciPy's spline interpolation I could get my derivatives by also maintaining a lower-order prefiltered volume and querying that a couple of times per derivative. This means adding a fair amount of memory (maybe competition with the "main" volume for cache), but presumably evaluation of the lower order splines is faster, so it's not obvious to me whether it would be faster or not overall than the central differencing using small offsets I'm doing currently. Haven't tried it yet.

like image 72
timday Avatar answered Sep 30 '22 12:09

timday