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/splin... (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 -- scipy-dev@python.org To unsubscribe send an email to scipy-dev-leave@python.org https://mail.python.org/mailman3/lists/scipy-dev.python.org/ Member address: mail@nilsgeib.de