AlgFit

Files
file	AlgBSpline.c
	Provides functions for fitting and evaluating B-Splines. This software is based on netlib Dierckx Fortran subroutines. In most cases the original comments have been preserved with little change and are inculded in the Doxygen documentation verbatim.
file	AlgLinearFit.c
	Provides functions for fitting linear models to data, ie linear regression.
file	AlgPolyLSQ.c
	Provides functions for fitting a polynomial using least squares.

Functions
void	AlgBSplineBspl (double *t, int k, double x, int l, double *h)
	Evaluates the (k+1) non-zero B-Splines of degree k at t[l] using the de Boor Cox recurrence. More...
AlgError	AlgBSplineNDFit (int iopt, int ipar, int idim, int m, double *u, int mx, double *x, double *w, double ub, double ue, int k, double s, int nest, int *n, double *t, int *nc, double *c, double *fp, double *wrk, int *iwrk)
	Computes a smooth approximating b-spline in multi-dimensional space. More...
AlgError	AlgBSplinePerFit (int iopt, int m, double *x, double *y, double *w, int k, double s, int nest, int *n, double *t, double *c, double *fp, double *wrk, int *iwrk)
	Computes a smooth periodic B-spline approximation through the given data. More...
AlgError	AlgBSplineFit (int iopt, int m, double *x, double *y, double *w, double xb, double xe, int k, double s, int nest, int *n, double *t, double *c, double *fp, double *wrk, int *iwrk)
	Computes a smooth B-spline approximation through the given data. More...
AlgError	AlgBSplineEval (double *t, int n, double *c, int k, double *x, double *y, int m)
	Evaluates a b-spline at a number of points. More...
AlgError	AlgBSplineDer (double *t, int n, double *c, int k, int nu, double *x, double *y, int m, double *wrk)
	Evaluates the derivatives of order nu of the B-spline of degree k. More...
AlgError	AlgLinearFit1D (int datSz, double *datXA, double *datYA, double *dstA, double *dstB, double *dstSigA, double *dstSigB, double *dstQ)
	Computes the least squares best fit straight line (y = a + bx) through the given data, ie linear regression. This function is based on the function fit(): Press W. H., Teukolsky S. A., Vetterling W. T. and Flannery B. P, Numerical Recipies in C, 1992, CUP. More...
AlgError	AlgLinearFitIdx1D (double *datXA, double *datYA, int *idxXA, int *idxYA, int idxASz, double *dstA, double *dstB, double *dstSigA, double *dstSigB, double *dstQ)
	Computes the least squares best fit straight line (y = a + bx) through the given data, ie linear regression. More...
AlgError	AlgPolynomialLSq (double *xVec, double *yVec, int vecSz, int polyDeg, double *cVec)
	Attempts to fit a polynomial to the given data using a least squares approach. More...

Detailed Description

Function Documentation

AlgBSplineBspl()

void AlgBSplineBspl	(	double *	t,
		int	k,
		double	x,
		int	l,
		double *	h
	)

Evaluates the (k+1) non-zero B-Splines of degree k at t[l] using the de Boor Cox recurrence.

This function has been derived from the netlib Dierckx function fpbspl(). The original fortran coments are:

 Subroutine fpbspl evaluates the (k+1) non-zero b-splines of
 degree k at t(l) <= x < t(l+1) using the stable recurrence
 relation of de Boor and Cox.

Parameters

t	Array of knot positions.
k	Degree of the B-Spline.
x	Point at which the spline is to be evaluated.
l	Knot index st k at t[l] <= x < t[l+1].
h	Array with length >= k in which to evaluate the spline.

Referenced by AlgBSplineDer(), and AlgBSplineEval().

AlgBSplineNDFit()

AlgError AlgBSplineNDFit	(	int	iopt,
		int	ipar,
		int	idim,
		int	m,
		double *	u,
		int	mx,
		double *	x,
		double *	w,
		double	ub,
		double	ue,
		int	k,
		double	s,
		int	nest,
		int *	n,
		double *	t,
		int *	nc,
		double *	c,
		double *	fp,
		double *	wrk,
		int *	iwrk
	)

Computes a smooth approximating b-spline in multi-dimensional space.

Returns

Alg error code.

