Initial commit

1 files changed, 293 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.