Biomedical Imaging Group Department of Microengineering

Spline Interpolation

Author: Philippe Thévenaz

• I. Download
• II. Related work
• III. Explanations
• IV. Underlying theory
  1. Interpolation coefficients
  2. Several dimensions
  3. Border conditions
• V. Examples

This distribution is dated January 5, 2009. It includes the complete set of source files, along with one data file.

• UNIX distribution [133,120 bytes]
• Macintosh distribution [83,879 bytes]
• PC distribution [71,028 bytes]

This C program is based on the following paper:

• P. Thévenaz, T. Blu, M. Unser, "Interpolation Revisited," IEEE Transactions on Medical Imaging, vol. 19, no. 7, pp. 739-758, July 2000.

The purpose of this program is to be a practical and didactic introduction on how to perform spline interpolation. It is a self-contained application that will apply a rigid-body transformation to an image (rotation and translation). The source code (ANSI-C) is divided into 4 principal files ('.c') and 3 associated header files ('.h'). A standard data file (the famous Lena image) is provided so that you can test the program as soon as it has been downloaded and compiled.

1. The file demo.c contains the topmost procedure main(). It calls the routines that perform the I/O-associated tasks, calls the routine that converts the image samples into B-spline coefficients, and cares about geometric issues.
2. The files coeff.c and coeff.h implement the algorithm that converts the image samples into B-spline coefficients. This efficient procedure essentially relies on the paper cited above; data are processed in-place. Even though this algorithm is robust with respect to quantization, we advocate the use of a floating-point format for the data.
3. The files interpol.c and interpol.h perform the bidimensional interpolation. Given an array of spline coefficients, they return the value of the underlying continuous spline model, sampled at the location (x, y). The model degree can range from 2 to 9.
4. The files io.c and io.h provide a rudimentary user interface that is supposed to work on any compiler/system/OS that is in compliance with the ANSI-C norm. The image to be processed is read from file; the file format is binary raw data, with 1 unsigned byte per pixel, organized in raster fashion. The image resulting from executing the program is also stored on file, with the same conventions. You will need to use your own software to display it. A platform-independent candidate is perhaps ImageJ, a public-domain image-processing package written in Java (use the menu File:Import:Raw... and supply the requested information about the image).
5. The data file lena.img stores a 256 x 256 image that, along the years, has become the de facto standard in image processing. The format corresponds to that described above. The visualization convention is to map low numeric values to dark gray-levels (black) and high numeric values to light gray-levels (white).

To simplify the explanations, let us examine the one-dimensional case first and let us defer to later the considerations on additional dimensions. Given a set of samples s[j] with j integer, the fundamental task of interpolation is to compute the following equation for any desired value of x:

ƒ(x) = ∑_k c[k] g(x - k),

where x is a real number and where k takes on every possible integer value. For a specific integer value j, it is also desired that the interpolated function ƒ(j) takes the same value as the sample s[j] of the set of available samples we want to interpolate

s[j] = ƒ(j) = ∑_k c[k] g(j - k).

We discourage the traditional approach one usually follows in order to satisfy this equation. Traditionally, one makes sure that two conditions are met. Firstly, set c[k] = s[k] for all k. Secondly, select a function g() that takes value 0 for all integer arguments, except at the origin where it must equal unity (for any non-integer arguments, g() is free to take any value it may like). This would do the trick, and would be the classical approach to interpolation.

Now, B-splines have very attractive properties, which makes them very good candidates in the role of the function g(). Unfortunately, they don't have all the qualities in the world; in particular, they may take a non-zero value for integer arguments, and they may differ from unity at the origin. Thus, one must depart from traditional interpolation so as to use B-splines for interpolating data. Instead, one must determine the suitable set of coefficients c[k] ≠ s[k] = ƒ(k) such that ƒ(j) = ∑ c[k] g(j - k) still holds true for any integer j. This time, setting c[k] = s[k] doesn't do the trick. We need instead some routine that produces a set of coefficients c[] out of a set of samples s[]. This is the role of the routine available in the file coeff.c. The core of this routine consists in a forward-backward recursive filter (for more details, see the reference [1]).

To process an image, one must consider more than a single dimension. The traditional extension--that we follow here--goes by the name of tensor product. It consists in weighting the interpolation in one direction by that of the remaining directions (this is also called a separable approach). More precisely, in two dimensions we compute

ƒ(x, y) = ∑_j ∑_i c[i, j] g(x - i) g(y - j),

So far, we have taken for granted that the number of samples s[] was infinite, same with the number of coefficients c[]. Obviously, the amount of RAM in your computer is only so much, which asks for a finite number of data. Reconciliation of these contradictory requirements is the object of the present section.

Let us examine the question in one dimension, and let us suppose you have N samples at your disposal. Given only those samples you know of, if we find a way to predict (extrapolate) the value of any coefficient outside the range of your known data, then we are back to the first (infinite) formulation, without its cost. This is known as an implicit scheme. Extrapolating the samples s[] is less relevant than extrapolating the coefficients c[], for the former do not explicitly appear in what we compute in practice, that is, ∑ c[k] g(x - k).

There is no limit to the number of ways one may select to perform extrapolation. After all, any single one is but a wild guess, even though some are more informed than others. We advocate here an extrapolation scheme where coefficients are mirrored around their extremities, namely, around c[0] and around c[N-1]. This double mirroring ensures that every c[k] takes a definite value, no matter how far k is from the interval [0, N-1]. Moreover, it imposes the same doubly-mirrored structure on the samples s[k], and also on the whole curve ƒ(x). An important benefit of this scheme is that there is no need to "invent" special values (e.g., zero); all extrapolated values are replicated from those of the known interval. Moreover, this scheme ensures that no discontinuities (or strong changes of value) of ƒ() arise at the extremities of the known interval.

The following pair of images is typical of the kind of results you can obtain with the present program. The image has been rotated by some angle around an arbitrary origin, while an additional translation has been thrown in for good measure. Please observe the influence of mirroring the data (the program allows for the masking out of the extrapolated data, if desired).

Left: 256x256 Lena input image. Right: image resulting from an arbitrary transformation (no masking)

Here is a list of additional examples that illustrate the use of our interpolation routine. You may want to compare your own results to the ones below so as to ensure that your implementation is correct. In terms of ease of check, we have ordered the examples from most easy to more difficult. The images that we provide here are subject to lossless coding.

1. Let us check that specifying an identity transformation indeed results in identity;
2. Check that identity does not depend on where the center of coordinates is;
3. Perform a simple, integer translation;
4. If we go back, the image should be partially restored. Note that a translation does not depend on where the center of coordinates is;
5. We can also perform non-integer translations. In addition, we are now unraveling the part of the image that was previously hidden by masking. The resulting image clearly shows the consequence of using mirror boundaries;
6. A translation by a non-trivial (dx, dy) can still result in identity, provided {dx mod (2 Nx - 2)} = 0 and {dy mod (2 Ny - 2)} = 0;
7. Because of mirroring, the upper-left ghost image is indeed a rotated version of the original image;
8. In fact, the ghost image above shares only one pixel with the original image, as shown by masking out the irrelevant pixels (only the upper-left corner pixel remains);
9. Let us perform a true rotation;
10. For rotations, the center of coordinates plays an important role;
11. Let us rotate the image around its exact center, by a half turn;
12. Finally, let us perform a more arbitrary transformation, this time with a higher spline degree.