This function has been derived from the netlib Dierckx function curfit(). The original fortran coments are:

 Given the ordered set of m points x(i) in the idim-dimensional space
 and given also a corresponding set of strictly increasing values u(i)
 and the set of positive numbers w(i),i=1,2,...,m, subroutine parcur
 determines a smooth approximating spline curve s(u), i.e.
     x1 = s1(u)
     x2 = s2(u)       ub <= u <= ue
     .........
     xidim = sidim(u)
 with sj(u),j=1,2,...,idim spline functions of degree k with common
 knots t(j),j=1,2,...,n.
 if ipar=1 the values ub,ue and u(i),i=1,2,...,m must be supplied by
 the user. if ipar=0 these values are chosen automatically by parcur
 as  v(1) = 0
     v(i) = v(i-1) + dist(x(i),x(i-1)) ,i=2,3,...,m
     u(i) = v(i)/v(m) ,i=1,2,...,m
     ub = u(1) = 0, ue = u(m) = 1.
 if iopt=-1 parcur calculates the weighted least-squares spline curve
 according to a given set of knots.
 if iopt>=0 the number of knots of the splines sj(u) and the position
 t(j),j=1,2,...,n is chosen automatically by the routine. the smooth-
 ness of s(u) is then achieved by minimalizing the discontinuity
 jumps of the k-th derivative of s(u) at the knots t(j),j=k+2,k+3,...,
 n-k-1. the amount of smoothness is determined by the condition that
 f(p)=sum((w(i)*dist(x(i),s(u(i))))**2) be <= s, with s a given non-
 negative constant, called the smoothing factor.
 the fit s(u) is given in the b-spline representation and can be
 evaluated by means of subroutine curev.

 calling sequence:
    call parcur(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,nc,c,
   * fp,wrk,lwrk,iwrk,ier)

 parameters:
  iopt  : integer flag. on entry iopt must specify whether a weighted
          least-squares spline curve (iopt=-1) or a smoothing spline
          curve (iopt=0 or 1) must be determined.if iopt=0 the routine
          will start with an initial set of knots t(i)=ub,t(i+k+1)=ue,
          i=1,2,...,k+1. if iopt=1 the routine will continue with the
          knots found at the last call of the routine.
          attention: a call with iopt=1 must always be immediately
          preceded by another call with iopt=1 or iopt=0.
          unchanged on exit.
  ipar  : integer flag. on entry ipar must specify whether (ipar=1)
          the user will supply the parameter values u(i),ub and ue
          or whether (ipar=0) these values are to be calculated by
          parcur. unchanged on exit.
  idim  : integer. on entry idim must specify the dimension of the
          curve. 0 < idim < 11.
          unchanged on exit.
  m     : integer. on entry m must specify the number of data points.
          m > k. unchanged on exit.
  u     : real array of dimension at least (m). in case ipar=1,before
          entry, u(i) must be set to the i-th value of the parameter
          variable u for i=1,2,...,m. these values must then be
          supplied in strictly ascending order and will be unchanged
          on exit. in case ipar=0, on exit,array u will contain the
          values u(i) as determined by parcur.
  mx    : integer. on entry mx must specify the actual dimension of
          the array x as declared in the calling (sub)program. mx must
          not be too small (see x). unchanged on exit.
  x     : real array of dimension at least idim*m.
          before entry, x(idim*(i-1)+j) must contain the j-th coord-
          inate of the i-th data point for i=1,2,...,m and j=1,2,...,
          idim. unchanged on exit.
  w     : real array of dimension at least (m). before entry, w(i)
          must be set to the i-th value in the set of weights. the
          w(i) must be strictly positive. unchanged on exit.
          see also further comments.
  ub,ue : real values. on entry (in case ipar=1) ub and ue must
          contain the lower and upper bound for the parameter u.
          ub <=u(1), ue>= u(m). if ipar = 0 these values will
          automatically be set to 0 and 1 by parcur.
  k     : integer. on entry k must specify the degree of the splines.
          1<=k<=5. it is recommended to use cubic splines (k=3).
          the user is strongly dissuaded from choosing k even,together
          with a small s-value. unchanged on exit.
  s     : real.on entry (in case iopt>=0) s must specify the smoothing
          factor. s >=0. unchanged on exit.
          for advice on the choice of s see further comments.
  nest  : integer. on entry nest must contain an over-estimate of the
          total number of knots of the splines returned, to indicate
          the storage space available to the routine. nest >=2*k+2.
          in most practical situation nest=m/2 will be sufficient.
          always large enough is nest=m+k+1, the number of knots
          needed for interpolation (s=0). unchanged on exit.
  n     : integer.
          unless ier = 10 (in case iopt >=0), n will contain the
          total number of knots of the smoothing spline curve returned
          if the computation mode iopt=1 is used this value of n
          should be left unchanged between subsequent calls.
          in case iopt=-1, the value of n must be specified on entry.
  t     : real array of dimension at least (nest).
          on succesful exit, this array will contain the knots of the
          spline curve,i.e. the position of the interior knots t(k+2),
          t(k+3),..,t(n-k-1) as well as the position of the additional
          t(1)=t(2)=...=t(k+1)=ub and t(n-k)=...=t(n)=ue needed for
          the b-spline representation.
          if the computation mode iopt=1 is used, the values of t(1),
          t(2),...,t(n) should be left unchanged between subsequent
          calls. if the computation mode iopt=-1 is used, the values
          t(k+2),...,t(n-k-1) must be supplied by the user, before
          entry. see also the restrictions (ier=10).
  nc    : integer. on entry nc must specify the actual dimension of
          the array c as declared in the calling (sub)program. nc
          must not be too small (see c). unchanged on exit.
  c     : real array of dimension at least (nest*idim).
          on succesful exit, this array will contain the coefficients
          in the b-spline representation of the spline curve s(u),i.e.
          the b-spline coefficients of the spline sj(u) will be given
          in c(n*(j-1)+i),i=1,2,...,n-k-1 for j=1,2,...,idim.
  fp    : real. unless ier = 10, fp contains the weighted sum of
          squared residuals of the spline curve returned.
  wrk   : real array of dimension at least m*(k+1)+nest*(6+idim+3*k).
          used as working space. if the computation mode iopt=1 is
          used, the values wrk(1),...,wrk(n) should be left unchanged
          between subsequent calls.
  lwrk  : integer. on entry,lwrk must specify the actual dimension of
          the array wrk as declared in the calling (sub)program. lwrk
          must not be too small (see wrk). unchanged on exit.
  iwrk  : integer array of dimension at least (nest).
          used as working space. if the computation mode iopt=1 is
          used,the values iwrk(1),...,iwrk(n) should be left unchanged
          between subsequent calls.
  ier   : integer. unless the routine detects an error, ier contains a
          non-positive value on exit, i.e.
   ier=0  : normal return. the curve returned has a residual sum of
            squares fp such that abs(fp-s)/s <= tol with tol a relat-
            ive tolerance set to 0.001 by the program.
   ier=-1 : normal return. the curve returned is an interpolating
            spline curve (fp=0).
   ier=-2 : normal return. the curve returned is the weighted least-
            squares polynomial curve of degree k.in this extreme case
            fp gives the upper bound fp0 for the smoothing factor s.
   ier=1  : error. the required storage space exceeds the available
            storage space, as specified by the parameter nest.
            probably causes : nest too small. if nest is already
            large (say nest > m/2), it may also indicate that s is
            too small
            the approximation returned is the least-squares spline
            curve according to the knots t(1),t(2),...,t(n). (n=nest)
            the parameter fp gives the corresponding weighted sum of
            squared residuals (fp>s).
   ier=2  : error. a theoretically impossible result was found during
            the iteration proces for finding a smoothing spline curve
            with fp = s. probably causes : s too small.
            there is an approximation returned but the corresponding
            weighted sum of squared residuals does not satisfy the
            condition abs(fp-s)/s < tol.
   ier=3  : error. the maximal number of iterations maxit (set to 20
            by the program) allowed for finding a smoothing curve
            with fp=s has been reached. probably causes : s too small
            there is an approximation returned but the corresponding
            weighted sum of squared residuals does not satisfy the
            condition abs(fp-s)/s < tol.
   ier=10 : error. on entry, the input data are controlled on validity
            the following restrictions must be satisfied.
            -1<=iopt<=1, 1<=k<=5, m>k, nest>2*k+2, w(i)>0,i=1,2,...,m
            0<=ipar<=1, 0<idim<=10, lwrk>=(k+1)*m+nest*(6+idim+3*k),
            nc>=nest*idim
            if ipar=0: sum j=1,idim (x(idim*i+j)-x(idim*(i-1)+j))**2>0
                       i=1,2,...,m-1.
            if ipar=1: ub<=u(1)<u(2)<...<u(m)<=ue
            if iopt=-1: 2*k+2<=n<=min(nest,m+k+1)
                        ub<t(k+2)<t(k+3)<...<t(n-k-1)<ue
                           (ub=0 and ue=1 in case ipar=0)
                      the schoenberg-whitney conditions, i.e. there
                      must be a subset of data points uu(j) such that
                        t(j) < uu(j) < t(j+k+1), j=1,2,...,n-k-1
            if iopt>=0: s>=0
                        if s=0 : nest >= m+k+1
            if one of these conditions is found to be violated,control
            is immediately repassed to the calling program. in that
            case there is no approximation returned.

 further comments:
  by means of the parameter s, the user can control the tradeoff
  between closeness of fit and smoothness of fit of the approximation.
  if s is too large, the curve will be too smooth and signal will be
  lost ; if s is too small the curve will pick up too much noise. in
  the extreme cases the program will return an interpolating curve if
  s=0 and the least-squares polynomial curve of degree k if s is
  very large. between these extremes, a properly chosen s will result
  in a good compromise between closeness of fit and smoothness of fit.
  to decide whether an approximation, corresponding to a certain s is
  satisfactory the user is highly recommended to inspect the fits
  graphically.
  recommended values for s depend on the weights w(i). if these are
  taken as 1/d(i) with d(i) an estimate of the standard deviation of
  x(i), a good s-value should be found in the range (m-sqrt(2*m),m+
  sqrt(2*m)). if nothing is known about the statistical error in x(i)
  each w(i) can be set equal to one and s determined by trial and
  error, taking account of the comments above. the best is then to
  start with a very large value of s ( to determine the least-squares
  polynomial curve and the upper bound fp0 for s) and then to
  progressively decrease the value of s ( say by a factor 10 in the
  beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
  approximating curve shows more detail) to obtain closer fits.
  to economize the search for a good s-value the program provides with
  different modes of computation. at the first call of the routine, or
  whenever he wants to restart with the initial set of knots the user
  must set iopt=0.
  if iopt=1 the program will continue with the set of knots found at
  the last call of the routine. this will save a lot of computation
  time if parcur is called repeatedly for different values of s.
  the number of knots of the spline returned and their location will
  depend on the value of s and on the complexity of the shape of the
  curve underlying the data. but, if the computation mode iopt=1 is
  used, the knots returned may also depend on the s-values at previous
  calls (if these were smaller). therefore, if after a number of
  trials with different s-values and iopt=1, the user can finally
  accept a fit as satisfactory, it may be worthwhile for him to call
  parcur once more with the selected value for s but now with iopt=0.
  indeed, parcur may then return an approximation of the same quality
  of fit but with fewer knots and therefore better if data reduction
  is also an important objective for the user.

  the form of the approximating curve can strongly be affected by
  the choice of the parameter values u(i). if there is no physical
  reason for choosing a particular parameter u, often good results
  will be obtained with the choice of parcur (in case ipar=0), i.e.
       v(1)=0, v(i)=v(i-1)+q(i), i=2,...,m, u(i)=v(i)/v(m), i=1,..,m
  where
       q(i)= sqrt(sum j=1,idim (xj(i)-xj(i-1))**2 )
  other possibilities for q(i) are
       q(i)= sum j=1,idim (xj(i)-xj(i-1))**2
       q(i)= sum j=1,idim abs(xj(i)-xj(i-1))
       q(i)= max j=1,idim abs(xj(i)-xj(i-1))
       q(i)= 1

 other subroutines required:
   fpback,fpbspl,fpchec,fppara,fpdisc,fpgivs,fpknot,fprati,fprota

 references:
  dierckx p. : algorithms for smoothing data with periodic and
               parametric splines, computer graphics and image
               processing 20 (1982) 171-184.
  dierckx p. : algorithms for smoothing data with periodic and param-
               etric splines, report tw55, dept. computer science,
               k.u.leuven, 1981.
  dierckx p. : curve and surface fitting with splines, monographs on
               numerical analysis, oxford university press, 1993.

 author:
   p.dierckx
   dept. computer science, k.u. leuven
   celestijnenlaan 200a, b-3001 heverlee, belgium.
   e-mail : Paul.Dierckx@cs.kuleuven.ac.be

 creation date : may 1979
 latest update : march 1987

