diff options
Diffstat (limited to 'math/iminterp/doc/im1dinterp.spc')
| -rw-r--r-- | math/iminterp/doc/im1dinterp.spc | 525 |
1 files changed, 525 insertions, 0 deletions
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 <iminterp.h> 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 <iminterp.h> + ... + 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 <iminterp.h> + ... + 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 <iminterp.h> + + 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 |