How to get uneven interpolation gradient directly?

Question

I have a fairly large array of 3D numbers with numeric values ​​(if you do, click OK). I want to interpolate a smooth scalar field over this in a sequence of irregular, not all previously known non-integer xyz coordinates.

Now, Scipy support for this is simply excellent: I filter the volume with

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

and call

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

for (x, y, z), which are of interest for obtaining apparently beautifully executed (smooth, etc.) interpolated values.

So far so good. However, my application also needs local derivatives of the interpolated field. Currently, I get them according to the central difference: I also tried the volume at 6 additional points (this can be done with at least one call to map_coordinates ) and calculate, for example, the derivative of x from (i(x+h,y,z)-i(xh,y,z))/(2*h) . (Yes, I know that I could reduce the number of extra taps to 3 and make “one-way” differences, but the asymmetry annoys me.)

My instinct is that there should be a more direct way to get the gradient, but I don’t know enough spline math (for now) to understand this, or to understand what is going on in the guts of the Scipy implementation: scipy/scipy/ndimage/src/ni_interpolation.c .

Is there a better way to get my gradients "more straightforward" than the central difference? Preferably one that allows them to get using existing functionality rather than hacking inside Scipy shells.

Answer 1

Aha: according to the classical work on splines indicated in 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) 

Thus, using SciPy spline interpolation, I could get my derivatives, also supporting a pre-filtered lower order volume and querying a couple of times for the derivative. This means adding enough memory (perhaps competing with the “main” cache volume), but it seems that lower-order splines are estimated faster, so it’s not clear to me whether this will be faster or not overall than the central difference using the small offsets I'm doing now. I have not tried it yet.