• Parameters

  iopt	Option for how spline is computed.
  ipar	If parameter values u[], ub and ue are supplied the must be 1, else 0.
  idim	Curve dimension (0 < idim < 11).
  m	Number of data points.
  u	Array of parameter variables.
  mx	Dimension of the array x.
  x	Independent data.
  w	Weights for data points.
  ub	Lower bound of parameter values if ipar == 1.
  ue	Upper bound of parameter values if ipar == 1.
  k	Degree of spline.
  s	Smoothing parameter.
  nest	Over estimate of the number of knots to be returned (eg)
  n	Destination pointer for the number of knots returned.
  t	Array for returned computed knots.
  nc
  c	Array for returned spline coefficients.
  fp	Weighted sum of squared residuals of the spline approximation returned.
  wrk	Workspace with minimum size m * (k + 1) + nest * (idim + 6 + k * 3).
  iwrk	Workspace with minimum size nest * idim.

References ALG_ERR_FUNC, and ALG_ERR_NONE.

Referenced by WlzBSplineFromVertices().

AlgBSplinePerFit()

AlgError AlgBSplinePerFit	(	int	iopt,
		int	m,
		double *	x,
		double *	y,
		double *	w,
		int	k,
		double	s,
		int	nest,
		int *	n,
		double *	t,
		double *	c,
		double *	fp,
		double *	wrk,
		int *	iwrk
	)

Computes a smooth periodic B-spline approximation through the given data.

Returns

Alg error code.

