Re: Interpolating in three dimensions

Hello,

I may be able to offer some comments as I recently had to implement the 
same method by Lekien and Marsden:

- It is a method to explicitly construct the spline coefficients, i.e., 
64 coefficients for every voxel. This is useful if one requires the 
repeated evaluation of the same three-dimensional function and/or 
several of its derivatives, e.g., at different possibly irregular 
positions for every function call. It may come with a performance and 
memory penalty compared to other interpolation methods, which construct 
the interpolant _or_ a _single_ derivative directly from the data 
without explicitly forming the coefficients. In other words, it may be 
not the best choice when the interpolant or a single derivative is only 
evaluated once.

- The interpolant is in C1, i.e., the first derivatives are continuous 
while its second derivative may be discontinuous. This is in contrast to 
a global tricubic spline interpolation which is in C2. As far as I know 
the method implemented in interpn() in scipy is the repeated application 
of a 1d cubic interpolation in each dimension, which should produce a C2 
interpolant. Matlab's interp3() makes this distinction clear with the 
'cubic' and 'spline' modes (mathworks.com/help/matlab/ref/interp3.html).

- The method is used in localization microscopy to model point spread 
functions when they are required in an optimization process. In this 
case the interpolant and its derivatives are evaluated repeatedly. This 
has spawned several implementations, e.g.:
     - https://github.com/TuragaLab/SplinePSF (a python package with the 
actual implementation in C++ and CUDA, usable and differentiable by pytorch)
     - 
https://github.com/ZhuangLab/storm-analysis/tree/master/storm_analysis/spliner 
(in C and Python, afaik)

- Other implementations are:
     - https://github.com/danielguterding/pytricubic (this is on pypi)
     - https://github.com/nbigaouette/libtricubic/ (the original C code 
by the authors, albeit GPL licensed here)

Two comments on Fabian's implementation:

- Using a sparse matrix format for the Binv matrix should improve 
performance.
- Using Horner's method to evaluate the interpolant and its derivatives 
should improve performance and accuracy.

IMHO, scipy would profit most if this interpolator is included together 
with an update of the documentation to be clear and explicit on the 
advantages and drawbacks of the different interpolation methods (e.g., 
continuity class, relative memory and computation requirements, 
requirements on the grid points (number, regularity), and suitability 
for repeated evaluation).
Also, the big advantage of Lekien and Marsden's method realizes if you 
want the derivatives in addition to the function value itself. As far as 
I see this is currently not exposed in the multivariate interpolator 
classes and functions - even if the interpolation methods themselves may 
be able to provide this information. I think this requires a discussion 
about how this should be included.

I would be interested in helping bringing an tricubic interpolator 
suitable for repeated evaluation to scipy. I think it is in scope, used 
in my field and currently there is no high-quality implementation that 
can be used in standard environments (no cuda etc).

Cheers
Nils

On 06.07.23 11:01, Fabian Gittins via SciPy-Dev wrote:
> Hi both,
>
> Thanks for the warm welcome and the useful feedback! I am grateful for the references.
>
> Stéfan: Thank you for the suggestions. I will compare my implementation with some of SciPy's offerings. If derivatives can be estimated using RBFInterpolator that would be a useful benchmark (and potentially better than my current method).
>
> I was able to find the repository for the quadcubic interpolator you referenced (https://github.com/DurhamDecLab/ARBInterp). While obviously very useful, it requires a regular grid.
>
> Best wishes,
>
> Fabian
> _______________________________________________
> SciPy-Dev mailing list -- [email protected]
> To unsubscribe send an email to [email protected]
> https://mail.python.org/mailman3/lists/scipy-dev.python.org/
> Member address: [email protected]
_______________________________________________
SciPy-Dev mailing list -- [email protected]
To unsubscribe send an email to [email protected]
https://mail.python.org/mailman3/lists/scipy-dev.python.org/
Member address: [email protected]