diff options
Diffstat (limited to 'math/surfit/doc')
| -rw-r--r-- | math/surfit/doc/iscoeff.hlp | 63 | ||||
| -rw-r--r-- | math/surfit/doc/iseval.hlp | 34 | ||||
| -rw-r--r-- | math/surfit/doc/isfree.hlp | 27 | ||||
| -rw-r--r-- | math/surfit/doc/isinit.hlp | 61 | ||||
| -rw-r--r-- | math/surfit/doc/islaccum.hlp | 60 | ||||
| -rw-r--r-- | math/surfit/doc/islfit.hlp | 67 | ||||
| -rw-r--r-- | math/surfit/doc/islrefit.hlp | 51 | ||||
| -rw-r--r-- | math/surfit/doc/islsolve.hlp | 45 | ||||
| -rw-r--r-- | math/surfit/doc/islzero.hlp | 32 | ||||
| -rw-r--r-- | math/surfit/doc/isreplace.hlp | 33 | ||||
| -rw-r--r-- | math/surfit/doc/isresolve.hlp | 45 | ||||
| -rw-r--r-- | math/surfit/doc/issave.hlp | 55 | ||||
| -rw-r--r-- | math/surfit/doc/issolve.hlp | 49 | ||||
| -rw-r--r-- | math/surfit/doc/isvector.hlp | 41 | ||||
| -rw-r--r-- | math/surfit/doc/iszero.hlp | 26 | ||||
| -rw-r--r-- | math/surfit/doc/surfit.hd | 19 | ||||
| -rw-r--r-- | math/surfit/doc/surfit.hlp | 157 | ||||
| -rw-r--r-- | math/surfit/doc/surfit.men | 15 | ||||
| -rw-r--r-- | math/surfit/doc/surfit.spc | 500 |
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 |