This function has been derived from the netlib Dierckx function percur(). The original fortran coments are:

Given the set of data points (x(i),y(i)) and the set of positive 
numbers w(i),i=1,2,...,m-1, subroutine percur determines a smooth 
periodic spline approximation of degree k with period per=x(m)-x(1). 
if iopt=-1 percur calculates the weighted least-squares periodic 
spline according to a given set of knots. 
if iopt>=0 the number of knots of the spline s(x) and the position 
t(j),j=1,2,...,n is chosen automatically by the routine. the smooth- 
ness of s(x) is then achieved by minimalizing the discontinuity 
jumps of the k-th derivative of s(x) at the knots t(j),j=k+2,k+3,..., 
n-k-1. the amount of smoothness is determined by the condition that 
f(p)=sum((w(i)*(y(i)-s(x(i))))**2) be <= s, with s a given non- 
negative constant, called the smoothing factor. 
the fit s(x) is given in the b-spline representation (b-spline coef- 
ficients c(j),j=1,2,...,n-k-1) and can be evaluated by means of 
subroutine splev. 

calling sequence: 
   call percur(iopt,m,x,y,w,k,s,nest,n,t,c,fp,wrk, lwrk,iwrk,ier) 

parameters: 
 iopt  : integer flag. on entry iopt must specify whether a weighted 
         least-squares spline (iopt=-1) or a smoothing spline (iopt= 
         0 or 1) must be determined. if iopt=0 the routine will start 
         with an initial set of knots t(i)=x(1)+(x(m)-x(1))*(i-k-1), 
         i=1,2,...,2*k+2. if iopt=1 the routine will continue with 
         the knots found at the last call of the routine. 
         attention: a call with iopt=1 must always be immediately 
         preceded by another call with iopt=1 or iopt=0. 
         unchanged on exit. 
 m     : integer. on entry m must specify the number of data points. 
         m > 1. unchanged on exit. 
 x     : real array of dimension at least (m). before entry, x(i) 
         must be set to the i-th value of the independent variable x, 
         for i=1,2,...,m. these values must be supplied in strictly 
         ascending order. x(m) only indicates the length of the 
         period of the spline, i.e per=x(m)-x(1). 
         unchanged on exit. 
 y     : real array of dimension at least (m). before entry, y(i) 
         must be set to the i-th value of the dependent variable y, 
         for i=1,2,...,m-1. the element y(m) is not used. 
         unchanged on exit. 
 w     : real array of dimension at least (m). before entry, w(i) 
         must be set to the i-th value in the set of weights. the 
         w(i) must be strictly positive. w(m) is not used. 
         see also further comments. unchanged on exit. 
 k     : integer. on entry k must specify the degree of the spline. 
         1<=k<=5. it is recommended to use cubic splines (k=3). 
         the user is strongly dissuaded from choosing k even,together 
         with a small s-value. unchanged on exit. 
 s     : real.on entry (in case iopt>=0) s must specify the smoothing 
         factor. s >=0. unchanged on exit. 
         for advice on the choice of s see further comments. 
 nest  : integer. on entry nest must contain an over-estimate of the 
         total number of knots of the spline returned, to indicate 
         the storage space available to the routine. nest >=2*k+2. 
         in most practical situation nest=m/2 will be sufficient. 
         always large enough is nest=m+2*k,the number of knots needed 
         for interpolation (s=0). unchanged on exit. 
 n     : integer. 
         unless ier = 10 (in case iopt >=0), n will contain the 
         total number of knots of the spline approximation returned. 
         if the computation mode iopt=1 is used this value of n 
         should be left unchanged between subsequent calls. 
         in case iopt=-1, the value of n must be specified on entry. 
 t     : real array of dimension at least (nest). 
         on succesful exit, this array will contain the knots of the 
         spline,i.e. the position of the interior knots t(k+2),t(k+3) 
         ...,t(n-k-1) as well as the position of the additional knots 
         t(1),t(2),...,t(k+1)=x(1) and t(n-k)=x(m),..,t(n) needed for 
         the b-spline representation. 
         if the computation mode iopt=1 is used, the values of t(1), 
         t(2),...,t(n) should be left unchanged between subsequent 
         calls. if the computation mode iopt=-1 is used, the values 
         t(k+2),...,t(n-k-1) must be supplied by the user, before 
         entry. see also the restrictions (ier=10). 
 c     : real array of dimension at least (nest). 
         on succesful exit, this array will contain the coefficients 
         c(1),c(2),..,c(n-k-1) in the b-spline representation of s(x) 
 fp    : real. unless ier = 10, fp contains the weighted sum of 
         squared residuals of the spline approximation returned. 
 wrk   : real array of dimension at least (m*(k+1)+nest*(8+5*k)). 
         used as working space. if the computation mode iopt=1 is 
         used, the values wrk(1),...,wrk(n) should be left unchanged 
         between subsequent calls. 
 lwrk  : integer. on entry,lwrk must specify the actual dimension of 
         the array wrk as declared in the calling (sub)program. lwrk 
         must not be too small (see wrk). unchanged on exit. 
 iwrk  : integer array of dimension at least (nest). 
         used as working space. if the computation mode iopt=1 is 
         used,the values iwrk(1),...,iwrk(n) should be left unchanged 
         between subsequent calls. 
 ier   : integer. unless the routine detects an error, ier contains a 
         non-positive value on exit, i.e. 
  ier=0  : normal return. the spline returned has a residual sum of 
           squares fp such that abs(fp-s)/s <= tol with tol a relat- 
           ive tolerance set to 0.001 by the program. 
  ier=-1 : normal return. the spline returned is an interpolating 
           periodic spline (fp=0). 
  ier=-2 : normal return. the spline returned is the weighted least- 
           squares constant. in this extreme case fp gives the upper 
           bound fp0 for the smoothing factor s. 
  ier=1  : error. the required storage space exceeds the available 
           storage space, as specified by the parameter nest. 
           probably causes : nest too small. if nest is already 
           large (say nest > m/2), it may also indicate that s is 
           too small 
           the approximation returned is the least-squares periodic 
           spline according to the knots t(1),t(2),...,t(n). (n=nest) 
           the parameter fp gives the corresponding weighted sum of 
           squared residuals (fp>s). 
  ier=2  : error. a theoretically impossible result was found during 
           the iteration proces for finding a smoothing spline with 
           fp = s. probably causes : s too small. 
           there is an approximation returned but the corresponding 
           weighted sum of squared residuals does not satisfy the 
           condition abs(fp-s)/s < tol. 
  ier=3  : error. the maximal number of iterations maxit (set to 20 
           by the program) allowed for finding a smoothing spline 
           with fp=s has been reached. probably causes : s too small 
           there is an approximation returned but the corresponding 
           weighted sum of squared residuals does not satisfy the 
           condition abs(fp-s)/s < tol. 
  ier=10 : error. on entry, the input data are controlled on validity 
           the following restrictions must be satisfied. 
           -1<=iopt<=1, 1<=k<=5, m>1, nest>2*k+2, w(i)>0,i=1,...,m-1 
           x(1)<x(2)<...<x(m), lwrk>=(k+1)*m+nest*(8+5*k) 
           if iopt=-1: 2*k+2<=n<=min(nest,m+2*k) 
                       x(1)<t(k+2)<t(k+3)<...<t(n-k-1)<x(m) 
                     the schoenberg-whitney conditions, i.e. there 
                     must be a subset of data points xx(j) with 
                     xx(j) = x(i) or x(i)+(x(m)-x(1)) such that 
                       t(j) < xx(j) < t(j+k+1), j=k+1,...,n-k-1 
           if iopt>=0: s>=0 
                       if s=0 : nest >= m+2*k 
           if one of these conditions is found to be violated,control 
           is immediately repassed to the calling program. in that 
           case there is no approximation returned. 

