diff options
Diffstat (limited to 'math/interp')
| -rw-r--r-- | math/interp/Iinterp.hlp | 293 | ||||
| -rw-r--r-- | math/interp/README | 7 | ||||
| -rw-r--r-- | math/interp/arbpix.x | 203 | ||||
| -rw-r--r-- | math/interp/arider.x | 214 | ||||
| -rw-r--r-- | math/interp/arival.x | 124 | ||||
| -rw-r--r-- | math/interp/asidef.h | 16 | ||||
| -rw-r--r-- | math/interp/asider.x | 121 | ||||
| -rw-r--r-- | math/interp/asieva.x | 38 | ||||
| -rw-r--r-- | math/interp/asifit.x | 75 | ||||
| -rw-r--r-- | math/interp/asigrl.x | 201 | ||||
| -rw-r--r-- | math/interp/asiset.x | 20 | ||||
| -rw-r--r-- | math/interp/asival.x | 49 | ||||
| -rw-r--r-- | math/interp/bench.x | 55 | ||||
| -rw-r--r-- | math/interp/cubspl.f | 119 | ||||
| -rw-r--r-- | math/interp/iieval.x | 137 | ||||
| -rw-r--r-- | math/interp/iif_spline.x | 67 | ||||
| -rw-r--r-- | math/interp/iimail | 213 | ||||
| -rw-r--r-- | math/interp/iipol_terp.x | 41 | ||||
| -rw-r--r-- | math/interp/interp.h | 19 | ||||
| -rw-r--r-- | math/interp/interpdef.h | 19 | ||||
| -rw-r--r-- | math/interp/mkpkg | 21 | ||||
| -rw-r--r-- | math/interp/noteari | 64 | ||||
| -rw-r--r-- | math/interp/noteasi | 101 | ||||
| -rw-r--r-- | math/interp/notes3 | 216 | ||||
| -rw-r--r-- | math/interp/overview | 43 | ||||
| -rw-r--r-- | math/interp/usernote | 118 |
26 files changed, 2594 insertions, 0 deletions
diff --git a/math/interp/Iinterp.hlp b/math/interp/Iinterp.hlp new file mode 100644 index 00000000..527f0375 --- /dev/null +++ b/math/interp/Iinterp.hlp @@ -0,0 +1,293 @@ + +.help INTERP Jan83 "Image Interpolation Package" +.sh +Introduction + + One of the most common operations in image processing is interpolation. +Due to the very large amount of data involved, efficiency is quite important, +and the general purpose routines available from many sources cannot be used. + + +The task of interpolating image data is simplified by the following +considerations: +.ls 4 +.ls o +The pixels (data points) are equally spaced along a line or on a rectangular +grid. +.le +.ls o +There is no need for coordinate transformations. The coordinates of a pixel +or point to be interpolated are the same as the subscript of the pixel in the +data array. The coordinate of the first pixel in the array may be assumed +to be 1.0, and the spacing between pixels is also 1.0. +.le +.le + + +The task of interpolating image data is complicated by the following +considerations: +.ls +.ls o +We would like the "quality" of the interpolator to be something that can +be set at run time by the user. This is done by permitting the user to +select the type of interpolator to be used, and the order of interpolation. +The following interpolators will be provided: + +.nf + - nearest neighbor + - linear + - natural cubic spline + - divided differences, either 3rd or 5th order. +.fi +.le +.ls o +The interpolator must be able to deal with indefinite valued pixels. +In the IRAF, an indefinite valued pixel has the special value INDEF. +If the interpolator cannot return a meaningful interpolated value +(due to the presence of indefinite pixels in the data, or because the +coordinates of the point to be interpolated are out of bounds) the value +INDEF should be returned instead. +.le +.le + +.sh +Sequential Interpolation of One Dimensional Arrays + + Ideally, we should be able to define a single set of procedures +which can perform interpolation by any of the techniques mentioned above. +The procedures for interpolating one dimensional arrays are optimized +for the case where the entire array is to be sampled, probably in a sequential +fashion. + +The following calls should suffice. The package prefix is "asi", meaning +array sequential interpolate, and the last three letters of each procedure +name are used to describe the action performed by the procedure. The +interpolator is restricted to single precision real data, so there is no +need for a type suffix. + +.ks +.nf + asiset (coeff, interpolator_type, boundary_type) + asifit (data, npts, coeff) + y = asival (x, coeff) # evaluate y(x) + asider (x, derivs, nderiv, coeff) # derivatives + v = asigrl (x1, x2, coeff) # take integral +.fi +.ke + +Here, COEFF is a structure used to describe the interpolant, and +the interpolator type is chosen from the following set of options: + +.ks +.nf + NEAREST return value of nearest neighbor + LINEAR linear interpolation + DIVDIF3 third order divided differences + DIVDIF5 fifth order divided differences + SPLINE3 cubic spline +.fi +.ke + +When the VAL procedure is called to evaluate the interpolant at a point X, +the code must check to see if X is out of bounds. If so, the type of value +returned may be specified by the third and final argument to the SET procedure. + +.ks +.nf + B_NEAREST return nearest boundary (nb) value + B_INDEF (default; return indefinite) + B_REFLECT interpolate at abs(x) + B_WRAP wrap around to opposite boundary + B_PROJECT return y(nb) - (y(abs(x)) - y(nb)) +.fi +.ke + +For example, the following subset preprocessor code fragment uses +the array interpolation routines to shift a data array by a constant +amount: + + +.ks +.nf + real shift, coeff[NPIX+SZ_ASIHDR] + real inrow[NPIX], outrow[NPIX], asival() + int i, interpolator_type, boundary_type, clgeti() + + begin + # get parameters from CL, setup interpolator + interpolator_type = clgeti ("interpolator") + boundary_type = clgeti ("boundary_type") + call asiset (coeff, interpolator_type, boundary_type) + + [code to read pixels into inrow] + call asifit (inrow, NPIX, coeff) + + do i = 1, NPIX # do the interpolation + outrow[i] = asival (i + shift, coeff) +.fi +.ke + + +The array COEFF is actually a simple structure. An initial call to ASISET +is required to save the interpolator parameters in the COEFF structure. +ASIFIT computes the coefficients of the interpolator curve, leaving the +coefficients in the appropriate part of the COEFF structure. Thereafter, +COEFF is read only. ASIVAL evaluates the zeroth derivative of the curve +at the given X coordinate. ASIDER evaluates the first NDERIV derivatives +at X, writing the values of these derivatives into the DERIVS output array. +ASIGRL integrates the interpolant over the indicated range (returning INDEF +if there are any indefinite pixels in that range). + +.ks +.nf + element usage + + 1 type of interpolator + 2 type of boundary condition + 3 number of pixels in x + 4 number of pixels in y + 5-9 reserved for future use + 10-N coefficients + + The COEFF structure +.fi +.ke + + +In the case of the nearest neighbor and linear interpolation algorithms, +the coefficient array will contain the same data as the input data array. +The cubic spline algorithm requires space for NPIX+2 b-spline coefficients. +The space requirements of the divided differences algorithm are similar. +The constant SZ_ASIHDR should allow enough space for the worst case. + +Note that both 3rd and 5th order divided differences interpolators are +provided, but only the 3rd order (cubic) spline. This is because the +5th order spline is considerably harder than the cubic spline. We can easily +add a 5th order spline, or indeed a completely new algorithm such as the sinc +function interpolator, without changing the calling sequences or data +structures. + +.sh +Random Interpolation of One Dimensional Arrays + + If we wish only to interpolate a small region of an array, or a few +points scattered along the array at random, or if memory space is more +precious than execution speed, then a different sort of interpolator is +desirable. We can dispense with the coefficient array in this case, and +operate directly on the data array. + +.ks +.nf + ariset (interpolator_type, boundary_type) + arifit (data, npts) + y = arival (x, data) + arider (x, derivs, nderiv, data) +.fi +.ke + +Since the number of points required to interpolate at a given X value +is fixed, these routines are internally quite different than the sequential +procedures. + +.sh +Sequential Interpolation of Two Dimensional Arrays + + None of the interpolation algorithms mentioned above are difficult +to generalize to two dimensions. The calling sequences should be kept as +similar as possible. The sequential two dimensional procedures are limited +by the size of the two dimensional array which can be kept in memory, and +therefore are most useful for subrasters. The COEFF array should be +declared as a one dimensional array, even though it used to store a two +dimensional array of coefficients. + +.ks +.nf + msiset (coeff, interpolator_type), boundary_type) + msifit (data, nx, ny, coeff) + y = msival (x, y, coeff) + msider (x, y, derivs, nderiv, coeff) +.fi +.ke + +Note that the coefficients of the bicubic spline are obtained by first +doing a one dimensional spline fit to the rows of the data array, leaving +the resultant coefficients in the rows of the coeff array, followed by +a series of one dimensional fits to the coefficients in the coeff array. +The row coefficients are overwritten by the coefficients of the tensor +product spline. Since all the work is done by the one dimensional spline +fitting routine, little extra code is required to handle the two dimensional +case. + +The bicubic spline is evaluated by summing the product C(i,j)*Bi(x)*Bj(y) +at each of the sixteen coefficients contributing to the point (x,y). +Once again, the one dimensional routines for evaluating the b-spline +(the functions Bi(x) and Bj(y)) do most of the work. + +Although it is not required by very many applications, it would also +be useful to be able to compute the surface integral of the interpolant +over an arbitrary surface (this is useful for surface photometry +applications). Never having done this, I do not have any recommendations +on how to set up the calling sequence. + +.sh +Random Interpolation of Two Dimensional Arrays + + The random interpolation procedures are particularly useful for two +dimensional arrays due to the increased size of the arrays. Also, if an +application generally finds linear interpolation to be sufficient, and rarely +uses a more expensive interpolator, the random routines will be more +efficient overall. + +.ks +.nf + mriset (interpolator_type, boundary_type) + msifit (data, nx, ny) + y = mrival (x, y, data) + mrider (x, y, derivs, nderiv, data) +.fi +.ke + +.sh +Cardinal B-Splines + + There are many different kinds of splines. One of the best +for the case where the data points are equally spaced is the cardinal spline. +The so called "natural" end point conditions are probably the best one can +do with noisy data. Other end point conditions, such as not-a-knot, may be +best for approximating analytic functions, but they are unreliable with +noisy data. The natural cubic cardinal spline is fitted by solving the +following tridiagonal system of equations (Prenter, 1975). + + +.ks +.nf + / : \ + | 6 -12 6 : 0 | + | 1 4 1 : y(1) | + | 1 4 1 : y(2) | + | ... : | + | 1 4 1 : y(n-1) | + | 1 4 1 : y(n) | + | 6 -12 6 : 0 | + \ : / +.fi +.ke + + +This matrix can most efficiently be solved by hardwiring the appropriate +recursion relations into a procedure. The problem of what to do when some +of the Y(i) are indefinite has not yet been solved. + +The b-spline basis function used above is defined by + +.ks +.nf + B(x) = (x-x1)**3 x:[x1,x2] + B(x) = 1 + 3(x-x2) + 3(x-x2)**2 - 3(x-x2)**3 x:[x2,x3] + B(x) = 1 + 3(x3-x) + 3(x3-x)**2 - 3(x3-x)**3 x:[x3,x4] + B(x) = (x4-x)**3 x:[x4,x5] +.fi +.ke + +where X1 through X4 are the X coordinates of the four b-spline +coefficients contributing to the b-spline function at X. diff --git a/math/interp/README b/math/interp/README new file mode 100644 index 00000000..915afd6e --- /dev/null +++ b/math/interp/README @@ -0,0 +1,7 @@ +Image interpolation package. Contains interpolator routines specially +designed to interpolate image data. This interpolator task is characterized +by data arrays in which the data points are equally spaced. Data points +and interpolation points are referenced simply by array index. Indefinite +valued pixels may be present in the data. Since the interpolation routines +will be applied to very large data arrays, efficiency is emphasized in these +routines. diff --git a/math/interp/arbpix.x b/math/interp/arbpix.x new file mode 100644 index 00000000..0b61f3d2 --- /dev/null +++ b/math/interp/arbpix.x @@ -0,0 +1,203 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +.help + arbpix -- fix bad data + +Takes care of bad pixels by replacing the indef's with interpolated values. + +In order to replace bad points, the spline interpolator uses a limited +data array, the maximum total length is given by SPLPTS. + +The process is divided as follows: + 1. Take care of points below first good point and above last good + good point. + 2. Load an array with only good points. + 3. Interpolate that array. + +.endhelp + +procedure arbpix(datain,n,dataout,terptype) +include "interpdef.h" +include "asidef.h" + +real datain[ARB] # data in array +int n # no. of data points +real dataout[ARB] # data out array - cannot be same as data in +int terptype # interpolator type - see interpdef.h + +int i, badnc +int k, ka, kb + +real iirbfval() + +begin + + # count bad points + badnc = 0 + do i = 1,n + if (IS_INDEFR (datain[i])) + badnc = badnc + 1 + + # return if all bad or all good + if(badnc == n || badnc == 0) + return + + + # find first good point + for (ka = 1; IS_INDEFR (datain[ka]); ka = ka + 1) + ; + + # bad points below first good point are set at first value + do k = 1,ka-1 + dataout[k] = datain[ka] + + # find last good point + for (kb = n; IS_INDEFR (datain[kb]); kb = kb - 1) + ; + + # bad points beyond last good point get set at last value + do k = n, kb+1, -1 + dataout[k] = datain[kb] + + # load the other points interpolating the bad points as needed + do k = ka, kb { + if (!IS_INDEFR (datain[k])) # good point + dataout[k] = datain[k] + + else # bad point -- generate interpolated value + dataout[k] = iirbfval(datain[ka],kb-ka+1,k-ka+1,terptype) + } + +end + +# This part fills a temporary array with good points that bracket the +# bad point and calls the interpolating routine. + +real procedure iirbfval(y, n, k, terptype) + +real y[ARB] # data_in array, y[1] and y[n] guaranteed to be good. +int n # length of data_in array. +int k # index of bad point to replace. +int terptype + +int j, jj, pd, pu, tk, pns +real td[SPLPTS], tx[SPLPTS] # temporary arrays for interpolation + +real iirbint() + +begin + # The following test is done to improve speed. + # This code will work only if subroutines are implemented by + # using static storage - i.e. the old internal values survive. + # This avoids reloading of temporary arrays if + # there are consequetive bad points. + + if (!IS_INDEFR (y[k-1])) { + # set number of good points needed on each side of bad point + switch (terptype) { + case IT_NEAREST : + pns = 1 + case IT_LINEAR : + pns = 1 + case IT_POLY3 : + pns = 2 + case IT_POLY5 : + pns = 3 + case IT_SPLINE3 : + pns = SPLPTS / 2 + } + + # search down + pd = 0 + for (j = k-1; j >= 1 && pd < pns; j = j-1) + if (!IS_INDEFR (y[j])) + pd = pd + 1 + + # load temp. arrays for values below our indef. + tk = 0 + for(jj = j + 1; jj < k; jj = jj + 1) + if (!IS_INDEFR (y[jj])) { + tk = tk + 1 + td[tk] = y[jj] + tx[tk] = jj + } + + # search and load up from indef. + pu = 0 + for (j = k + 1; j <= n && pu < pns; j = j + 1) + if (!IS_INDEFR (y[j])) { + pu = pu + 1 + tk = tk + 1 + td[tk] = y[j] + tx[tk] = j + } + } + + # return value interpolated from these arrays. + return(iirbint(real(k), tx, td, tk, pd, terptype)) + +end + + +# This part interpolates the temporary arrays. +# It does not represent a general purpose routine because the +# previous part has determined the proper indices etc. so that +# effort is not duplicated here. + +real procedure iirbint (x, tx, td, tk, pd, terptype) + +real x # point to interpolate +real tx[ARB] # xvalues +real td[ARB] # data values +int tk # size of data array +int pd # index such that tx[pd] < x < tx[pd+1] +int terptype + +int i, ks, tpol +real cc[4,SPLPTS] +real h + +real iipol_terp() + +begin + switch (terptype) { + + case IT_NEAREST : + if (x - tx[1] > tx[2] - x) + return(td[2]) + else + return(td[1]) + + case IT_LINEAR : + return(td[1] + (x - tx[1]) * + (td[2] - td[1]) / (tx[2] - tx[1])) + + case IT_SPLINE3 : + do i = 1,tk + cc[1,i] = td[i] + cc[2,1] = 0. + cc[2,tk] = 0. + + # use spline routine from C. de Boor's book + # A Practical Guide to Splines + call cubspl(tx,cc,tk,2,2) + h = x - tx[pd] + return(cc[1,pd] + h * (cc[2,pd] + h * + (cc[3,pd] + h * cc[4,pd]/3.)/2.)) + + default : # one of the polynomial types + # allow lower order if not enough points on one side + tpol = tk + ks = 1 + if (tk - pd < pd) { + tpol = 2 * (tk - pd) + ks = 2 * pd - tk + 1 + } + if (tk - pd > pd) + tpol = 2 * pd + + # finally polynomial interpolate + return(iipol_terp(tx[ks], td[ks], tpol, x)) + } + +end diff --git a/math/interp/arider.x b/math/interp/arider.x new file mode 100644 index 00000000..fef55e0e --- /dev/null +++ b/math/interp/arider.x @@ -0,0 +1,214 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +.help + arider -- random interpolator returns derivatives + +First looks at type of interpolator requested, then calls routine +to evaluate. +No error traps are included -- unreasonable values for the number of +data points or the x position will either produce hard errors or garbage. + +.endhelp + +procedure arider(x, datain, n, derivs, nderiv, interptype) +include "interpdef.h" + +real x # need 1 <= x <= n +real datain[ARB] # data values +int n # number of data values +real derivs[ARB] # derivatives out -- beware derivs[1] is + # function value +int nderiv # total number of values returned in derivs +int interptype + +int i, j, k, nt, nd, nx +real pc[6], ac, s + +begin + if(nderiv <= 0) + return + + # zero out derivs array + do i = 1, nderiv + derivs[i] = 0. + + switch (interptype) { + case IT_NEAREST : + derivs[1] = datain[int(x + 0.5)] + return + case IT_LINEAR : + nx = x + if( nx >= n ) + nx = nx - 1 + derivs[1] = (x - nx) * datain[nx+1] + (nx + 1 - x) * datain[nx] + if(nderiv >= 2) + derivs[2] = datain[nx+1] - datain[nx] + return + # The other cases call subroutines to generate polynomial coeff. + case IT_POLY3 : + call iidr_poly3(x, datain, n, pc) + nt = 4 + case IT_POLY5 : + call iidr_poly5(x, datain, n, pc) + nt = 6 + case IT_SPLINE3 : + call iidr_spline3(x, datain, n, pc) + nt = 4 + } + + nx = x + s = x - nx + nd = nderiv + if (nderiv > nt) + nd = nt + + do k = 1,nd { + ac = pc[nt - k + 1] # evluate using nested multiplication + do j = nt - k, 1, -1 + ac = pc[j] + s * ac + + derivs[k] = ac + + do j = 1, nt - k # differentiate + pc[j] = j * pc[j + 1] + } +end + +procedure iidr_poly3(x, datain, n, pc) + +real x +real datain[ARB] +int n +real pc[ARB] + +int i, k, nx, nt +real a[4] + +begin + + nx = x + + # The major complication is that near the edge interior polynomial + # must somehow be defined. + k = 0 + for(i = nx - 1; i <= nx + 2; i = i + 1){ + k = k + 1 + # project data points into temporary array + if ( i < 1 ) + a[k] = 2. * datain[1] - datain[2 - i] + else if ( i > n ) + a[k] = 2. * datain[n] - datain[2 * n - i] + else + a[k] = datain[i] + } + + nt = 4 + + # generate diffrence table for Newton's form + do k = 1, nt-1 + do i = 1, nt-k + a[i] = (a[i+1] - a[i]) / k + + # shift to generate polynomial coefficients + do k = nt,2,-1 + do i = 2,k + a[i] = a[i] + a[i-1] * (k - i - nt/2) + + do i = 1,nt + pc[i] = a[nt+1-i] + + return + +end + +procedure iidr_poly5(x, datain, n, pc) + +real x +real datain[ARB] +int n +real pc[ARB] + +int i, k, nx, nt +real a[6] + +begin + + nx = x + + # The major complication is that near the edge interior polynomial + # must somehow be defined. + k = 0 + for(i = nx - 2; i <= nx + 3; i = i + 1){ + k = k + 1 + # project data points into temporary array + if ( i < 1 ) + a[k] = 2. * datain[1] - datain[2 - i] + else if ( i > n ) + a[k] = 2. * datain[n] - datain[2 * n - i] + else + a[k] = datain[i] + } + + nt = 6 + + # generate diffrence table for Newton's form + do k = 1, nt-1 + do i = 1, nt-k + a[i] = (a[i+1] - a[i]) / k + + # shift to generate polynomial coefficients + do k = nt,2,-1 + do i = 2,k + a[i] = a[i] + a[i-1] * (k - i - nt/2) + + do i = 1,nt + pc[i] = a[nt+1-i] + + return + +end + +procedure iidr_spline3(x, datain, n, pc) + +real x +real datain[ARB] +int n +real pc[ARB] + +int i, k, nx, px +real temp[SPLPTS+2], bcoeff[SPLPTS+2], h + +begin + + nx = x + + h = x - nx + k = 0 + # maximum number of points used is SPLPTS + for(i = nx - SPLPTS/2 + 1; i <= nx + SPLPTS/2; i = i + 1){ + if(i < 1 || i > n) + ; + else { + k = k + 1 + if(k == 1) + px = nx - i + 1 + bcoeff[k+1] = datain[i] + } + } + + bcoeff[1] = 0. + bcoeff[k+2] = 0. + + # Use special routine for cardinal splines. + call iif_spline(bcoeff, temp, k) + + px = px + 1 + + pc[1] = bcoeff[px-1] + 4. * bcoeff[px] + bcoeff[px+1] + pc[2] = 3. * (bcoeff[px+1] - bcoeff[px-1]) + pc[3] = 3. * (bcoeff[px-1] - 2. * bcoeff[px] + bcoeff[px+1]) + pc[4] = -bcoeff[px-1] + 3. * bcoeff[px] - 3. * bcoeff[px+1] + + bcoeff[px+2] + + return +end diff --git a/math/interp/arival.x b/math/interp/arival.x new file mode 100644 index 00000000..50ee1c3e --- /dev/null +++ b/math/interp/arival.x @@ -0,0 +1,124 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +.help + arival -- random interpolator returns value + +No error traps are included -- unreasonable values for the number of +data points or the x position will either produce hard errors or garbage. +This version is written in-line for speed. + +.endhelp + +real procedure arival(x, datain, n, interptype) +include "interpdef.h" + +real x # need 1 <= x <= n +real datain[ARB] # data values +int n # number of data values +int interptype + +int i, k, nx, px +real a[6], cd20, cd21, cd40, cd41, s, t, h +real bcoeff[SPLPTS+2], temp[SPLPTS+2], pc[4] + +begin + switch (interptype) { + case IT_NEAREST : + return(datain[int(x + 0.5)]) + case IT_LINEAR : + nx = x + # protect against x = n case + # else it will reference datain[n + 1] + if(nx >= n) + nx = nx - 1 + return((x - nx) * datain[nx+1] + (nx + 1 - x) * datain[nx]) + case IT_POLY3 : + nx = x + + # The major complication is that near the edge interior polynomial + # must somehow be defined. + k = 0 + for(i = nx - 1; i <= nx + 2; i = i + 1){ + k = k + 1 + # project data points into temporary array + if ( i < 1 ) + a[k] = 2. * datain[1] - datain[2 - i] + else if ( i > n ) + a[k] = 2. * datain[n] - datain[2 * n - i] + else + a[k] = datain[i] + } + + s = x - nx + t = 1. - s + + # second central differences + cd20 = 1./6. * (a[3] - 2. * a[2] + a[1]) + cd21 = 1./6. * (a[4] - 2. * a[3] + a[2]) + + return( s * (a[3] + (s * s - 1.) * cd21) + + t * (a[2] + (t * t - 1.) * cd20) ) + case IT_POLY5 : + nx = x + + # The major complication is that near the edge interior polynomial + # must somehow be defined. + k = 0 + for(i = nx - 2; i <= nx + 3; i = i + 1){ + k = k + 1 + # project data points into temporary array + if ( i < 1 ) + a[k] = 2. * datain[1] - datain[2 - i] + else if ( i > n ) + a[k] = 2. * datain[n] - datain[2 * n - i] + else + a[k] = datain[i] + } + + s = x - nx + t = 1. - s + + # second central differences + cd20 = 1./6. * (a[4] - 2. * a[3] + a[2]) + cd21 = 1./6. * (a[5] - 2. * a[4] + a[3]) + + # fourth central differences + cd40 = 1./120. * (a[1] - 4. * a[2] + 6. * a[3] - 4. * a[4] + a[5]) + cd41 = 1./120. * (a[2] - 4. * a[3] + 6. * a[4] - 4. * a[5] + a[6]) + + return( s * (a[4] + (s * s - 1.) * (cd21 + (s * s - 4.) * cd41)) + + t * (a[3] + (t * t - 1.) * (cd20 + (t * t - 4.) * cd40)) ) + case IT_SPLINE3 : + nx = x + + h = x - nx + k = 0 + # maximum number of points used is SPLPTS + for(i = nx - SPLPTS/2 + 1; i <= nx + SPLPTS/2; i = i + 1){ + if(i < 1 || i > n) + ; + else { + k = k + 1 + if(k == 1) + px = nx - i + 1 + bcoeff[k+1] = datain[i] + } + } + + bcoeff[1] = 0. + bcoeff[k+2] = 0. + + # Use special routine for cardinal splines. + call iif_spline(bcoeff, temp, k) + + px = px + 1 + + pc[1] = bcoeff[px-1] + 4. * bcoeff[px] + bcoeff[px+1] + pc[2] = 3. * (bcoeff[px+1] - bcoeff[px-1]) + pc[3] = 3. * (bcoeff[px-1] - 2. * bcoeff[px] + bcoeff[px+1]) + pc[4] = -bcoeff[px-1] + 3. * bcoeff[px] - 3. * bcoeff[px+1] + + bcoeff[px+2] + + return(pc[1] + h * (pc[2] + h * (pc[3] + h * pc[4]))) + } +end diff --git a/math/interp/asidef.h b/math/interp/asidef.h new file mode 100644 index 00000000..9fd06084 --- /dev/null +++ b/math/interp/asidef.h @@ -0,0 +1,16 @@ +# defines for use with interpolator package +# intended for internal consumption only +# used to set up storage in coeff array + +define TYPEI coeff[1] # code for interpolator type +define NPTS coeff[2] # no. of data points + +define ITYPEI int(coeff[1]) +define INPTS int(coeff[2]) + +define ASITYPEERR 1 +define ASIFITERR 2 + +define COFF 10 # offset into coeff array + # beware coeff[COFF - 2] is first + # location used !! diff --git a/math/interp/asider.x b/math/interp/asider.x new file mode 100644 index 00000000..2c9a2819 --- /dev/null +++ b/math/interp/asider.x @@ -0,0 +1,121 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +.help + procedure asider + + This procedure finds derivatives assuming that x lands in +array i.e. 1 <= x <= npts. +It must be called by a routine that checks for out of bound references +and takes care of bad points and checks to make sure that the number +of derivatives requested is reasonable (max. is 5 derivatives or nder +can be 6 at most and min. nder is 1) + + This version assumes interior polynomial interpolants are stored as +the data points with corections for bad points, and the spline interpolant +is stored as a series of b-spline coefficients. + + The storage and evaluation for nearest neighbor and linear interpolation +are simpler so they are evaluated separately from other piecewise +polynomials. +.endhelp + +procedure asider(x,der,nder,coeff) +include "interpdef.h" +include "asidef.h" + +real x +real coeff[ARB] +int nder # no. items returned = 1 + no. of deriv. + +real der[ARB] # der[1] is value der[2] is f prime etc. + +int nx,n0,i,j,k,nt,nd +real s,ac,pc[6],d[6] + +begin + do i = 1, nder # return zero for derivatives that are zero + der[i] = 0. + + nt = 0 # nt is number of terms in case polynomial type + + switch (ITYPEI) { # switch on interpolator type + + case IT_NEAREST : + der[1] = (coeff[COFF + nint(x)]) + return + + case IT_LINEAR : + nx = x + der[1] = (x - nx) * coeff[COFF + nx + 1] + + (nx + 1 - x) * coeff[COFF + nx] + if ( nder > 1 ) # try to return exact number requested + der[2] = coeff[COFF + nx + 1] - coeff[COFF + nx] + return + + case IT_POLY3 : + nt = 4 + + case IT_POLY5 : + nt = 6 + + case IT_SPLINE3 : + nt = 4 + + } + + # falls through to here if interpolant is one of + # the higher order polynomial types or third order spline + + nx = x + n0 = COFF + nx + s = x - nx + + nd = nder # no. of derivatives needed + if ( nder > nt ) + nd = nt + + # generate polynomial coefficients, first for spline. + if (ITYPEI == IT_SPLINE3) { + pc[1] = coeff[n0-1] + 4. * coeff[n0] + coeff[n0+1] + pc[2] = 3. * (coeff[n0+1] - coeff[n0-1]) + pc[3] = 3. * (coeff[n0-1] - 2. * coeff[n0] + coeff[n0+1]) + pc[4] = -coeff[n0-1] + 3. * coeff[n0] - 3. * coeff[n0+1] + + coeff[n0+2] + + } else { + + # Newton's form written in line to get polynomial from data + # load data + do i = 1,nt + d[i] = coeff[n0 - nt/2 + i] + + # generate difference table + do k = 1, nt-1 + do i = 1,nt-k + d[i] = (d[i+1] - d[i]) / k + + # shift to generate polynomial coefficients of (x - n0) + do k = nt,2,-1 + do i = 2,k + d[i] = d[i] + d[i-1] * (k - i - nt/2) + + do i = 1,nt + pc[i] = d[nt + 1 - i] + } + + do k = 1,nd { # as loop progresses pc contains coefficients of + # higher and higher derivatives + + ac = pc[nt - k + 1] + do j = nt - k, 1, -1 # evaluate using nested mult. + ac = pc[j] + s * ac + + der[k] = ac + + do j = 1,nt - k # differentiate polynomial + pc[j] = j * pc[j + 1] + } + + return + +end diff --git a/math/interp/asieva.x b/math/interp/asieva.x new file mode 100644 index 00000000..bc5fcbb3 --- /dev/null +++ b/math/interp/asieva.x @@ -0,0 +1,38 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +# procedure to evaluate interpolator at a sequence of ordered points +# aasumes that all x_values land within [1,n]. + +procedure asieva(x, y, n, coeff) +include "interpdef.h" +include "asidef.h" + +real x[ARB] # ordered x array +int n # no. of points in x +real coeff[ARB] + +real y[ARB] # interpolated values + +begin + switch (ITYPEI) { # switch on interpolator type + + case IT_NEAREST : + call iievne(x, y, n, coeff[COFF + 1]) + + case IT_LINEAR : + call iievli(x, y, n, coeff[COFF + 1]) + + case IT_POLY3 : + call iievp3(x, y, n, coeff[COFF + 1]) + + case IT_POLY5 : + call iievp5(x, y, n, coeff[COFF + 1]) + + case IT_SPLINE3 : + call iievs3(x, y, n, coeff[COFF + 1]) + + } + + return + +end diff --git a/math/interp/asifit.x b/math/interp/asifit.x new file mode 100644 index 00000000..12b16b0f --- /dev/null +++ b/math/interp/asifit.x @@ -0,0 +1,75 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +.help + asifit -- fit interpolator to data + +This version uses the "basis-spline" representation for the spline. +It stores the data array itself for the polynomial interpolants. + +.endhelp + +procedure asifit(datain,n,coeff) +include "interpdef.h" +include "asidef.h" + +real datain[ARB] # data array +int n # no. of data points +real coeff[ARB] + +int i + +begin + NPTS = n + + # error trap - check data array for size + switch (ITYPEI) { + case IT_SPLINE3: + if(n < 4) + call error(ASIFITERR,"too few points for spline") + case IT_POLY5: + if(n < 6) + call error(ASIFITERR,"too few points for poly5") + case IT_POLY3: + if(n < 4) + call error(ASIFITERR,"too few points for poly3") + case IT_LINEAR: + if(n < 2) + call error(ASIFITERR,"too few points for linear") + default: + if(n < 1) + call error(ASIFITERR,"too few points for interpolation") + + } + + do i = 1,n + coeff[COFF + i] = datain[i] + + + if (ITYPEI == IT_SPLINE3) { + + # specify natural end conditions - second deriv. zero + coeff[COFF] = 0. + coeff[COFF+n+1] = 0. + + # fit spline - generate b-spline coefficients. + call iif_spline(coeff[COFF],coeff[COFF + n + 2],n) + + + } else { # not the spline + # We extend array to take care of values near edge. + + # Assign arbitrary values in case there are only a few points. + for(i = n + 1; i < 5; i = i + 1) + coeff[COFF + i] = coeff[COFF + n] + + # Extend for worst case - poly 5 may need 3 extra points. + do i = 1,3 { # same recipe as project + coeff[COFF+1 - i] = 2. * coeff[COFF + 1] - + coeff[COFF + 1 + i] + coeff[COFF+n + i] = 2. * coeff[COFF + n] - + coeff[COFF + n - i] + } + } + + return +end diff --git a/math/interp/asigrl.x b/math/interp/asigrl.x new file mode 100644 index 00000000..d528d0be --- /dev/null +++ b/math/interp/asigrl.x @@ -0,0 +1,201 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +.help + procedure asigrl + + This procedure finds the integral of the interpolant from a to b +assuming both a and b land in the array. +.endhelp + +real procedure asigrl(a,b,coeff) # returns value of integral +include "interpdef.h" +include "asidef.h" + +real a,b # integral limits +real coeff[ARB] + +int na,nb,i,j,n0,nt +real s,t,ac,xa,xb,pc[6] + +begin + xa = a + xb = b + if ( a > b ) { # flip order and sign at end + xa = b + xb = a + } + + na = xa + nb = xb + ac = 0. # zero accumulator + + # set number of terms + switch (ITYPEI) { # switch on interpolator type + case IT_NEAREST : + nt = 0 + case IT_LINEAR : + nt = 1 + case IT_POLY3 : + nt = 4 + case IT_POLY5 : + nt = 6 + case IT_SPLINE3 : + nt = 4 + } + + # NEAREST_NEIGHBOR and LINEAR are handled differently because of + # storage. Also probably good for speed. + + if (nt == 0) { # NEAREST_NEIGHBOR + # reset segment to center values + na = xa + 0.5 + nb = xb + 0.5 + + # set up for first segment + s = xa - na + + # for clarity one segment case is handled separately + if ( nb == na ) { # only one segment involved + t = xb - nb + n0 = COFF + na + ac = ac + (t - s) * coeff[n0] + } else { # more than one segment + + # first segment + n0 = COFF + na + ac = ac + (0.5 - s) * coeff[n0] + + # middle segments + do j = na+1, nb-1 { + n0 = COFF + j + ac = ac + coeff[n0] + } + + # last segment + n0 = COFF + nb + t = xb - nb + ac = ac + (t + 0.5) * coeff[n0] + } + + } else if (nt == 1) { # LINEAR + # set up for first segment + s = xa - na + + # for clarity one segment case is handled separately + if ( nb == na ) { # only one segment involved + t = xb - nb + n0 = COFF + na + ac = ac + (t - s) * coeff[n0] + + 0.5 * (coeff[n0+1] - coeff[n0]) * (t*t - s*s) + } else { # more than one segment + + # first segment + n0 = COFF + na + ac = ac + (1. - s) * coeff[n0] + + 0.5 * (coeff[n0+1] - coeff[n0]) * (1. - s*s) + + # middle segments + do j = na+1, nb-1 { + n0 = COFF + j + ac = ac + 0.5 * (coeff[n0+1] + coeff[n0]) + } + + # last segment + n0 = COFF + nb + t = xb - nb + ac = ac + coeff[n0] * t + 0.5 * + (coeff[n0+1] - coeff[n0]) * t * t + } + + } else { # A higher order interpolant + + # set up for first segment + s = xa - na + + # for clarity one segment case is handled separately + if ( nb == na ) { # only one segment involved + t = xb - nb + n0 = COFF + na + call iigetpc(n0,pc,coeff) + do i = 1,nt + ac = ac + (1./i) * pc[i] * (t ** i - s ** i) + } else { # more than one segment + + # first segment + n0 = COFF + na + call iigetpc(n0,pc,coeff) + do i = 1,nt + ac = ac + (1./i) * pc[i] * (1. - s ** i) + + # middle segments + do j = na+1, nb-1 { + n0 = COFF + j + call iigetpc(n0,pc,coeff) + do i = 1,nt + ac = ac + (1./i) * pc[i] + } + + # last segment + n0 = COFF + nb + t = xb - nb + call iigetpc(n0,pc,coeff) + do i = 1,nt + ac = ac + (1./i) * pc[i] * t ** i + } + } + + if ( a < b ) + return(ac) + else + return(-ac) +end + + +procedure iigetpc(n0, pc, coeff) # generates polynomial coefficients + # if spline or poly3 or poly5 + +int n0 # coefficients wanted for n0 < x n0 + 1 +real coeff[ARB] + +real pc[ARB] + +int i,k,nt +real d[6] + +begin + # generate polynomial coefficients, first for spline. + if (ITYPEI == IT_SPLINE3) { + pc[1] = coeff[n0-1] + 4. * coeff[n0] + coeff[n0+1] + pc[2] = 3. * (coeff[n0+1] - coeff[n0-1]) + pc[3] = 3. * (coeff[n0-1] - 2. * coeff[n0] + coeff[n0+1]) + pc[4] = -coeff[n0-1] + 3. * coeff[n0] - 3. * coeff[n0+1] + + coeff[n0+2] + + } else { + if (ITYPEI == IT_POLY5) + nt = 6 + else # must be POLY3 + nt = 4 + + # Newton's form written in line to get polynomial from data + + # load data + do i = 1,nt + d[i] = coeff[n0 - nt/2 + i] + + # generate difference table + do k = 1, nt-1 + do i = 1,nt-k + d[i] = (d[i+1] - d[i]) / k + + # shift to generate polynomial coefficients of (x - n0) + do k = nt,2,-1 + do i = 2,k + d[i] = d[i] + d[i-1] * (k - i - 2) + + do i = 1,nt + pc[i] = d[nt + 1 - i] + } + + return +end diff --git a/math/interp/asiset.x b/math/interp/asiset.x new file mode 100644 index 00000000..126f4e39 --- /dev/null +++ b/math/interp/asiset.x @@ -0,0 +1,20 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +# sets sequential interpolator type + +procedure asiset(coeff,interptype) +include "interpdef.h" +include "asidef.h" + +int interptype +real coeff[ARB] + +begin + + if(interptype <= 0 || interptype > ITNIT) + call error(ASITYPEERR,"illegal interpolator type") + else + TYPEI = interptype + + return +end diff --git a/math/interp/asival.x b/math/interp/asival.x new file mode 100644 index 00000000..56f18938 --- /dev/null +++ b/math/interp/asival.x @@ -0,0 +1,49 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +.help + procedure asival + + This procedure finds interpolated value assuming that x lands in +array i.e. 1 <= x < npts. +It must be called by a routine that checks for out of bound references +and takes care of bad pixels. + This version assumes that the data are stored for polynomials and +the spline appears as basis-spline coefficients. + Also the sequential evaluators are used to obtain the values in +order to reduce the amount of duplicated code. +.endhelp + +real procedure asival(x,coeff) +include "interpdef.h" +include "asidef.h" + +real x +real coeff[ARB] + +real t + +begin + switch (ITYPEI) { # switch on interpolator type + + case IT_NEAREST : + call iievne(x,t,1,coeff[COFF+1]) + return(t) + + case IT_LINEAR : + call iievli(x,t,1,coeff[COFF+1]) + return(t) + + case IT_POLY3 : + call iievp3(x,t,1,coeff[COFF+1]) + return(t) + + case IT_POLY5 : + call iievp5(x,t,1,coeff[COFF+1]) + return(t) + + case IT_SPLINE3 : + call iievs3(x,t,1,coeff[COFF+1]) + return(t) + + } +end diff --git a/math/interp/bench.x b/math/interp/bench.x new file mode 100644 index 00000000..539ddf67 --- /dev/null +++ b/math/interp/bench.x @@ -0,0 +1,55 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + + +include "/usr/vesa/spec/terpdef.h" + +# Quickie benchmark, 1-d sequential image interpolators. + +define MAX_PIXELS 2048 +define SZ_NAME 20 + + +task interp + + +procedure interp() + +real data[MAX_PIXELS], x[MAX_PIXELS], y[MAX_PIXELS] +real coeff[2 * MAX_PIXELS + SZ_ASI] +char interp_name[SZ_NAME] +int i, j, npix, nlines, interpolant +bool streq() +int clgeti() + +begin + npix = min (MAX_PIXELS, clgeti ("npix")) + nlines = npix # square image + call clgstr ("interpolant", interp_name, SZ_NAME) + + if (streq (interp_name, "nearest")) + interpolant = IT_NEAR + else if (streq (interp_name, "linear")) + interpolant = IT_LINEAR + else if (streq (interp_name, "poly3")) + interpolant = IT_POLY3 + else if (streq (interp_name, "poly5")) + interpolant = IT_POLY5 + else if (streq (interp_name, "spline3")) + interpolant = IT_SPLINE3 + else + call error (0, "Unknown interpolant keyword") + + call asiset (coeff, II_INIT, 0) + call asiset (coeff, II_INTERPOLANT, interpolant) + call asiset (coeff, II_ALLGOODPIX, YES) + + do i = 1, npix { # initialize data, x arrays + data[i] = max(1.0, min(real(npix), real(i) + 0.15151)) + x[i] = i + } + + do j = 1, nlines { # interpolate npix ** 2 + call asifit (data, npix, coeff) + call asieva (x, y, npix, coeff) + } +end diff --git a/math/interp/cubspl.f b/math/interp/cubspl.f new file mode 100644 index 00000000..4bb54964 --- /dev/null +++ b/math/interp/cubspl.f @@ -0,0 +1,119 @@ + subroutine cubspl ( tau, c, n, ibcbeg, ibcend ) +c from * a practical guide to splines * by c. de boor +c ************************ input *************************** +c n = number of data points. assumed to be .ge. 2. +c (tau(i), c(1,i), i=1,...,n) = abscissae and ordinates of the +c data points. tau is assumed to be strictly increasing. +c ibcbeg, ibcend = boundary condition indicators, and +c c(2,1), c(2,n) = boundary condition information. specifically, +c ibcbeg = 0 means no boundary condition at tau(1) is given. +c in this case, the not-a-knot condition is used, i.e. the +c jump in the third derivative across tau(2) is forced to +c zero, thus the first and the second cubic polynomial pieces +c are made to coincide.) +c ibcbeg = 1 means that the slope at tau(1) is made to equal +c c(2,1), supplied by input. +c ibcbeg = 2 means that the second derivative at tau(1) is +c made to equal c(2,1), supplied by input. +c ibcend = 0, 1, or 2 has analogous meaning concerning the +c boundary condition at tau(n), with the additional infor- +c mation taken from c(2,n). +c *********************** output ************************** +c c(j,i), j=1,...,4; i=1,...,l (= n-1) = the polynomial coefficients +c of the cubic interpolating spline with interior knots (or +c joints) tau(2), ..., tau(n-1). precisely, in the interval +c interval (tau(i), tau(i+1)), the spline f is given by +c f(x) = c(1,i)+h*(c(2,i)+h*(c(3,i)+h*c(4,i)/3.)/2.) +c where h = x - tau(i). the function program *ppvalu* may be +c used to evaluate f or its derivatives from tau,c, l = n-1, +c and k=4. + integer ibcbeg,ibcend,n, i,j,l,m + real c(4,n),tau(n), divdf1,divdf3,dtau,g +c****** a tridiagonal linear system for the unknown slopes s(i) of +c f at tau(i), i=1,...,n, is generated and then solved by gauss elim- +c ination, with s(i) ending up in c(2,i), all i. +c c(3,.) and c(4,.) are used initially for temporary storage. + l = n-1 +compute first differences of tau sequence and store in c(3,.). also, +compute first divided difference of data and store in c(4,.). + do 10 m=2,n + c(3,m) = tau(m) - tau(m-1) + 10 c(4,m) = (c(1,m) - c(1,m-1))/c(3,m) +construct first equation from the boundary condition, of the form +c c(4,1)*s(1) + c(3,1)*s(2) = c(2,1) + if (ibcbeg-1) 11,15,16 + 11 if (n .gt. 2) go to 12 +c no condition at left end and n = 2. + c(4,1) = 1. + c(3,1) = 1. + c(2,1) = 2.*c(4,2) + go to 25 +c not-a-knot condition at left end and n .gt. 2. + 12 c(4,1) = c(3,3) + c(3,1) = c(3,2) + c(3,3) + c(2,1) =((c(3,2)+2.*c(3,1))*c(4,2)*c(3,3)+c(3,2)**2*c(4,3))/c(3,1) + go to 19 +c slope prescribed at left end. + 15 c(4,1) = 1. + c(3,1) = 0. + go to 18 +c second derivative prescribed at left end. + 16 c(4,1) = 2. + c(3,1) = 1. + c(2,1) = 3.*c(4,2) - c(3,2)/2.*c(2,1) + 18 if(n .eq. 2) go to 25 +c if there are interior knots, generate the corresp. equations and car- +c ry out the forward pass of gauss elimination, after which the m-th +c equation reads c(4,m)*s(m) + c(3,m)*s(m+1) = c(2,m). + 19 do 20 m=2,l + g = -c(3,m+1)/c(4,m-1) + c(2,m) = g*c(2,m-1) + 3.*(c(3,m)*c(4,m+1)+c(3,m+1)*c(4,m)) + 20 c(4,m) = g*c(3,m-1) + 2.*(c(3,m) + c(3,m+1)) +construct last equation from the second boundary condition, of the form +c (-g*c(4,n-1))*s(n-1) + c(4,n)*s(n) = c(2,n) +c if slope is prescribed at right end, one can go directly to back- +c substitution, since c array happens to be set up just right for it +c at this point. + if (ibcend-1) 21,30,24 + 21 if (n .eq. 3 .and. ibcbeg .eq. 0) go to 22 +c not-a-knot and n .ge. 3, and either n.gt.3 or also not-a-knot at +c left end point. + g = c(3,n-1) + c(3,n) + c(2,n) = ((c(3,n)+2.*g)*c(4,n)*c(3,n-1) + * + c(3,n)**2*(c(1,n-1)-c(1,n-2))/c(3,n-1))/g + g = -g/c(4,n-1) + c(4,n) = c(3,n-1) + go to 29 +c either (n=3 and not-a-knot also at left) or (n=2 and not not-a- +c knot at left end point). + 22 c(2,n) = 2.*c(4,n) + c(4,n) = 1. + go to 28 +c second derivative prescribed at right endpoint. + 24 c(2,n) = 3.*c(4,n) + c(3,n)/2.*c(2,n) + c(4,n) = 2. + go to 28 + 25 if (ibcend-1) 26,30,24 + 26 if (ibcbeg .gt. 0) go to 22 +c not-a-knot at right endpoint and at left endpoint and n = 2. + c(2,n) = c(4,n) + go to 30 + 28 g = -1./c(4,n-1) +complete forward pass of gauss elimination. + 29 c(4,n) = g*c(3,n-1) + c(4,n) + c(2,n) = (g*c(2,n-1) + c(2,n))/c(4,n) +carry out back substitution + 30 j = l + 40 c(2,j) = (c(2,j) - c(3,j)*c(2,j+1))/c(4,j) + j = j - 1 + if (j .gt. 0) go to 40 +c****** generate cubic coefficients in each interval, i.e., the deriv.s +c at its left endpoint, from value and slope at its endpoints. + do 50 i=2,n + dtau = c(3,i) + divdf1 = (c(1,i) - c(1,i-1))/dtau + divdf3 = c(2,i-1) + c(2,i) - 2.*divdf1 + c(3,i-1) = 2.*(divdf1 - c(2,i-1) - divdf3)/dtau + 50 c(4,i-1) = (divdf3/dtau)*(6./dtau) + return + end diff --git a/math/interp/iieval.x b/math/interp/iieval.x new file mode 100644 index 00000000..87538525 --- /dev/null +++ b/math/interp/iieval.x @@ -0,0 +1,137 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +.help + procedures to evaluate interpolants + + These procedures are seperated from the other programs in order +to make it easier to optimise these in case it becomes necessary to +obtain the fastest possible interpolants. The interpolation package +spends most of its time in these routines. +.endhelp + +procedure iievne(x,y,n,a) # nearest neighbor + +real x[ARB] # x values, must be within [1,n] +real y[ARB] # interpolated values returned to user +int n # number of x values +real a[ARB] # data to be interpolated + +int i + +begin + do i = 1, n + y[i] = a[int(x[i] + 0.5)] + return +end + +procedure iievli(x,y,n,a) # linear + +real x[ARB] # x values, must be within [1,n] +real y[ARB] # interpolated values returned to user +int n # number of x values +real a[ARB] # data to be interpolated + +int i,nx + +begin + + do i = 1,n { + nx = x[i] + y[i] = (x[i] - nx) * a[nx + 1] + (nx + 1 - x[i]) * a[nx] + } + return +end + +procedure iievp3(x,y,n,a) # interior third order polynomial + +real x[ARB] # x values, must be within [1,n] +real y[ARB] # interpolated values returned to user +int n # number of x values +real a[ARB] # data to be interpolated from a[0] to a[n+2] + +int i,nx,nxold +real s,t,cd20,cd21 + +begin + nxold = -1 + do i = 1,n { + nx = x[i] + s = x[i] - nx + t = 1. - s + if (nx != nxold) { + + # second central differences: + cd20 = 1./6. * (a[nx+1] - 2. * a[nx] + a[nx-1]) + cd21 = 1./6. * (a[nx+2] - 2. * a[nx+1] + a[nx]) + nxold = nx + } + y[i] = s * (a[nx+1] + (s * s - 1.) * cd21) + + t * (a[nx] + (t * t - 1.) * cd20) + + } + return +end + +procedure iievp5(x,y,n,a) # interior fifth order polynomial + +real x[ARB] # x values, must be within [1,n] +real y[ARB] # interpolated values returned to user +int n # number of x values +real a[ARB] # data to be interpolated - from a[-1] to a[n+3] + +int i,nx,nxold +real s,t,cd20,cd21,cd40,cd41 + +begin + nxold = -1 + do i = 1,n { + nx = x[i] + s = x[i] - nx + t = 1. - s + if (nx != nxold) { + cd20 = 1./6. * (a[nx+1] - 2. * a[nx] + a[nx-1]) + cd21 = 1./6. * (a[nx+2] - 2. * a[nx+1] + a[nx]) + + # fourth central differences + cd40 = 1./120. * (a[nx-2] - 4. * a[nx-1] + + 6. * a[nx] - 4. * a[nx+1] + a[nx+2]) + cd41 = 1./120. * (a[nx-1] - 4. * a[nx] + + 6. * a[nx+1] - 4. * a[nx+2] + a[nx+3]) + nxold = nx + } + y[i] = s * (a[nx+1] + (s * s - 1.) * + (cd21 + (s * s - 4.) * cd41)) + + t * (a[nx] + (t * t - 1.) * + (cd20 + (t * t - 4.) * cd40)) + } + return +end + +procedure iievs3(x,y,n,a) # cubic spline evaluator + +real x[ARB] # x values, must be within [1,n] +real y[ARB] # interpolated values returned to user +int n # number of x values +real a[ARB] # basis spline coefficients - from a[0] to a[n+1] + +int i,nx,nxold +real s,c0,c1,c2,c3 + +begin + nxold = -1 + do i = 1,n { + nx = x[i] + s = x[i] - nx + if (nx != nxold) { + + # convert b-spline coeff's to poly. coeff's + c0 = a[nx-1] + 4. * a[nx] + a[nx+1] + c1 = 3. * (a[nx+1] - a[nx-1]) + c2 = 3. * (a[nx-1] - 2. * a[nx] + a[nx+1]) + c3 = -a[nx-1] + 3. * a[nx] - 3. * a[nx+1] + a[nx+2] + nxold = nx + } + y[i] = c0 + s * (c1 + s * (c2 + s * c3) ) + } + return +end diff --git a/math/interp/iif_spline.x b/math/interp/iif_spline.x new file mode 100644 index 00000000..a8574e0b --- /dev/null +++ b/math/interp/iif_spline.x @@ -0,0 +1,67 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +.help + procedure iif_spline + + This fits uniformly spaced data with a cubic spline. The spline is +given as basis-spline coefficients that replace the data values. + +Storage at call time: + +b[1] = second derivative at x = 1 +b[2] = first data point y[1] +b[3] = y[2] +. +. +. +b[n+1] = y[n] +b[n+2] = second derivative at x = n + +Storage after call: + +b[1] ... b[n+2] = the n + 2 basis-spline coefficients for the basis splines + as defined in P. M. Prenter's book Splines and Variational Methods, + Wiley, 1975. + +.endhelp + +procedure iif_spline(b,d,n) + +real b[ARB] # data in and also bspline coefficients out +real d[ARB] # needed for offdiagnol matrix elements +int n # number of data points + +int i + +begin + + d[1] = -2. + b[1] = b[1] / 6. + + d[2] = 0. + b[2] = (b[2] - b[1]) / 6. + + # Gaussian elimination - diagnol below main is made zero + # and main diagnol is made all 1's + do i = 3,n+1 { + d[i] = 1. / (4. - d[i-1]) + b[i] = d[i] * (b[i] - b[i-1]) + } + + # Find last b spline coefficient first - overlay r.h.s.'s + b[n+2] = ( (d[n] + 2.) * b[n+1] - b[n] + b[n+2]/6.) / (1. + + d[n+1] * (d[n] + 2.)) + + # back substitute filling in coefficients for b splines + do i = n+1, 3, -1 + b[i] = b[i] - d[i] * b[i+1] + + # note b[n+1] is evaluated correctly as can be checked + + # b[2] is already set since offdiagnol is 0. + + # evaluate b[1] + b[1] = b[1] + 2. * b[2] - b[3] + + return +end diff --git a/math/interp/iimail b/math/interp/iimail new file mode 100644 index 00000000..69c19e9f --- /dev/null +++ b/math/interp/iimail @@ -0,0 +1,213 @@ +From tody Wed Apr 20 18:50:45 1983 +Date: 20 Apr 1983 18:50-MST +To: vesa +Subject: ii comments +Cc: tody + +I read through "tnotes1", and felt that enough information was conveyed +to explain how the asi___ routines work, and how to use them. It is +good documentation for the package, though some sections are probably +to technical for the general user, and should probably be moved to an +appendix in the final documentation. + +I only have one minor suggestion, and that is that there seems little +reason to abbreviate NEAR, REFL, and PROJ, since only three characters +are saved in each case, and the other parameter names are all spelled out. + +I also looked at the code for ASIFIT. I noticed that the recursion +relations for the cubic spline are coded inline. I suspect that this +is actually LESS efficient than using a separate procedure, due to +the complex coeff offset expressions appearing in the inner loops +(unless the Fortran compiler does a very good job of optimizing loops). + +By using a separate procedure, for example, the statement + + coeff[COFF+i] = coeff[off2+i] * (coeff[COFF+i] - coeff[COFF+i-1]) + +can be converted to something like + + bcoeff[i] = acoeff[i] * (bcoeff[i] - bcoeff[i-1]) + +where in the main routine, there is a procedure call something like the +following: + + call fit_bspline3 (coeff[COFF], coeff[off2], ...) + +In other words, the routines that actually do the work do not need to +know that everything is saved in a single coeff array. This simplifies +the code, as well as making it run faster on some machines. The technique +has an additional advantage: by isolating the compute bound code into +a separate (small) procedure, it becomes easy to hand optimize that procedure +in assembler, should we wish to do so. + +From tody Wed Apr 20 19:48:51 1983 +Date: 20 Apr 1983 19:48-MST +To: vesa +Subject: first benchmarks +Cc: tody + +I ran some preliminary benchmarks on the asi code, assuming no bad pixels. +The benchmark routine is in math/interp. Results: + +512 by 512 image interpolation, no bad pixels, using ASIEVA: + + nearest 47 (cpu seconds) + linear 22 + poly3 50 + poly5 89 + spline3 94 + +This is quite encouraging. I expect that these timings can be decreased by +a factor of from 3 to 5, given a good Fortran compiler (the current UNIX +compiler is very poor), and possibly hand optimising critical sections of +the code in assembler. + +The nearest neighbor timings don't make any sense (that should be the +fastest interpolant). Are you using sqrt, or something like that? +The spline timings are worse than I expected, compared to poly3. May +be due to complex offset expressions, which F77 does not optimize. + +From tody Fri Feb 18 14:10:36 1983 +To: vesa +Subject: x code +Cc: tody + + +We want to make IRAF source code as readable as possible, therefore there +is a (currently unofficial) standard for ".x" code. Your code is pretty +close to the standard, but here are some things to consider: + + (1) It is preferable to enclose large blocks of documentation in + .help .. .endhelp sections, rather than making each line a comment. + This has two advantages: (1) it is awkward to edit text set off with + a # on each line, and (2) the help utilities can only extract + source code documentation which is set in help blocks. It is not + necessary to use text formatter commands in help blocks, i.e., + you can simply type it in as it will appear on the terminal. + + In the same vein, it is better to put the discussion of the input + and output parameters in the help text, rather than in the + declarations section of the procedure. Placing each input and + output parameter on a single line, followed at the right by a short + comment, gets the job done without cluttering up the declarations + section. More extensive comments should be placed in the help text. + + + (2) It is mostly up to the programmer to decide where to put whitespace + to make the structure of a program clearest. Whitespace should + be inserted in a statement to make the important logical components + of the statement stand out. + + In particular, there should be a space after every keyword, before + the left parenthesis. If the statement is compound, the left brace + should be set off from the right paren by a space. In short + expressions, operators should be set off by whitespace. + + Adding whitespace between a function name and the left paren is + optional, as is whitespace between the elements of an argument list + (after commas). + + Blank lines should be used to set off sections of code, just as + blank lines are used to separate paragraphs of text. + + It is possible to overcomment code, making it hard to read. If + more than a third of the lines in a procedure are comments, the + code is probably overcommented. Place detailed discussion of the + algorithm in the help text. + + + (3) Indenting structures by full tab stops is acceptable, but 4 space + indents are preferred (the code rapidly runs off the right side of + the screen if full tab indents are used, encouraging use of too + short identifiers and omission of whitespace). The autoindent + feature of VI makes this easy (in .login file: 'set ai sw=4'). + + + (4) The standard form of the if .. else construct is + + if (expr) { + stmt + } else { + stmt + + or (if "stmt" is several lines, and whitespace is desired) + + if (expr) { + stmt + + } else { + stmt + + rather than + + if (expr) { + stmt + } + else { + stmt + +From tody Fri Feb 18 15:21:33 1983 +Date: 18 Feb 1983 15:21-MST +To: vesa +Subject: spec changes +Cc: tody + + +Now that Garth, Harvey, Steve and others have had a chance to comment on the +specs for the interpolator routines, it is time to make the final revisions +to the specifications. Often ones perspective on a problem changes when one +actually begins coding, and even the most carefully prepared specifications +need to be changed. Please feel free to suggest further changes. + + +Modifications + + (1) Everyone wanted the interpolators to be able to generate estimated + values to replace indefinites (as an option). + + (2) Steve wanted sinc function interpolation added to the set of + interpolators. I would also like to try this for undersampled + CCD data. A fairly high "order" (basis function size) seems + justified for this interpolant. + + (3) A surface integration routine is desired (MSIGRL). + + +SET Routines + + I think we should generalize the form of the "set" routines to permit +options such as (1), and others we may wish to add in the future (such +as telling the interpolator not to check for indefinites in the data). + +Set routines are also used in the FIO, IMIO, VSIO, and GIO packages in +the program interface. For example, in FIO, the number of buffers for a +file could be set by a call of the form + + call fset (fd, NBUFFERS, 2) + +The form of a set call is therefore + + call ___set (object, parameter, value) + +i.e., + + call asiset (coeff, II_INTERPOLANT, II_DIVDIF3) + call asiset (coeff, II_TYBNDRY, II_NEAREST) + call asiset (coeff, II_REPLACE_INDEFS, YES) + +Note that the prefix "II_" must be added to all Image Interpolation +defined parameters and parameter values. This is unattractive, but necessary +to avoid redefinitions of identifiers in other packages, such as IMIO (IMIO +has the same boundary extension options as II). + + +Surface Integration + + I suggest the following calling sequence for the surface integration +routine: + + integral = msigrl (x, y, npts, coeff) + +where the surface to be integrated is defined by the closed curve {x[i],y[i]}, +where 1 <= i <= npts, defined by the caller. + diff --git a/math/interp/iipol_terp.x b/math/interp/iipol_terp.x new file mode 100644 index 00000000..030d1964 --- /dev/null +++ b/math/interp/iipol_terp.x @@ -0,0 +1,41 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +.help + Procedure iipol_terp + +A polynomial interpolator with x and y arrays given. +This algorithm is based on the Newton form as described in +de Boor's book, A Practical Guide to Splines, 1978. +There is no error checking - this is meant to be used only by calls +from more complete routines that take care of such things. + +Maximum number of terms is 6. +.endhelp + +real procedure iipol_terp(x,y,n,x0) + +real x[ARB],y[ARB] # x and y array +real x0 # desired x +int n # number of points in x and y = number of + # terms in polynomial = order + 1 + +int k,i +real d[6] + +begin + + # Fill in entries for divided difference table. + do i = 1,n + d[i] = y[i] + + # Generate divided differences + do k = 1,n-1 + do i = 1,n-k + d[i] = (d[i+1] - d[i])/(x[i+k] - x[i]) + + # Shift divided difference table to center on x0 + do i = 2,n + d[i] = d[i] + d[i-1] * (x0 - x[i]) + + return(d[n]) +end diff --git a/math/interp/interp.h b/math/interp/interp.h new file mode 100644 index 00000000..e54fc1b5 --- /dev/null +++ b/math/interp/interp.h @@ -0,0 +1,19 @@ +# File of definitions for interpolator package. + +# Interpolator Types: + +define IT_NEAREST 1 +define IT_LINEAR 2 +define IT_POLY3 3 +define IT_POLY5 4 +define IT_SPLINE3 5 +define ITNIT 5 # number of types + +# Size of header etc. used for part of coeff dimension + +define SZ_ASI 20 + +# Total number of points used in spline interpolation of of bad pixels by +# subroutine arbpix. + +define SPLPTS 16 diff --git a/math/interp/interpdef.h b/math/interp/interpdef.h new file mode 100644 index 00000000..e54fc1b5 --- /dev/null +++ b/math/interp/interpdef.h @@ -0,0 +1,19 @@ +# File of definitions for interpolator package. + +# Interpolator Types: + +define IT_NEAREST 1 +define IT_LINEAR 2 +define IT_POLY3 3 +define IT_POLY5 4 +define IT_SPLINE3 5 +define ITNIT 5 # number of types + +# Size of header etc. used for part of coeff dimension + +define SZ_ASI 20 + +# Total number of points used in spline interpolation of of bad pixels by +# subroutine arbpix. + +define SPLPTS 16 diff --git a/math/interp/mkpkg b/math/interp/mkpkg new file mode 100644 index 00000000..55ab29c6 --- /dev/null +++ b/math/interp/mkpkg @@ -0,0 +1,21 @@ +# Old image interpolator tools library. + +$checkout libinterp.a lib$ +$update libinterp.a +$checkin libinterp.a lib$ +$exit + +libinterp.a: + arbpix.x interpdef.h asidef.h + arider.x interpdef.h + arival.x interpdef.h + asider.x interpdef.h asidef.h + asieva.x interpdef.h asidef.h + asifit.x interpdef.h asidef.h + asigrl.x interpdef.h asidef.h + asiset.x interpdef.h asidef.h + asival.x interpdef.h asidef.h + iieval.x + iif_spline.x + iipol_terp.x + ; diff --git a/math/interp/noteari b/math/interp/noteari new file mode 100644 index 00000000..d42ed48b --- /dev/null +++ b/math/interp/noteari @@ -0,0 +1,64 @@ + Using the Random One Dimensional Interpolation Routines + A Quick Note for Programmers + + +I. GENERAL NOTES: + + 1. Defines are found in file interpdef.h. The routines + are in library interplib.a. + + 2. All pixels are assumed to be good. Except for routine + arbpix. + + 3. This is for uniformly spaced data -- thus for the i-th + data value, y[i], the corresponding x[i] = i by assumption. + + 4. All x references are assumed to land in the closed interval + 1 <= x <= NPTS, where NPTS is the number of points in the + data array. If x is outside this range, the result is + not specified. + + +II. PROCEDURES: + + + **** Subroutine to replace INDEF's with interpolated values. + * + * arbpix (datain, n, dataout, interpolator_type) + * + * real datain[ARB] # data_in array + * int n # number of points in data_in + * real dataout[ARB] # array out, may not be same as + * # data_in + * int interpolator_type + * + * The interpolator type can be set to: + * + * IT_NEAREST nearest neighbor + * IT_LINEAR linear interpolation + * IT_POLY3 interior polynomial 3rd order + * IT_POLY5 interior polynomial 5th order + **** IT_SPLINE3 cubic natural spline + + + **** Subroutine to return interpolated value + * + * real arival (x, datain, n , interpolator_type) + * + * real x + * real datain[ARB] # data array + * int n # number of points in data array + **** int interpolator_type # see above + + + **** Subroutine to evaluate derivatives at a point. + * + * arider (x, datain, n, derivs, nder, interpolator_type) + * + * real x + * real datain[ARB] # data array + * int n # number of points in data array + * real derivs[ARB] # output array containing derivatives + * # derivs[1] is function value + * int nder # input: nder - 1 derivs are evaluated + **** int interpolator_type # see above diff --git a/math/interp/noteasi b/math/interp/noteasi new file mode 100644 index 00000000..ae22df83 --- /dev/null +++ b/math/interp/noteasi @@ -0,0 +1,101 @@ + Using the Sequential One Dimensional Interpolation Routines + A Quick Note for Programmers + + +I. GENERAL NOTES: + + 1. Defines are found in file interpdef.h. The routines + are in library interplib.a. + + 2. All pixels are assumed to be good. Except for routine + arbpix. + + 3. This is for uniformly spaced data -- thus for the i-th + data value, y[i], the corresponding x[i] = i by assumption. + + 4. All x references are assumed to land in the closed interval + 1 <= x <= NPTS, where NPTS is the number of points in the + data array. If x is outside this range, the result is + not specified. + + 5. A storage area and work space array must be provided. + It is of type real and the size is 2 * NPTS + SZ_ASI + + + +II. PROCEDURES: + + + **** Subroutine to replace INDEF's with interpolated values. + * + * arbpix (datain, n, dataout, interpolator_type) + * + * real datain[ARB] # data_in array + * int n # number of points in data_in + * real dataout[ARB] # array out, may not be same as + * # data_in + * int interpolator_type + * + * The interpolator type can be set to: + * + * IT_NEAREST nearest neighbor + * IT_LINEAR linear interpolation + * IT_POLY3 interior polynomial 3rd order + * IT_POLY5 interior polynomial 5th order + **** IT_SPLINE3 cubic natural spline + + + **** Subroutine to set interpolator type. + * + * asiset(coeff,interpolator_type) + * + * real coeff[ARB] # work + storage array, dimension + * # 2 * NPTS + SZ_ASI + **** int interpolator_type # see above + + + **** Subroutine to fit interpolator to data. + * + * asifit (datain, npts, coeff) + * + * real datain[ARB] # data array + * int npts # number of points in data array + **** real coeff[ARB] # work+storage area dim 2*npts + SZ_ASI + + + **** Subroutine to return interpolated value after fitting. + * + * real asival (x, coeff) + * + * real x + **** real coeff[ARB] + + + **** Subroutine to evaluate an ordered sequence of values. + * + * asieva (x, y, n, coeff) + * + * real x[ARB] # input array of x values + * real y[ARB] # output array of interpolated values + * int n # number of x values + **** real coeff[ARB] + + + **** Subroutine to evaluate derivatives at a point. + * + * asider (x, derivs, nder, coeff) + * + * real x + * real derivs[ARB] # output array containing derivatives + * # derivs[1] is function value + * int nder # input: nder - 1 derivatives are evaluated + **** real coeff[ARB] + + + **** Subroutine to integrate the interpolated data. + * + * real asigrl (a, b, coeff) + * + * real a # lower limit + * real b # upper limit + **** real coeff[ARB] diff --git a/math/interp/notes3 b/math/interp/notes3 new file mode 100644 index 00000000..e7121c99 --- /dev/null +++ b/math/interp/notes3 @@ -0,0 +1,216 @@ + Notes on the Interpolator Package for IRAF Version 2 + + Some parts of the image interpolator package are available. +This is a reminder of changes to the document TRP prepared by +Doug Tody. Also it is a place to keep track of why things were +done the way they were. The individual procedures (subroutines) +will be discussed; making it easy to add, subtract and change things. + +************************************************************************** +** ** +** ** + + GENERAL NOTES FOR 1-D PACKAGE: + +1. Handling of errors. + + Some routines call the error routines of the IRAF package. + Logically it is reasonable to trap silly errors at certain points. + For instance if the interpolator type is set to a nonexistant type + an error call is generated. + + To the programmer, the two most important error features are: + + 1. All pixels are assumed to be good. + + 2. All x references are assumed to land in the closed interval + 1 <= x <= NPTS. + + Point 2 probably needs some explanation, because at first hand + it seems rather easy to test for x < 1 || x > NPTS and take + appropriate action. There are two objections to including this + test. First it is often possible to generate x values in such + a way that they are guaranteed to lie in the given interval. + Testing again within the interpolator package is a waste of + computer time. Second the interpolator package becomes too large + if code is included to handle out of bounds references and such + code duplicates other parts of the IRAF package that handles + out of bound values. The alternative is to return INDEF when + x is out of bounds. This is not very appealing for many + routine applications such as picture shifting. + +2. Size of "coeff" array: + + coeff is a real 1-d array of size -- 2 * NPTS + SZ_ASI + + where SZ_ASI is currently 30. It, like the other defines needed for + the interp package, is available in the file -- terpdef.h . + + This coeff array serves as work area for asifit, a structure to save + interpolator type, flags etc., an array to save data points or b-spline + coefficients as needed. + + The polynomial interpolators are saved as the data points and the +interpolation is done "on the fly". Rather than compute polynomial +coefficients before-hand, the interpolated value is calculated from +the data values at time of request. (The actual calculation is done +by adding and multiplying things in an order suggested by Everett's +interpolation formula rather than the polynomial coefficient form. +The resulting values are the same up to differences due to numerical +roundoff. As a practical matter, all forms work except polynomials +shifted far from the point of interest. Polynomials shifted to a +center far away suffer from bad roundoff errors.) + + The splines must be calculated at a single pass because (formally) +they are global (i.e use all of the data points). The b-splines (b for +basis) are just certain piecewise polynomials for solving the interpo- +lation problem, it is convenient to calculate the coefficients of these +b-splines and to save these in the array. + +3. Method of solution: + +************************************************************************** +** ** +** ** +** ** + PROCEDURES: + +*************************************************************************** + + ***** asiset(coeff,interptype) ***** + +The interpolator type can be set to: + + IT_NEAREST nearest neighbor + IT_LINEAR linear interpolation + IT_POLY3 interior polynomial 3rd order + IT_POLY5 interior polynomial 5th order + IT_SPLINE3 cubic natural spline + +Subroutines needed: + + NONE + + + +************************************************************************* + + ******* asifit (data, npts, coeff) ***** + + This fits the interpolator to the data. Takes care of bad points by +replacing them with interpolated values to generate a +uniformly spaced array. For polynomial interpolators, this array is +sent to the evaluation routines. For the spline interpolator, the uniform +array is fitted and the basis-spline coefficients are saved. + +Interior polynomial near boundary: + Array is extended according to choice of boundary conditions. +Out of bounds values and values near the edge depend on choice of +boundary conditions. + +Spline near boundary: + Natural (second derivative equals zero at edge) spline is fitted to +data only. The resulting function is reflected, wrapped etc. to generate +out of bounds values. Out of bounds values depend on choice of boundary +conditions but values within the array near the edges do not. + For the wrap boundary condition, the one segment connecting point +number n to point 1 is not generated by the spline fit. Instead the +polynomial that is continuous and has a continuous first and second +derivative n is inserted. In this one case, the first and second derivative +are not continuous at 1. + +Replacing bad points: + The same kind of interpolator is used to replace bad points as that +used to generate interpolated values. The boundary conditions are not +handled in the same manner. Bad points near the edge are replaced by +using lower order polynomial interpolators if there are not enough +points on both sides of the bad point. Bad points that are not +straddled by good points are assigned the value of the nearest good +point. This is done in the temporary array before the uniformly spaced +interpolator is applied. + +subroutines needed: + + real iipint(x,y,nterm,x0) + Returns interpolated value at x0 using polynomial with + <nterm> terms to connect the <nterm> data pairs + in x,y. + + cubspl(tau,c,n,ibcbeg,ibcend) + Unmodified cubic spline interpolator for non-uniformly + spaced data taken from C. de Boor's book: + A Practical Guide to Splines, (1978 + Springer-Verlag). + + iifits(b,d,n) + Fits cubic spline to uniformly spaced array of y values + using a basis-spline representation for + the spline. + + + +***************************************************************************** + + ***** y = asival (x, coeff) ***** + + Subroutine to return interpolated value after fitting. + +Subroutines needed: + + real iivalu (x, coeff) + Returns interpolated value assuming x is in array. Also + has no code for bad point propagation. + + + + +*************************************************************************** + + ***** asieva (x, y, n, coeff) ***** + + Given an ordered array of x values returns interpolated values in y. +This is needed for speed. Reduces number of subroutine calls per point, +reduces the number of tests for out of bounds, and bad pixel +propagation can be made more efficient. + +Subroutines needed: + + real iivalu (x, coeff) + + + +************************************************************************** + + ***** asider (x, derivs, nderiv, coeff) ***** + + This evaluates derivatives at point x. Returns the first +nderiv - 1 derivatives. derivs[1] is function value, derivs[2] is +first derivative, derivs[3] is second derivative etc. + There is no check that nderiv is a reasonable value. Polynomials +have a limited number of non-zero derivatives, if more are asked for +they come back as zero. + +Subroutines needed: + + iivalu (x, coeff) + + iiderv (x, derivs, nderiv, coeff) + Assumes x lands in array returns derivatives. No handling + of bad points. + +**************************************************************************** + + ****** integral = asigrl (a, b, coeff) ******* + + This integrates the array from to a to b. The boundary conditions are +incorporated into the code so if function values can be returned from all +points between a and b, the integral will be defined. Thus for WRAP, +REFLECT and PROJECT boundary conditions, the distance from a to b can be +as large as 3 times the array length if a and b are chosen appropriately. +For these boundary conditions, attempts to extend beyond +/- one cycle +produce INDEF. + +Subroutines needed: + + iiigrl (a, b, coeff) + returns integral when a and b land in array. diff --git a/math/interp/overview b/math/interp/overview new file mode 100644 index 00000000..4ad813ca --- /dev/null +++ b/math/interp/overview @@ -0,0 +1,43 @@ +Interpolation Package + + The routines here are intended to allow programmers to interpolate +one dimensional arrays and two dimensional +images in the IRAF system. The routines are +mathematical in nature and are not heavily tied into the IRAF package. +The advantages of having locally written interpolators include the +ability to optimize for uniformly spaced data and the possibility of +adding features that are useful for the final application. The +one major feature that has been implemented is the ability to change +the kind of interpolator used at run time by carrying along a variable +that gives the interpolator type in the higher level code. + + The kinds of interpolators include: + 1. Nearest neighbor + 2. Linear + 3. Interior polynomial order = 3 + 4 Interior polynomial order = 5 + 5 Spline order = 3 + + The package is divided into array interpolators and image +interpolators. Within the array interpolators, there is further +division into sequential and random interpolators. The sequential +interpolators are optimized for returning many values as is the +case when an array is shifted or when an array is 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. + + The original specification called for the interpolator package +to deal with indefinite valued pixels and references out of bounds. +A version of the array sequential interpolators was written that handled +both situations allowing some options to the user. The extension of +these features to the two dimensional case becomes more dependent on +details within IRAF and less mathematical in nature. It would also +duplicate other code within IRAF that deals with out of bound references. +Finally the need to test for the various situations slows down the operation. +So all features referring to indefinite valued pixels and out of bound +references have been removed. The user is responsible for ensuring that +"x" values land within the array. The code does not test to see if x +is in bounds. A routine has been added that interpolates +an array with indefinites and returns an array with all legitimate +real numbers. diff --git a/math/interp/usernote b/math/interp/usernote new file mode 100644 index 00000000..55483022 --- /dev/null +++ b/math/interp/usernote @@ -0,0 +1,118 @@ + Using the Sequential One Dimensional Interpolation Routines + A Quick Note for Programmers + + +I. GENERAL NOTES: + + 1. Defines are found in file interpdef.h. The routines + are in library interplib.a. + + 2. All pixels are assumed to be good. Except for routine + arbpix. + + 3. This is for uniformly spaced data -- thus for the i-th + data value, y[i], the corresponding x[i] = i by assumption. + + 4. All x references are assumed to land in the closed interval + 1 <= x <= NPTS, where NPTS is the number of points in the + data array. If x is outside this range, the result is + not specified. + + 5. A storage area and work space array must be provided. + It is of type real and the size is 2 * NPTS + SZ_ASI + + + +II. PROCEDURES: + + + **** Subroutine to select the interpolator type + * + * interpolator_type = clginterp (param) + * + * char param[ARB] + * + * The parameter prompt string is sent to the CL and the + * routine evaluates the returned string to select one of the + * interpolator types given in interpdef.h. The input string + * may be one of the following with possible abbreviations: + * nearest + * linear + * 3rd order poly + * 5th order poly + **** cubic spline + + + **** Subroutine to replace INDEF's with interpolated values. + * + * arbpix (datain, n, dataout, interpolator_type) + * + * real datain[ARB] # data_in array + * int n # number of points in data_in + * real dataout[ARB] # array out, may not be same as + * # data_in + * int interpolator_type + * + * The interpolator type can be set to: + * + * IT_NEAREST nearest neighbor + * IT_LINEAR linear interpolation + * IT_POLY3 interior polynomial 3rd order + * IT_POLY5 interior polynomial 5th order + **** IT_SPLINE3 cubic natural spline + + + **** Subroutine to set interpolator type. + * + * asiset(coeff,interpolator_type) + * + * real coeff[ARB] # work + storage array, dimension + * # 2 * NPTS + SZ_ASI + **** int interpolator_type # see above + + + **** Subroutine to fit interpolator to data. + * + * asifit (datain, npts, coeff) + * + * real datain[ARB] # data array + * int npts # number of points in data array + **** real coeff[ARB] # work+storage area dim 2*npts + SZ_ASI + + + **** Subroutine to return interpolated value after fitting. + * + * real asival (x, coeff) + * + * real x + **** real coeff[ARB] + + + **** Subroutine to evaluate an ordered sequence of values. + * + * asieva (x, y, n, coeff) + * + * real x[ARB] # input array of x values + * real y[ARB] # output array of interpolated values + * int n # number of x values + **** real coeff[ARB] + + + **** Subroutine to evaluate derivatives at a point. + * + * asider (x, derivs, nder, coeff) + * + * real x + * real derivs[ARB] # output array containing derivatives + * # derivs[1] is function value + * int nder # input: nder - 1 derivatives are evaluated + **** real coeff[ARB] + + + **** Subroutine to integrate the interpolated data. + * + * real asigrl (a, b, coeff) + * + * real a # lower limit + * real b # upper limit + **** real coeff[ARB] |