diff options
Diffstat (limited to 'math/interp/notes3')
| -rw-r--r-- | math/interp/notes3 | 216 |
1 files changed, 216 insertions, 0 deletions
diff --git a/math/interp/notes3 b/math/interp/notes3 new file mode 100644 index 00000000..e7121c99 --- /dev/null +++ b/math/interp/notes3 @@ -0,0 +1,216 @@ + 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 + <nterm> terms to connect the <nterm> 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. |