further comments: 
 by means of the parameter s, the user can control the tradeoff 
 between closeness of fit and smoothness of fit of the approximation. 
 if s is too large, the spline will be too smooth and signal will be 
 lost ; if s is too small the spline will pick up too much noise. in 
 the extreme cases the program will return an interpolating periodic 
 spline if s=0 and the weighted least-squares constant if s is very 
 large. between these extremes, a properly chosen s will result in 
 a good compromise between closeness of fit and smoothness of fit. 
 to decide whether an approximation, corresponding to a certain s is 
 satisfactory the user is highly recommended to inspect the fits 
 graphically. 
 recommended values for s depend on the weights w(i). if these are 
 taken as 1/d(i) with d(i) an estimate of the standard deviation of 
 y(i), a good s-value should be found in the range (m-sqrt(2*m),m+ 
 sqrt(2*m)). if nothing is known about the statistical error in y(i) 
 each w(i) can be set equal to one and s determined by trial and 
 error, taking account of the comments above. the best is then to 
 start with a very large value of s ( to determine the least-squares 
 constant and the corresponding upper bound fp0 for s) and then to 
 progressively decrease the value of s ( say by a factor 10 in the 
 beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the 
 approximation shows more detail) to obtain closer fits. 
 to economize the search for a good s-value the program provides with 
 different modes of computation. at the first call of the routine, or 
 whenever he wants to restart with the initial set of knots the user 
 must set iopt=0. 
 if iopt=1 the program will continue with the set of knots found at 
 the last call of the routine. this will save a lot of computation 
 time if percur is called repeatedly for different values of s. 
 the number of knots of the spline returned and their location will 
 depend on the value of s and on the complexity of the shape of the 
 function underlying the data. but, if the computation mode iopt=1 
 is used, the knots returned may also depend on the s-values at 
 previous calls (if these were smaller). therefore, if after a number 
 of trials with different s-values and iopt=1, the user can finally 
 accept a fit as satisfactory, it may be worthwhile for him to call 
 percur once more with the selected value for s but now with iopt=0. 
 indeed, percur may then return an approximation of the same quality 
 of fit but with fewer knots and therefore better if data reduction 
 is also an important objective for the user. 

other subroutines required: 
  fpbacp,fpbspl,fpchep,fpperi,fpdisc,fpgivs,fpknot,fprati,fprota 

references: 
 dierckx p. : algorithms for smoothing data with periodic and 
              parametric splines, computer graphics and image 
              processing 20 (1982) 171-184. 
 dierckx p. : algorithms for smoothing data with periodic and param- 
              etric splines, report tw55, dept. computer science, 
              k.u.leuven, 1981. 
 dierckx p. : curve and surface fitting with splines, monographs on 
              numerical analysis, oxford university press, 1993. 

author: 
  p.dierckx 
  dept. computer science, k.u. leuven 
  celestijnenlaan 200a, b-3001 heverlee, belgium. 
  e-mail : Paul.Dierckx@cs.kuleuven.ac.be 

creation date : may 1979 
latest update : march 1987 

Parameters

iopt	Option for how spline is computed.
m	Number of data points.
x	Independent data.
y	Dependent data.
w	Weights for data points.
k	Degree of spline.
s	Smoothing parameter.
nest	Over estimate of the number of knots to be returned (eg m + 2 * (k + 1)).
n	Destination pointer for the number of knots returned.
t	Array for returned computed knots.
c	Array for returned spline coefficients.
fp	Weighted sum of squared residuals of the spline approximation returned.
wrk	Workspace with minimum size (m * (k + 1) + nest * (8 + 5 * k)).
iwrk	Workspace with minimum size nest.

References ALG_ERR_FUNC, and ALG_ERR_NONE.

AlgBSplineFit()

AlgError AlgBSplineFit	(	int	iopt,
		int	m,
		double *	x,
		double *	y,
		double *	w,
		double	xb,
		double	xe,
		int	k,
		double	s,
		int	nest,
		int *	n,
		double *	t,
		double *	c,
		double *	fp,
		double *	wrk,
		int *	iwrk
	)

Computes a smooth B-spline approximation through the given data.

Returns

Alg error code.

