Initial commit

19 files changed, 1380 insertions, 0 deletions

diff --git a/math/surfit/doc/iscoeff.hlp b/math/surfit/doc/iscoeff.hlp
new file mode 100644
index 00000000..ea0b292c
--- /dev/null
+++ b/math/surfit/doc/iscoeff.hlp
@@ -0,0 +1,63 @@
+.help iscoeff Apr85 "Surfit Package"
+.ih
+NAME
+iscoeff - get the number and values of the surface coefficients
+.ih
+SYNOPSIS
+iscoeff (sf, coeff, ncoeff)
+
+.nf
+pointer	sf		# surface descriptor
+real	coeff[ncoeff]	# coefficient array
+int	ncoeff		# number of coefficients
+.fi
+.ih
+ARGUMENTS
+.ls sf   
+Pointer to the surface descriptor.
+.le
+.ls coeff
+Array of coefficients.
+.le
+.ls ncoeff
+The number of coefficients.
+.le
+
+.ih
+DESCRIPTION
+ISCOEFF fetches the coefficient array and the number of coefficients from
+the surface descriptor structure. A 1-D array ncoeff elements long
+is required to hold the coefficients where ncoeff is defined as follows:
+
+.nf
+	SF_LEGENDRE:		nxcoeff = xorder
+				nycoeff = yorder
+				ncoeff = xorder * yorder        xterms = yes
+				ncoeff = nycoeff + nxcoeff - 1  xterms = no
+
+	SF_CHEBYSHEV:		nxcoeff = xorder
+				nycoeff = yorder
+				ncoeff = xorder * yorder        xterms = yes
+				ncoeff = nycoeff + nxcoeff - 1  xterms = no
+
+	SF_SPLINE3:		nxcoeff = xorder + 3
+				nycoeff = yorder + 3
+				ncoeff = (xorder + 3) * (yorder + 3)
+
+	SF_SPLINE1:		nxcoeff = xorder + 1
+				nycoeff = yorder + 1
+				ncoeff = (xorder + 1) * (yorder + 1)
+
+.fi
+
+The coefficient of basis function B(i,x) * B(j,y) will be stored in
+element (i - 1) *
+nycoeff + j of the array coeff if the xterms parameter was set
+to yes by ISINIT. Otherwise the nycoeff y-term coefficients will be output
+first followed by the (nxcoeff - 1) x-term coefficients.
+
+.ih
+NOTES
+.ih
+SEE ALSO
+.endhelp
diff --git a/math/surfit/doc/iseval.hlp b/math/surfit/doc/iseval.hlp
new file mode 100644
index 00000000..53853cfe
--- /dev/null
+++ b/math/surfit/doc/iseval.hlp
@@ -0,0 +1,34 @@
+.help iseval Apr85 "Surfit Package"
+.ih
+NAME
+iseval -- evaluate the fitted surface at x and y
+.ih
+SYNOPSIS
+y = iseval (sf, x, y)
+
+.nf
+pointer	sf		# surface descriptor
+real	x		# x value, 1 <= x <= ncols
+real	y		# y value, 1 <= y <= nlines
+.fi
+.ih
+ARGUMENTS
+.ls sf      
+Pointer to the surface descriptor structure.
+.le
+.ls x, y
+X and y values at which the surface is to be evaluated.
+.le
+.ih
+DESCRIPTION
+Evaluate the surface at the specified value of x and y. ISEVAL is a real
+valued function which returns the fitted value.
+.ih
+NOTES
+ISEVAL uses the coefficient array stored in the surface descriptor structure.
+Checking for out of bounds x and y values is the responsibility of the calling
+program.
+.ih
+SEE ALSO
+isvector
+.endhelp
diff --git a/math/surfit/doc/isfree.hlp b/math/surfit/doc/isfree.hlp
new file mode 100644
index 00000000..04f52b5f
--- /dev/null
+++ b/math/surfit/doc/isfree.hlp
@@ -0,0 +1,27 @@
+.help isfree Apr85 "Surfit Package"
+.ih
+NAME
+isfree -- free the surface grid descriptor structure
+.ih
+SYNOPSIS
+isfree (sf)
+
+.nf
+pointer	sf		# surface descriptor
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface descriptor structure.
+.le
+.ih
+DESCRIPTION
+Frees the surface descriptor structure and all the vectors and arrays used
+by the numerical routines.
+.ih
+NOTES
+ISFREE should be called after the surface fit is complete.
+.ih
+SEE ALSO
+isinit
+.endhelp
diff --git a/math/surfit/doc/isinit.hlp b/math/surfit/doc/isinit.hlp
new file mode 100644
index 00000000..e0f2123c
--- /dev/null
+++ b/math/surfit/doc/isinit.hlp
@@ -0,0 +1,61 @@
+.help isinit Apr85 "Surfit Package"
+.ih
+NAME
+isinit -- initialize surface grid descriptor
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+.nf
+isinit (sf, surface_type, xorder, yorder, xterms, ncols, nlines)
+.fi
+
+.nf
+pointer	sf		# surface descriptor
+int	surface_type	# surface function
+int	xorder		# order of function in x
+int	yorder		# order of function in y
+int	xterms		# include cross-terms? (YES/NO)
+int	ncols		# number of columns in the surface grid
+int	nlines		# number of lines in the surface grid
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface descriptor structure.
+.le
+.ls surface_type
+Fitting function. Permitted values are SF_LEGENDRE and SF_CHEBYSHEV for
+the Legendre and Chebyshev polynomials and SF_SPLINE1 and SF_SPLINE3
+for the linear and bicubic splines.
+.le
+.ls xorder, yorder
+Order of the polynomial to be fit in x and y or the number of spline pieces to 
+be fit in x and y. The orders must be greater than or equal to 1.
+.le
+.ls xterms
+Include cross-terms? If xterms = YES coefficients are fit to terms containing
+the cross products of x and y polynomials. Xterms defaults to YES for the
+spline functions.
+.le
+.ls ncols
+The number of columns in the surface grid. The surface is assumed to lie
+on a rectangular grid such that 1 <= x <= ncols.
+.le
+.ls nlines
+The number of lines in the surface to be grid. The surface is assumed to lie
+on a rectangular grid such that 1 <= y <= nlines.
+.le
+.ih
+DESCRIPTION
+ISINIT allocates space for the surface descriptor and the arrays and vectors
+used by the numerical routines. It initializes all arrays and vectors to zero,
+calculates and stores the basis functions in x and y
+and returns the surface descriptor to the calling routine.
+.ih
+NOTES
+ISINIT must be the first SURFIT routine called.
+.ih
+SEE ALSO
+isfree
+.endhelp
diff --git a/math/surfit/doc/islaccum.hlp b/math/surfit/doc/islaccum.hlp
new file mode 100644
index 00000000..48481cf0
--- /dev/null
+++ b/math/surfit/doc/islaccum.hlp
@@ -0,0 +1,60 @@
+.help islaccum Apr85 "Surfit Package"
+.ih
+NAME
+islaccum -- accumulate a surface grid line into the fit
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+islaccum (sf, cols, lineno, z, w, npts, wtflag)
+
+.nf
+pointer	sf		# surface grid descriptor
+int	cols[npts]	# column values, 1 <= col[i] <= ncols
+int	lineno		# number of line, 1 <= lineno <= nlines 
+real	z[npts]		# surface values on lineno at cols
+real	w[npts]		# array of weights
+int	npts		# number of surface values, npts <= ncols
+int	wtflag		# type of weighting desired
+.fi
+.ih
+ARGUMENTS
+.ls sf      
+Pointer to the surface grid descriptor structure.
+.le
+.ls cols
+The column numbers of surface grid points on line lineno to be added to the
+dataset.
+.le
+.ls lineno
+The line number of the surface grid line to be added to the data set.
+.le
+.ls z      
+The data values.
+.le
+.ls w
+The array of weights for the data points.
+.le
+.ls npts
+The number of surface grid values on line number lineno to be added to the
+data set.
+.le
+.ls wtflag
+Type of weighting. The options are SF_USER and SF_UNIFORM. If wtflag
+equals SF_USER the weight for each data point is supplied by the user
+and points with zero-valued weights are not included in the fit.
+If wtflag equals SF_UNIFORM the routine sets the weights to 1.
+.le
+.ih
+DESCRIPTION
+ISLACCUM computes the contribution of each data point to the normal
+equations in x, sums that contribution into the appropriate arrays and
+vectors and stores these intermediate results for use by ISLSOLVE.
+.ih
+NOTES
+Checking for out of bounds col values and INDEF valued data points is
+the responsibility of the calling program.
+.ih
+SEE ALSO
+isinit, islfit, islrefit, islsolve, islzero, issolve, isfree
+.endhelp
diff --git a/math/surfit/doc/islfit.hlp b/math/surfit/doc/islfit.hlp
new file mode 100644
index 00000000..5adf0a59
--- /dev/null
+++ b/math/surfit/doc/islfit.hlp
@@ -0,0 +1,67 @@
+.help islfit Apr85 "Surfit Package"
+.ih
+NAME
+islfit -- fit a surface grid line
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+islfit (sf, cols, lineno, z, w, npts, wtflag, ier)
+
+.nf
+pointer	sf		# surface grid descriptor
+real	cols[npts]	# array of column numbers, 1 <= cols[i] <= ncols
+int	lineno		# number of surface grid line to be added
+real	z[npts]		# data values
+real	w[npts]		# weight array
+int	npts		# number of data points, npts <= ncols
+int	wtflag		# type of weighting
+int	ier		# error code
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface grid descriptor structure.
+.le
+.ls cols
+The column numbers of surface grid points to be added to the dataset.
+.le
+.ls lineno
+The line number of the surface grid line to be added to the dataset.
+.le
+.ls z      
+Array of data values.
+.le
+.ls w
+Array of weights.
+.le
+.ls npts
+Number of data points.
+.le
+.ls wtflag
+Type of weighting. The options are SF_USER and SF_UNIFORM. If wtflag =
+SF_USER individual weights for each data point are supplied by the calling
+program and points with zero-valued weights are not included in the fit.
+If wtflag = SF_UNIFORM, all weights are assigned values of 1.
+.le
+.ls ier     
+Error code for the fit. The options are OK, SINGULAR and NO_DEG_FREEDOM.
+If ier = SINGULAR, the numerical routines will compute a solution but one
+or more of the coefficients will be zero. If ier = NO_DEG_FREEDOM there
+were too few data points to solve the matrix equations and the routine
+returns without fitting the data.
+.le
+.ih
+DESCRIPTION
+ISLFIT zeroes the appropriate arrays and vectors,
+computes the contribution of each data point to the normal equations
+in x and accumulates it into the appropriate array and vector elements.  The
+x coefficients are stored for later use by ISSOLVE.
+.ih
+NOTES
+Checking for out of bounds col values and INDEF values pixels is the
+responsibility of the user.
+.ih
+SEE ALSO
+isinit, islrefit, islaccum, islsolve, islzero, issolve, isfree
+.endhelp
diff --git a/math/surfit/doc/islrefit.hlp b/math/surfit/doc/islrefit.hlp
new file mode 100644
index 00000000..c34476d6
--- /dev/null
+++ b/math/surfit/doc/islrefit.hlp
@@ -0,0 +1,51 @@
+.help islrefit Aug85 "Gsurfit Package"
+.ih
+NAME
+islrefit -- refit surface grid line  assuming old cols, w vector
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+islrefit (sf, cols, lineno, z, w)
+
+.nf
+pointer	sf			# surface grid descriptor
+int	cols[ARB]		# columns to be fit, 1 <= cols[i] <= ncols
+int	lineno			# line number of surface grid line to be fit
+real	z[ARB]			# array of data values
+real	w[ARB]			# array of weights
+.fi
+.ih
+ARGUMENTS
+.ls sf          
+Pointer to the surface grid descriptor structure.
+.le
+.ls cols
+The column numbers of the surface grid line to be added to the dataset.
+.le
+.ls lineno
+The number of the surface grid line to be added to the dataset.
+.le
+.ls z    
+Array of data values.
+.le
+.ls w    
+Array of weights.
+.le
+.ih
+DESCRIPTION
+In some applications the cols, and w values remain unchanged from fit
+to fit and only the z values vary. In this case it is redundant to
+reaccumulate the matrix and perform the Cholesky factorization for each
+succeeding surface grid line. ISLREFIT
+zeros and reaccumulates the vector on the right hand side of the matrix
+equation and performs the forward and back substitution phase to fit for
+a new coefficient vector.
+.ih
+NOTES
+It is the responsibility of the calling program to check for out of bounds
+column values and INDEF valued pixels.
+.ih
+SEE ALSO
+isinit, islfit, islaccum, islsolve, islzero, issolve, isfree
+.endhelp
diff --git a/math/surfit/doc/islsolve.hlp b/math/surfit/doc/islsolve.hlp
new file mode 100644
index 00000000..752cbc19
--- /dev/null
+++ b/math/surfit/doc/islsolve.hlp
@@ -0,0 +1,45 @@
+.help islsolve Apr85 "Surfit Package"
+.ih
+NAME
+islsolve -- solve the equations for a surface grid line
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+islsolve (sf, lineno, ier)
+
+.nf
+pointer	sf		# surface grid descriptor
+int	lineno		# surface grid line number
+int	ier		# error code
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface grid descriptor structure.
+.le
+.ls lineno
+The number of the surface grid line to be added to the dataset.
+.le
+.ls ier      
+Error code returned by the fitting routines. The options are OK, SINGULAR,
+and NO_DEG_FREEDOM. If ier = SINGULAR the matrix is singular, ISLSOLVE
+will compute a solution to the normal equations but one or more of the
+coefficients will be zero. If ier = NO_DEG_FREEDOM, too few data points
+exist for a reasonable solution to be computed and ISLSOLVE returns without
+fitting the data.
+.le
+.ih
+DESCRIPTION
+ISLSOLVE computes the Cholesky factorization of the data matrix and
+solves for the coefficients
+of the x fitting function by forward and back substitution.
+An error code is
+returned by ISLSOLVE if it is unable to solve the normal equations as
+formulated.
+.ih
+NOTES
+.ih
+SEE ALSO
+isinit, islfit, islrefit, islaccum, islzero, issolve, isfree
+.endhelp
diff --git a/math/surfit/doc/islzero.hlp b/math/surfit/doc/islzero.hlp
new file mode 100644
index 00000000..cbdb9515
--- /dev/null
+++ b/math/surfit/doc/islzero.hlp
@@ -0,0 +1,32 @@
+.help islzero Apr85 "Surfit Package"
+.ih
+NAME
+islzero -- set up for a new surface grid line fit
+.ih
+SYNOPSIS
+islzero (sf, lineno)
+
+.nf
+pointer	sf		# surface grid descriptor
+int	lineno		# surface grid line number
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface grid descriptor structure.
+.le
+.ls lineno
+Line number of the surface grid line fit to be zeroed.
+.le
+.ih
+DESCRIPTION
+ISLZERO zeros the appropriate arrays and vectors for the surface grid line
+number line number. The line can then be refit by calls to ISLACCUM and
+ISLSOLVE. Alternatively a single call to ISLFIT can perform the combined
+functions of ISLZERO, ISLACCUM and ISLSOLVE.
+.ih
+NOTES
+.ih
+SEE ALSO
+isinit, islfit, islrefit, islaccum, islsolve, issolve
+.endhelp
diff --git a/math/surfit/doc/isreplace.hlp b/math/surfit/doc/isreplace.hlp
new file mode 100644
index 00000000..d248e5f8
--- /dev/null
+++ b/math/surfit/doc/isreplace.hlp
@@ -0,0 +1,33 @@
+.help isreplace Apr1985 "Surfit Package"
+.ih
+NAME
+isreplace -- restore a surface fit saved by issave
+.ih
+SYNOPSIS
+isreplace (sf, fit)
+
+.nf
+pointer	sf		# surface grid descriptor
+real	fit[ARB]	# fit array
+.fi
+
+.ih
+ARGUMENTS
+
+.ls sf
+The pointer to the surface grid descriptor.
+.le
+.ls fit
+The array containing the fit parameters and coefficients.
+.le
+
+.ih
+DESCRIPTION
+
+ISREPLACE restores the fit saved by ISSAVE for use by ISEVAL or ISVECTOR.
+.ih
+NOTES
+.ih
+SEE ALSO
+issave
+.endhelp
diff --git a/math/surfit/doc/isresolve.hlp b/math/surfit/doc/isresolve.hlp
new file mode 100644
index 00000000..afed8ed8
--- /dev/null
+++ b/math/surfit/doc/isresolve.hlp
@@ -0,0 +1,45 @@
+.help isresolve Aug85 "Surfit Package"
+.ih
+NAME
+isresolve -- resolve the surface, same lines array
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+isresolve (sf, lines, ier)
+
+.nf
+pointer	sf		# surface grid descriptor
+int	lines[nlines]	# surface grid line numbers, 1 <= lines[i] <= nlines
+int	ier		# error code
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface grid descriptor structure.
+.le
+.ls nlines
+The number of surface grid lines to be fit.
+.le
+.ls ier      
+Error code returned by the fitting routines. The options are OK, SINGULAR,
+and NO_DEG_FREEDOM. If ier = SINGULAR the matrix is singular, ISSOLVE
+will compute a solution to the normal equations but one or more of the
+coefficients will be zero. If ier = NO_DEG_FREEDOM, too few data points
+exist for a reasonable solution to be computed. ISSOLVE returns without
+fitting the data.
+.le
+.ih
+DESCRIPTION
+In some applications it is necessary to refit the same surface several
+times with the same lines array. In this case it is inefficient to reaccumulate
+the matrices and calculate the Cholesky factorization for each fit.
+ISERESOLVE reaccumulates only the right side of the equation.
+.ih
+NOTES
+It is the responsibility of the calling program to check for out of bounds
+lines values.
+.ih
+SEE ALSO
+issolve, isinit, slzero, islfit, islrefit, islaccum, islresolve, isfree
+.endhelp
diff --git a/math/surfit/doc/issave.hlp b/math/surfit/doc/issave.hlp
new file mode 100644
index 00000000..6ec97cd6
--- /dev/null
+++ b/math/surfit/doc/issave.hlp
@@ -0,0 +1,55 @@
+.help issave Apr85 "Surfit Package"
+.ih
+NAME
+issave - save the surface fit
+.ih
+SYNOPSIS
+issave (sf, fit)
+
+.nf
+pointer	sf		# surface descriptor
+real	fit[ARB]	# fit array
+.fi
+.ih
+ARGUMENTS
+.ls sf   
+Pointer to the surface grid descriptor.
+.le
+.ls fit
+Array for storing the fit. The fit array must be at least SAVE_COEFF + ncoeff
+elements long.
+.le
+
+.ih
+DESCRIPTION
+ISSAVE stores the surface parameters and coefficient array
+in the 1-D reall array fit. Fit must be at least ncoeff + SAVE_COEFF
+elements long where ncoeff is defined as follows.
+
+.nf
+	SF_LEGENDRE:		nxcoeff = xorder
+				nycoeff = yorder
+				ncoeff = xorder * yorder        xterms = yes
+				ncoeff = nycoeff + nxcoeff - 1  xterms = no
+
+	SF_CHEBYSHEV:		nxcoeff = xorder
+				nycoeff = yorder
+				ncoeff = xorder * yorder        xterms = yes
+				ncoeff = nycoeff + nxcoeff - 1  xterms = no
+
+	SF_SPLINE3:		nxcoeff = xorder + 3
+				nycoeff = yorder + 3
+				ncoeff = (xorder + 3) * (yorder + 3)
+
+	SF_SPLINE1:		nxcoeff = xorder + 1
+				nycoeff = yorder + 1
+				ncoeff = (xorder + 1) * (yorder + 1)
+
+.fi
+
+.ih
+NOTES
+.ih
+SEE ALSO
+isreplace
+.endhelp
diff --git a/math/surfit/doc/issolve.hlp b/math/surfit/doc/issolve.hlp
new file mode 100644
index 00000000..312d69a3
--- /dev/null
+++ b/math/surfit/doc/issolve.hlp
@@ -0,0 +1,49 @@
+.help issolve Aug85 "Surfit Package"
+.ih
+NAME
+issolve -- solve the surface
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+issolve (sf, lines, nlines, ier)
+
+.nf
+pointer	sf		# surface grid descriptor
+int	lines[nlines]	# surface grid line numbers, 1 <= lines[i] <= nlines
+int	nlines		# number of lines
+int	ier		# error code
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface grid descriptor structure.
+.le
+.ls lines
+The line number of surface grid lines to be included in the fit.
+.le
+.ls nlines
+The number of surface grid lines to be fit.
+.le
+.ls ier      
+Error code returned by the fitting routines. The options are OK, SINGULAR,
+and NO_DEG_FREEDOM. If ier = SINGULAR the matrix is singular, ISSOLVE
+will compute a solution to the normal equations but one or more of the
+coefficients will be zero. If ier = NO_DEG_FREEDOM, too few data points
+exist for a reasonable solution to be computed. ISSOLVE returns without
+fitting the data.
+.le
+.ih
+DESCRIPTION
+ISSOLVE solves for the surface coefficients.
+An error code is
+returned by ISSOLVE if it is unable to solve the normal equations as
+formulated.
+.ih
+NOTES
+It is the responsibility of the calling program to check for out of bounds
+lines values.
+.ih
+SEE ALSO
+isinit, slzero, islfit, islrefit, islaccum, islresolve, isfree
+.endhelp
diff --git a/math/surfit/doc/isvector.hlp b/math/surfit/doc/isvector.hlp
new file mode 100644
index 00000000..fabbdcc7
--- /dev/null
+++ b/math/surfit/doc/isvector.hlp
@@ -0,0 +1,41 @@
+.help isvector Aug85 "Surfit Package"
+.ih
+NAME
+isvector -- evaluate the fitted surface at a set of points
+.ih
+SYNOPSIS
+isvector (sf, x, y, zfit, npts)
+
+.nf
+pointer	sf		# surface grid descriptor
+real	x[npts]		# x array, 1 <= x[i] <= ncols
+real	y[npts]		# y array, 1 <= y[i] <= nlines
+real	zfit[npts]	# data values
+int	npts		# number of data points
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface grid descriptor structure.
+.le
+.ls x, y
+Arrays of x and y values.
+.le
+.ls zfit
+Array of fitted values.
+.le
+.ls npts
+The number of points to be fit.
+.le
+.ih
+DESCRIPTION
+Fit the surface to an array of x and y values. ISVECTOR uses the coefficients
+stored in the surface descriptor structure.
+.ih
+NOTES
+Checking for out of bounds x and y values is the responsibility of the
+calling program.
+.ih
+SEE ALSO
+iseval
+.endhelp
diff --git a/math/surfit/doc/iszero.hlp b/math/surfit/doc/iszero.hlp
new file mode 100644
index 00000000..4c338807
--- /dev/null
+++ b/math/surfit/doc/iszero.hlp
@@ -0,0 +1,26 @@
+.help iszero Aug85 "Surfit Package"
+.ih
+NAME
+iszero -- set up for a new surface fit
+.ih
+SYNOPSIS
+iszero (sf)
+
+.nf
+pointer	sf		# surface grid descriptor
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface grid descriptor structure.
+.le
+.ih
+DESCRIPTION
+ISZERO zeros all matrices and vectors used by the fitting program. It
+should be used only if one wishes to refit the surface from scratch.
+.ih
+NOTES
+.ih
+SEE ALSO
+isinit, islzero, islaccum, islfit, islrefit, islsolve, issolve, isfree
+.endhelp
diff --git a/math/surfit/doc/surfit.hd b/math/surfit/doc/surfit.hd
new file mode 100644
index 00000000..18e7ec1a
--- /dev/null
+++ b/math/surfit/doc/surfit.hd
@@ -0,0 +1,19 @@
+# HELP Directory for SURFIT package
+
+$surfit = "math$surfit/"
+
+iscoeff		hlp = iscoeff.hlp,		src = surfit$iscoeff.x
+iseval		hlp = iseval.hlp,		src = surfit$iseval.x
+isfree		hlp = isfree.hlp,		src = surfit$isfree.x
+isinit		hlp = isinit.hlp,		src = surfit$isinit.x
+islaccum	hlp = islaccum.hlp,		src = surfit$islaccum.x
+islfit		hlp = islfit.hlp,		src = surfit$islfit.x
+islrefit	hlp = islrefit.hlp,		src = surfit$islrefit.x
+islsolve	hlp = islsolve.hlp,		src = surfit$islsolve.x
+islzero		hlp = islzero.hlp,		src = surfit$islzero.x
+isreplace	hlp = isreplace.hlp,		src = surfit$isreplace.x
+isresolve	hlp = isresolve.hlp,		src = surfit$isresolve.x
+issave		hlp = issave.hlp,		src = surfit$issave.x
+issolve		hlp = issolve.hlp,		src = surfit$issolve.x
+isvector	hlp = isvector.hlp,		src = surfit$isvector.x
+iszero		hlp = iszero.hlp,		src = surfit$iszero.x
diff --git a/math/surfit/doc/surfit.hlp b/math/surfit/doc/surfit.hlp
new file mode 100644
index 00000000..a7c347cc
--- /dev/null
+++ b/math/surfit/doc/surfit.hlp
@@ -0,0 +1,157 @@
+.help surfit Aug84 "Math Package"
+
+.ih
+NAME
+surfit -- set of routines for fitting gridded surface data
+
+
+.ih
+SYNOPSIS
+
+The surface is defined in the following way. The x and y values are
+assumed to be integers and lie in the range 0 <= x <= ncols + 1 and
+0 <= y <= nlines + 1. The package prefix is for image surface fit.
+Package routines with a prefix followed by l operate on individual
+image lines only.
+
+.nf
+        isinit	(sf, surf_type, xorder, yorder, xterms, ncols, nlines)
+	iszero	(sf)
+       islzero	(sf, lineno)
+      islaccum	(sf, cols, lineno, z, w, ncols, wtflag)
+      islsolve	(sf, lineno, ier)
+        islfit	(sf, cols, lineno, z, w, ncols, wtflag, ier)
+      islrefit	(sf, cols, lineno, z, w, ier)
+       issolve	(sf, lines, nlines, ier)
+     isresolve	(sf, lines, ier)
+    z = iseval	(sf, x, y)
+      isvector	(sf, x, y, zfit, npts)
+        issave	(sf, surface)
+     isrestore	(sf, surface)
+        isfree	(sf)
+.fi
+
+.ih
+DESCRIPTION
+
+The SURFIT package provides a set of routines for fitting surfaces to functions
+linear in their coefficients using the tensor product method and least squares
+techniques. The basic numerical
+technique employed is the solution of the normal equations by the Cholesky
+method.
+
+.ih
+NOTES
+
+The following series of steps illustrates the use of the package.
+
+.ls 4
+.ls (1)
+Include an include <math/surfit.h> statement in the calling program to make the
+SURFIT package definitions available to the user program.
+.le
+.ls (2)
+Call ISINIT to initialize the surface fitting parameters.
+.le
+.ls (3)
+Call ISLACCUM to select a weighting function and accumulate data points
+for each line into the appropriate arrays and vectors. Call ISLSOLVE
+to compute the coefficients in x for each line. The coefficients are
+stored inside SURFIT. Repeat this procedure for each line. ISLACCUM
+and SFLSOLVE can be combined by a call to SFLFIT. If the x values
+and weights remain the same from line to line ISLREFIT can be called.
+.le
+.ls (4)
+Call ISSOLVE to solve for the surface coefficients.
+.le
+.ls (5)
+Call ISEVAL or ISVECTOR to evaluate the fitted surface at the x and y
+value(s) of interest.
+.le
+.ls (6)
+Call ISFREE to release the space allocated for the fit.
+.le
+.le
+
+.ih
+EXAMPLES
+
+.nf
+Example 1: Fit a 2nd order Lengendre polynomial in x and y to an image
+and output the fitted image
+
+include <imhdr.h>
+include <math/surfit.h>
+
+ 	
+	old = immap (old_image, READ_ONLY, 0)
+	new = immap (new_image, NEW_COPY, 0)
+	
+	ncols = IM_LEN(old, 1)
+	nlines = IM_LEN(old, 2)
+	
+	# initialize surface fit
+	call isinit (sf, LEGENDRE, 3, 3, YES, ncols, nlines)
+        
+	# allocate space for lines, columns and weight arrays
+	call malloc (cols, ncols, TY_INT)
+	call malloc (lines, nlines, TY_INT)
+	call malloc (weight, ncols, TY_REAL)
+
+	# initialize lines and columns arrays
+	do i = 1, ncols
+	    Memi[cols - 1 + i] = i
+	do i = 1, nlines
+	    Memi[lines - 1 + i] = i
+
+	# fit the surface in x line by line
+	call amovkl (long(1), v, IM_MAXDIM)
+	do i = 1, nlines {
+	    if (imgnlr (old, inbuf, v) == EOF)
+		call error (0, "Error reading image.")
+	    if (i == 1) {
+		call islfit (sf, Memi[cols], i, Memr[inbuf], Memr[weight],
+				ncols, SF_WTSUNIFORM, ier)
+		if (ier != OK)
+		    ...
+	    } else
+		call islrefit (sf, Memi[cols], i, Memr[inbuf], Memr[weight],
+				ier)
+	}
+
+	# solve for surface coefficients
+	call issolve (sf, Memi[lines], nlines, ier)
+
+	# free space used in fitting arrays
+	call mfree (cols, TY_INT)
+	call mfree (lines, TY_INT)
+	call mfree (weight, TY_REAL)
+
+	# allocate space for x and y arrays
+	call malloc (x, ncols, TY_REAL)
+	call malloc (y, ncols, TY_REAL)
+
+	# intialize z array
+	do i = 1, ncols
+	    Memr[x - 1 + i] = real (i)
+
+	# create fitted image
+	call amovkl (long(10, v, IM_MAXDIM)
+	do i = 1, nlines {
+	    if (impnlr (new, outbuf, v) == EOF)
+		call error (0, "Error writing image.")
+	    call amovkr (real (i), Memr[y], ncols)
+	    call isvector (sf, Memr[x], Memr[y], Memr[outbuf], ncols)
+	}
+
+	# close files and cleanup
+	call mfree (x, TY_REAL)
+	call mfree (y, TY_REAL
+
+	call isfree (sf)
+
+	call imunmap (old)
+	call imunmap (new)
+
+.fi
+.endhelp
diff --git a/math/surfit/doc/surfit.men b/math/surfit/doc/surfit.men
new file mode 100644
index 00000000..4a3b3258
--- /dev/null
+++ b/math/surfit/doc/surfit.men
@@ -0,0 +1,15 @@
+         isinit - initialize surface descriptor
+        iscoeff - get coefficients of surface
+         iseval - evaluate the surface at a point
+         isfree - free surface descriptor
+       islaccum - accumulate surface at a given line number
+         islfit - fit the surface at a given line number
+       islrefit - refit the surface at a given line number, same columns
+       islsolve - solve the surface at a given line number
+        islzero - zero the fit for a given line number
+        issolve - solve the surface
+      isreplace - restore the surface fit
+      isresolve - resolve surface, same lines
+         issave - save the surface fit
+       isvector - fit surface at a set of (x,y) values
+         iszero - zero the surface fit
diff --git a/math/surfit/doc/surfit.spc b/math/surfit/doc/surfit.spc
new file mode 100644
index 00000000..07f2b934
--- /dev/null
+++ b/math/surfit/doc/surfit.spc
@@ -0,0 +1,500 @@
+.help surfit Aug84 "Math Package"
+.CE
+Specifications for the Surface Fitting Package
+.CE
+Lindsey Davis
+.CE
+August 1984
+.CE
+First Draft
+
+.sh
+1. Introduction
+
+The SURFIT package provides a set of routines for fitting surfaces to functions
+linear in their coefficients using least squares techniques. The basic numerical
+technique employed is the solution of the normal equations by the Cholesky
+method.
+
+.sh
+2. Requirements
+
+.ls (1)
+The package shall take as input a surface or section of
+a surface.
+The data are assumed to be on a rectangular grid and to have the following
+ranges, 0 <= x <= ncols + 1 and 0 <= y <= nlines + 1.
+The package routines assume that data values equal to INDEF have
+been removed from the data set or replaced with appropriate interpolated
+values prior to entering the package routines. It is not necessary that
+all the data be present to perform the fit.
+.le
+.ls (2)
+The package shall perform the following operations:
+.ls o
+Determine the coefficients of the fitting function by solving the normal
+equations. The surface fitting function is selected at run time from the
+following list: (1) SF_LEGENDRE, Legendre polynomials in x and y,
+(2) SF_CHEBYSHEV, Chebyshev polynomials in x and y, (3) SF_SPLINE3, bicubic
+spline with break points evenly spaced in x and y. The calling sequence
+must be invariant to the form of the fitting function.
+.le
+.ls o
+Set an error code if the numerical routines are unable to fit the specified
+function.
+.le
+.ls o
+Optionally output the values of the coefficients. The coefficients will be
+stored internal to the SURFIT package. However in some applications it is
+the coefficients which are of primary interest. A package routine shall
+exist to extract the coefficients from the curve descriptor structure.
+.le
+.ls o
+Evaluate the surface at arbitrary value(s) of x and y. The evaluating
+routines shall use the calculated coefficients and the user supplied
+x values(s).
+.le
+.ls o
+Be capable of storing the fitted surface parameters and coefficients in a
+user supplied array and restoring the saved fit for later use by the
+evaluating routines.
+.le
+.ls o
+Input surface parameters and surface coefficients derived external to
+the surfit package into the surfit format for use by the evaluate
+routines.
+.le
+.le
+.ls (3)
+The program shall perform a weighted fit to the surface using a user
+supplied weight array and weight flag. The weighting options will be
+SF_WTSUSER and SF_WTSUNIFORM. In SF_WTSUNIFORM mode the package routines
+set all the weights to 1. In SF_WTSUSER mode the package routines apply
+user supplied weights to the individual data points.
+.le
+.ls (4)
+The input data set, output coefficients, error and fitted z arrays are
+single precision real quantities. All package arithmetic shall be performed
+in single precision.
+.le
+
+.sh
+3. Specifications
+
+.sh
+3.1. List of Routines
+
+
+The surface is defined in the following way. The x and y values are
+assumed to be integers and lie in the range 0 <= x <= ncols + 1 and
+0 <= y <= nlines + 1. The package prefix is is for image surface fit.
+Package routines with a prefix followed by l operate on individual
+image lines only.
+
+.nf
+	z = f (col, line)
+.fi
+
+.nf
+        isinit	(sf, surf_type, xorder, yorder, xterms, ncols, nlines)
+	iszero	(sf)
+       islzero	(sf, lineno)
+      islaccum	(sf, cols, lineno, z, w, ncols, wtflag)
+      islsolve	(sf, lineno, ier)
+        islfit	(sf, cols, lineno, z, w, ncols, wtflag, ier)
+      islrefit	(sf, cols, lineno, z, w, ier)
+       issolve	(sf, lines, nlines, ier)
+     isresolve	(sf, lines, ier)
+    z = iseval	(sf, x, y)
+      isvector	(sf, x, y, zfit, npts)
+        issave	(sf, surface)
+     isrestore	(sf, surface)
+        isfree	(sf)
+.fi
+
+.sh
+3.2. Algorithms
+
+.sh
+3.2.1. Polynomial Basis Functions
+
+The approximating function is assumed to be of the following form,
+
+.nf
+    f(x,y) = sum (a[i,j] * F(i,x) * F(j,y))    i = 1...nxcoeff  j = 1...nycoeff
+.fi
+
+where the F(i,x) are the polynomial basis functions containing terms of order
+x**(i-1) and the a(i,j) are the coefficients. In order to avoid a very
+ill-conditioned linear system for moderate or large i, the Legendre and
+Chebyshev polynomials were chosen for the basis functions. The Chebyshev
+and Legendre polynomials are orthogonal over -1. <= x,y <= 1. The data x and
+y values are normalized to this regime using the number of lines and columns
+supplied by the user. For each data point the basis
+functions are calculated using the following recursion relations. The
+cross terms are optional.
+
+Legendre series:
+.nf
+
+    F(1,x) = 1.
+    F(2,x) = x
+    F(i,x) = [(2*i-1)*x*F(i-1,x)-(i-2)*F(i-2,x)]/(i-1)
+    F(1,y) = 1.
+    F(2,y) = y
+    F(j,y) = [(2*j-1)*y*F(j-1,y)-(j-2)*F(j-2,y)]/(j-1)
+
+.fi
+Chebyshev series:
+.nf
+
+    F(1,x) = 1.
+    F(2,x) = x
+    F(i,x) = 2.*x*F(i-1,x)-F(i-2,x)
+    F(1,y) = 1.
+    F(2,y) = y
+    F(j,y) = 2.*y*F(j-1,y)-F(j-2,y)
+.fi
+
+
+.sh
+3.2.2. Bicubic Cardinal B-spline
+
+The approximating function is of the form
+
+.nf
+    f(x,y) =  sum (x(i,j) * F(i,x) * F(j,y))   i=1...nxcoeff j=1...nycoeff
+.fi
+
+where the basis functions, F(i,x), are the cubic cardinal B-splines
+(Prenter 1975). The user supplies the number of columns and lines and the
+number of polynomial pieces in x and y, nxpieces and nypieces, to be fit
+to the data set. The number of bicubic spline coefficients, ncoeff, will be
+
+.nf
+    nxcoeff = (nxpieces + 3)
+    nycoeff = (nypieces + 3)
+    ncoeff = nxcoeff * nycoeff
+.fi
+
+The cardinal B-spline is stored in a lookup table. For each x and y the
+appropriate break point is selected and the four non-zero B-splines are
+calculated by nearest neighbour interpolation in the lookup table.
+
+.sh
+3.2.3. Method of Solution
+
+.sh
+3.2.3.1. Full Surface Fit
+
+The normal equations are accumulated in the following form.
+
+.nf
+    c * x = b
+    c[i,j] = (B(i,x,y), B(j,x,y))
+    b[j] = (B(j,x,y), S(x,y))
+.fi
+
+B(i,x,y) is the ith basis function at x and y, S(x,y) is the surface to
+be approximated and the inner product of two functions G and H is
+given by 
+
+.nf
+    (G,H) = sum (w[i] * G(x[i],y[i]) * H(x[i],y[i]))  i = 1...npts
+.fi
+
+Since the matrix c is symmetric and positive semi-definite it may
+be solved by the Cholesky method.
+The coefficient matrix c  can be written as
+
+.nf
+    c = l * d * l-transpose
+.fi
+
+where l is a unit lower triangular matrix and d is the diagonal of c
+(de Boor 1978). Near zero pivots are handled in the following way.
+At the nth elimination step the current value of the nth diagonal
+element is compared with the original nth diagonal element. If the
+diagonal element has been reduced by one computer word length, the
+entire nth row is declared linearly dependent on the previous n-1
+rows and x(n) = 0.
+
+The triangular system
+
+.nf
+    l * w = b
+.fi
+
+is solved for w (forward substitution), the vector d  ** (-1) * w is
+computed and the triangular system
+
+.nf
+    l-transpose * x = d ** (-1) * w
+.fi
+
+solved for the coefficients, x (backward substitution).
+
+
+.sh
+3.2.3.2. Line by Line Fit
+
+Each line of the surface can be represented by the following equation.
+
+.nf
+    S(x,yo) = a[i](yo) * B[i](x)            i = 1...nxcoeff for each yo
+.fi
+
+The normal equations for each image line are formed
+
+.nf
+    c * a = b
+    c[i,j] = (B(i,x), B(j,x))
+    b[j] = (B(j,x), S(x,yo))
+.fi
+
+and solved for a as described in the previous section.
+After fitting the entire image matrix a has nxcoeff columns and
+nlines rows.
+
+For each column i in a the normal equations
+are formed and solved for the c[i,j]
+
+.nf
+    c * x  = b
+    c[j,k] = (B(j,y), B(k,x))
+    b[i,j] = (a[i](y), B(j,y))
+.fi
+
+.sh
+3.2.4. Number of Operations
+
+It is worth while to calculate the number of operations reqired to compute
+the surface fit by the full surface method versus the line by line method.
+The number of operations required to accumulate the normal equations
+and solve the set is the following
+
+.nf
+    nop = npts * order ** 2 / 2 + order ** 3 /6
+.fi
+
+where the first term is the number of operations required to accumulate the
+matrix and the second is the number required to solve the system.
+
+.nf
+Example 1: full surface, 3 by 3 polynomial no cross terms, 512 by 512 image
+
+    norder = 5
+    npts = 512 * 512
+    nops(accum) = 3,276,800
+    nops(solve) = 21
+
+Example 2: full surface, 30 by 30 spline, 512 by 512 image
+
+    norder = 900
+    npts = 512 * 512
+    nops(accum) = 1.062 * 10 ** 11
+    nops(solve) = 1.216 * 10 ** 8
+
+Example 3: line by line method, 3 by 3 polynomial, 512 by 512 image
+
+    step 1 solve for a coefficients
+    norder = 3
+    npts = 512
+    nlines = 512
+    nops(accum) = 1,179,648
+    nops(solve) = 2304
+
+    step 2 solve for c coefficients
+    norder = 3
+    npts = 512
+    nlines = 3
+    nops(accum) = 6912
+    nops(solve) = 14
+
+Example 4: line by line method, 30 by 30 spline, 512 by 512 image
+
+    step 1 solve coefficients
+    norder = 30
+    npts = 512
+    nlines = 512
+    nops(accum) = 117,964,800
+    nops(solve) = 2,304,000
+
+    step 2 solve for c coefficients
+    norder = 30
+    npts = 512
+    nlines = 30
+    nops(accum) = 6,912,000
+    nops(solve) = 135,000
+.fi
+
+The figures for the line by line method are worst case numbers. If the
+x values remain the same from line to line then the coefficient matrix
+only has to be accumulated and inverted twice.
+For the bicubic spline function the number of operations is significantly
+reduced by taking advantage of the banded nature of the matrices.
+
+.sh
+4. Usage
+
+.sh
+4.1. User Notes
+
+The following series of steps illustrates the use of the package.
+
+.ls 4
+.ls (1)
+Include an include <surfit.h> statement in the calling program to make the
+SURFIT package definitions available to the user program.
+.le
+.ls (2)
+Call SFINIT to initialize the surface fitting parameters.
+.le
+.ls (3)
+Call SFLACCUM to select a weighting function and accumulate data points
+for each line into the appropriate arrays and vectors. Call SFLSOLVE
+to compute the coefficients in x for each line. The coefficents are
+stored inside SURFIT. Repeat this procedure for each line. SFLACCUM
+and SFLSOLVE can be combined by a call to SFLFIT. If the x values
+and weights remain the same from line to line SFLREFIT can be called.
+.le
+.ls (4)
+Call SFSOLVE to solve for the surface coefficients.
+.le
+.ls (5)
+Call SFEVAL or SFVECTOR to evaluate the fitted surface at the x and y
+value(s) of interest.
+.le
+.ls (6)
+Call SFFREE to release the space allocated for the fit.
+.le
+.le
+
+.sh
+4.2. Examples
+
+.nf
+Example 1: Fit a 2nd order Lengendre polynomial in x and y to an image
+and output the fitted image
+
+include <imhdr.h>
+include <math/surfit.h>
+
+ 	
+	old = immap (old_image, READ_ONLY, 0)
+	new = immap (new_image, NEW_COPY, 0)
+	
+	ncols = IM_LEN(old, 1)
+	nlines = IM_LEN(old, 2)
+	
+	# initialize surface fit
+	call isinit (sf, LEGENDRE, 3, 3, YES, ncols, nlines)
+        
+	# allocate space for lines, columns and weight arrays
+	call malloc (cols, ncols, TY_INT)
+	call malloc (lines, nlines, TY_INT)
+	call malloc (weight, ncols, TY_REAL)
+
+	# initialize lines and columns arrays
+	do i = 1, ncols
+	    Memi[cols - 1 + i] = i
+	do i = 1, nlines
+	    Memi[lines - 1 + i] = i
+
+	# fit the surface in x line by line
+	call amovkl (long(1), v, IM_MAXDIM)
+	do i = 1, nlines {
+	    if (imgnlr (old, inbuf, v) == EOF)
+		call error (0, "Error reading image.")
+	    if (i == 1) {
+		call islfit (sf, Memi[cols], i, Memr[inbuf], Memr[weight],
+				ncols, SF_WTSUNIFORM, ier)
+		if (ier != OK)
+		    ...
+	    } else
+		call sflrefit (sf, Memi[cols], i, Memr[inbuf], Memr[weight],
+				ier)
+	}
+
+	# solve for surface coefficients
+	call issolve (sf, Memi[lines], nlines, ier)
+
+	# free space used in fitting arrays
+	call mfree (cols, TY_INT)
+	call mfree (lines, TY_INT)
+	call mfree (weight, TY_REAL)
+
+	# allocate space for x and y arrays
+	call malloc (x, ncols, TY_REAL)
+	call malloc (y, ncols, TY_REAL)
+
+	# intialize z array
+	do i = 1, ncols
+	    Memr[x - 1 + i] = real (i)
+
+	# create fitted image
+	call amovkl (long(10, v, IM_MAXDIM)
+	do i = 1, nlines {
+	    if (impnlr (new, outbuf, v) == EOF)
+		call error (0, "Error writing image.")
+	    call amovkr (real (i), Memr[y], ncols)
+	    call isvector (sf, Memr[x], Memr[y], Memr[outbuf], ncols)
+	}
+
+	# close files and cleanup
+	call mfree (x, TY_REAL)
+	call mfree (y, TY_REAL
+
+	call isfree (sf)
+
+	call imunmap (old)
+	call imunmap (new)
+
+.fi
+.sh
+5. Detailed Design
+
+.sh
+5.1. Surface Descriptor Structure
+
+To be written when specifications are finalised.
+
+.sh
+5.2. Storage Requirements
+
+The minimum storage requirements for the full surface fit method
+assuming that points are to be rejected from the fit without revaluating
+the matrix are the following.
+
+.nf
+    (nxcoeff*nycoeff) ** 2    -- coefficient array plus triangularized matrix
+    nxcoeff*nycoeff           -- right side
+    nxcoeff*nycoeff           -- solution vector
+.fi
+
+For a 30 by 30 spline roughly 811,800 storage units are required.
+For a 3 by 3 polynomial the requirements are 675 storage units. The
+requirements are roughly half when rejection is disabled.
+
+The minimum storage requirements for the line by line method are
+
+.nf
+    nx*nx*2   -- The x matrix and its factorization
+    nx*nlines -- The x coefficient matrix
+    ny*ny*2 -- The y matrix and its factorization
+    nx*ny -- The solution
+.fi
+
+
+.sh
+6. References
+
+.ls (1)
+Carl de Boor, "A Practical Guide to Splines", 1978, Springer-Verlag New York
+Inc.
+.le
+.ls (2)
+P.M. Prenter, "Splines and Variational Methods", 1975, John Wiley and Sons
+Inc.
+.le
+.endhelp