Multidimensional B-Spline Interpolation

I just started a new modern Fortran software library called bspline-fortran, which is for multidimensional (multivariate) b-spline interpolation of data defined on a regular grid. It is available on GitHub, and released under a permissive BSD-style license.

It seems impossible to find code for higher than 3D spline interpolation on the internet. Einspline only has 1D-3D, as do the NIST Core Math Library DBSPLIN and DTENSBS routines. Python’s SciPy stops at 2D (Bivariate splines). At first glance, the SciPy routine map_coordinates seems to do what I want, but not really since there doesn’t seem to be any way to specify the abscissa vectors (x, y, z, etc.).  Note: maybe there is a way, but looking at the less-than-great documentation, the minimally-commented C code, and then reading several confusing Stack Overflow posts didn’t help me much. Apparently, I’m not the only one confused by this routine.

After much searching, I finally just took the 2D and 3D NIST routines (written in 1982), refactored them significantly so I could better understand what they were doing, and then extended the algorithm into higher dimensions (which is actually quite easy to do for b-splines). The original FORTRAN 77 code was somewhat hard to follow, since it was stuffing matrices into a giant WORK vector, which required a lot of obscure bookkeeping of indices. For example, I replaced this bit of code:

   IZM1 = IZ - 1
   KCOLY = LEFTY - KY + 1
   DO 60 K=1,KZ
      I = (K-1)*KY + 1
      J = IZM1 + K
      WORK(J) = DBVALU(TY(KCOLY),WORK(I),KY,KY,IDY,YVAL,INBV1,WORK(IW))
60 CONTINUE

with this:

    kcoly = lefty - ky + 1
    do k=1,kz
        temp2(k) = dbvalu(ty(kcoly:),temp1(:,k),ky,ky,idy,yval,inbv1,work)
    end do

which makes it a lot easier to deal with and add the extra dimensions as necessary.

The new library contains routines for 2D, 3D, 4D, 5D, and 6D interpolation. It seems like someone else may find them useful, I don’t know why they don’t seem to be available anywhere else. Eventually, I’ll add object-oriented wrappers for the core routines, to make them easier to use.

Posted in Programming Tagged with: , , , ,
3 comments on “Multidimensional B-Spline Interpolation
  1. Stefano says:

    Good work!

    Very helpful, I will surely read it!

    Thank you for sharing it.

  2. Manuel says:

    Very powerful! How hard would it be to extend your code to unstructured interpolation points?

  3. Jacob says:

    So, scattered data? I’m not sure. I’ve never needed to interpolate scattered data before and am not that familiar with the methods for doing it.

Leave a Reply

Your email address will not be published. Required fields are marked *

*