Cubic splines on a grid

Question

I trying to work out how to interpolate on a grid with cubic splines. Let the point at which I'm trying to interpolate be at {xp,yp}. At the moment I am fitting splines across the rows and then interpolating along these splines at xp. I then have a set of interpolated values (one for each row) which I can fit another vertical spline through and evaluate at yp. Now my trouble is that each such evaluation requires another spline to be fitted (I can save the horizontal splines). Is there a quicker way to do this that would ensure C^2? Some way of averaging spline coefficients? Please help.

Answer 1

The usual way to get good performance is is to store enough information about the spline, so that evaluation can be a simple local calculation. You certainly don't want to recalculate the spline coefficients every time you want to calculate a point.

So, in the curve case, the spline calculation gives you the first derivative at each data point, and you store these derivatives. Then, when you need to calculate a point on the spline, it's easy: you figure out which segment it lies on, you retrieve the points and first derivatives at the end-points of this segment, and do a simple Hermite cubic interpolation.

Alternatively, you can compute and store the b-spline coefficients of the curve, and use these to evaluate points. In the case of $C_2$ curves, the b-spline coeffients will occupy roughly 50% less space (compared to the "Hermite" approach of storing points and first derivatives).

Now, on to the surface case ...

The principle is the same -- after computing the surface, store enough information to make subsequent evaluations easy. There is a direct analog of the Hermite scheme for curves outlined above. At each data point $(x_i, y_j)$ you compute and store $$ S(x_i, y_j) \quad ; \quad \frac{\partial S}{\partial x}(x_i, y_j) \quad ; \quad \frac{\partial S}{\partial y}(x_i, y_j) \quad ; \quad \frac{\partial^2 S}{\partial x \partial y}(x_i, y_j) $$ The computation is done by separately "splining" in the $x$ and $y$ directions. It sounds like you know how to do this, so I won't go into details unless you ask. Then, to compute the surface value at given coordinates $(x,y)$, you proceed as follows. First, find out which "patch" the given point lies in. In other words, find indices $i$ and $j$ such that $x_i \le x \le x_{i+1}$ and $y_j \le y \le y_{j+1}$. Then retrieve the four pieces of stored data at the four corners $(x_i,y_j)$, $(x_i,y_{j+1})$, $(x_{i+1},y_j)$, $(x_{i+1},y_{j+1})$ of this patch. You now have 16 pieces of data that will allow you to do a simple Hermite bicubic patch calculation to get $S(x,y)$, as explained on page 9 of these notes.

As in the curve case, you can compute and store b-spline coefficients rather than partial derivatives, and this will give you a roughly 4x reduction in data.