.help INTERP Jan83 "Image Interpolation Package" .sh Introduction One of the most common operations in image processing is interpolation. Due to the very large amount of data involved, efficiency is quite important, and the general purpose routines available from many sources cannot be used. The task of interpolating image data is simplified by the following considerations: .ls 4 .ls o The pixels (data points) are equally spaced along a line or on a rectangular grid. .le .ls o There is no need for coordinate transformations. The coordinates of a pixel or point to be interpolated are the same as the subscript of the pixel in the data array. The coordinate of the first pixel in the array may be assumed to be 1.0, and the spacing between pixels is also 1.0. .le .le The task of interpolating image data is complicated by the following considerations: .ls .ls o We would like the "quality" of the interpolator to be something that can be set at run time by the user. This is done by permitting the user to select the type of interpolator to be used, and the order of interpolation. The following interpolators will be provided: .nf - nearest neighbor - linear - natural cubic spline - divided differences, either 3rd or 5th order. .fi .le .ls o The interpolator must be able to deal with indefinite valued pixels. In the IRAF, an indefinite valued pixel has the special value INDEF. If the interpolator cannot return a meaningful interpolated value (due to the presence of indefinite pixels in the data, or because the coordinates of the point to be interpolated are out of bounds) the value INDEF should be returned instead. .le .le .sh Sequential Interpolation of One Dimensional Arrays Ideally, we should be able to define a single set of procedures which can perform interpolation by any of the techniques mentioned above. The procedures for interpolating one dimensional arrays are optimized for the case where the entire array is to be sampled, probably in a sequential fashion. The following calls should suffice. The package prefix is "asi", meaning array sequential interpolate, and the last three letters of each procedure name are used to describe the action performed by the procedure. The interpolator is restricted to single precision real data, so there is no need for a type suffix. .ks .nf asiset (coeff, interpolator_type, boundary_type) asifit (data, npts, coeff) y = asival (x, coeff) # evaluate y(x) asider (x, derivs, nderiv, coeff) # derivatives v = asigrl (x1, x2, coeff) # take integral .fi .ke Here, COEFF is a structure used to describe the interpolant, and the interpolator type is chosen from the following set of options: .ks .nf NEAREST return value of nearest neighbor LINEAR linear interpolation DIVDIF3 third order divided differences DIVDIF5 fifth order divided differences SPLINE3 cubic spline .fi .ke When the VAL procedure is called to evaluate the interpolant at a point X, the code must check to see if X is out of bounds. If so, the type of value returned may be specified by the third and final argument to the SET procedure. .ks .nf B_NEAREST return nearest boundary (nb) value B_INDEF (default; return indefinite) B_REFLECT interpolate at abs(x) B_WRAP wrap around to opposite boundary B_PROJECT return y(nb) - (y(abs(x)) - y(nb)) .fi .ke For example, the following subset preprocessor code fragment uses the array interpolation routines to shift a data array by a constant amount: .ks .nf real shift, coeff[NPIX+SZ_ASIHDR] real inrow[NPIX], outrow[NPIX], asival() int i, interpolator_type, boundary_type, clgeti() begin # get parameters from CL, setup interpolator interpolator_type = clgeti ("interpolator") boundary_type = clgeti ("boundary_type") call asiset (coeff, interpolator_type, boundary_type) [code to read pixels into inrow] call asifit (inrow, NPIX, coeff) do i = 1, NPIX # do the interpolation outrow[i] = asival (i + shift, coeff) .fi .ke The array COEFF is actually a simple structure. An initial call to ASISET is required to save the interpolator parameters in the COEFF structure. ASIFIT computes the coefficients of the interpolator curve, leaving the coefficients in the appropriate part of the COEFF structure. Thereafter, COEFF is read only. ASIVAL evaluates the zeroth derivative of the curve at the given X coordinate. ASIDER evaluates the first NDERIV derivatives at X, writing the values of these derivatives into the DERIVS output array. ASIGRL integrates the interpolant over the indicated range (returning INDEF if there are any indefinite pixels in that range). .ks .nf element usage 1 type of interpolator 2 type of boundary condition 3 number of pixels in x 4 number of pixels in y 5-9 reserved for future use 10-N coefficients The COEFF structure .fi .ke In the case of the nearest neighbor and linear interpolation algorithms, the coefficient array will contain the same data as the input data array. The cubic spline algorithm requires space for NPIX+2 b-spline coefficients. The space requirements of the divided differences algorithm are similar. The constant SZ_ASIHDR should allow enough space for the worst case. Note that both 3rd and 5th order divided differences interpolators are provided, but only the 3rd order (cubic) spline. This is because the 5th order spline is considerably harder than the cubic spline. We can easily add a 5th order spline, or indeed a completely new algorithm such as the sinc function interpolator, without changing the calling sequences or data structures. .sh Random Interpolation of One Dimensional Arrays If we wish only to interpolate a small region of an array, or a few points scattered along the array at random, or if memory space is more precious than execution speed, then a different sort of interpolator is desirable. We can dispense with the coefficient array in this case, and operate directly on the data array. .ks .nf ariset (interpolator_type, boundary_type) arifit (data, npts) y = arival (x, data) arider (x, derivs, nderiv, data) .fi .ke Since the number of points required to interpolate at a given X value is fixed, these routines are internally quite different than the sequential procedures. .sh Sequential Interpolation of Two Dimensional Arrays None of the interpolation algorithms mentioned above are difficult to generalize to two dimensions. The calling sequences should be kept as similar as possible. The sequential two dimensional procedures are limited by the size of the two dimensional array which can be kept in memory, and therefore are most useful for subrasters. The COEFF array should be declared as a one dimensional array, even though it used to store a two dimensional array of coefficients. .ks .nf msiset (coeff, interpolator_type), boundary_type) msifit (data, nx, ny, coeff) y = msival (x, y, coeff) msider (x, y, derivs, nderiv, coeff) .fi .ke Note that the coefficients of the bicubic spline are obtained by first doing a one dimensional spline fit to the rows of the data array, leaving the resultant coefficients in the rows of the coeff array, followed by a series of one dimensional fits to the coefficients in the coeff array. The row coefficients are overwritten by the coefficients of the tensor product spline. Since all the work is done by the one dimensional spline fitting routine, little extra code is required to handle the two dimensional case. The bicubic spline is evaluated by summing the product C(i,j)*Bi(x)*Bj(y) at each of the sixteen coefficients contributing to the point (x,y). Once again, the one dimensional routines for evaluating the b-spline (the functions Bi(x) and Bj(y)) do most of the work. Although it is not required by very many applications, it would also be useful to be able to compute the surface integral of the interpolant over an arbitrary surface (this is useful for surface photometry applications). Never having done this, I do not have any recommendations on how to set up the calling sequence. .sh Random Interpolation of Two Dimensional Arrays The random interpolation procedures are particularly useful for two dimensional arrays due to the increased size of the arrays. Also, if an application generally finds linear interpolation to be sufficient, and rarely uses a more expensive interpolator, the random routines will be more efficient overall. .ks .nf mriset (interpolator_type, boundary_type) msifit (data, nx, ny) y = mrival (x, y, data) mrider (x, y, derivs, nderiv, data) .fi .ke .sh Cardinal B-Splines There are many different kinds of splines. One of the best for the case where the data points are equally spaced is the cardinal spline. The so called "natural" end point conditions are probably the best one can do with noisy data. Other end point conditions, such as not-a-knot, may be best for approximating analytic functions, but they are unreliable with noisy data. The natural cubic cardinal spline is fitted by solving the following tridiagonal system of equations (Prenter, 1975). .ks .nf / : \ | 6 -12 6 : 0 | | 1 4 1 : y(1) | | 1 4 1 : y(2) | | ... : | | 1 4 1 : y(n-1) | | 1 4 1 : y(n) | | 6 -12 6 : 0 | \ : / .fi .ke This matrix can most efficiently be solved by hardwiring the appropriate recursion relations into a procedure. The problem of what to do when some of the Y(i) are indefinite has not yet been solved. The b-spline basis function used above is defined by .ks .nf B(x) = (x-x1)**3 x:[x1,x2] B(x) = 1 + 3(x-x2) + 3(x-x2)**2 - 3(x-x2)**3 x:[x2,x3] B(x) = 1 + 3(x3-x) + 3(x3-x)**2 - 3(x3-x)**3 x:[x3,x4] B(x) = (x4-x)**3 x:[x4,x5] .fi .ke where X1 through X4 are the X coordinates of the four b-spline coefficients contributing to the b-spline function at X.