This function has been derived from the netlib Dierckx function curfit(). The original fortran coments are:

 given the set of data points (x(i),y(i)) and the set of positive
 numbers w(i),i=1,2,...,m,subroutine curfit determines a smooth spline
 approximation of degree k on the interval xb <= x <= xe.
 if iopt=-1 curfit calculates the weighted least-squares spline
 according to a given set of knots.
 if iopt>=0 the number of knots of the spline s(x) and the position
 t(j),j=1,2,...,n is chosen automatically by the routine. the smooth-
 ness of s(x) is then achieved by minimalizing the discontinuity
 jumps of the k-th derivative of s(x) at the knots t(j),j=k+2,k+3,...,
 n-k-1. the amount of smoothness is determined by the condition that
 f(p)=sum((w(i)*(y(i)-s(x(i))))**2) be <= s, with s a given non-
 negative constant, called the smoothing factor.
 the fit s(x) is given in the b-spline representation (b-spline coef-
 ficients c(j),j=1,2,...,n-k-1) and can be evaluated by means of
 subroutine splev.

 calling sequence:
    call curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,wrk,
   * lwrk,iwrk,ier)

 parameters:
  iopt  : integer flag. on entry iopt must specify whether a weighted
          least-squares spline (iopt=-1) or a smoothing spline (iopt=
          0 or 1) must be determined. if iopt=0 the routine will start
          with an initial set of knots t(i)=xb, t(i+k+1)=xe, i=1,2,...
          k+1. if iopt=1 the routine will continue with the knots
          found at the last call of the routine.
          attention: a call with iopt=1 must always be immediately
          preceded by another call with iopt=1 or iopt=0.
          unchanged on exit.
  m     : integer. on entry m must specify the number of data points.
          m > k. unchanged on exit.
  x     : real array of dimension at least (m). before entry, x(i)
          must be set to the i-th value of the independent variable x,
          for i=1,2,...,m. these values must be supplied in strictly
          ascending order. unchanged on exit.
  y     : real array of dimension at least (m). before entry, y(i)
          must be set to the i-th value of the dependent variable y,
          for i=1,2,...,m. unchanged on exit.
  w     : real array of dimension at least (m). before entry, w(i)
          must be set to the i-th value in the set of weights. the
          w(i) must be strictly positive. unchanged on exit.
          see also further comments.
  xb,xe : real values. on entry xb and xe must specify the boundaries
          of the approximation interval. xb<=x(1), xe>=x(m).
          unchanged on exit.
  k     : integer. on entry k must specify the degree of the spline.
          1<=k<=5. it is recommended to use cubic splines (k=3).
          the user is strongly dissuaded from choosing k even,together
          with a small s-value. unchanged on exit.
  s     : real.on entry (in case iopt>=0) s must specify the smoothing
          factor. s >=0. unchanged on exit.
          for advice on the choice of s see further comments.
  nest  : integer. on entry nest must contain an over-estimate of the
          total number of knots of the spline returned, to indicate
          the storage space available to the routine. nest >=2*k+2.
          in most practical situation nest=m/2 will be sufficient.
          always large enough is  nest=m+k+1, the number of knots
          needed for interpolation (s=0). unchanged on exit.
  n     : integer.
          unless ier =10 (in case iopt >=0), n will contain the
          total number of knots of the spline approximation returned.
          if the computation mode iopt=1 is used this value of n
          should be left unchanged between subsequent calls.
          in case iopt=-1, the value of n must be specified on entry.
  t     : real array of dimension at least (nest).
          on succesful exit, this array will contain the knots of the
          spline,i.e. the position of the interior knots t(k+2),t(k+3)
          ...,t(n-k-1) as well as the position of the additional knots
          t(1)=t(2)=...=t(k+1)=xb and t(n-k)=...=t(n)=xe needed for
          the b-spline representation.
          if the computation mode iopt=1 is used, the values of t(1),
          t(2),...,t(n) should be left unchanged between subsequent
          calls. if the computation mode iopt=-1 is used, the values
          t(k+2),...,t(n-k-1) must be supplied by the user, before
          entry. see also the restrictions (ier=10).
  c     : real array of dimension at least (nest).
          on succesful exit, this array will contain the coefficients
          c(1),c(2),..,c(n-k-1) in the b-spline representation of s(x)
  fp    : real. unless ier=10, fp contains the weighted sum of
          squared residuals of the spline approximation returned.
  wrk   : real array of dimension at least (m*(k+1)+nest*(7+3*k)).
          used as working space. if the computation mode iopt=1 is
          used, the values wrk(1),...,wrk(n) should be left unchanged
          between subsequent calls.
  lwrk  : integer. on entry,lwrk must specify the actual dimension of
          the array wrk as declared in the calling (sub)program.lwrk
          must not be too small (see wrk). unchanged on exit.
  iwrk  : integer array of dimension at least (nest).
          used as working space. if the computation mode iopt=1 is
          used,the values iwrk(1),...,iwrk(n) should be left unchanged
          between subsequent calls.
  ier   : integer. unless the routine detects an error, ier contains a
          non-positive value on exit, i.e.
   ier=0  : normal return. the spline returned has a residual sum of
            squares fp such that abs(fp-s)/s <= tol with tol a relat-
            ive tolerance set to 0.001 by the program.
   ier=-1 : normal return. the spline returned is an interpolating
            spline (fp=0).
   ier=-2 : normal return. the spline returned is the weighted least-
            squares polynomial of degree k. in this extreme case fp
            gives the upper bound fp0 for the smoothing factor s.
   ier=1  : error. the required storage space exceeds the available
            storage space, as specified by the parameter nest.
            probably causes : nest too small. if nest is already
            large (say nest > m/2), it may also indicate that s is
            too small
            the approximation returned is the weighted least-squares
            spline according to the knots t(1),t(2),...,t(n). (n=nest)
            the parameter fp gives the corresponding weighted sum of
            squared residuals (fp>s).
   ier=2  : error. a theoretically impossible result was found during
            the iteration proces for finding a smoothing spline with
            fp = s. probably causes : s too small.
            there is an approximation returned but the corresponding
            weighted sum of squared residuals does not satisfy the
            condition abs(fp-s)/s < tol.
   ier=3  : error. the maximal number of iterations maxit (set to 20
            by the program) allowed for finding a smoothing spline
            with fp=s has been reached. probably causes : s too small
            there is an approximation returned but the corresponding
            weighted sum of squared residuals does not satisfy the
            condition abs(fp-s)/s < tol.
   ier=10 : error. on entry, the input data are controlled on validity
            the following restrictions must be satisfied.
            -1<=iopt<=1, 1<=k<=5, m>k, nest>2*k+2, w(i)>0,i=1,2,...,m
            xb<=x(1)<x(2)<...<x(m)<=xe, lwrk>=(k+1)*m+nest*(7+3*k)
            if iopt=-1: 2*k+2<=n<=min(nest,m+k+1)
                        xb<t(k+2)<t(k+3)<...<t(n-k-1)<xe
                      the schoenberg-whitney conditions, i.e. there
                      must be a subset of data points xx(j) such that
                        t(j) < xx(j) < t(j+k+1), j=1,2,...,n-k-1
            if iopt>=0: s>=0
                        if s=0 : nest >= m+k+1
            if one of these conditions is found to be violated,control
            is immediately repassed to the calling program. in that
            case there is no approximation returned.

 further comments:
  by means of the parameter s, the user can control the tradeoff
  between closeness of fit and smoothness of fit of the approximation.
  if s is too large, the spline will be too smooth and signal will be
  lost ; if s is too small the spline will pick up too much noise. in
  the extreme cases the program will return an interpolating spline if
  s=0 and the weighted least-squares polynomial of degree k if s is
  very large. between these extremes, a properly chosen s will result
  in a good compromise between closeness of fit and smoothness of fit.
  to decide whether an approximation, corresponding to a certain s is
  satisfactory the user is highly recommended to inspect the fits
  graphically.
  recommended values for s depend on the weights w(i). if these are
  taken as 1/d(i) with d(i) an estimate of the standard deviation of
  y(i), a good s-value should be found in the range (m-sqrt(2*m),m+
  sqrt(2*m)). if nothing is known about the statistical error in y(i)
  each w(i) can be set equal to one and s determined by trial and
  error, taking account of the comments above. the best is then to
  start with a very large value of s ( to determine the least-squares
  polynomial and the corresponding upper bound fp0 for s) and then to
  progressively decrease the value of s ( say by a factor 10 in the
  beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
  approximation shows more detail) to obtain closer fits.
  to economize the search for a good s-value the program provides with
  different modes of computation. at the first call of the routine, or
  whenever he wants to restart with the initial set of knots the user
  must set iopt=0.
  if iopt=1 the program will continue with the set of knots found at
  the last call of the routine. this will save a lot of computation
  time if curfit is called repeatedly for different values of s.
  the number of knots of the spline returned and their location will
  depend on the value of s and on the complexity of the shape of the
  function underlying the data. but, if the computation mode iopt=1
  is used, the knots returned may also depend on the s-values at
  previous calls (if these were smaller). therefore, if after a number
  of trials with different s-values and iopt=1, the user can finally
  accept a fit as satisfactory, it may be worthwhile for him to call
  curfit once more with the selected value for s but now with iopt=0.
  indeed, curfit may then return an approximation of the same quality
  of fit but with fewer knots and therefore better if data reduction
  is also an important objective for the user.

 other subroutines required:
   fpback,fpbspl,fpchec,fpcurf,fpdisc,fpgivs,fpknot,fprati,fprota

 references:
  dierckx p. : an algorithm for smoothing, differentiation and integ-
               ration of experimental data using spline functions,
               j.comp.appl.maths 1 (1975) 165-184.
  dierckx p. : a fast algorithm for smoothing data on a rectangular
               grid while using spline functions, siam j.numer.anal.
               19 (1982) 1286-1304.
  dierckx p. : an improved algorithm for curve fitting with spline
               functions, report tw54, dept. computer science,k.u.
               leuven, 1981.
  dierckx p. : curve and surface fitting with splines, monographs on
               numerical analysis, oxford university press, 1993.

 author:
   p.dierckx
   dept. computer science, k.u. leuven
   celestijnenlaan 200a, b-3001 heverlee, belgium.
   e-mail : Paul.Dierckx@cs.kuleuven.ac.be

 creation date : may 1979
 latest update : march 1987

