Description:
Given the set of data points (x(i),y(i),z(i)) and the set of positive
numbers w(i),i=1,...,m, subroutine pda_surfit determines a smooth
bivariate spline approximation s(x,y) of degrees kx and ky on the
rectangle xb <= x <= xe, yb <= y <= ye.
If iopt = -1 pda_surfit calculates the weighted least-squares spline
according to a given set of knots.
If iopt >= 0 the total numbers nx and ny of these knots and their
position tx(j),j=1,...,nx and ty(j),j=1,...,ny are chosen automatically
by the routine. The smoothness of s(x,y) is then achieved by minimalizing
the discontinuity jumps in the derivatives of s(x,y) across the
boundaries of the sub-panels (tx(i),tx(i+1))*(ty(j),ty(j+1). The amount of
smoothness is determined by the condition that f(p) = sum
((w(i)*(z(i)-s(x(i),y(i))))**2) be <= s, with s a given non-negative
constant, called the smoothing factor.
The fit is given in the b-spline representation (b-spline coefficients
c((ny-ky-1)*(i-1)+j),i=1,...,nx-kx-1;j=1,...,ny-ky-1) and can be
evaluated by means of subroutine bispev.
Calling Sequence:
call pda_surfit( iopt, m, x, y, z, w, xb, xe, yb, ye, kx, ky, s, nxest,
nyest, nmax, eps, nx, tx, ny, ty, c, fp, wrk1, lwrk1,
wrk2, lwrk2, iwrk, kwrk, 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:
tx(i)=xb,tx(i+kx+1)=xe,i=1,...,kx+1;ty(i)=yb,ty(i+ky+1)=ye,i=1,...,ky+1.
If iopt=1 the routine will continue with the set of 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 >= (kx+1)*(ky+1). Unchanged on exit.
x : Real array of dimension at least (m).
y : Real array of dimension at least (m).
z : Real array of dimension at least (m).
Before entry, x(i),y(i),z(i) must be set to the co-ordinates
of the i-th data point, for i=1,...,m. The order of the data
points is immaterial. 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.
xb,xe : Real values. On entry xb,xe,yb and ye must specify the
yb,ye boundaries of the rectangular approximation domain.
xb<=x(i)<=xe,yb<=y(i)<=ye,i=1,...,m. Unchanged on exit.
kx,ky : Integer values. On entry kx and ky must specify the degrees
of the spline. 1<=kx,ky<=5. It is recommended to use bi-cubic
(kx=ky=3) splines. 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
nxest : Integer. Unchanged on exit.
nyest : Integer. Unchanged on exit.
On entry, nxest and nyest must specify an upper bound for the
number of knots required in the x- and y-directions respect.
These numbers will also determine the storage space needed by
the routine. nxest >= 2*(kx+1), nyest >= 2*(ky+1). In most
practical situation nxest = kx+1+sqrt(m/2),
nyest = ky+1+sqrt(m/2) will be sufficient. See also further
comments.
nmax : Integer. On entry nmax must specify the actual dimension of
the arrays tx and ty. Nmax >= nxest, nmax >=nyest.
unchanged on exit.
eps : Real. On entry, eps must specify a threshold for determining the
effective rank of an over-determined linear system of equations.
0 < eps < 1. If the number of decimal digits in the computer
representation of a real number is q, then 10**(-q) is a
suitable value for eps in most practical applications. Unchanged
on exit.
nx : Integer. Unless ier=10 (in case iopt >=0), nx will contain the
total number of knots with respect to the x-variable, of the spline
approximation returned. If the computation mode iopt=1 is used,
the value of nx should be left unchanged between subsequent calls.
In case iopt=-1, the value of nx should be specified on entry
tx : Real array of dimension nmax. On successful exit, this array will
contain the knots of the spline with respect to the x-variable,
i.e. The position of the interior knots tx(kx+2),...,tx(nx-kx-1)
as well as the position of the additional knots
tx(1)=...=tx(kx+1)=xb and tx(nx-kx)=...=tx(nx)=xe needed for
the b-spline representation. If the computation mode iopt=1 is
used, the values of tx(1), ...,tx(nx) should be left unchanged
between subsequent calls. If the computation mode iopt=-1 is
used, the values tx(kx+2), ...tx(nx-kx-1) must be supplied by
the user, before entry. See also the restrictions (ier=10).
ny : Integer. Unless ier=10 (in case iopt >=0), ny will contain the
total number of knots with respect to the y-variable, of the
spline approximation returned. If the computation mode iopt=1
is used, the value of ny should be left unchanged between
subsequent calls. In case iopt=-1, the value of ny should be
specified on entry
ty : Real array of dimension nmax. On successful exit, this array
will contain the knots of the spline with respect to the
y-variable, i.e. The position of the interior knots
ty(ky+2),...,ty(ny-ky-1) as well as the position of the
additional knots ty(1)=...=ty(ky+1)=yb and
ty(ny-ky)=...=ty(ny)=ye needed for the b-spline representation.
If the computation mode iopt=1 is used, the values of ty(1),
...,ty(ny) should be left unchanged between subsequent calls.
If the computation mode iopt=-1 is used, the values ty(ky+2),
...ty(ny-ky-1) must be supplied by the user, before entry. See
also the restrictions (ier=10).
c : Real array of dimension at least (nxest-kx-1)*(nyest-ky-1).
On successful exit, c contains the coefficients of the spline
approximation s(x,y)
fp : Real. Unless ier=10, fp contains the weighted sum of
squared residuals of the spline approximation returned.
wrk1 : Real array of dimension (lwrk1). Used as workspace. If the
computation mode iopt=1 is used the value of wrk1(1) should be
left unchanged between subsequent calls. On exit
wrk1(2),wrk1(3),...,wrk1(1+(nx-kx-1)*(ny-ky-1)) will contain
the values d(i)/max(d(i)),i=1,...,(nx-kx-1)*(ny-ky-1) with d(i)
the i-th diagonal element of the reduced triangular matrix for
calculating the b-spline coefficients. It includes those
elements whose square is less than eps, which are treated as 0
in the case of presumed rank deficiency (ier<-2).
lwrk1 : Integer. On entry lwrk1 must specify the actual dimension of
the array wrk1 as declared in the calling (sub)program.
lwrk1 must not be too small. Let:
u = nxest-kx-1, v = nyest-ky-1, km = max(kx,ky)+1,
ne = max(nxest,nyest), bx = kx*v+ky+1, by = ky*u+kx+1,
if(bx.le.by) b1 = bx, b2 = b1+v-ky
if(bx.gt.by) b1 = by, b2 = b1+u-kx
then
lwrk1 >= u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1
wrk2 : Real array of dimension (lwrk2). Used as workspace, but
only in the case a rank deficient system is encountered.
lwrk2 : Integer. On entry lwrk2 must specify the actual dimension of
the array wrk2 as declared in the calling (sub)program.
lwrk2 > 0 . A safe upper bound for lwrk2 = u*v*(b2+1)+b2
where u,v and b2 are as above. If there are enough data
points, scattered uniformly over the approximation domain
and if the smoothing factor s is not too small, there is a
good chance that this extra workspace is not needed. A lot
of memory might therefore be saved by setting lwrk2=1.
(see also ier > 10).
iwrk : integer array of dimension (kwrk). Used as workspace.
kwrk : Integer. On entry kwrk must specify the actual dimension of
the array iwrk as declared in the calling (sub)program.
kwrk >= m+(nxest-2*kx-1)*(nyest-2*ky-1).
ier : Integer. Unless the routine detects an error, ier contains a
non-positive value on exit:
0 : Normal return. The spline returned has a residual sum of
squares fp such that abs(fp-s)/s <= tol with tol a
relative tolerance set to 0.001 by the program.
-1 : Normal return. The spline returned is an interpolating
spline (fp=0).
-2 : Normal return. The spline returned is the weighted least
squares polynomial of degrees kx and ky. In this extreme
case fp gives the upper bound for the smoothing factor s.
< -2 : Warning. The coefficients of the spline returned have been
computed as the minimal norm least-squares solution of a
(numerically) rank deficient system. (-ier) gives the rank.
Especially if the rank deficiency which can be computed as
(nx-kx-1)*(ny-ky-1)+ier, is large the results may be
inaccurate. They could also seriously depend on the value of
eps.
1 : Error. The required storage space exceeds the available
storage space, as specified by the parameters nxest and
nyest. Probable causes - nxest or nyest too small. If these
parameters are already large, it may also indicate that s is
too small. The approximation returned is the weighted
least-squares spline according to the current set of knots.
The parameter fp gives the corresponding weighted sum of
squared residuals (fp>s).
2 : Error. A theoretically impossible result was found during
the iteration process for finding a smoothing spline with
fp = s. Probable causes - s too small or badly chosen eps.
There is an approximation returned but the corresponding
weighted sum of squared residuals does not satisfy the
condition abs(fp-s)/s < tol.
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. Probable 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.
4 : Error. No more knots can be added because the number of
b-spline coefficients (nx-kx-1)*(ny-ky-1) already exceeds
the number of data points m. Probable causes - either s or
m too small. The approximation returned is the weighted
least-squares spline according to the current set of knots.
The parameter fp gives the corresponding weighted sum of
squared residuals (fp>s).
5 : Error. No more knots can be added because the additional
knot would (quasi) coincide with an old one. Probable
causes - s too small or too large a weight to an inaccurate
data point. The approximation returned is the weighted
least-squares spline according to the current set of knots.
The parameter fp gives the corresponding weighted sum of
squared residuals (fp>s).
10 : Error. On entry, the input data are controlled on validity.
The following restrictions must be satisfied:
-1<=iopt<=1, 1<=kx,ky<=5, m>=(kx+1)*(ky+1), nxest>=2*kx+2,
nyest>=2*ky+2, 0<eps<1, nmax>=nxest, nmax>=nyest,
xb<=x(i)<=xe, yb<=y(i)<=ye, w(i)>0, i=1,...,m
lwrk1 >= u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1
kwrk >= m+(nxest-2*kx-1)*(nyest-2*ky-1)
if iopt=-1: 2*kx+2<=nx<=nxest
xb<tx(kx+2)<tx(kx+3)<...<tx(nx-kx-1)<xe
2*ky+2<=ny<=nyest
yb<ty(ky+2)<ty(ky+3)<...<ty(ny-ky-1)<ye
if iopt>=0: s>=0
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.
> 10 : Error. lwrk2 is too small, i.e. There is not enough work
space for computing the minimal least-squares solution of
a rank deficient system of linear equations. Ier gives the
requested value for lwrk2. There is no approximation
returned but, having saved the information contained in nx,
ny,tx,ty,wrk1, and having adjusted the value of lwrk2 and
the dimension of the array wrk2 accordingly, the user can
continue at the point the program was left, by calling
pda_surfit with iopt=1.
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 (degrees kx,ky)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 z(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 z(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 choose s very small is strongly discouraged.
This considerably increases computation time and memory requirements.
It may also cause rank-deficiency (ier<-2) and endanger numerical
stability.
To economize the search for a good s-value the program provides
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
pda_surfit 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. 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 pda_surfit once more with the
selected value for s but now with iopt=0. Indeed, pda_surfit 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 number of knots may also depend on the upper bounds nxest and
nyest. Indeed, if at a certain stage in pda_surfit the number of knots
in one direction (say nx) has reached the value of its upper bound
(nxest), then from that moment on all subsequent knots are added in the
other (y) direction. This may indicate that the value of nxest is too
small. On the other hand, it gives the user the option of limiting the
number of knots the routine locates in any direction for example, by
setting nxest=2*kx+2 (the lowest allowable value for nxest), the user
can indicate that he wants an approximation which is a simple
polynomial of degree kx in the variable x.
Other Subroutines Required:
pda_fpback, pda_fpbspl, pda_fpsurf, pda_fpdisc, pda_fpgivs, pda_fprank,
pda_fprati, pda_fprota, pda_fporde
References:
Dierckx P. : "An algorithm for surface fitting with spline functions"
Ima J. Numer. Anal. 1 (1981) 267-283.
Dierckx P. : "An algorithm for surface fitting with spline functions"
Report tw50, dept. Computer science,k.u.leuven, 1980.
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