From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001 From: Joe Hunkeler Date: Tue, 11 Aug 2015 16:51:37 -0400 Subject: Repatch (from linux) of OSX IRAF --- math/iminterp/doc/arbpix.hlp | 57 +++++ math/iminterp/doc/arider.hlp | 59 +++++ math/iminterp/doc/arieval.hlp | 48 ++++ math/iminterp/doc/asider.hlp | 52 ++++ math/iminterp/doc/asieval.hlp | 44 ++++ math/iminterp/doc/asifit.hlp | 40 +++ math/iminterp/doc/asifree.hlp | 25 ++ math/iminterp/doc/asigeti.hlp | 36 +++ math/iminterp/doc/asigetr.hlp | 36 +++ math/iminterp/doc/asigrl.hlp | 40 +++ math/iminterp/doc/asiinit.hlp | 39 +++ math/iminterp/doc/asirestore.hlp | 36 +++ math/iminterp/doc/asisave.hlp | 39 +++ math/iminterp/doc/asisinit.hlp | 60 +++++ math/iminterp/doc/asitype.hlp | 95 +++++++ math/iminterp/doc/asivector.hlp | 52 ++++ math/iminterp/doc/im1dinterp.spc | 525 +++++++++++++++++++++++++++++++++++++++ math/iminterp/doc/im2dinterp.spc | 432 ++++++++++++++++++++++++++++++++ math/iminterp/doc/iminterp.hd | 37 +++ math/iminterp/doc/iminterp.hlp | 234 +++++++++++++++++ math/iminterp/doc/iminterp.men | 32 +++ math/iminterp/doc/iminterp.spc | 525 +++++++++++++++++++++++++++++++++++++++ math/iminterp/doc/mrider.hlp | 79 ++++++ math/iminterp/doc/mrieval.hlp | 57 +++++ math/iminterp/doc/msider.hlp | 52 ++++ math/iminterp/doc/msieval.hlp | 46 ++++ math/iminterp/doc/msifit.hlp | 45 ++++ math/iminterp/doc/msifree.hlp | 26 ++ math/iminterp/doc/msigeti.hlp | 35 +++ math/iminterp/doc/msigetr.hlp | 37 +++ math/iminterp/doc/msigrid.hlp | 51 ++++ math/iminterp/doc/msigrl.hlp | 43 ++++ math/iminterp/doc/msiinit.hlp | 41 +++ math/iminterp/doc/msirestore.hlp | 36 +++ math/iminterp/doc/msisave.hlp | 38 +++ math/iminterp/doc/msisinit.hlp | 61 +++++ math/iminterp/doc/msisqgrl.hlp | 38 +++ math/iminterp/doc/msitype.hlp | 95 +++++++ math/iminterp/doc/msivector.hlp | 54 ++++ 39 files changed, 3377 insertions(+) create mode 100644 math/iminterp/doc/arbpix.hlp create mode 100644 math/iminterp/doc/arider.hlp create mode 100644 math/iminterp/doc/arieval.hlp create mode 100644 math/iminterp/doc/asider.hlp create mode 100644 math/iminterp/doc/asieval.hlp create mode 100644 math/iminterp/doc/asifit.hlp create mode 100644 math/iminterp/doc/asifree.hlp create mode 100644 math/iminterp/doc/asigeti.hlp create mode 100644 math/iminterp/doc/asigetr.hlp create mode 100644 math/iminterp/doc/asigrl.hlp create mode 100644 math/iminterp/doc/asiinit.hlp create mode 100644 math/iminterp/doc/asirestore.hlp create mode 100644 math/iminterp/doc/asisave.hlp create mode 100644 math/iminterp/doc/asisinit.hlp create mode 100644 math/iminterp/doc/asitype.hlp create mode 100644 math/iminterp/doc/asivector.hlp create mode 100644 math/iminterp/doc/im1dinterp.spc create mode 100644 math/iminterp/doc/im2dinterp.spc create mode 100644 math/iminterp/doc/iminterp.hd create mode 100644 math/iminterp/doc/iminterp.hlp create mode 100644 math/iminterp/doc/iminterp.men create mode 100644 math/iminterp/doc/iminterp.spc create mode 100644 math/iminterp/doc/mrider.hlp create mode 100644 math/iminterp/doc/mrieval.hlp create mode 100644 math/iminterp/doc/msider.hlp create mode 100644 math/iminterp/doc/msieval.hlp create mode 100644 math/iminterp/doc/msifit.hlp create mode 100644 math/iminterp/doc/msifree.hlp create mode 100644 math/iminterp/doc/msigeti.hlp create mode 100644 math/iminterp/doc/msigetr.hlp create mode 100644 math/iminterp/doc/msigrid.hlp create mode 100644 math/iminterp/doc/msigrl.hlp create mode 100644 math/iminterp/doc/msiinit.hlp create mode 100644 math/iminterp/doc/msirestore.hlp create mode 100644 math/iminterp/doc/msisave.hlp create mode 100644 math/iminterp/doc/msisinit.hlp create mode 100644 math/iminterp/doc/msisqgrl.hlp create mode 100644 math/iminterp/doc/msitype.hlp create mode 100644 math/iminterp/doc/msivector.hlp (limited to 'math/iminterp/doc') diff --git a/math/iminterp/doc/arbpix.hlp b/math/iminterp/doc/arbpix.hlp new file mode 100644 index 00000000..0d7ec9ec --- /dev/null +++ b/math/iminterp/doc/arbpix.hlp @@ -0,0 +1,57 @@ +.help arbpix Dec98 "Image Interpolator Package" +.ih +NAME +arbpix -- replace INDEF valued pixels with interpolated values +.ih +SYNOPSIS +include + +arbpix (datain, dataout, npix, interp_type, boundary_type) + +.nf + real datain[npix] #I input data + real dataout[npix] #O output array, dataout != datain + int npix #I number of data points + int interp_type #I type of interpolant + int boundary_type #I type of boundary condition +.fi +.ih +ARGUMENTS +.ls datain +Array of input data containing 0 or more INDEF valued pixels. +.le +.ls dataout +Array of output data with INDEFS replaced by interpolated values. +The dataout array must be different from the datain array. +.le +.ls npix +Number of data points. +.le +.ls interp_type +Type of interpolant. Options are II_NEAREST, II_LINEAR, II_POLY3, II_POLY5, +II_SPLINE3, II_SINC / II_LSINC, and II_DRIZZLE. The look-up table sinc +interpolant is not supported, and defaults to the sinc interpolant. +The sinc interpolant width is 31 pixels. The drizzle interpolant is not +supported and defaults to the linear interpolant. The interpolant type +definitions are stored in the file math/iminterp.h. +.le +.ls boundary_type +Type of boundary extension. The only supported option is II_BOUNDARYEXT. +Polynomial interpolants of lower order are used if there are not enough +good pixels to define the requested interpolant. Nearest neighbor boundary +extension is used if there are not enough good points to define the sinc +interpolant. The boundary type definitions are stored in the header file +math/iminterp.h. +.le +.ih +DESCRIPTION +If there are no good points in datain, ARBPIX returns INDEFS in dataout. +Points below and above the first and last good point are replaced by the +first and last good point values respectively. +.ih +NOTES +The sinc function actually evaluates the interpolant by computing the +average of two interpolations at +-0.05 pixels about the bad pixel since +the interpolant is undefined exactly at a pixel. +.ih +SEE ALSO diff --git a/math/iminterp/doc/arider.hlp b/math/iminterp/doc/arider.hlp new file mode 100644 index 00000000..b8631eaa --- /dev/null +++ b/math/iminterp/doc/arider.hlp @@ -0,0 +1,59 @@ +.help arider Dec98 "Image Interpolator Package" +.ih +NAME +arider -- calculate the interpolant derivatives at x +.ih +SYNOPSIS +include + +arider (x, datain, npix, der, nder, interp_type) + +.nf + real x[2] #I x value, 1 <= x[1-2] <= npts + real datain[npix] #I array of data points + int npix #I number of data points + real der[nder] #O derivatives, der[1] = function value + int nder #I number of derivatives, 1 + max order + int interp_type #I interpolant type +.fi +.ih +ARGUMENTS +.ls x +Single X value, or pair of X values defining a range in the case of the +drizzle interpolant. +.le +.ls datain +Array of data values. +.le +.ls npix +Number of data points. +.le +.ls der +Array of derivatives. Der[1] contains the function value, der[2] the +first derivative, and so on. +.le +.ls nder +Number of derivatives. ARIDER checks that the requested number of derivatives +is sensible. The sinc interpolant returns the function value and the first +two derivatives. The drizzle interpolant returns the function and the first +derivative. +.le +.ls interp_type +Interpolant type. The options are II_NEAREST, II_LINEAR, II_POLY3, II_POLY5, +II_SPLINE3, II_SINC / II_LSINC, and II_DRIZZLE. The look-up table sinc +is not supported and defaults to sinc. The sinc interpolant width is 31 pixels. +The drizzle pixel fraction is 1.0. The interpolant type definitions are found +in the package header file math/iminterp.h. +.le +.ih +DESCRIPTION +ARIDER permits the evaluation of the interpolant at a few randomly spaced +points within datain without the storage requirements of the sequential +version. +.ih +NOTES +Checking for INDEF valued or out of bounds pixels is the responsibility +of the user. +.ih +SEE ALSO +asider diff --git a/math/iminterp/doc/arieval.hlp b/math/iminterp/doc/arieval.hlp new file mode 100644 index 00000000..52c62148 --- /dev/null +++ b/math/iminterp/doc/arieval.hlp @@ -0,0 +1,48 @@ +.help arieval Dec98 "Image Interpolator Package" +.ih +NAME +arieval -- evaluate the interpolant at x +.ih +SYNOPSIS +include + +y = arieval (x, datain, npix, interp_type) + +.nf + real x[2] #I x value, 1 <= x[1-2] <= npix + real datain[npix] #I data values + int npix #I number of data values + int interp_type #I interpolant type +.fi +.ih +ARGUMENTS +.ls x +Single X value, or a pair of X values specifying a range in the case +of the drizzle interpolant. +.le +.ls datain +Array of input data. +.le +.ls npix +Number of data points. +.le +.ls interp_type +Interpolant type. Options are II_NEAREST, II_LINEAR, II_POLY3, II_POLY5, +II_SPLINE3, II_SINC / II_LSINC, and II_DRIZZLE, for nearest neighbor, +linear, 3rd and fifth order polynomials, cubic spline, sinc, look-up +table sinc, and drizzle interpolants respectively. The look-up table sinc +interpolant is not supported and defaults to the sinc interpolant. The sinc +width is 31 pixels. The drizzle pixel fraction is 1.0. The interpolant +type definitions are contained in the package header file math/iminterp.h +.le +.ih +DESCRIPTION +ARIEVAL allows the evaluation of a few interpolated points without the +storage required for the sequential interpolant. +.ih +NOTES +Checking for out of bounds and INDEF valued pixels is the responsibility of +the user. +.ih +SEE ALSO +arider, asieval, asivector diff --git a/math/iminterp/doc/asider.hlp b/math/iminterp/doc/asider.hlp new file mode 100644 index 00000000..0c27ffbc --- /dev/null +++ b/math/iminterp/doc/asider.hlp @@ -0,0 +1,52 @@ +.help asider Dec98 "Image Interpolator Package" +.ih +NAME +asider -- evaluate the interpolant derivatives at x +.ih +SYNOPSIS +asider (asi, x, der, nder) + +.nf + pointer asi #I interpolant descriptor + real x[2] #I x value, 1 <= x[1-2] <= npix + real der[] #O der[1] = interpolant, der[2] = 1st derivative + int nder #I number of derivatives +.fi +.ih +ARGUMENTS +.ls asi +Pointer to the sequential interpolant descriptor. +.le +.ls x +Single X value, or pair of X values defining a range in the case of +the drizzle interpolant. +.le +.ls der +Array containing the derivatives. Der[1] = interpolant at x, der[2] the +first derivative of the interpolant at x and so on. +.le +.ls nder +Number of derivatives. Nder = 1 + order of the maximum desired derivative. +ASIDER checks that nder is reasonable. The sinc interpolant returns the +interpolant value and first two derivatives. The drizzle interpolant returns +the interpolant value and the first derivative. +.le +.ih +DESCRIPTION +The polynomial coefficients are evaluated directly from the data points +for the polynomial interpolants and from the B-spline coefficients +for the cubic spline interpolant. The derivatives are evaluated from +the polynomial coefficients using nested multiplication. The sinc +derivatives are analytic but are defined only for the first two derivatives. +The drizzle derivative is an approximation defined for the first derivative +only. +.ih +NOTES +ASIDER checks that the number of derivatives requested is reasonable. +Checking for out of bounds and INDEF valued pixels is the responsibility +of the user. ASIINIT or ASISINIT and ASIFIT must be called before ASIDER +is called. +.ih +SEE ALSO +asieval, asivector, arieval, arider +.endhelp diff --git a/math/iminterp/doc/asieval.hlp b/math/iminterp/doc/asieval.hlp new file mode 100644 index 00000000..20f70abe --- /dev/null +++ b/math/iminterp/doc/asieval.hlp @@ -0,0 +1,44 @@ +.help asieval Dec98 "Image Interpolator Package" +.ih +NAME +asieval -- procedure to evaluate interpolant at x +.ih +SYNOPSIS +y = asieval (asi, x) + +.nf + pointer asi #I interpolant descriptor + real x[2] #I x value, 1 <= x[1-2] <= npts +.fi +.ih +ARGUMENTS +.ls asi +Pointer to the sequential interpolant descriptor structure. +.le +.ls x +Single X value, or pair of X values defining a range in the case of the +drizzle interpolant. +.le +.ih +DESCRIPTION +The polynomial coefficients are calculated directly from the data points +for the polynomial interpolants, and from the B-spline coefficients for +the cubic spline interpolant. The actual calculation is done by adding and +multiplying terms according to Everett's central difference interpolation +formula. The boundary extension algorithm is projection. + +The sinc interpolant is computed using a range of data points around +the desired position. Look-up table sinc interpolation is computed +using the most appropriate entry in a precomputed look-up table. +The boundary extension algorithm is nearest neighbor. + +The drizzle interpolant is computed by summing the data over the user +supplied X interval. +.ih +NOTES +Checking for out of bounds and INDEF valued pixels is the responsibility of +the user. ASIINIT or ASISINIT and ASIFIT must be called before using ASIEVAL. +.ih +SEE ALSO +asivector, arieval +.endhelp diff --git a/math/iminterp/doc/asifit.hlp b/math/iminterp/doc/asifit.hlp new file mode 100644 index 00000000..bcd1fdc8 --- /dev/null +++ b/math/iminterp/doc/asifit.hlp @@ -0,0 +1,40 @@ +.help asifit Dec98 "Image Interpolator Package" +.ih +NAME +asifit - fit the interpolant to data +.ih +SYNOPSIS +asifit (asi, datain, npix) + +.nf + pointer asi #I interpolant descriptor + real datain[npix] #I input data + int npix #I the number of data points +.fi +.ih +ARGUMENTS +.ls asi +Pointer to sequential interpolant descriptor structure. +.le +.ls datain +Array of input data. +.le +.ls npix +Number of data points. +.le +.ih +DESCRIPTION +The datain array is checked for size, memory is allocated for the coefficient +array, and the end conditions are specified. The interior polynomial, sinc and +drizzle interpolants are saved as the data points. The polynomial coefficients +are calculated directly from the data points in the evaluation stage. The +B-spline coefficients are calculated in ASIFIT as they depend on the entire +data array. +.ih +NOTES +Checking for INDEF valued and out of bounds pixels is the responsibility +of the user. ASIINIT or ASISINIT and ASIFIT must be called before using +ASIEVAL, ASIVECTOR, ASIDER or ASIGRL. +.ih +SEE ALSO +.endhelp diff --git a/math/iminterp/doc/asifree.hlp b/math/iminterp/doc/asifree.hlp new file mode 100644 index 00000000..c61f2ce0 --- /dev/null +++ b/math/iminterp/doc/asifree.hlp @@ -0,0 +1,25 @@ +.help asifree Dec98 "Image Interpolator Package" +.ih +NAME +asifree - free sequential interpolant descriptor +.ih +SYNOPSIS +asifree (asi) + +.nf +pointer asi #U interpolant descriptor +.fi +.ih +ARGUMENTS +.ls asi +Pointer to the sequential interpolant descriptor structure. +.le +.ih +DESCRIPTION +ASIFREE frees the sequential interpolant descriptor structure allocated by +ASIINIT or ASISINIT. ASIFREE should be called when interpolation is complete. +.ih +NOTES +.ih +SEE ALSO +asiinit, asisinit diff --git a/math/iminterp/doc/asigeti.hlp b/math/iminterp/doc/asigeti.hlp new file mode 100644 index 00000000..0e9d1964 --- /dev/null +++ b/math/iminterp/doc/asigeti.hlp @@ -0,0 +1,36 @@ +.help asigeti Dec98 asigeti.hlp +.ih +NAME +asigeti -- fetch an asi integer parameter +.ih +SYNOPSIS +include + +value = asigeti (asi, param) + +.nf + pointer asi #I interpolant descriptor + int param #I parameter +.fi +.ih +ARGUMENTS +.ls asi +Pointer to the sequential interpolant descriptor structure. +.le +.ls param +The parameter to be fetched. The choices are: II_ASITYPE the interpolant +type, II_ASINSAVE the length of the saved coefficient array, and +II_ASINSINC the half-width of the sinc interpolant. The parameter +definitions are contained in the package header file math/iminterp.h. +.le +.ih +DESCRIPTION +ASIGETI is used to determine the size of the coefficient array that +must be allocated to save the sequential interpolant description structure, +and to fetch selected sequential interpolant parameters. +.ih +NOTES +.ih +SEE ALSO +asiinit, asisinit, asigetr +.endhelp diff --git a/math/iminterp/doc/asigetr.hlp b/math/iminterp/doc/asigetr.hlp new file mode 100644 index 00000000..8a31deef --- /dev/null +++ b/math/iminterp/doc/asigetr.hlp @@ -0,0 +1,36 @@ +.help asigetr Dec98 asigetr.hlp +.ih +NAME +asigetr -- fetch an asi integer parameter +.ih +SYNOPSIS +include + +value = asigetr (asi, param) + +.nf + pointer asi #I interpolant descriptor + int param #I parameter +.fi +.ih +ARGUMENTS +.ls asi +Pointer to the sequential interpolant descriptor structure. +.le +.ls param +The parameter to be fetched. The choices are: II_ASIBADVAL the undefined +pixel value for the drizzle interpolant. The parameter definitions are +contained in the package header file math/iminterp.h. +.le +.ih +DESCRIPTION +ASIGETR is used to set the value of undefined drizzle interpolant pixels. +Undefined pixels are those for which the interpolation coordinates do not +overlap the input coordinates, but are still, within the boundaries of the input +image, a situation which may occur when the pixel fraction is < 1.0. +.ih +NOTES +.ih +SEE ALSO +asiinit, asisinit, asigeti +.endhelp diff --git a/math/iminterp/doc/asigrl.hlp b/math/iminterp/doc/asigrl.hlp new file mode 100644 index 00000000..4c3087bb --- /dev/null +++ b/math/iminterp/doc/asigrl.hlp @@ -0,0 +1,40 @@ +.help asigrl Dec98 "Image Interpolator Package" +.ih +NAME +asigrl -- integrate interpolant from a to b +.ih +SYNOPSIS +integral = asigrl (asi, a, b) + +.nf + pointer asi #I interpolant descriptor + real a #I lower limit for integral, 1 <= a <= npix + real b #I upper limit for integral, 1 <= b <= npix +.fi +.ih +ARGUMENTS +.ls asi +Pointer to the sequential interpolant descriptor structure. +.le +.ls a +Lower limit to the integral, where 1 <= a <= npix. +.le +.ls b +Upper limit to the integral, where 1 <= b <= npix. +.le +.ih +DESCRIPTION +The integral is calculated assuming that the interior polynomial, sinc, and +drizzle interpolants are stored as the data points, and that the spline +interpolant is stored as an array of B-spline coefficients. + +The integral of the sinc interpolant is computed by dividing the integration +interval into a number of equal size subintervals which are at most one pixel +wide. The integral of each subinterval is the central value times the interval +width. The look-up table sinc interpolant is not supported and defaults to +the sinc interpolant. +.ih +NOTES +ASIINIT or ASISINIT and ASIFIT must be called before using ASIGRL. +.ih +SEE ALSO diff --git a/math/iminterp/doc/asiinit.hlp b/math/iminterp/doc/asiinit.hlp new file mode 100644 index 00000000..711d7969 --- /dev/null +++ b/math/iminterp/doc/asiinit.hlp @@ -0,0 +1,39 @@ +.help asiinit Dec98 "Image Interpolator Package" +.ih +NAME +asiinit -- initialize the sequential interpolant descriptor using default parameters +.ih +SYNOPSIS +include + +asiinit (asi, interp_type) + +.nf + pointer asi #O interpolant descriptor + int interp_type #I interpolant type +.fi + +.ih +ARGUMENTS +.ls asi +Pointer to sequential interpolant descriptor. +.le +.ls interp_type +Interpolant type. The options are II_NEAREST, II_LINEAR, II_POLY3, II_POLY5, +II_SPLINE3, II_SINC, II_LSINC, and II_DRIZZLE for nearest neighbor, linear, +3rd order polynomial, 5th order polynomial, cubic spline, sinc, look-up +table sinc, and drizzle respectively. The interpolant type definitions are +found in the package header file math/iminterp.h. +.le +.ih +DESCRIPTION +The interpolant descriptor is allocated and initialized. The pointer asi is +returned by ASIINIT. The sinc interpolant width defaults to 31 pixels. The +look-up table sinc interpolant resolution defaults to 20 intervals or +0.05 pixels. The drizzle pixel fraction defaults to 1.0. +.ih +NOTES +ASIINIT or ASISINIT must be called before using any other ASI routines. +.ih +SEE ALSO +asisinit, asifree diff --git a/math/iminterp/doc/asirestore.hlp b/math/iminterp/doc/asirestore.hlp new file mode 100644 index 00000000..6dd70262 --- /dev/null +++ b/math/iminterp/doc/asirestore.hlp @@ -0,0 +1,36 @@ +.help asirestore Dec98 "Image Interpolator Package" +.ih +NAME +asirestore -- restore interpolant +.ih +SYNOPSIS +asirestore (asi, interpolant) + +.nf + pointer asi #O interpolant descriptor + real interpolant[] #I array containing interpolant +.fi +.ih +ARGUMENTS +.ls asi +Pointer to the interpolant descriptor structure. +.le +.ls interpolant +Array containing the interpolant. The length of interpolant can be determined +by a call to ASIGETI. +.le + +.nf + len_interpolant = asigeti (asi, II_ASINSAVE) +.fi +.ih +DESCRIPTION +ASIRESTORE allocates space for the interpolant descriptor and restores the +parameters and coefficients stored in the interpolant array to a structure +for use with ASIEVAL, ASIVECTOR, ASIDER and ASIGRL. +.ih +NOTES +.ih +SEE ALSO +asisave +.endhelp diff --git a/math/iminterp/doc/asisave.hlp b/math/iminterp/doc/asisave.hlp new file mode 100644 index 00000000..7c8ff37a --- /dev/null +++ b/math/iminterp/doc/asisave.hlp @@ -0,0 +1,39 @@ +.help asisave Dec98 "Image Interpolator Package" +.ih +NAME +asisave -- save interpolant +.ih +SYNOPSIS +asisave (asi, interpolant) + +.nf + pointer asi #I interpolant descriptor + real interpolant[] #O array containing the interpolant +.fi +.ih +ARGUMENTS +.ls asi +Pointer to the interpolant descriptor structure. +.le +.ls interpolant +Array where the interpolant is stored. The size of interpolant can be +determined by a call to asigeti. +.le + +.nf + len_interpolant = asigeti (asi, II_ASINSAVE) +.fi +.ih +DESCRIPTION +The interpolant type, number of coefficients and the position of +the first data point in the coefficient array, along with various +parameters such as the sinc interpolant width, sinc look-up table +resolution, and drizzle pixel fraction, are stored in the first +7 elements of the interpolant array. The remaining elements contain +the coefficients calculated by ASIFIT. +.ih +NOTES +.ih +SEE ALSO +asirestore +.endhelp diff --git a/math/iminterp/doc/asisinit.hlp b/math/iminterp/doc/asisinit.hlp new file mode 100644 index 00000000..1ec9bc4e --- /dev/null +++ b/math/iminterp/doc/asisinit.hlp @@ -0,0 +1,60 @@ +.help asisinit Dec98 "Image Interpolator Package" +.ih +NAME +asisinit -- initialize the sequential interpolant descriptor using user parameters +.ih +SYNOPSIS +include + +asisinit (asi, interp_type, nsinc, nincr, rparam, badval) + +.nf + pointer asi #O interpolant descriptor + int interp_type #I interpolant type + int nsinc #I sinc interpolant width in pixels + int nincr #I sinc look-up table resolution + real pixfrac #I sinc shift or drizzle pixel fraction + real badval #I drizzle undefined pixel value +.fi + +.ih +ARGUMENTS +.ls asi +Pointer to sequential interpolant descriptor. +.le +.ls interp_type +Interpolant type. The options are II_NEAREST, II_LINEAR, II_POLY3, II_POLY5, +II_SPLINE3, II_SINC, II_LSINC, and II_DRIZZLE for the nearest neighbour, linear, +3rd order polynomial, 5th order polynomial, cubic spline, sinc, look-up +table sinc, and drizzle interpolants respectively. The interpolant type +definitions are found in the package header file math/iminterp.h. +.le +.ls nsinc +The sinc and look-up table sinc interpolant width in pixels. Nsinc is +rounded up internally to the nearest odd number. +.le +.ls nincr +The look-up table sinc resolution in number of entries. Nincr = 10 implies +a pixel resolution of 0.1 pixels, nincr = 20 a pixel resolution of 0.05 +pixels, etc. +.le +.ls pixfrac +The look-up table sinc fractional pixel shift if nincr = 1 in which case +-0.5 <= pixfrac <= 0.5 , or the drizzle pixel fraction in which case +0.0 <= pixfrac <= 1.0. A minimum value of 0.001 is imposed on pixfrac. +.le +.ls badval +The undefined pixel value for the drizzle interpolant. Pixels within +the boundaries of the input image which do not overlap the input image +pixels are assigned a value of badval. +.le +.ih +DESCRIPTION +The interpolant descriptor is allocated and initialized. The pointer asi is +returned by ASISINIT. +.ih +NOTES +ASIINIT or ASISINIT must be called before using any other ASI routines. +.ih +SEE ALSO +asisinit, asifree diff --git a/math/iminterp/doc/asitype.hlp b/math/iminterp/doc/asitype.hlp new file mode 100644 index 00000000..d8c78b44 --- /dev/null +++ b/math/iminterp/doc/asitype.hlp @@ -0,0 +1,95 @@ +.help asitype Dec98 "Image Interpolator Package" +.ih +NAME +asitype -- decode an interpolant string +.ih +SYNOPSIS +include + +asitype (interpstr, interp_type, nsinc, nincr, rparam) + +.nf + char interpstr #I the input interpolant string + int interp_type #O interpolant type + int nsinc #O sinc interpolant width in pixels + int nincr #O sinc look-up table resolution + real rparam #O sinc or drizzle pixel fraction +.fi + +.ih +ARGUMENTS +.ls interpstr +The user supplied interpolant string to be decoded. The options are +.ls nearest +Nearest neighbor interpolation. +.le +.ls linear +Linear interpolation +.le +.ls poly3 +Cubic polynomial interpolation. +.le +.ls poly5 +Quintic polynomial interpolation. +.le +.ls spline3 +Cubic spline interpolation. +.le +.ls sinc +Sinc interpolation. Users can specify the sinc interpolant width by +appending a width value to the interpolant string, e.g. sinc51 specifies +a 51 pixel wide sinc interpolant. The sinc width will be rounded up to +the nearest odd number. The default sinc width is 31. +.le +.ls lsinc +Look-up table sinc interpolation. Users can specify the look-up table sinc +interpolant width by appending a width value to the interpolant string, e.g. +lsinc51 specifies a 51 pixel wide look-up table sinc interpolant. The user +supplied sinc width will be rounded up to the nearest odd number. The default +sinc width is 31 pixels. Users can specify the resolution of the sinc lookup +up table by appending the look-up table size in square brackets to the +interpolant string, e.g. lsinc51[20] specifies a 20 element sinc look-up +table interpolant with a pixel resolution of 0.05 pixels. The default +look-up table size and resolution are 20 and 0.05 pixels respectively. +The fractional pixel shift for a 1 element look-up table sinc may be +specified by replacing the number of lookup-table elements with the fractional +shift, e.g. sinc51[0.5] will precompute a lookup table for a 0.5 pixel shift. +.le +.ls drizzle +Drizzle interpolation. Users can specify the drizzle pixel fraction by +appending the pixel fraction value to the interpolant string in square +brackets, e.g. drizzle[0.5] specifies a pixel fraction of 0.5. +The default pixel fraction is 1.0. +.le +.le +.ls interp_type +The output interpolant type. The options are II_NEAREST, II_LINEAR, II_POLY3, +II_POLY5, II_SPLINE3, II_SINC, II_LSINC, and II_DRIZZLE for the nearest +neighbor, linear, 3rd order polynomial, 5th order polynomial, cubic spline, +sinc, look-up table sinc, and drizzle interpolants respectively. The +interpolant type definitions are found in the package header file +math/iminterp.h. +.le +.ls nsinc +The output sinc and look-up table sinc interpolant width in pixels. The +default value is 31 pixels. +.le +.ls nincr +The output sinc look-up table size. Nincr = 10 implies a pixel resolution +of 0.1 pixels, nincr = 20 a pixel resolution of 0.05 pixels, etc. The +default value of nincr is 20. +.le +.ls rparam +The output look-up table sinc fractional pixel shift if nincr = 1 in which case +-0.5 <= rparam <= 0.5 , or the drizzle pixel fraction in which case +0.0 <= rparam <= 1.0. +.le +.ih +DESCRIPTION +The interpolant string is decoded into values suitable for the ASISINIT +or ASIINIT routines. +.ih +NOTES +.ih +SEE ALSO +asinit, asisinit, asifree diff --git a/math/iminterp/doc/asivector.hlp b/math/iminterp/doc/asivector.hlp new file mode 100644 index 00000000..95bac138 --- /dev/null +++ b/math/iminterp/doc/asivector.hlp @@ -0,0 +1,52 @@ +.help asivector Dec98 "Image Interpolator Package" +.ih +NAME +asivector -- evaluate the interpolant +.ih +SYNOPSIS +asivector (asi, x, y, npix) + +.nf + pointer asi #I interpolator descriptor + real x[npix/2*npix] #I x array, 1 <= x[i] <= npix + real y[npix] #O array of interpolated values + int npix #I number of x values +.fi +.ih +ARGUMENTS +.ls asi +Pointer to the sequential interpolator descriptor structure. +.le +.ls x +Array of npix x values, or array of npix x ranges if the interpolant is +drizzle. +.le +.ls y +Array of interpolated values. +.le +.ls npix +The number of x values. +.le +.ih +DESCRIPTION +The polynomial coefficients are calculated directly from the data points +for the polynomial interpolants, and from the B-spline coefficients for +the cubic spline interpolant. The actual calculation is done by adding and +multiplying terms according to Everett's central difference interpolation +formula. The boundary extension algorithm is projection. + +The sinc interpolant is computed using a range of data points around +the desired position. Look-up table sinc interpolation is computed +using the most appropriate entry in a precomputed look-up table. +The boundary extension algorithm is nearest neighbor. + +The drizzle interpolant is computed by summing the data over the user +supplied X intervals. +.ih +NOTES +Checking for out of bounds and INDEF valued pixels is the responsibility of the +user. ASIINIT or ASISINIT and ASIFIT must be called before calling ASIVECTOR. +.ih +SEE ALSO +asieval, asider, arieval, arider +.endhelp diff --git a/math/iminterp/doc/im1dinterp.spc b/math/iminterp/doc/im1dinterp.spc new file mode 100644 index 00000000..ce5b8680 --- /dev/null +++ b/math/iminterp/doc/im1dinterp.spc @@ -0,0 +1,525 @@ +.help iminterp Jul84 "Math Package" + +.ce +Specifications for the Image Interpolator Package +.ce +Lindsey Davis +.ce +Vesa Junkkarinen +.ce +August 1984 + +.sh +1. Introduction + + One of the most common operations in image processing is +interpolation in a data array. Due to the large amount of data involved, +efficiency is highly important. The advantage of having locally written +interpolators, includes the ability to optimize for uniformly spaced data +and the possibility of adding features that are useful to the final +application. + +.sh +2. Requirements + +.ls (1) +The package shall take as input a one-dimensional array containing image +data. The pixels are assumed to be equally spaced along a line. +The coordinates of a pixel are assumed to be +the same as the subscript of the pixel in the data array. +The coordinate of the first pixel in the array and the spacing between pixels +is assumed to be 1.0. All pixels are assumed to be good. +Checking for INDEF valued and out of bounds pixels is the responsibility of the +user. A routine to remove INDEF valued pixels from a data array shall be +included in the package. +.le +.ls (2) +The package is divided into array sequential interpolators and array +random interpolators. The sequential interpolators have been optimized +for returning many values as is the case when an array is shifted, or +oversampled at many points in order to produce a +smooth plot. +The random interpolators allow the evaluation of a few interpolated +points without the computing time and storage overhead required for +setting up the sequential version. +.le +.ls (3) +The quality of the interpolant will be set at run time. The options are: + +.nf + II_NEAREST - nearest neighbour + II_LINEAR - linear interpolation + II_POLY3 - 3rd order divided differences + II_POLY5 - 5th order divided differences + II_SPLINE3 - cubic spline +.fi + +The calling sequences shall be invariant to the interpolant selected. +Routines should be designed so that new interpolants can be added +with minimal changes to the code and no change to the calling sequences. +.le +.ls (4) +The interpolant parameters and the arrays necessary to store the coefficients +are stored in a structure referenced by a pointer. The pointer is returned +to the user program by the initial call to ASIINIT or ASIRESTORE and freed +by a call to ASIFREE (see section 3.1). +.le +.ls (5) +The package routines shall be able to: +.ls o +Calculate the coefficients of the interpolant and store these coefficients in +the appropriate part of the interpolant descriptor structure. +.le +.ls o +Evaluate the interplant at a given x(s) coordinate(s). +.le +.ls o +Calculate the derivatives of the interpolant at a given value of x. +.le +.ls o +Integrate the interpolant over a specified x interval. +.le +.ls o +Store the interpolant in a user supplied array. Restore the saved interpolant +to the interpolant descriptor structure for later use by ASIEVAL, ASIVECTOR, +ASIDER and ASIGRL. +.le + +.sh +3. Specifications + +.sh +3.1. The Array Sequential Interpolator Routines + + The package prefix is asi and the package routines are: + +.nf + asiinit (asi, interp_type) + asifit (asi, datain, npix) + y = asieval (asi, x) + asivector (asi, x, yfit, npix) + asider (asi, x, der, nder) + v = asigrl (asi, a, b) + asisave (asi, interpolant) + asirestore (asi, interpolant) + asifree (asi) +.fi + +.sh +3.2. The Array Random Interpolator Routines + + The package prefix is ari and the package routines are: + +.nf + y = arieval (x, datain, npix, interp_type) + arider (x, datain, npix, der, nder, interp_type) +.fi + +.sh +3.3. Miscellaneous + + A routine has been included in the package to remove INDEF valued +pixels from an array. + +.nf + arbpix (datain, dataout, npix, interp_type, boundary_type) +.fi + +.sh +3.4. Algorithms + +.sh +3.4.1. Coefficients + + The coefficient array used by the asi routines is calculated by ASIFIT. +ASIFIT accepts an array of data, checks that the number +of data points is appropriate for the interpolant selected, allocates +space for the interpolant, and calculates the coefficients. +Boundary coefficient values are calculated +using boundary projection. With the exception of the cubic spline interpolant, +the coefficients are stored as the data points. +The B-spline coefficients are +calculated using natural end conditions (Prenter 1975). +After a call to ASIFIT the coefficient array contains the following. + +.nf + case II_NEAREST: + + # no boundary conditions necessary + coeff[1] = datain[1] + . + . + . + coeff[npts] = datain[npix] + + case II_LINEAR: + + # coeff[npxix+1] required if x = npix + coeff[1] = datain[1] + . + . + . + coeff[npix] = datain[npix] + coeff[npix+1] = 2. * datain[npix] - datain[npix-1] + + case II_POLY3: + + # coeff[0] required if x = 1 + # coeff[npix+1], coeff[npix+2] required if x = npix + coeff[0] = 2. * datain[1] - datain[2] + coeff[1] = datain[1] + . + . + . + coeff[npix] = datain[npix] + coeff[npix+1] = 2. * datain[npix] - datain[npix-1] + coeff[npix+2] = 2. * datain[npix] - datain[npix-2] + + case II_POLY5: + + # coeff[1], coeff[0] reqired if x = 1 + # coeff[npix+1], coeff[npix+2], coeff[npix=3] + # required if x = npix + + coeff[-1] = 2. * datain[1] - datain[3] + coeff[0] = 2. * datain[1] - datain[2] + coeff[1] = datain[1] + . + . + . + coeff[npix] = datain[npix] + coeff[npix+1] = 2. * datain[npix] - datain[npix-1] + coeff[npix+2] = 2. * datain[npix] - datain[npix-2] + coeff[npix+3] = 2. * datain[npix] - datain[npix-3] + + case SPLINE3: + + # coeff[0] = 2nd der at x = 1, coeff[0] = 0. + # coeff[npix+1] = 2nd der at x = npts, coeff[npix+1] = 0. + # coeff[npix+2] = 0., required if x = npix + coeff[0] = b[1] + coeff[1] = b[2] + . + . + . + coeff[npix] = b[npix+1] + coeff[npix+1] = b[npix+2] + coeff[npix+2] = 0. +.fi + +.sh +3.4.2. Evaluation + + The ASIEVAL and ASIVECTOR routines have been optimized to be as efficient +as possible. The values of the II_NEAREST and II_LINEAR interpolants +are calculated directly. The II_SPLINE3 interpolant is evaluated using +polynomial coefficients calculated directly from the B-spline coefficients +(de Boor 1978). Values of the higher order polynomial interpolants +are calculated using central differences. The equations for each case are +listed below. + +.nf +case II_NEAREST: + + y = coeff[int (x + 0.5)] + +case II_LINEAR: + + nx = x + y = (x - nx) * coeff[nx+1] + (nx + 1 - x) * coeff[nx] + +case II_POLY3: + + nx = x + s = x - nx + t = 1. - s + + # second central differences + cd20 = 1./6. * (coeff[nx+1] - 2. * coeff[nx] + coeff[nx-1]) + cd21 = 1./6. * (coeff[nx+2] - 2. * coeff[nx+1] + coeff[nx]) + + y = s * (coeff[nx+1] + (s * s - 1.) * cd21) + t * (coeff[nx] + + (t * t - 1.) * cd20) + +case II_POLY5: + + nx = x + s = x - nx + t = 1. - s + + # second central differences + cd20 = 1./6. * (coeff[nx+1] - 2. * coeff[nx] + coeff[nx-1]) + cd21 = 1./6. * (coeff[nx+2] - 2. * coeff[nx+1] + coeff[nx]) + + # fourth central diffreences + cd40 = 1./120. * (coeff[nx-2] - 4. * coeff[nx-1] + 6. * coeff[nx] - 4. * + coeff[nx+1] + a[nx+2]) + cd41 = 1./120. * (coeff[nx-1] - 4. * coeff[nx] + 6. * coeff[nx+1] - 4. * + coeff[nx+2] + coeff[nx+3] + + y = s * (coeff[nx+1] + (s * s - 1.) * (cd21 + (s * s - 4.) * cd41)) + + t * (coeff[nx] + (t * t - 1.) * (cd20 + (t * t - 4.) * cd40)) + +case II_SPLINE3: + + nx = x + s = x - nx + + pc[1] = coeff[nx-1] + 4. * coeff[nx] + coeff[nx+1] + pc[2] = 3. * (coeff[nx+1] - coeff[nx-1]) + pc[3] = 3. * (coeff[nx-1] - 2. * coeff[nx] + coeff[nx+1]) + pc[4] = -coeff[nx-1] + 3. * coeff[nx] - 3. * coeff[nx+1] + coeff[nx+2] + + y = pc[1] + s * (pc[2] + s * (pc[3] + s * pc[4])) +.fi + + + The ARIEVAL routine uses the expressions above to evaluate the +interpolant. However unlike ASIEVAL, ARIEVAL does not use a previously +calculated coefficient array. Instead ARIEVAL selects the appropriate +portion of the data array, calculates the coefficients and boundary +coefficients if necessary, and evaluates the interpolant at the time it +is called. The cubic spline interpolant uses at most SPLTS (currently 16) +data points to calculate the B-spline coefficients. + +.sh +3.4.3. Derivatives + + Derivatives of the interpolant are calculated by evaluating the +derivatives of the interpolating polynomial. For all interpolants der[1] +equals the value of the interpolant at x. +For the sake of efficiency the derivatives +of the II_NEAREST and II_LINEAR interpolants are calculated directly. + +.nf + case II_NEAREST: + + der[1] = coeff[int (x+0.5)] + + case II_LINEAR: + + der[1] = (x - nx) * coeff [nx+1] + (nx + 1 - x) * coeff[nx] + der[2] = coeff[nx+1] - coeff[nx] +.fi + + In order to calculate the derivatives of the cubic spline and +polynomial interpolants +the coefficients of the interpolating polynomial must be calculated. +The polynomial +coefficients for the cubic spline interpolant are computed directly from the +B-spline coefficients (see 3.4.2.). The higher order polynomial +interpolant coefficients are calculated as follows. + +First the appropriate portion of the coefficient array is loaded. + +.nf + do i = 1, nterms + d[i] = coeff[nx - nterms/2 + i] +.fi + +Next the divided differences are calculated (Conte and de Boor 1972). + +.nf + do k = 1, nterms - 1 + do i = 1, nterms - k + d[i] = (d[i+1] - d[i]) / k +.fi + +The d[i] are the coefficients of an interpolating polynomial of the +following form. The x[i] are the nterms data points surrounding the +point of interest. + +.nf + p(x) = d[1] * (x-x[1]) * ... * (x-x[nterms-1) + + d[2] * (x-x[2]) * ... * (x-x[nterms-1]) + ... + d[nterms] +.fi + +Next a new set of polynomial coefficients are calculated +(Conte and de Boor 1972). + +.nf + do k = nterms, 2, -1 + do i = 2, k + d[i] = d[i] + d[i-1] * (k - i - nterms/2) +.fi + +The new d[i] are the coefficients of the follwoing polynomial. + +.nf + nx = x + p(x) = d[1] * (x-nx) ** (nterms-1) + d[2] * (x-nx) ** (nterms-2) + ... + d[nterms] +.fi + +The d[i] array is flipped. The value and derivatives +of the interpolant are then calculated using the d[i] array and +nested multiplication. + +.nf + s = x - nx + + do k = 1, nder { + + accum = d[nterms-k+1] + + do j = nterms - k, 1, -1 + accum = d[j] + s * accum + + der[k] = accum + + # differnetiate + do j = 1, nterms - k + d[j] = j * d[j + 1] + } +.fi + + ARIDER calculates the derivatives of the interpolant using the same +technique ASIDER. However ARIDER does not use a previously calculated +coefficient array like ASIDER. Instead ARIDER selects the appropriate portion +of the data array, calculates the coefficients and boundary coefficients, +and computes the derivatives at the time it is called. + +.sh +3.4.5. Integration + + ASIGRL calculates the integral of the interpolant between fixed limits +by integrating the interpolating polynomial. The coefficients of the +interpolating polynomial are calculated as discussed in section 3.4.4. + +.sh +4. Usage + +.sh +4.1. User Notes + +The following series of steps illustrates the use of the package. + +.ls 4 +.ls (1) +Insert an include statement in the calling program to make +the IINTERP definitions available to the user program. +.le +.ls (2) +Remove INDEF valued pixels from the data using ARBPIX. +.le +.ls (3) +Call ASIINIT to initialize the interpolant parameters. +.le +.ls (4) +Call ASIFIT to calculate the coefficients of the interpolant. +.le +.ls (5) +Evaluate the interpolant at a given value of x(s) using ASIEVAL or +ASIVECTOR. +.le +.ls (6) +Calculate the derivatives and integral or the interpolant using +ASIDER and ASIGRL. +.le +.ls (7) +Free the interpolator structure by calling ASIFREE. +.le +.le + + The interpolant can be saved and restored using ASISAVE and ASIRESTORE. +If the values and derivatives of only a few points in an array are desired +ARIEVAL and ARIDER can be called. + +.sh +4.2. Examples + +.nf +Example 1: Shift a data array by a constant amount + + include + ... + call asiinit (asi, II_POLY5) + call asifit (asi, inrow, npix) + + do i = 1, npix + outrow[i] = asieval (asi, i + shift) + + call asifree (asi) + ... + +Example 2: Calculate the integral under the data array + + include + ... + call asiinit (asi, II_POLY5) + call asifit (asi, datain, npix) + + integral = asigrl (asi, 1. real (npix)) + + call asifree (asi) + ... + +Example 2: Store interpolant for later use by ASIEVAL + LEN_INTERP must be at least npix + 8 units long where npix is + is defined in the call to ASIFIT. + + include + + real interpolant[LEN_INTERP] + ... + call asiinit (asi, II_POLY3) + call asifit (asi, datain, npix) + call asisave (asi, interpolant) + call asifree (asi) + ... + call asirestore (asi, interpolant) + do i = 1, npts + yfit[i] = asieval (asi, x[i]) + call asifree (asi) + ... +.fi +.sh +5. Detailed Design + +.sh +5.1. Interpolator Descriptor + + The interpolant parameters and coefficients are stored in a +structure listed below. + +.nf + define LEN_ASISTRUCT 4 # Length in structure units of + # interpolant descriptor + + define ASI_TYPE Memi[$1] # Interpolant type + define ASI_NCOEFF Memi[$1+1] # No. of coefficients + define ASI_OFFSET Memi[$1+2] # First "data" point in + # coefficient array + define ASI_COEFF Memi[$1+3] # Pointer to coefficient array +.fi + +.sh +5.2. Storage Requirements + + The interpolant descriptor requires LEN_ASISTRUCT storage units. The +coefficient array requires ASI_NCOEFF(asi) real storage elements, where +ASI_NCOEFF(asi) is defined as follows. + +.nf + ASI_NCOEFF(asi) = npix # II_NEAREST + ASI_NCOEFF(asi) = npix+1 # II_LINEAR + ASI_NCOEFF(asi) = npix+3 # II_POLY3 + ASI_NCOEFF(asi) = npix+5 # II_POLY5 + ASI_NCOEFF(asi) = npix+3 # II_SPLINE3 +.fi + +.sh +6. References + +.ls (1) +Carl de Boor, "A Practical Guide to Splines", 1978, Springer-Verlag New +York Inc. +.le +.ls (2) +S.D. Conte and C. de Boor, "Elementary Numerical Analysis", 1972, McGraw-Hill, +Inc. +.le +.ls (3) +P.M. Prenter, "Splines and Variational Methods", 1975, John Wiley and Sons Inc. +.le +.endhelp diff --git a/math/iminterp/doc/im2dinterp.spc b/math/iminterp/doc/im2dinterp.spc new file mode 100644 index 00000000..e4d88be3 --- /dev/null +++ b/math/iminterp/doc/im2dinterp.spc @@ -0,0 +1,432 @@ +.help iinterp Dec84 "Math Package" +.ce +Specifications for the 2D-Image Interpolator Package +.ce +Lindsey Davis +.ce +December 1984 + +.sh +1. Introduction + +One of the most common operations required in image processing is +two-dimensional interpolation in a data array. Any image operator which +physically moves pixels around requires image interpolation to calculate +the new gray levels of the pixels. The advantage +of having a locally written image interpolation package includes the ability to +optimize for uniformly spaced data and the possibility of adding features +that are useful to the final application. + +.sh +2. Requirements + +.ls (1) +The package shall take as input a 2-D array containing an image or image +section. The pixels are assumed to be equally spaced on a rectangular grid. +The coordinates of the pixels +are assumed to be the same as the subscripts of the pixel in the data array. +Therefore the coordinates of the first pixel in the array are assumed +to be (1,1). For operations on image sections the calling program must +keep track of any transformations between image and section coordinates. +All pixels are assumed to be good. +Checking for INDEF and out of bounds pixels is +the responsibility of the user. A routine to remove INDEF valued pixels +ARBPIX is available in the 1-D package. +.le +.ls (2) +The package is divided into array sequential interpolants and array random +interpolants. The sequential interpolants have been optimized for returning +many values as when an array is shifted or when an array is oversampled +at many points in order to produce a smooth surface plot. +The random interpolants +allow the evaluation of a few interpolated points without the computing +time and storage overhead required for setting up the sequential version. +.le +.ls (3) +The quality of the interpolant will be set at run time. The options are: + +.nf + II_BINEAREST # nearest neighbour + II_BILINEAR # bilinear + II_BIPOLY3 # bicubic interior polynomial + II_BIPOLY5 # biquintic interior polynomial + II_BISPLINE3 # bicubic spline +.fi + +The calling sequences shall be invariant to the interpolant selected. +Routines should be designed so that new interpolants can be added with +minimal changes to the code and no changes to the calling sequences. +.le +.ls (4) +The interpolant parameters and the arrays necessary to store the +coefficients are stored in a structure referenced by a pointer. The pointer +is returned to the user program by the initial call to MSIINIT or +MSIRESTORE and freed by a call to MSIFREE. +.le +.ls (5) +The package routines should be able to: +.ls o +Calculate the coefficients of the interpolant and store these coefficients +in the appropriate part of the interpolant descriptor structure. +.le +.ls o +Evaluate the interpolant at a given x-y(s) coordinate. +.le +.ls o +Evaluate the first nder derivatives at a given value of x-y. +.le +.ls o +Integrate the interpolant over a specified interval in x-y. +.le + +.sh +3. Specifications + +.sh +3.1. The Matrix Sequential Interpolator Routines + +The package prefix is msi and the package routines are: + +.nf + msiinit (msi, interp_type) + msifit (msi, datain, nxpix, nypix, len_datain) + y = msieval (msi, x, y) + msivector (msi, x, y, zfit, npix) + msigrid (msi, x, y, zfit, nx, ny, len_zfit) + msider (msi, x, y, der, nxder, nyder, len_der) + v = msigrl (msi, x, y, npts) + v = msisqgrl (msi, x1, x2, y1, y2) + msisave (msi, interpolant) + msirestore (msi, interpolant) + msifree (msi) +.fi + +.sh +3.2. The Matrix Random Interpolator Routines + +The package prefix is mri and the package routines are: + +.nf + y = mrieval (x, y, datain, nx, ny, len_datain, interp_type) + mrider (x, y, datain, nx, ny, len_datain, der, nxder, nyder, + len_der, interp_type) +.fi + +.sh +3.3. Miscellaneous + +A routine ARBPIX will remove INDEF valued pixels from an +array and is available in the 1-D package. + +.nf + arbpix (datain, dataout, npix, interp_type, boundary_type) +.fi + +.sh +3.4. Algorithms + +.sh +3.4.1. Coefficients + +The coefficients for the msi routines are calculated by MSIFIT. MSIFIT +accepts the input data, checks that the number of data pixels is +appropriate for the interpolant selected, and allocates space for the +interpolant. Boundary coefficient values are calculated using boundary +projection. With the exception of the II_BISPLINE3 option, the interpolant +is stored as the data points. The B-spline coefficients are calculated +using the "natural" end conditions in two steps. +First the B-spline coefficients in x at each value of y are calculated. +The B-x coefficients are then solved for the B-spline coefficients in x-y. +After a call to MSIFIT the coefficient +array contains the following. + +.nf + CASE II_BINEAREST: + + # no boundary extension required, coeff[1:nx,1:ny] + + coeff[i,j] = data[i,j] i = 1,...,nx + j = 1,...,ny + + case II_BILINEAR: + + # must extend nx by 1 and ny by 1, coeff[1:nx+1,1:ny+1] + + coeff[i,j] = data[i,j] i = 1,...,nx + j = 1,...,ny + + coeff[nx+1,j] = 2 * data[nx] - data[nx-1] j = 1,...,ny + coeff[i,ny+1] = 2 * data[ny] - data[ny-1] i = 1,...,nx + + coeff[nx+1,ny+1] = 2 * coeff[nx+1,ny] - data[nx,ny] + + case II_BIPOLY3: + + # must extend nx by -1 and 2 and ny by -1 and 2, coeff[0:nx+2,0:ny+2] + + coeff[i,j] = data[i,j] i = 1,...,nx + j = 1,...,ny + + coeff[0,j] = 2 * data[1,j] - data[2,j] j = 1,...,ny + coeff[nx+1,j] = 2 * data[nx,j] - data[nx-1,j] j = 1,...,ny + coeff[nx+2,j] = 2 * data[nx,j] - data[nx-2,j] j = 1,...,ny + + coeff[i,0] = 2 * data[i,1] - data[i,2] i = 1,...,ny + coeff[i,ny+1] = 2 * data[i,ny] - data[i,ny-1] i = 1,...,nx + coeff[i,ny+2] = 2 * data[i,ny] - data[i,ny-2] i = 1,...,nx + + # plus remaining points + + case II_BIPOLY5: + + # extend -2 and 3 in nx and -2 and 3 in ny, coeff[-1:nx+3,-1:ny+3] + + coeff[i,j] = data[i,j] i = 1,...,nx + j = 1,...,ny + + coeff[-1,j] = 2 * data[1,j] - data[3,j] j = 1,...,ny + coeff[0,j] = 2 * data[1,j] - data[2,j] j = 1,...,ny + coeff[nx+1,j] = 2 * data[nx,j] - data[nx-1,j] j = 1,...,ny + coeff[nx+2,j] = 2 * data[nx,j] - data[nx-2,j] j = 1,...,ny + coeff[nx+3,j] = 2 * data[nx,j] - data[nx-3,j] j = 1,...,ny + + coeff[i,-1] = 2 * data[i,1] - data[i,3] i = 1,...,nx + coeff[i,0] = 2 * data[i,1] - data[i,2] i = 1,...,nx + coeff[i,ny+1] = 2 * data[i,ny] - data[i,ny-1] i = 1,...,nx + coeff[i,ny+2] = 2 * data[i,ny] - data[i,ny-2] i = 1,...,nx + coeff[i,ny+3] = 2 * data[i,ny] - data[i,ny-3] i = 1,...,nx + + # plus remaining conditions + + case II_BISPLINE3: + + # the natural boundary conditions are used, coeff[0:nx+2,0:ny+2] + + coeff[i,j] = data[i,j] i = 1,...,nx + j = 1,...,ny + + coeff[i,0] = 0. i = 0,...,nx+2 + coeff[i,ny+1] = 0. i = 0,...,nx+2 + coeff[i,ny+2] = 0. i = 0,...,nx+2 + + coeff[0,j] = 0. j = 1,...,ny + coeff[nx+1,j] = 0. j = 1,...,ny + coeff[nx+2,j] = 0. j = 1,...,ny + + # plus remaining coefficients + +.fi + +.sh +3.4.2. Evaluation + +The MSIEVAL and MSIVECTOR routines will be optimized to be as efficient +as possible for evaluating arbitrarily spaced data. A special function +MSIGRID is included for evaluating closely spaced points on a rectangular grid. +For the options II_BINEAREST and II_BILINEAR the value of the +interpolant is calculated directly. The II_BISPLINE3 interpolant is +evaluated using polynomial coefficients calculated directly from the +B-spline coefficients. Values of the higher order polynomial interpolants +are calculated using Everett's central difference formula. The equations +are listed below. + +.nf +case II_BINEAREST: + + z = coeff[int (x + 0.5), int (y + 0.5)) + +case II_BILINEAR: + + sx = x - nx sy = y - ny + tx = 1. - sx ty = 1. - sy + + z = tx * ty * coeff[nx,ny] + sx * ty * coeff[nx+1,ny] + + sy * tx * coeff[nx,ny+1] + sx * sy * coeff[nx+1,ny+1] + + +case II_BIPOLY3: + + nx = x + sx = x - nx + tx = 1. - sx + + # interpolate in x + + i = 1 + + do j = ny-1, ny+2 { + + cd20[i] = 1./6. * (coeff[nx+1,j] - 2. * coeff[nx,j] + coeff[nx-1,j]) + cd21[i] = 1./6. * (coeff[nx+2,j] - 2. * coeff[nx+1,j] + coeff[nx,j]) + + z[i] = sx * (coeff[nx+1,j] + (sx * sx - 1.) * cd21[i]) + + tx * (coeff[nx,j] + (tx * tx - 1.) * cd20[i]) + + i = i + 1 + } + + ny = y + sy = y - ny + ty = 1. - sy + + # interpolate in y + + cd20y = 1./6. * (z[3] - 2. * z[2] + z[1]) + cd21y = 1./6. * (z[4] - 2. * z[3] + z[2]) + + value = sy * (z[3] + (sy * sy - 1.) * cd21y) + + ty * (z[2] + (ty * ty - 1.) *cd20y) + + +case II_BIPOLY5: + + nx = x + sx = x - nx + sx2 = sx * sx + tx = 1. - sx + tx2 = tx * tx + + # interpolate in x + + i = 1 + + do j = ny-2, ny+3 { + + cd20[i] = 1./6. * (coeff[nx+1,j] - 2. * coeff[nx,j] + coeff[nx-1,j]) + cd21[i] = 1./6. * (coeff[nx+2,j] - 2. * coeff[nx+1,j] + coeff[nx,j]) + cd40[i] = 1./120. * (coeff[nx-2,j] - 4. * coeff[nx-1,j] + + 6. * coeff[nx,j] - 4. * coeff[nx+1,j] + coeff[nx+2,j]) + cd41[i] = 1./120. * (coeff[nx-1,j] - 4. * coeff[nx,j] + + 6. * coeff[nx+1,j] - 4. * coeff[nx+2,j] + coeff[nx+3,j]) + + z[i] = sx * (coeff[nx+1,j] + (sx2 - 1.) * (cd21[j] + (sx2 - 4.) * + cd41[j])) + tx * (coeff[nx,j] + (tx2 - 1.) * + (cd20[j] + (tx2 - 4.) * cd40[j])) + + i = i + 1 + } + + ny = y + sy = y - ny + sy2 = sy * sy + ty = 1. - sy + ty2 = ty * ty + + # interpolate in y + + cd20y = 1./6. * (z[3] - 2. * z[2] + z[1]) + cd21y = 1./6. * (z[4] - 2. * z[3] + z[2]) + cd40y = 1./120. * (z[1] - 4. * z[2] + 6. * z[3] - 4. * z[4] + z[5]) + cd41y = 1./120. * (z[2] - 4. * z[3] + 6. * z[4] - 4. * z[5] + z[6]) + + value = sy * (z[4] + (sy2 - 1.) * (cd21y + (sy2 - 4.) * cd41y)) + + ty * (z[3] + (ty2 - 1.) * (cd20y + (ty2 - 4.) * cd40y)) + +case II_BISPLINE3: + + # use the B-spline representation + + value = coeff[i,j] * B[i](x) * B[j](y) + + # the B-splines are the following + + B1(x) = (x - x1) ** 3 + B2(x) = 1 + 3 * (x - x2) + 3 * (x - x2) ** 2 - 3 * (x - x2) ** 3 + B3(x) = 1 + 3 * (x3 - x) + 3 * (x3 - x) ** 2 - 3 * (x3 - x) ** 3 + B4(x) = (x4 - x) ** 3 + + # the y B-splines are identical + +.fi + +.sh +3.4.3. Derivatives + +The derivatives are calculated by evaluating the derivatives of the +interpolating polynomial. The 0-th derivative is the value of the +interpolant. The 1-st derivatives are d/dx and d/dy. The second are +d2/dx2, d2/dy2 and d2/dxdy. The derivatives in the case II_BINEAREST +and II_BILINEAR are calculated directly. + +.nf +case II_BINEAREST: + + der[1,1] = value # see previous section + +case II_BILINEAR: + + der[1] = value # see previous section + + # d/dx + der[2,1] = -ty * coeff[nx,ny] + ty * coeff[nx+1,ny] - + sy * coeff[nx,ny+1] + sy * coeff[nx+1,ny+1] + + # d/dy + der[1,2] = -tx * coeff[nx,ny] + sx * coeff[nx+1,ny] + + tx * coeff[nx,ny+1] + sx * coeff[nx+1,ny+1] + + # d2/dxdy + der[2,2] = coeff[nx,ny] - coeff[nx+1,ny] - coeff[nx,ny+1] + coeff[nx+1,ny+1] + +case II_BIPOLY3, II_BIPOLY5, II_BISPLINE3: +.fi + +For the higher order interpolants the coefficients of the interpolating +polynomial in x and y are calculated. In the case of II_BIPOLY3 and II_BIPOLY5 +this is the Everett polynomial of the 3-rd and 5-th order respectively. +In the case of II_BISPLINE3 the pp-representation is calculated for +the B-spline coefficients. The value of the interpolant and its +derivatives is calculated using nested multiplication. + +.sh +3.4.5. Integration + +Integration is most easily accomplished by integrating the interpolant in +x for each value of y. The resulting function of y can then be integrated +in y to derive the 2-D integral. The limits of the integral are assumed +to be the corners of a polygon. At minimum of three points which define +a triangular region in x-y are required. The integral will be an approximation. +A special function for integrating over a rectangular region is also +provided. + +.sh +4. Detailed Design + +.sh +4.1. Interpolant Descriptor + +The interpolant parameters and coefficients will be stored in a structure +as listed below. + +.nf + define LEN_MSISTRUCT 10 + + struct { + + int msi_type # interpolant type + int msi_nxcoeff # size of coefficient array in x + int msi_nycoeff # size of coefficient array in y + int msi_fstpnt # first datapoint in coefficient array + int asi_coeff # pointer to 1-D coefficient array + + } msistruct + +.fi + +.sh +4.2. Storage Requirements + +The interpolant descriptor requires LEN_MSISTRUCT storage units. +The coefficient array storage is dynamically allocated and requires +msi_nxcoeff * msi_nycoeff real storage elements. The requirements for +each interpolant type are listed below. + +.nf + II_BINEAREST nxpix * nypix + II_BILINEAR (nxpix + 1) * (nypix + 1) + II_BIPOLY3 (nxpix + 3) * (nypix + 3) + II_BIPOLY5 (nxpix + 5) * (nypix + 5) + II_BISPLINE3 (nxpix + 3) * (nypix + 3) +.fi + +.endhelp diff --git a/math/iminterp/doc/iminterp.hd b/math/iminterp/doc/iminterp.hd new file mode 100644 index 00000000..de836c3e --- /dev/null +++ b/math/iminterp/doc/iminterp.hd @@ -0,0 +1,37 @@ +# Help directory for the IMINTERP (image interpolator) package. + +$iminterp = "math$iminterp/" + +arbpix hlp = arbpix.hlp, src = iminterp$arbpix.x +arider hlp = arider.hlp, src = iminterp$arider.x +arieval hlp = arieval.hlp, src = iminterp$arieval.x +asider hlp = asider.hlp, src = iminterp$asider.x +asieval hlp = asieval.hlp, src = iminterp$asieval.x +asifit hlp = asifit.hlp, src = iminterp$asifit.x +asifree hlp = asifree.hlp, src = iminterp$asifree.x +asigeti hlp = asigeti.hlp, src = iminterp$asigeti.x +asigetr hlp = asigetr.hlp, src = iminterp$asigetr.x +asigrl hlp = asigrl.hlp, src = iminterp$asigrl.x +asiinit hlp = asiinit.hlp, src = iminterp$asiinit.x +asirestore hlp = asirestore.hlp, src = iminterp$asirestore.x +asisave hlp = asisave.hlp, src = iminterp$asisave.x +asisinit hlp = asisinit.hlp, src = iminterp$asisinit.x +asivector hlp = asivector.hlp, src = iminterp$asivector.x +asitype hlp = asitype.hlp, src = iminterp$asitype.x + +mrider hlp = mrider.hlp, src = iminterp$mrider.x +mrieval hlp = mrieval.hlp, src = iminterp$mrieval.x +msider hlp = msider.hlp, src = iminterp$msider.x +msieval hlp = msieval.hlp, src = iminterp$msieval.x +msifit hlp = msifit.hlp, src = iminterp$msifit.x +msifree hlp = msifree.hlp, src = iminterp$msifree.x +msigeti hlp = msigeti.hlp, src = iminterp$msigeti.x +msigetr hlp = msigetr.hlp, src = iminterp$msigetr.x +msigrid hlp = msigrid.hlp, src = iminterp$msigrid.x +msigrl hlp = msigrl.hlp, src = iminterp$msigrl.x +msiinit hlp = msiinit.hlp, src = iminterp$msiinit.x +msirestore hlp = msirestore.hlp, src = iminterp$msirestore.x +msisave hlp = msisave.hlp, src = iminterp$msisave.x +msisinit hlp = msisinit.hlp, src = iminterp$msisinit.x +msitype hlp = msitype.hlp, src = iminterp$msitype.x +msivector hlp = msivector.hlp, src = iminterp$msivector.x diff --git a/math/iminterp/doc/iminterp.hlp b/math/iminterp/doc/iminterp.hlp new file mode 100644 index 00000000..c602b43e --- /dev/null +++ b/math/iminterp/doc/iminterp.hlp @@ -0,0 +1,234 @@ +.help iminterp Dec98 "Math Package" +.ih +NAME +iminterp -- image interpolator package +.ih +SYNOPSIS + +.nf + asitype (interpstr, interp_type, nsinc, nincr, rparam) + asiinit (asi, interp_type) + asisinit (asi, interp_type, nsinc, nincr, rparam) + asifit (asi, datain, npix) +ivalue = asigeti (asi, param) +rvalue = asigetr (asi, param) + y = asieval (asi, x) + asivector (asi, x, yfit, npix) + asider (asi, x, der, nder) + v = asigrl (asi, a, b) + asisave (asi, interpolant) + asirestore (asi, interpolant) + asifree (asi) + + y = arieval (x, datain, npix, interp_type) + arider (x, datain, npix, der, nder, interp_type) + + arbpix (datain, dataout, npix, interp_type, boundary_type) +.fi + +.nf + msitype (interpstr, interp_type, nsinc, nincr, rparam) + msisinit (msi, interp_type, nsinc, nxincr, nyincr, rparam1, rparam2) + msiinit (msi, interp_type) + msifit (msi, datain, nxpix, nypix, len_datain) +ivalue = msigeti (msi, param) +rvalue = msigetr (msi, param) + y = msieval (msi, x, y) + msivector (msi, x, y, zfit, npts) + msider (msi, x, y, der, nxder, nyder, len_der) + v = msigrl (msi, x, y, npts) + v = msisqgrl (msi, x1, x2, y1, y2) + msisave (msi, interpolant) + msirestore (msi, interpolant) + msifree (msi) + + y = mrieval (x, y, datain, nxpix, nypix, len_dataina, interp_type) + mrider (x, y, datain, nxpix, nypix, len_datain, der, nxder, nyder, + len_der, interp_type) +.fi + +.ih +DESCRIPTION +The iminterp package provides a set of routines for interpolating uniformly +spaced data assuming that the spacing between data points is 1.0. The +package is divided into 1D and 2D array sequential interpolants, +prefixes asi and msi, and 1D and 2D +array random interpolants, prefixes ari and mri. +The sequential interpolants have +been optimized for returning many values as is the case when an array is +shifted. The random interpolants allow evaluation of a few interpolated +points without the computing time and storage overhead required for +setting up the sequential version. +.ih +NOTES +The interpolant is chosen at run time from the following list. + +.nf + II_NEAREST # nearest neighbour in x + II_LINEAR # linear interpolation in x + II_POLY3 # 3rd order interior polynomial in x + II_POLY5 # fifth order interior polynomial in x + II_SPLINE3 # cubic spline in x + II_SINC # sinc interpolation in x + II_LSINC # look-up table sinc interpolation in x + II_DRIZZLE # drizzle interpolation in x + + II_BINEAREST # nearest neighbour in x and y + II_BILINEAR # bilinear interpolation + II_BIPOLY3 # 3rd order interior polynomial in x and y + II_BIPOLY5 # 5th order interior polynomial in x and y + II_BISPLINE3 # bicubic spline + II_BISINC # sinc interpolation in x and y + II_BILSINC # look-up table sinc interpolation in x and y + II_BIDRIZZLE # drizzle interpolation in x and y +.fi + +The routines assume that all x (1D, 2D) and y (2D) values of interest lie in +the region 1 <= x <= nxpix, 1 <= y <= nypix. +Checking for out of bounds x and/or y values is the responsibility +of the calling program. The asi, ari, msi, and mri routines assume that INDEF +valued pixels have been removed from the data prior to entering the +package. The routine ARBPIX has been added to the package to facilitate +INDEF valued pixel removal. + +In order to make the package definitions available to the calling program +an include statement must appear in the calling program. +Either ASIINIT, ASISINIT or ASIRESTORE must be called before using the +asi routines. ASIFREE frees the space used by the asi routines. For the +msi routines the corresponding examples are MSIINIT, MSISINIT, MSIRESTORE +and MSIFREE. +.ih +EXAMPLES +.nf +Example 1: Shift a 1D data array by a constant amount using a 5th order +polynomial interpolant and the drizzle routine respectively. Note that +in this example the drizzle interpolant is equivalent to the linear +interpolant since the default drizzle pixel fraction is 1.0 and there +is no scale change. Out-of-bounds pixels are set to 0.0 + + include + ... + call asiinit (asi, II_POLY5) + call asifit (asi, inrow, npix) + + do i = 1, npix + if (i + xshift < 1.0 || i + xshift > npix) + outrow[i] = 0.0 + else + outrow[i] = asieval (asi, i + xshift) + + call asifree (asi) + ... + + + include + + real tmpx[2] + ... + call asiinit (asi, II_DRIZZLE) + call asifit (asi, inrow, npix) + + do i = 1, npix + tmpx[1] = i + xshift - 0.5 + tmpx[2] = i + xshift + 0.5 + if (tmpx[1] < 1 || tmpx[2] > npix) + outrow[i] = 0.0 + else + outrow[i] = asieval (asi, tmpx) + + call asifree (asi) + ... + + +Example 2: Shift a 2D array by a constant amount using a 3rd order polynomial +interpolant and the drizzle interpolant respectively. Note that +in this example the drizzle interpolant is equivalent to the linear +interpolant since the default drizzle pixel fraction is 1.0 and there +is no scale change. Out-of-bounds pixels are set to 0.0. + + include + ... + call msiinit (msi, II_BIPOLY3) + call msifit (msi, insection, nxpix, nypix, nxpix) + + do j = 1, nypix + if (j + yshift < 1 || j + yshift > nypix) + do i = 1, nxpix + outsection[i,j] = 0.0 + else + do i = 1, nxpix + if (i + xshift < 1 || i + xshift > nxpix) + outsection[i,j] = 0.0 + else + outsection[i,j] = msieval (msi, i + xshift, j + yshift) + + call msifree (msi) + ... + + + include + ... + real tmpx[4], tmpy[4] + ... + call msiinit (msi, II_BIDRIZZLE) + call msifit (msi, insection, nxpix, nypix, nxpix) + + do j = 1, nypix { + tmpy[1] = j + yshift - 0.5 + tmpy[2] = j + yshift - 0.5 + tmpy[3] = j + yshift + 0.5 + tmpy[4] = j + yshift + 0.5 + if (tmpy[1] < 1 || tmpy[4] > nypix) + do i = 1, nxpix + outsection[i,j] = 0.0 + else + do i = 1, nxpix + tmpx[1] = i + xshift - 0.5 + tmpx[2] = i + xshift + 0.5 + tmpx[3] = i + xshift + 0.5 + tmpx[4] = i + xshift - 0.5 + if (tmpx[1] < 1 || tmpx[2] > nxpix) + outsection[i,j] = 0.0 + else + outsection[i,j] = msieval (msi, tmpx, tmpy) + } + + call msifree (msi) + ... + + +Example 3: Calculate the integral under a 1D data array + + include + ... + call asiinit (asi, II_POLY5) + call asifit (asi, datain, npix) + + integral = asigrl (asi, 1. real (npix)) + + call asifree (asi) + ... + +Example 4: Store a 1D interpolant for later use by ASIEVAL + + include + + ... + call asiinit (asi, II_POLY3) + call asifit (asi, datain, npix) + + len_interpolant = asigeti (asi, ASINSAVE) + call salloc (interpolant, len_interpolant, TY_REAL) + call asisave (asi, Memr[interpolant]) + + call asifree (asi) + ... + call asirestore (asi, Memr[interpolant]) + + do i = 1, npts + yfit[i] = asieval (asi, x[i]) + + call asifree (asi) + ... +.fi +.endhelp diff --git a/math/iminterp/doc/iminterp.men b/math/iminterp/doc/iminterp.men new file mode 100644 index 00000000..9daf5883 --- /dev/null +++ b/math/iminterp/doc/iminterp.men @@ -0,0 +1,32 @@ + arbpix - Replace bad pixels in an array + arider - Calculate derivatives for a few points in a 1-D array + arieval - Evaluate the interpolant at a few points in a 1-D array + asider - Evaluate derivatives at x + asieval - Evaluate the interpolant at x + asifit - Fit the 1-D interpolant + asifree - Free space allocated by asiinit or asisinit + asigeti - Fetch an integer parameter + asigrl - Calculate the integral under an array + asiinit - Initialize 1-D interpolant using default parameters + asirestore - Restore interpolant parameters and coefficients + asisave - Save interpolant parameters and coefficients + asisinit - Initialize 1-D interpolant using user parameters + asitype - Decode string into interpolant type and parameters + asivector - Evaluate interpolant at an array of x + + mrider - Calculate derivatives for a few points in a 2-D array + mrieval - Evaluate the interpolant at a few points in a 2-D array + msider - Evaluate derivatives at x and y + msieval - Evaluate the interpolant at x and y + msifit - Fit the 2-D interpolant + msifree - Free space allocated by msiinit or msisinit + msigeti - Fetch an integer parameter + msigetr - Fetch a real parameter + msigrl - Calculate the integral inside a polygon + msiinit - Initialize the 2-D interpolant using default parameters + msirestore - Restore interpolant parameters and coefficients + msisave - Save 2-D interpolant parameters and coefficients + msisinit - Initialize the 2-D interpolant using user parameters + msisqgrl - Calculate the integral inside a rectangle + msitype - Decode string into interpolant type and parameters + msivector - Evaluate interpolant at an array of x and y diff --git a/math/iminterp/doc/iminterp.spc b/math/iminterp/doc/iminterp.spc new file mode 100644 index 00000000..ce5b8680 --- /dev/null +++ b/math/iminterp/doc/iminterp.spc @@ -0,0 +1,525 @@ +.help iminterp Jul84 "Math Package" + +.ce +Specifications for the Image Interpolator Package +.ce +Lindsey Davis +.ce +Vesa Junkkarinen +.ce +August 1984 + +.sh +1. Introduction + + One of the most common operations in image processing is +interpolation in a data array. Due to the large amount of data involved, +efficiency is highly important. The advantage of having locally written +interpolators, includes the ability to optimize for uniformly spaced data +and the possibility of adding features that are useful to the final +application. + +.sh +2. Requirements + +.ls (1) +The package shall take as input a one-dimensional array containing image +data. The pixels are assumed to be equally spaced along a line. +The coordinates of a pixel are assumed to be +the same as the subscript of the pixel in the data array. +The coordinate of the first pixel in the array and the spacing between pixels +is assumed to be 1.0. All pixels are assumed to be good. +Checking for INDEF valued and out of bounds pixels is the responsibility of the +user. A routine to remove INDEF valued pixels from a data array shall be +included in the package. +.le +.ls (2) +The package is divided into array sequential interpolators and array +random interpolators. The sequential interpolators have been optimized +for returning many values as is the case when an array is shifted, or +oversampled at many points in order to produce a +smooth plot. +The random interpolators allow the evaluation of a few interpolated +points without the computing time and storage overhead required for +setting up the sequential version. +.le +.ls (3) +The quality of the interpolant will be set at run time. The options are: + +.nf + II_NEAREST - nearest neighbour + II_LINEAR - linear interpolation + II_POLY3 - 3rd order divided differences + II_POLY5 - 5th order divided differences + II_SPLINE3 - cubic spline +.fi + +The calling sequences shall be invariant to the interpolant selected. +Routines should be designed so that new interpolants can be added +with minimal changes to the code and no change to the calling sequences. +.le +.ls (4) +The interpolant parameters and the arrays necessary to store the coefficients +are stored in a structure referenced by a pointer. The pointer is returned +to the user program by the initial call to ASIINIT or ASIRESTORE and freed +by a call to ASIFREE (see section 3.1). +.le +.ls (5) +The package routines shall be able to: +.ls o +Calculate the coefficients of the interpolant and store these coefficients in +the appropriate part of the interpolant descriptor structure. +.le +.ls o +Evaluate the interplant at a given x(s) coordinate(s). +.le +.ls o +Calculate the derivatives of the interpolant at a given value of x. +.le +.ls o +Integrate the interpolant over a specified x interval. +.le +.ls o +Store the interpolant in a user supplied array. Restore the saved interpolant +to the interpolant descriptor structure for later use by ASIEVAL, ASIVECTOR, +ASIDER and ASIGRL. +.le + +.sh +3. Specifications + +.sh +3.1. The Array Sequential Interpolator Routines + + The package prefix is asi and the package routines are: + +.nf + asiinit (asi, interp_type) + asifit (asi, datain, npix) + y = asieval (asi, x) + asivector (asi, x, yfit, npix) + asider (asi, x, der, nder) + v = asigrl (asi, a, b) + asisave (asi, interpolant) + asirestore (asi, interpolant) + asifree (asi) +.fi + +.sh +3.2. The Array Random Interpolator Routines + + The package prefix is ari and the package routines are: + +.nf + y = arieval (x, datain, npix, interp_type) + arider (x, datain, npix, der, nder, interp_type) +.fi + +.sh +3.3. Miscellaneous + + A routine has been included in the package to remove INDEF valued +pixels from an array. + +.nf + arbpix (datain, dataout, npix, interp_type, boundary_type) +.fi + +.sh +3.4. Algorithms + +.sh +3.4.1. Coefficients + + The coefficient array used by the asi routines is calculated by ASIFIT. +ASIFIT accepts an array of data, checks that the number +of data points is appropriate for the interpolant selected, allocates +space for the interpolant, and calculates the coefficients. +Boundary coefficient values are calculated +using boundary projection. With the exception of the cubic spline interpolant, +the coefficients are stored as the data points. +The B-spline coefficients are +calculated using natural end conditions (Prenter 1975). +After a call to ASIFIT the coefficient array contains the following. + +.nf + case II_NEAREST: + + # no boundary conditions necessary + coeff[1] = datain[1] + . + . + . + coeff[npts] = datain[npix] + + case II_LINEAR: + + # coeff[npxix+1] required if x = npix + coeff[1] = datain[1] + . + . + . + coeff[npix] = datain[npix] + coeff[npix+1] = 2. * datain[npix] - datain[npix-1] + + case II_POLY3: + + # coeff[0] required if x = 1 + # coeff[npix+1], coeff[npix+2] required if x = npix + coeff[0] = 2. * datain[1] - datain[2] + coeff[1] = datain[1] + . + . + . + coeff[npix] = datain[npix] + coeff[npix+1] = 2. * datain[npix] - datain[npix-1] + coeff[npix+2] = 2. * datain[npix] - datain[npix-2] + + case II_POLY5: + + # coeff[1], coeff[0] reqired if x = 1 + # coeff[npix+1], coeff[npix+2], coeff[npix=3] + # required if x = npix + + coeff[-1] = 2. * datain[1] - datain[3] + coeff[0] = 2. * datain[1] - datain[2] + coeff[1] = datain[1] + . + . + . + coeff[npix] = datain[npix] + coeff[npix+1] = 2. * datain[npix] - datain[npix-1] + coeff[npix+2] = 2. * datain[npix] - datain[npix-2] + coeff[npix+3] = 2. * datain[npix] - datain[npix-3] + + case SPLINE3: + + # coeff[0] = 2nd der at x = 1, coeff[0] = 0. + # coeff[npix+1] = 2nd der at x = npts, coeff[npix+1] = 0. + # coeff[npix+2] = 0., required if x = npix + coeff[0] = b[1] + coeff[1] = b[2] + . + . + . + coeff[npix] = b[npix+1] + coeff[npix+1] = b[npix+2] + coeff[npix+2] = 0. +.fi + +.sh +3.4.2. Evaluation + + The ASIEVAL and ASIVECTOR routines have been optimized to be as efficient +as possible. The values of the II_NEAREST and II_LINEAR interpolants +are calculated directly. The II_SPLINE3 interpolant is evaluated using +polynomial coefficients calculated directly from the B-spline coefficients +(de Boor 1978). Values of the higher order polynomial interpolants +are calculated using central differences. The equations for each case are +listed below. + +.nf +case II_NEAREST: + + y = coeff[int (x + 0.5)] + +case II_LINEAR: + + nx = x + y = (x - nx) * coeff[nx+1] + (nx + 1 - x) * coeff[nx] + +case II_POLY3: + + nx = x + s = x - nx + t = 1. - s + + # second central differences + cd20 = 1./6. * (coeff[nx+1] - 2. * coeff[nx] + coeff[nx-1]) + cd21 = 1./6. * (coeff[nx+2] - 2. * coeff[nx+1] + coeff[nx]) + + y = s * (coeff[nx+1] + (s * s - 1.) * cd21) + t * (coeff[nx] + + (t * t - 1.) * cd20) + +case II_POLY5: + + nx = x + s = x - nx + t = 1. - s + + # second central differences + cd20 = 1./6. * (coeff[nx+1] - 2. * coeff[nx] + coeff[nx-1]) + cd21 = 1./6. * (coeff[nx+2] - 2. * coeff[nx+1] + coeff[nx]) + + # fourth central diffreences + cd40 = 1./120. * (coeff[nx-2] - 4. * coeff[nx-1] + 6. * coeff[nx] - 4. * + coeff[nx+1] + a[nx+2]) + cd41 = 1./120. * (coeff[nx-1] - 4. * coeff[nx] + 6. * coeff[nx+1] - 4. * + coeff[nx+2] + coeff[nx+3] + + y = s * (coeff[nx+1] + (s * s - 1.) * (cd21 + (s * s - 4.) * cd41)) + + t * (coeff[nx] + (t * t - 1.) * (cd20 + (t * t - 4.) * cd40)) + +case II_SPLINE3: + + nx = x + s = x - nx + + pc[1] = coeff[nx-1] + 4. * coeff[nx] + coeff[nx+1] + pc[2] = 3. * (coeff[nx+1] - coeff[nx-1]) + pc[3] = 3. * (coeff[nx-1] - 2. * coeff[nx] + coeff[nx+1]) + pc[4] = -coeff[nx-1] + 3. * coeff[nx] - 3. * coeff[nx+1] + coeff[nx+2] + + y = pc[1] + s * (pc[2] + s * (pc[3] + s * pc[4])) +.fi + + + The ARIEVAL routine uses the expressions above to evaluate the +interpolant. However unlike ASIEVAL, ARIEVAL does not use a previously +calculated coefficient array. Instead ARIEVAL selects the appropriate +portion of the data array, calculates the coefficients and boundary +coefficients if necessary, and evaluates the interpolant at the time it +is called. The cubic spline interpolant uses at most SPLTS (currently 16) +data points to calculate the B-spline coefficients. + +.sh +3.4.3. Derivatives + + Derivatives of the interpolant are calculated by evaluating the +derivatives of the interpolating polynomial. For all interpolants der[1] +equals the value of the interpolant at x. +For the sake of efficiency the derivatives +of the II_NEAREST and II_LINEAR interpolants are calculated directly. + +.nf + case II_NEAREST: + + der[1] = coeff[int (x+0.5)] + + case II_LINEAR: + + der[1] = (x - nx) * coeff [nx+1] + (nx + 1 - x) * coeff[nx] + der[2] = coeff[nx+1] - coeff[nx] +.fi + + In order to calculate the derivatives of the cubic spline and +polynomial interpolants +the coefficients of the interpolating polynomial must be calculated. +The polynomial +coefficients for the cubic spline interpolant are computed directly from the +B-spline coefficients (see 3.4.2.). The higher order polynomial +interpolant coefficients are calculated as follows. + +First the appropriate portion of the coefficient array is loaded. + +.nf + do i = 1, nterms + d[i] = coeff[nx - nterms/2 + i] +.fi + +Next the divided differences are calculated (Conte and de Boor 1972). + +.nf + do k = 1, nterms - 1 + do i = 1, nterms - k + d[i] = (d[i+1] - d[i]) / k +.fi + +The d[i] are the coefficients of an interpolating polynomial of the +following form. The x[i] are the nterms data points surrounding the +point of interest. + +.nf + p(x) = d[1] * (x-x[1]) * ... * (x-x[nterms-1) + + d[2] * (x-x[2]) * ... * (x-x[nterms-1]) + ... + d[nterms] +.fi + +Next a new set of polynomial coefficients are calculated +(Conte and de Boor 1972). + +.nf + do k = nterms, 2, -1 + do i = 2, k + d[i] = d[i] + d[i-1] * (k - i - nterms/2) +.fi + +The new d[i] are the coefficients of the follwoing polynomial. + +.nf + nx = x + p(x) = d[1] * (x-nx) ** (nterms-1) + d[2] * (x-nx) ** (nterms-2) + ... + d[nterms] +.fi + +The d[i] array is flipped. The value and derivatives +of the interpolant are then calculated using the d[i] array and +nested multiplication. + +.nf + s = x - nx + + do k = 1, nder { + + accum = d[nterms-k+1] + + do j = nterms - k, 1, -1 + accum = d[j] + s * accum + + der[k] = accum + + # differnetiate + do j = 1, nterms - k + d[j] = j * d[j + 1] + } +.fi + + ARIDER calculates the derivatives of the interpolant using the same +technique ASIDER. However ARIDER does not use a previously calculated +coefficient array like ASIDER. Instead ARIDER selects the appropriate portion +of the data array, calculates the coefficients and boundary coefficients, +and computes the derivatives at the time it is called. + +.sh +3.4.5. Integration + + ASIGRL calculates the integral of the interpolant between fixed limits +by integrating the interpolating polynomial. The coefficients of the +interpolating polynomial are calculated as discussed in section 3.4.4. + +.sh +4. Usage + +.sh +4.1. User Notes + +The following series of steps illustrates the use of the package. + +.ls 4 +.ls (1) +Insert an include statement in the calling program to make +the IINTERP definitions available to the user program. +.le +.ls (2) +Remove INDEF valued pixels from the data using ARBPIX. +.le +.ls (3) +Call ASIINIT to initialize the interpolant parameters. +.le +.ls (4) +Call ASIFIT to calculate the coefficients of the interpolant. +.le +.ls (5) +Evaluate the interpolant at a given value of x(s) using ASIEVAL or +ASIVECTOR. +.le +.ls (6) +Calculate the derivatives and integral or the interpolant using +ASIDER and ASIGRL. +.le +.ls (7) +Free the interpolator structure by calling ASIFREE. +.le +.le + + The interpolant can be saved and restored using ASISAVE and ASIRESTORE. +If the values and derivatives of only a few points in an array are desired +ARIEVAL and ARIDER can be called. + +.sh +4.2. Examples + +.nf +Example 1: Shift a data array by a constant amount + + include + ... + call asiinit (asi, II_POLY5) + call asifit (asi, inrow, npix) + + do i = 1, npix + outrow[i] = asieval (asi, i + shift) + + call asifree (asi) + ... + +Example 2: Calculate the integral under the data array + + include + ... + call asiinit (asi, II_POLY5) + call asifit (asi, datain, npix) + + integral = asigrl (asi, 1. real (npix)) + + call asifree (asi) + ... + +Example 2: Store interpolant for later use by ASIEVAL + LEN_INTERP must be at least npix + 8 units long where npix is + is defined in the call to ASIFIT. + + include + + real interpolant[LEN_INTERP] + ... + call asiinit (asi, II_POLY3) + call asifit (asi, datain, npix) + call asisave (asi, interpolant) + call asifree (asi) + ... + call asirestore (asi, interpolant) + do i = 1, npts + yfit[i] = asieval (asi, x[i]) + call asifree (asi) + ... +.fi +.sh +5. Detailed Design + +.sh +5.1. Interpolator Descriptor + + The interpolant parameters and coefficients are stored in a +structure listed below. + +.nf + define LEN_ASISTRUCT 4 # Length in structure units of + # interpolant descriptor + + define ASI_TYPE Memi[$1] # Interpolant type + define ASI_NCOEFF Memi[$1+1] # No. of coefficients + define ASI_OFFSET Memi[$1+2] # First "data" point in + # coefficient array + define ASI_COEFF Memi[$1+3] # Pointer to coefficient array +.fi + +.sh +5.2. Storage Requirements + + The interpolant descriptor requires LEN_ASISTRUCT storage units. The +coefficient array requires ASI_NCOEFF(asi) real storage elements, where +ASI_NCOEFF(asi) is defined as follows. + +.nf + ASI_NCOEFF(asi) = npix # II_NEAREST + ASI_NCOEFF(asi) = npix+1 # II_LINEAR + ASI_NCOEFF(asi) = npix+3 # II_POLY3 + ASI_NCOEFF(asi) = npix+5 # II_POLY5 + ASI_NCOEFF(asi) = npix+3 # II_SPLINE3 +.fi + +.sh +6. References + +.ls (1) +Carl de Boor, "A Practical Guide to Splines", 1978, Springer-Verlag New +York Inc. +.le +.ls (2) +S.D. Conte and C. de Boor, "Elementary Numerical Analysis", 1972, McGraw-Hill, +Inc. +.le +.ls (3) +P.M. Prenter, "Splines and Variational Methods", 1975, John Wiley and Sons Inc. +.le +.endhelp diff --git a/math/iminterp/doc/mrider.hlp b/math/iminterp/doc/mrider.hlp new file mode 100644 index 00000000..47515c48 --- /dev/null +++ b/math/iminterp/doc/mrider.hlp @@ -0,0 +1,79 @@ +.help mrider Dec98 "Image Interpolation Package" +.ih +NAME +mrider -- calculate the derivatives at x and y +.ih +SYNOPSIS +include + +.nf +mrider (x, y, datain, nxpix, nypix, len_datain, der, nxder, nyder, len_der, + interp_type) +.fi + +.nf +real x[4] #I x value, 1. <= x[1-4] <= nxpix +real y[4] #I y value, 1. <= y[1-4] <= nypix +real datain[len_datain, ARB] #I data array +int nxpix #I number of data pixels in x +int nypix #I number of data pixels in y +int len_datain #I length of datain, len_datain >= nxpix +real der[len_der, ARB] #O derivative array +int nxder #I x order of the derivatives +int nyder #I y order of the derivatives +int len_der #I row length of der, len_der >= nxder +int interp_type #I interpolant type +.fi +.ih +ARGUMENTS +.ls x, y +The single x and y points or in the case of the drizzle interpolant the +single quadrilateral at / over which the derivatives are to be evaluated. +The quadrilateral vertices may be stored in clock-wise or counter-clockwise +order. +.le +.ls datain +Array of data values. +.le +.ls nxpix, nypix +The number of data values in the x and y directions +.le +.ls len_datain +The row length of the datain array. Len_datain must be >= nxpix. +.le +.ls der +The derivative array. Der[1,1] equals the function value at x and y and +der[2,1], der[1,2] are the first derivatives with respect to x and y +respectively. +.le +.ls nxder, nyder +The number of the derivatives in x and y to be returned. MRIDER checks +that the requested number of derivatives is sensible. The sinc interpolants +return the interpolant value and all the first and second order derivatives. +The drizzle interpolant returns the interpolant value and the first +derivative in x and y. +.le +.ls len_der +The row length of the derivative array. Len_der must be >= nxder. +.le +.ls interp_type +Interpolant type. The options are II_BINEAREST, II_BILINEAR, II_BIPOLY3, +II_BIPOLY5, II_BISPLINE3, II_SINC / II_LSINC, and II_DRIZZLE. The look-up +table sinc is not supported and defaults to the sinc interpolant. The +interpolant width is 31 pixels. The drizzle pixel fraction is 1.0. The +interpolant type definitions are found in the package header file +math/iminterp.h. +.le +.ih +DESCRIPTION +MRIDER is useful for evaluating the function and derivatives at a few +widely spaced points in a data array without the storage space required +by the sequential version. +.ih +NOTES +Checking for out of bounds and INDEF valued pixels is the +responsibility of the user. +.ih +SEE ALSO +msider +.endhelp diff --git a/math/iminterp/doc/mrieval.hlp b/math/iminterp/doc/mrieval.hlp new file mode 100644 index 00000000..35477614 --- /dev/null +++ b/math/iminterp/doc/mrieval.hlp @@ -0,0 +1,57 @@ +.help mrieval Dec98 "Image Interpolation Package" +.ih +NAME +mrieval -- evaluate the interpolant at x and y +.ih +SYNOPSIS +include + +y = mrieval (x, y, datain, nxpix, nypix, len_datain, interp_type) + +.nf + real x[4] #I x value, 1 <= x[1-4] <= nxpix + real y[4] #I y value, 1 <= y[1-4] <= nypix + real datain[len_datain, ARB] #I data array + int nxpix #I number of x values + int nypix #I number of y values + int len_datain #I length datain, len_datain >= nxpix + int interp_type #I interpolant type +.fi +.ih +ARGUMENTS +.ls x, y +The single x and y values or in the case of the drizzle interpolant the +single quadrilateral at / over which the interpolant is to be evaluated. +The vertices of the quadilateral must be defined in clock-wise or +counter-clockwise order. +.le +.ls datain +The array of data values. +.le +.ls nxpix, nypix +The number of data pixels in x and y. +.le +.ls len_datain +The row length of datain. Len_datain must be >= nxpix. +.le +.ls interp_type +Interpolant type. The options are II_BINEAREST, II_BILINEAR, II_BIPOLY3, +II_BIPOLY5, II_BISPLINE3, II_SINC / II_LSINC, and II_DRIZZLE. The look-up +table sinc interpolant is not supported and defaults to the sinc interpolant. +The sinc interpolant width is 31 pixels. The drizzle pixel fraction is 1.0. +The interpolant type definitions are found in the package header file +math/iminterp.h. +.le +.ih +DESCRIPTION +MRIEVAL is useful for evaluating the interpolant at a few selected points +in the datain array without the storage overhead required for the sequential +version. +.ih +NOTES +Checking for INDEF valued or out of bounds pixels is the +responsibility of the user. +.ih +SEE ALSO +msieval, msivector, mrider +.endhelp diff --git a/math/iminterp/doc/msider.hlp b/math/iminterp/doc/msider.hlp new file mode 100644 index 00000000..0139c0a0 --- /dev/null +++ b/math/iminterp/doc/msider.hlp @@ -0,0 +1,52 @@ +.help msider Dec98 "Image Interpolation Package" +.ih +NAME +msider -- evaluate the interpolant derivatives at x and y +.ih +SYNOPSIS +msider (msi, x, y, der, nxder, nyder, len_der) + +.nf + pointer msi #I interpolant descriptor + real x[4] #I x value, 1 <= x[1-4] <= nxpix + real y[4] #I y value, 1 <= y[1-4] <= nypix + real der[len_der, ARB] #O derivative array + int nxder #I number of x derivatives + int nyder #I number of y derivatives + int len_der #I row length of der, len_der >= nxder +.fi +.ih +ARGUMENTS +.ls msi +Pointer to the 2D sequential interpolant descriptor. +.le +.ls x, y +The single x and y values or in the case of the drizzle interpolant the +single quadrilateral at / over which the point is to be evaluated. +.le +.ls der +The array containing the derivatives. Der[1,1] contains the value of +the interpolant at x and y. Der[2,1] and der[1,2] contain the 1st +derivatives of x and y respectively. +.le +.ls nxder, nyder +The number derivatives in x and y. +.le +.ls len_der +The row length of der. Len_der must be >= nxder. +.le +.ih +DESCRIPTION +The polynomial and spline interpolants are evaluated using the polynomial +coefficients and nested multiplication. The polynomial interpolants are +stored as the data points. The spline interpolant is stored as a set of +B-spline coefficients. +.ih +NOTES +MRIDER checks that the number of derivatives requested is reasonable. +Checking for out of bounds and INDEF valued pixels is the responsibility of the +user. MSIINIT and MSIFIT must be called before using MSIDER. +.ih +SEE ALSO +mrider +.endhelp diff --git a/math/iminterp/doc/msieval.hlp b/math/iminterp/doc/msieval.hlp new file mode 100644 index 00000000..9a77c006 --- /dev/null +++ b/math/iminterp/doc/msieval.hlp @@ -0,0 +1,46 @@ +.help msieval Dec98 "Image Interpolation Package" +.ih +NAME +msieval -- procedure to evaluate the interpolant at x and y +.ih +SYNOPSIS +z = msieval (msi, x, y) + +.nf +pointer msi #I interpolant descriptor +real x[4] #I x value, 1 <= x[1-4] <= nxpix +real y[4] #I y value, 1 <= y[1-4] <= nypix +.fi +.ih +ARGUMENTS +.ls msi +The pointer to the sequential interpolant descriptor structure. +.le +.ls x, y +The single x and y values of or in the case of the drizzle interpolant +the single quadrilateral over which the point is to be evaluated. +.le +.ih +DESCRIPTION +The polynomial coefficients are calculated from the data points in the +case of the polynomial interpolants and the B-spline coefficients in +the case of the spline interpolant. The polynomial interpolants +are evaluated using Everett's central difference formula. The boundary +extension algorithm is projection. + +The sinc interpolant is evaluated using an array of data points around +the desired position. The look-up table sinc interpolant is computed +using an a pre-computed look--up table entry. The boundary extension +algorithm is nerest neighbor. + +The drizzle interpolant is computed by computing the mean value of the +data within the user supplied quadrilateral. +.ih +NOTES +Checking for out of bounds and INDEF valued pixels is the responsibility of +the user. MSIINIT or MSISINIT and MSIFIT must be called before calling +MSIEVAL. +.ih +SEE ALSO +msivector, mrieval, mrider +.endhelp diff --git a/math/iminterp/doc/msifit.hlp b/math/iminterp/doc/msifit.hlp new file mode 100644 index 00000000..9a63ae2d --- /dev/null +++ b/math/iminterp/doc/msifit.hlp @@ -0,0 +1,45 @@ +.help msifit Dec98 "Image Interpolation Package" +.ih +NAME +msifit - fit the interpolant to the data +.ih +SYNOPSIS +msifit (msi, datain, nxpix, nypix, len_datain) + +.nf + pointer msi #I interpolant descriptor + real datain[len_datain,ARB] #I data array + int nxpix #I number of x pixels + int nypix #I number of y pixels + int len_datain #I length of datain, len_datain >= nxpix +.fi +.ih +ARGUMENTS +.ls msi +Pointer to the sequential interpolant descriptor. +.le +.ls datain +Array containing the data. +.le +.ls nxpix, nypix +The number of pixels in x and y. +.le +.ls len_datain +The row length of the datain array. Len_datain must be >= nxpix. +.le +.ih +DESCRIPTION +The datain array is checked for size, memory is allocated for the coefficient +array and the end conditions are specified. The interior polynomial, sinc, +and drizzle interpolants are saved as the data points. The polynomial +coefficients are calculated from the data points in the evaluation stage. +The B-spline coefficients are calculated in MSIFIT as they depend on the +entire data array. +.ih +NOTES +Checking for INDEF valued pixels is the responsibility of the user. +MSIINIT or MSISINIT must be called before using MSIFIT. MSIFIT must be +called before using MSIEVAL, MSIVECTOR, MSIDER, MSIGRL or MSISQGRL. +.ih +SEE ALSO +.endhelp diff --git a/math/iminterp/doc/msifree.hlp b/math/iminterp/doc/msifree.hlp new file mode 100644 index 00000000..79b1966d --- /dev/null +++ b/math/iminterp/doc/msifree.hlp @@ -0,0 +1,26 @@ +.help msifree Dec98 "Image Interpolation Package" +.ih +NAME +msifree -- free sequential interpolant descriptor +.ih +SYNOPSIS +msifree (msi) + +.nf + pointer msi #U interpolant descriptor +.fi +.ih +ARGUMENTS +.ls msi +Pointer to the sequential interpolant descriptor structure. +.le +.ih +DESCRIPTION +MSIFREE frees the sequential interpolant descriptor structure. +MSIFREE should be called when interpolation is complete. +.ih +NOTES +.ih +SEE ALSO +msiinit, msisinit +.endhelp diff --git a/math/iminterp/doc/msigeti.hlp b/math/iminterp/doc/msigeti.hlp new file mode 100644 index 00000000..5ab46156 --- /dev/null +++ b/math/iminterp/doc/msigeti.hlp @@ -0,0 +1,35 @@ +.help msigeti Dec98 msigeti.hlp +.ih +NAME +msigeti -- fetch an msi integer parameter +.ih +SYNOPSIS +include + +ivalue = msigeti (msi, param) + +.nf + pointer msi #I interpolant descriptor + int param #I parameter +.fi +.ih +ARGUMENTS +.ls msi +Pointer to the sequential interpolant descriptor structure. +.le +.ls param +The parameter to be fetched. The choices are: II_MSITYPE, the interpolant +type, II_MSINSAVE, the length of the saved coefficient array, and +II_MSINSINC, the half-width of the sinc interpolant. +.le +.ih +DESCRIPTION +MSIGETI is used to determine the size of the coefficient array that +must be allocated to save the sequential interpolant description structure, +and to fetch selected sequential interpolant parameters. +.ih +NOTES +.ih +SEE ALSO +msiinit, msisinit, msigetr +.endhelp diff --git a/math/iminterp/doc/msigetr.hlp b/math/iminterp/doc/msigetr.hlp new file mode 100644 index 00000000..dadac5c2 --- /dev/null +++ b/math/iminterp/doc/msigetr.hlp @@ -0,0 +1,37 @@ +.help msigetr Dec98 msigetr.hlp +.ih +NAME +msigetr -- fetch an msi real parameter +.ih +SYNOPSIS +include + +rvalue = msigetr (msi, param) + +.nf + pointer msi #I interpolant descriptor + int param #I parameter +.fi +.ih +ARGUMENTS +.ls msi +Pointer to the sequential interpolant descriptor structure. +.le +.ls param +The parameter to be fetched. The choices are: II_MSIBADVAL, the undefined +pixel value for the drizzle interpolant. The parameter definitions are +contained in the package header file math/iminterp.h. +.le +.ih +DESCRIPTION +MSIGETR is used to set the value of undefined drizzle interpolant pixels. +Undefined pixels are those for which the interpolation coordinates do not +overlap the input coordinates, but are still, within the boundaries of the input +image, a situation which may occur when the pixel fraction is < 1.0. +.ih +.ih +NOTES +.ih +SEE ALSO +msiinit, msisinit, msigeti +.endhelp diff --git a/math/iminterp/doc/msigrid.hlp b/math/iminterp/doc/msigrid.hlp new file mode 100644 index 00000000..6bc110d5 --- /dev/null +++ b/math/iminterp/doc/msigrid.hlp @@ -0,0 +1,51 @@ +.help msigrid Dec98 "Image Interpolation Package" +.ih +NAME +msigrid -- evaluate the interpolant on a grid of points +.ih +SYNOPSIS +msigrid (msi, x, y, zfit, nx, ny, len_zfit) + +.nf + pointer msi #I interpolant descriptor + real x[2*nx] #I x values, 1 <= x[i] <= nx + real y[2*ny] #I y values, 1 <= y[i] <= ny + real zfit[len_zfit,ARB] #O grid of interpolated values + int nx #I number of x points + int ny #I number of y points + int len_zfit #I length zfit, len_zfit >= nx +.fi +.ih +ARGUMENTS +.ls msi +Pointer to the interpolant descriptor structure. +.le +.ls x, y +The x and y arrays of points to be evaluated, or in the case of the drizzle +interpolant the x and y ranges over which the points are to be evaluated. +The x and y arrays must be ordered in increasing values of x and y respectively. +.le +.ls zfit +The array of interpolated points. +.le +.ls nx, ny +The number of points in the x and y directions respectively. +.le +.ls len_zfit +The row length of the zfit array. Len_zfit >= nx. +.le +.ih +DESCRIPTION +MSIGRID evaluates the interpolant at a set of x and y values on a +rectangular grid or in the case of the drizzle interpolant within +rectangular regions. It is most efficient for evaluating the interpolant +at many values which are closely spaced in x and y. For widely spaced +points MSIVECTOR should be used. +.ih +NOTES +Checking for out of bounds and INDEF valued pixels is the responsibility +of the user. +.ih +SEE ALSO +msieval, msivector, msider, mrieval, mrider +.endhelp diff --git a/math/iminterp/doc/msigrl.hlp b/math/iminterp/doc/msigrl.hlp new file mode 100644 index 00000000..f6e2326a --- /dev/null +++ b/math/iminterp/doc/msigrl.hlp @@ -0,0 +1,43 @@ +.help msigrl Dec98 "Image Interpolation Package" +.ih +NAME +msigrl -- integrate the interpolant inside a polygon +.ih +SYNOPSIS +y = msigrl (msi, x, y, npts) + +.nf + pointer msi #I interpolant descriptor + real x[npts] #I x values, 1 <= x <= npts, x[1] = x[npts] + real y[npts] #I y values, 1 <= y <= npts, y[1] = y[npts] + int npts #I number of points +.fi +.ih +ARGUMENTS +.ls msi +Pointer to the sequential interpolant descriptor structure. +.le +.ls x, y +An array of x and y values describing a polygon, where x[1] = x[npts] and +y[1] = y[npts]. X and y describe a closed curve where any horizontal line +segment intersects the domain of integration at at most one point. +.le +.ls npts +The number of points describing the polygon. Npts must >= 4 (triangle). +.le +.ih +DESCRIPTION +MSIGRL integrates the interpolant exactly for rectangular domains +of integration. For more irregular regions of integration MSIGRL +returns an approximation whose accuracy depends on the size of the +integration region and the shape of the polygon. +.ih +NOTES +Checking for out of bound integration regimes is the responsibility of +the user. Non-rectangular partial pixel domains of integration default +to rectangular regions. MSIINIT or MSISINIT and MSIFIT must be called +before using MSIGRL. +.ih +SEE ALSO +msisqrgl +.endhelp diff --git a/math/iminterp/doc/msiinit.hlp b/math/iminterp/doc/msiinit.hlp new file mode 100644 index 00000000..68dba684 --- /dev/null +++ b/math/iminterp/doc/msiinit.hlp @@ -0,0 +1,41 @@ +.help msiinit Dec98 "Image Interpolation Package" +.ih +NAME +msiinit -- initialize the sequential interpolant descriptor +.ih +SYNOPSIS +include + +msiinit (msi, interp_type) + +.nf + pointer msi #U interpolant descriptor + int interp_type #I interpolant type +.fi +.ih +ARGUMENTS +.ls msi +Pointer to the sequential interpolant descriptor. +.le +.ls interp_type +Interpolant type. The options are II_BINEAREST, II_BILINEAR, II_BIPOLY3, +II_BIPOLY5, II_BISPLINE3, II_BISINC, II_BILSINC, and II_BIDRIZZLE, for +nearest neighbour, bilinear, 3rd and 5th order interior polynomials, bicubic +spline, sinc, look-up table sinc, and drizzle respectively. The interpolant +definitions are found in the package header file math/iminterp.h. +.le +.ih +DESCRIPTION +The interpolant type is allocated and initialized. The pointer msi is +returned by MSIINIT. The sinc interpolant width defaults to 31 pixels +in x and y. The look-up table sinc resolution defaults to 20 resolution +elements or 0.05 pixels in x and y. The drizzle pixel fraction defaults +to 1.0. +.ih +NOTES +MSIINIT, MSISINIT or MSIRESTORE must be called before using any other +MSI routines. +.ih +SEE ALSO +msirestore, msifree +.endhelp diff --git a/math/iminterp/doc/msirestore.hlp b/math/iminterp/doc/msirestore.hlp new file mode 100644 index 00000000..664c9b50 --- /dev/null +++ b/math/iminterp/doc/msirestore.hlp @@ -0,0 +1,36 @@ +.help msirestore Dec98 "Image Interpolator Package" +.ih +NAME +msirestore -- restore interpolant +.ih +SYNOPSIS +msirestore (msi, interpolant) + +.nf + pointer msi #U interpolant descriptor + real interpolant[] #I array containing interpolant +.fi +.ih +ARGUMENTS +.ls msi +Pointer to the interpolant descriptor structure. +.le +.ls interpolant +Array containing the interpolant. The amount of space required by interpolant +can be determined by a call to msigeti. +.le + +.nf + len_interpolant = msigeti (msi, II_MSINSAVE) +.fi +.ih +DESCRIPTION +MSIRESTORE allocates space for the interpolant descriptor and restores the +parameters and coefficients stored in the interpolant array to the +interpolant structure for use by MSIEVAL, MSIVECTOR, MSIDER and MSIGRL. +.ih +NOTES +.ih +SEE ALSO +msisave +.endhelp diff --git a/math/iminterp/doc/msisave.hlp b/math/iminterp/doc/msisave.hlp new file mode 100644 index 00000000..7d5f67ef --- /dev/null +++ b/math/iminterp/doc/msisave.hlp @@ -0,0 +1,38 @@ +.help msisave Dec98 "Image Interpolator Package" +.ih +NAME +msisave -- save interpolant +.ih +SYNOPSIS +msisave (msi, interpolant) + +.nf + pointer msi #I interpolant descriptor + real interpolant[] #O array containing the interpolant +.fi +.ih +ARGUMENTS +.ls msi +Pointer to the interpolant descriptor structure. +.le +.ls interpolant +Array where the interpolant is stored. The required interpolant array length +required can be determined by a call to msigeti. +.le + +.nf + len_interpolant = msigeti (msi, II_MSINSAVE) +.fi +.ih +DESCRIPTION +The interpolant type, number of coefficients in x and y, the position of +the first data point in the coefficient array, and the sinc and drizzle +interpolant parameters are stored in the first eleven elements of interpolant. +The remaining elements contain the coefficients and look-up tables +calculated by MSIFIT. +.ih +NOTES +.ih +SEE ALSO +msirestore +.endhelp diff --git a/math/iminterp/doc/msisinit.hlp b/math/iminterp/doc/msisinit.hlp new file mode 100644 index 00000000..0a05e7ab --- /dev/null +++ b/math/iminterp/doc/msisinit.hlp @@ -0,0 +1,61 @@ +.help msisinit Dec98 "Image Interpolator Package" +.ih +NAME +msisinit -- initialize the sequential interpolant descriptor using user parameters +.ih +SYNOPSIS +include + +msisinit (msi, interp_type, nsinc, nxincr, nyincr, pixfrac1, pixfrac2, badval) + +.nf + pointer msi #O interpolant descriptor + int interp_type #I interpolant type + int nsinc #I sinc interpolant width in pixels + int nxincr,nyincr #I sinc look-up table resolution + real pixfrac1,pixfrac2 #I sinc or drizzle pixel fractions + real badval #I drizzle undefined pixel value +.fi + +.ih +ARGUMENTS +.ls msi +Pointer to sequential interpolant descriptor. +.le +.ls interp_type +Interpolant type. The options are II_BINEAREST, II_BILINEAR, II_BIPOLY3, +II_BIPOLY5, II_BISPLINE3, II_BISINC, II_BILSINC, and II_BIDRIZZLE for the +nearest neighbour, linear, 3rd order polynomial, 5th order polynomial, +cubic spline, sinc, look-up table sinc, and drizzle interpolants respectively. +The interpolant type definitions are found in the package header file +math/iminterp.h. +.le +.ls nsinc +The sinc and look-up table sinc interpolant width in pixels. Nsinc is +rounded up internally to the nearest odd number. +.le +.ls nxincr, nyincr +The look-up table sinc resolution in x and y in number of entries. Nxincr = 10 +implies a pixel resolution of 0.1 pixels in x, nxincr = 20 a pixel resolution +of 0.05 pixels in x, etc. The default value of nxincr and nyincr are 20 and 20 +.le +.ls pixfrac1, pixfrac2 +The look-up table sinc fractional pixel shifts in x and y if nincr = 1 in +which case -0.5 <= rparam[1/2] <= 0.5 , or the drizzle pixel fractions in +which case 0.0 <= rparam[1/2] <= 1.0. +.le +.ls badval +The undefined pixel value for the drizzle interpolant. Pixels within +the boundaries of the input image which do not overlap the input image +pixels are assigned a value of badval. +.le +.ih +DESCRIPTION +The interpolant descriptor is allocated and initialized. The pointer msi is +returned by MSISINIT. +.ih +NOTES +MSIINIT or MSISINIT must be called before using any other MSI routines. +.ih +SEE ALSO +msisinit, msifree diff --git a/math/iminterp/doc/msisqgrl.hlp b/math/iminterp/doc/msisqgrl.hlp new file mode 100644 index 00000000..fc49514d --- /dev/null +++ b/math/iminterp/doc/msisqgrl.hlp @@ -0,0 +1,38 @@ +.help msisqgrl Dec98 "Image Interpolation Package" +.ih +NAME +msisqgrl -- integrate the interpolant over a rectangular region +.ih +SYNOPSIS +y = msisqgrl (msi, x1, x2, y1, y2) + +.nf + pointer msi #I interpolant descriptor + real x1 #I lower x limit, 1 <= x1 <= nxpix + real x2 #I upper x limit, 1 <= x2 <= nxpix + real y1 #I lower y limit, 1 <= y1 <= nypix + real y2 #I upper y limit, 1 <= y2 <= nypix +.fi +.ih +ARGUMENTS +.ls msi +Pointer to the sequential interpolant descriptor structure. +.le +.ls x1, x2 +The x limits of integration +.le +.ls y1, y2 +The y limits of integration. +.le +.ih +DESCRIPTION +MSISQGRL integrates the interpolant exactly for rectangular domains +of integration, including partial pixel regions. +.ih +NOTES +Checking for out of bound integration regimes is the responsibility of +the user. MSIFIT must be called before using MSISQGRL. +.ih +SEE ALSO +msigrl +.endhelp diff --git a/math/iminterp/doc/msitype.hlp b/math/iminterp/doc/msitype.hlp new file mode 100644 index 00000000..5da6dc6d --- /dev/null +++ b/math/iminterp/doc/msitype.hlp @@ -0,0 +1,95 @@ +.help msitype Dec98 "Image Interpolator Package" +.ih +NAME +msitype -- decode an interpolant string +.ih +SYNOPSIS +include + +msitype (interpstr, interp_type, nsinc, nincr, pixfrac) + +.nf + char interpstr #I the input interpolant string + int interp_type #O interpolant type + int nsinc #O sinc interpolant width in pixels + int nincr #O sinc look-up table resolution + real pixfrac #O sinc or drizzle pixel fraction +.fi + +.ih +ARGUMENTS +.ls interpstr +The user supplied interpolant string to be decoded. The options are +.ls nearest +Nearest neighbor interpolation. +.le +.ls linear +Bilinear interpolation +.le +.ls poly3 +Bicubic polynomial interpolation. +.le +.ls poly5 +Biquintic polynomial interpolation. +.le +.ls spline3 +Bicubic spline interpolation. +.le +.ls sinc +2D sinc interpolation. Users can specify the sinc interpolant width by +appending a width value to the interpolant string, e.g. sinc51 specifies +a 51 by 51 pixel wide sinc interpolant. The sinc width will be rounded up to +the nearest odd number. The default sinc width is 31 by 31. +.le +.ls lsinc +Look-up table sinc interpolation. Users can specify the look-up table sinc +interpolant width by appending a width value to the interpolant string, e.g. +lsinc51 specifies a 51 by 51 pixel wide look-up table sinc interpolant. The user +supplied sinc width will be rounded up to the nearest odd number. The default +sinc width is 31 by 31 pixels. Users can specify the resolution of the lookup +table sinc by appending the look-up table size in square brackets to the +interpolant string, e.g. lsinc51[20] specifies a 20 by 20 element sinc +look-up table interpolant with a pixel resolution of 0.05 pixels in x and y. +The default look-up table size and resolution are 20 by 20 and 0.05 pixels +in x and y respectively. +.le +.ls drizzle +Drizzle interpolation. Users can specify the drizzle pixel fraction by +appending the pixel fraction value to the interpolant string in square +brackets, e.g. drizzle[0.5] specifies a pixel fraction of 0.5 in x and y. +The default pixel fraction is 1.0. If either of the x or y pixel +fractions are 0.0, then both are set to 0.0. A minimum value of 0.001 +is imposed on the actual value of pixfrac. +.le +.le +.ls interp_type +The output interpolant type. The options are II_BINEAREST, II_BILINEAR, +II_BIPOLY3, II_BIPOLY5, II_BISPLINE3, II_BISINC, II_BILSINC, and II_BIDRIZZLE +for the nearest neighbor, linear, 3rd order polynomial, 5th order polynomial, +cubic spline, sinc, look-up table sinc, and drizzle interpolants respectively. +The interpolant type definitions are found in the package header file +math/iminterp.h. +.le +.ls nsinc +The output sinc and look-up table sinc interpolant width in pixels. The +default value is 31 pixels in x and y. +.le +.ls nincr +The output sinc look-up table size. Nincr = 10 implies a pixel resolution +of 0.1 pixels in x, nincr = 20 a pixel resolution of 0.05 pixels in y, etc. The +default value of nincr is 20. +.le +.ls pixfrac +The output look-up table sinc fractional pixel shift if nincr = 1 +in which case -0.5 <= pixfrac <= 0.5 , or the drizzle pixel +fraction in which case 0.0 <= pixfrac <= 1.0. +.le +.ih +DESCRIPTION +The interpolant string is decoded into values suitable for the MSISINIT +or MSIINIT routines. +.ih +NOTES +.ih +SEE ALSO +msinit, msisinit, msifree diff --git a/math/iminterp/doc/msivector.hlp b/math/iminterp/doc/msivector.hlp new file mode 100644 index 00000000..d6561be7 --- /dev/null +++ b/math/iminterp/doc/msivector.hlp @@ -0,0 +1,54 @@ +.help msivector Dec98 "Image Interpolation Package" +.ih +NAME +msivector -- evaluate the interpolant at an array of x and y points +.ih +SYNOPSIS +msivector (msi, x, y, zfit, npts) + +.nf + pointer msi #I interpolant descriptor + real x[npts/4*npts] #I x values, 1 <= x <= nxpix + real y[npts/4*npts] #I y values, 1 <= y <= nypix + real zfit[npts] #O interpolated values + int npts #I number of points +.fi +.ih +ARGUMENTS +.ls msi +The pointer to the sequential interpolant descriptor +.le +.ls x, y +The array of x and y values at or in the case tof the drizzle interpolant the +array of quadrilaterals over which to evaluate the interpolant. +.le +.ls zfit +The interpolated values. +.le +.ls npts +The number of points. +.le +.ih +DESCRIPTION +The polynomial coefficients are calculated directly from the data points, +The polynomial interpolants are evaluated using Everett's central difference +formula. The spline interpolant uses the B-spline coefficients +calculated using the MSIFIT routine. The boundary extension algorithm is +projection. + +The sinc interpolant is evaluated using a array of data points around +the point in question. The look-up table since is computed by convolving +the data with a pre-computed look-up table entry. The boundary extension +algorithm is nearest neighbor. + +The drizzle interpolant is evaluated by summing the data over the +list of user supplied quadrilaterals. +.ih +NOTES +Checking for out of bounds and INDEF valued pixels is the responsibility +of the user. MSIINIT or MSISINIT and MSIFIT must be called before using +MSIVECTOR. +.ih +SEE ALSO +msieval, mrieval +.endhelp -- cgit