Degenerate Conic

Algorithms • Modern Fortran Programming • Orbital Mechanics

Feb 25, 2015

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.