Notes on the Interpolator Package for IRAF Version 2 Some parts of the image interpolator package are available. This is a reminder of changes to the document TRP prepared by Doug Tody. Also it is a place to keep track of why things were done the way they were. The individual procedures (subroutines) will be discussed; making it easy to add, subtract and change things. ************************************************************************** ** ** ** ** GENERAL NOTES FOR 1-D PACKAGE: 1. Handling of errors. Some routines call the error routines of the IRAF package. Logically it is reasonable to trap silly errors at certain points. For instance if the interpolator type is set to a nonexistant type an error call is generated. To the programmer, the two most important error features are: 1. All pixels are assumed to be good. 2. All x references are assumed to land in the closed interval 1 <= x <= NPTS. Point 2 probably needs some explanation, because at first hand it seems rather easy to test for x < 1 || x > NPTS and take appropriate action. There are two objections to including this test. First it is often possible to generate x values in such a way that they are guaranteed to lie in the given interval. Testing again within the interpolator package is a waste of computer time. Second the interpolator package becomes too large if code is included to handle out of bounds references and such code duplicates other parts of the IRAF package that handles out of bound values. The alternative is to return INDEF when x is out of bounds. This is not very appealing for many routine applications such as picture shifting. 2. Size of "coeff" array: coeff is a real 1-d array of size -- 2 * NPTS + SZ_ASI where SZ_ASI is currently 30. It, like the other defines needed for the interp package, is available in the file -- terpdef.h . This coeff array serves as work area for asifit, a structure to save interpolator type, flags etc., an array to save data points or b-spline coefficients as needed. The polynomial interpolators are saved as the data points and the interpolation is done "on the fly". Rather than compute polynomial coefficients before-hand, the interpolated value is calculated from the data values at time of request. (The actual calculation is done by adding and multiplying things in an order suggested by Everett's interpolation formula rather than the polynomial coefficient form. The resulting values are the same up to differences due to numerical roundoff. As a practical matter, all forms work except polynomials shifted far from the point of interest. Polynomials shifted to a center far away suffer from bad roundoff errors.) The splines must be calculated at a single pass because (formally) they are global (i.e use all of the data points). The b-splines (b for basis) are just certain piecewise polynomials for solving the interpo- lation problem, it is convenient to calculate the coefficients of these b-splines and to save these in the array. 3. Method of solution: ************************************************************************** ** ** ** ** ** ** PROCEDURES: *************************************************************************** ***** asiset(coeff,interptype) ***** The interpolator type can be set to: IT_NEAREST nearest neighbor IT_LINEAR linear interpolation IT_POLY3 interior polynomial 3rd order IT_POLY5 interior polynomial 5th order IT_SPLINE3 cubic natural spline Subroutines needed: NONE ************************************************************************* ******* asifit (data, npts, coeff) ***** This fits the interpolator to the data. Takes care of bad points by replacing them with interpolated values to generate a uniformly spaced array. For polynomial interpolators, this array is sent to the evaluation routines. For the spline interpolator, the uniform array is fitted and the basis-spline coefficients are saved. Interior polynomial near boundary: Array is extended according to choice of boundary conditions. Out of bounds values and values near the edge depend on choice of boundary conditions. Spline near boundary: Natural (second derivative equals zero at edge) spline is fitted to data only. The resulting function is reflected, wrapped etc. to generate out of bounds values. Out of bounds values depend on choice of boundary conditions but values within the array near the edges do not. For the wrap boundary condition, the one segment connecting point number n to point 1 is not generated by the spline fit. Instead the polynomial that is continuous and has a continuous first and second derivative n is inserted. In this one case, the first and second derivative are not continuous at 1. Replacing bad points: The same kind of interpolator is used to replace bad points as that used to generate interpolated values. The boundary conditions are not handled in the same manner. Bad points near the edge are replaced by using lower order polynomial interpolators if there are not enough points on both sides of the bad point. Bad points that are not straddled by good points are assigned the value of the nearest good point. This is done in the temporary array before the uniformly spaced interpolator is applied. subroutines needed: real iipint(x,y,nterm,x0) Returns interpolated value at x0 using polynomial with terms to connect the data pairs in x,y. cubspl(tau,c,n,ibcbeg,ibcend) Unmodified cubic spline interpolator for non-uniformly spaced data taken from C. de Boor's book: A Practical Guide to Splines, (1978 Springer-Verlag). iifits(b,d,n) Fits cubic spline to uniformly spaced array of y values using a basis-spline representation for the spline. ***************************************************************************** ***** y = asival (x, coeff) ***** Subroutine to return interpolated value after fitting. Subroutines needed: real iivalu (x, coeff) Returns interpolated value assuming x is in array. Also has no code for bad point propagation. *************************************************************************** ***** asieva (x, y, n, coeff) ***** Given an ordered array of x values returns interpolated values in y. This is needed for speed. Reduces number of subroutine calls per point, reduces the number of tests for out of bounds, and bad pixel propagation can be made more efficient. Subroutines needed: real iivalu (x, coeff) ************************************************************************** ***** asider (x, derivs, nderiv, coeff) ***** This evaluates derivatives at point x. Returns the first nderiv - 1 derivatives. derivs[1] is function value, derivs[2] is first derivative, derivs[3] is second derivative etc. There is no check that nderiv is a reasonable value. Polynomials have a limited number of non-zero derivatives, if more are asked for they come back as zero. Subroutines needed: iivalu (x, coeff) iiderv (x, derivs, nderiv, coeff) Assumes x lands in array returns derivatives. No handling of bad points. **************************************************************************** ****** integral = asigrl (a, b, coeff) ******* This integrates the array from to a to b. The boundary conditions are incorporated into the code so if function values can be returned from all points between a and b, the integral will be defined. Thus for WRAP, REFLECT and PROJECT boundary conditions, the distance from a to b can be as large as 3 times the array length if a and b are chosen appropriately. For these boundary conditions, attempts to extend beyond +/- one cycle produce INDEF. Subroutines needed: iiigrl (a, b, coeff) returns integral when a and b land in array.