Parameters

iopt	Option for how spline is computed.
m	Number of data points.
x	Independent data.
y	Dependent data.
w	Weights for data points.
xb
xe
k	Degree of spline.
s	Smoothing parameter.
nest	Over estimate of the number of knots to be returned (eg m + k + 1).
n	Destination pointer for the number of knots returned.
t	Array for returned computed knots.
c	Array for returned spline coefficients.
fp	Weighted sum of squared residuals of the spline approximation returned.
wrk	Workspace with minimum size m * k1 + nest * (k * 3 + 7)
iwrk	Workspace with minimum size nest.

References ALG_ERR_FUNC, and ALG_ERR_NONE.

AlgBSplineEval()

AlgError AlgBSplineEval	(	double *	t,
		int	n,
		double *	c,
		int	k,
		double *	x,
		double *	y,
		int	m
	)

Evaluates a b-spline at a number of points.

Returns

Alg error code.

This function has been derived from the netlib Dierckx function splev(). The original fortran coments are:

 Subroutine splev evaluates in a number of points x(i),i=1,2,...,m
 a spline s(x) of degree k, given in its b-spline representation.

 calling sequence:
    call splev(t,n,c,k,x,y,m,ier)

 input parameters:
   t    : array,length n, which contains the position of the knots.
   n    : integer, giving the total number of knots of s(x).
   c    : array,length n, which contains the b-spline coefficients.
   k    : integer, giving the degree of s(x).
   x    : array,length m, which contains the points where s(x) must
          be evaluated.
   m    : integer, giving the number of points where s(x) must be
          evaluated.

 output parameter:
   y    : array,length m, giving the value of s(x) at the different
          points.
   ier  : error flag
     ier = 0 : normal return
     ier =10 : invalid input data (see restrictions)

 restrictions:
   m >= 1
   t(k+1) <= x(i) <= x(i+1) <= t(n-k) , i=1,2,...,m-1.

 other subroutines required: fpbspl.

 references :
   de boor c  : on calculating with b-splines, j. approximation theory
                6 (1972) 50-62.
   cox m.g.   : the numerical evaluation of b-splines, j. inst. maths
                applics 10 (1972) 134-149.
   dierckx p. : curve and surface fitting with splines, monographs on
                numerical analysis, oxford university press, 1993.

 author :
   p.dierckx
   dept. computer science, k.u.leuven
   celestijnenlaan 200a, b-3001 heverlee, belgium.
   e-mail : Paul.Dierckx@cs.kuleuven.ac.be

 latest update : march 1987

Parameters

