Repatch (from linux) of OSX IRAF

1 files changed, 500 insertions, 0 deletions

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