| Biomedical Imaging Group |
| School of Engineering |
 |
|
|
This distribution is dated June 19, 2008. It includes the complete set of source files, along with one data file.
This C program is based on the following paper:
The purpose of this program is to be a practical and didactic introduction on how to perform shifted-linear 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.
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, it so happens that the specific basis function g() that is considered for shifted-linear interpolation takes a non-zero value for g(1) = ½ (1 - √3 ⁄ 3), and also differs from unity at the origin since g(0) = ½ (1 + √3 ⁄ 3). Thus, one must depart from traditional interpolation when performing shifted-linear interpolation of data. 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 recursive filter (for more details, see 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 of space (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),
where the notations should be obvious. The only difficulty is that we now have to deal with a bidimensional array of coefficients c[i, j]. Fortunately, these coefficients can be obtained by processing to the bidimensional set of samples in one spatial direction first, and then by processing the remaining perpendicular direction. This separable approach is the most efficient way to process data with more than one dimension. Provided that the set of coefficients c[i, j] has already been obtained from coeff.c, the role of the part of the program available in the file interpol.c is to carry out the computation of the tensor product and of the function g(). More specifically, this latter is a translated hat function.
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. For shifted-linear interpolation, we advocate an extrapolation scheme where coefficients are given a constant value on one side of the signal. (Since the production of the coefficients c[] proceeds by a causal filter, what happens on the other side of the signal is not arbitrary but is determined by the recursion.) This scheme ensures that only weak discontinuities (or strong changes of value) of ƒ() will 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 how the program has masked out the extrapolated data.
 |  |
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.