t	Knot positions.
n	Number of knots.
c	B-spline coefficients.
k	Degree of the B-spline.
x	Points at which the B-spline must be evaluated.
y	Array for evaluated B-spline at given points.
m	Number of points at which B-spline is to be evaluated.

References ALG_CLAMP, ALG_ERR_FUNC, ALG_ERR_NONE, and AlgBSplineBspl().

Referenced by WlzBSplineEval().

AlgBSplineDer()

AlgError AlgBSplineDer	(	double *	t,
		int	n,
		double *	c,
		int	k,
		int	nu,
		double *	x,
		double *	y,
		int	m,
		double *	wrk
	)

Evaluates the derivatives of order nu of the B-spline of degree k.

Returns

Alg error code.

This function has been derived from the netlib Dierckx function splder(). The original fortran coments are:

 Subroutine splder evaluates in a number of points x(i),i=1,2,...,m
 the derivative of order nu of a spline s(x) of degree k,given in
 its b-spline representation.

 calling sequence:
    call splder(t,n,c,k,nu,x,y,m,wrk,ier)

 input parameters:
   t    : array,length n, which contains the position of the knots.
   n    : integer, giving the total number of knots of s(x).
   c    : array,length n, which contains the b-spline coefficients.
   k    : integer, giving the degree of s(x).
   nu   : integer, specifying the order of the derivative. 0<=nu<=k
   x    : array,length m, which contains the points where the deriv-
          ative of s(x) must be evaluated.
   m    : integer, giving the number of points where the derivative
          of s(x) must be evaluated
   wrk  : real array of dimension n. used as working space.

 output parameters:
   y    : array,length m, giving the value of the derivative of s(x)
          at the different points.
   ier  : error flag
     ier = 0 : normal return
     ier =10 : invalid input data (see restrictions)

 restrictions:
   0 <= nu <= k
   m >= 1
   t(k+1) <= x(i) <= x(i+1) <= t(n-k) , i=1,2,...,m-1.

 other subroutines required: fpbspl

 references :
   de boor c : on calculating with b-splines, j. approximation theory
               6 (1972) 50-62.
   cox m.g.  : the numerical evaluation of b-splines, j. inst. maths
               applics 10 (1972) 134-149.
  dierckx p. : curve and surface fitting with splines, monographs on
               numerical analysis, oxford university press, 1993.

 author :
   p.dierckx
   dept. computer science, k.u.leuven
   celestijnenlaan 200a, b-3001 heverlee, belgium.
   e-mail : Paul.Dierckx@cs.kuleuven.ac.be

 latest update : march 1987

Parameters

t	Array of knot positions.
n	Number of knots.
c	B-spline coefficients.
k	Degree of the B-Spline.
nu	Order of the derivative [0-k].
x	Positions at which the derivative of the spline is to be evaluated.
y	Destination for the evaluated derivatives.
m	Number of data points.
wrk	Workspace with minimum size n.

References ALG_CLAMP, ALG_ERR_FUNC, ALG_ERR_NONE, and AlgBSplineBspl().

Referenced by WlzBSplineEval().

AlgLinearFit1D()

AlgError AlgLinearFit1D	(	int	datSz,
		double *	datXA,
		double *	datYA,
		double *	dstA,
		double *	dstB,
		double *	dstSigA,
		double *	dstSigB,
		double *	dstQ
	)

Computes the least squares best fit straight line (y = a + bx) through the given data, ie linear regression. This function is based on the function fit(): Press W. H., Teukolsky S. A., Vetterling W. T. and Flannery B. P, Numerical Recipies in C, 1992, CUP.

Returns

Error code.

Parameters

datSz	Number of elements in given data arrays array.
datXA	Data array with 'x' values.
datYA	Data array with 'y' values.
dstA	Destination ptr for intercept 'a', may be NULL.
dstB	Destination ptr for gradient 'b', may be NULL.
dstSigA	Destination ptr for std dev of 'a', may be NULL.
dstSigB	Destination ptr for std dev of 'b', may be NULL.
dstQ	Destination ptr for goodness of fit, may be NULL.

References ALG_ERR_FUNC, ALG_ERR_NONE, and AlgGammaP().

AlgLinearFitIdx1D()

AlgError AlgLinearFitIdx1D	(	double *	datXA,
		double *	datYA,
		int *	idxXA,
		int *	idxYA,
		int	idxASz,
		double *	dstA,
		double *	dstB,
		double *	dstSigA,
		double *	dstSigB,
		double *	dstQ
	)

Computes the least squares best fit straight line (y = a + bx) through the given data, ie linear regression.

Returns

Error code.

Parameters

datXA	Data array with 'x' values.
datYA	Data array with 'y' values.
idxXA	Index array with indicies into the 'x' data buffer. for the values use examine.
idxYA	Index array with indicies into the 'y' data buffer. for the values use examine.
idxASz	Number of elements in each of the given index arrays.
dstA	Destination ptr for intercept 'a', may be NULL.
dstB	Destination ptr for gradient 'b', may be NULL.
dstSigA	Destination ptr for std dev of 'a', may be NULL.
dstSigB	Destination ptr for std dev of 'b', may be NULL.
dstQ	Destination ptr for goodness of fit, may be NULL.

References ALG_ERR_FUNC, ALG_ERR_NONE, and AlgGammaP().

AlgPolynomialLSq()

AlgError AlgPolynomialLSq	(	double *	xVec,
		double *	yVec,
		int	vecSz,
		int	polyDeg,
		double *	cVec
	)

Attempts to fit a polynomial to the given data using a least squares approach.

Returns

Error code.

Parameters

xVec	Data vector x of size vecSz.
yVec	Data vector y of size vecSz.
vecSz	Size of data vectors.
polyDeg	Degree of ploynomial.
cVec	Destination vector for the polynomial coefficients, which must have at least polyDeg + 1 elements.

References AlcFree(), AlcMalloc(), ALG_DBG, ALG_DBG_LVL_1, ALG_DBG_LVL_FN, ALG_ERR_FUNC, ALG_ERR_MALLOC, AlgMatrixGaussSolve(), AlgMatrixRectFree(), AlgMatrixRectNew(), _AlgMatrixRect::array, _AlgMatrix::core, and _AlgMatrix::rect.