Previous: fitpb3 Up: ../plot79_f.html Next: fitpc2
SUBROUTINE FITPC1 (N, X, Y, XPP, YPP, TEMP, ARCLEN,
X NSIGMA, SIGMA)
C$ (Fit Tensioned P-Spline - Closed - Parameter Evaluation)
C$ Fit a set of points forming a closed curve in the XY plane
C$ to a parametric function dependent upon a set of continuous
C$ tension parameters, SIGMA(*), which can be used to adjust
C$ the shape of the fit. The curve is closed in the sense
C$ that (X(N),Y(N)) is connected to (X(1),Y(1)); it is not
C$ necessary to duplicate the first point at the end of the
C$ arrays.
C$
C$ This routine computes the necessary parameters for
C$ subsequent interpolation at individual points by routine
C$ FITPC2.
C$
C$ The input arguments are:
C$
C$ X(*),Y(*)......Original data points. As with most fitting
C$ methods, duplicate adjacent points will lead
C$ to floating-point overflows and invalidate
C$ the fit. Duplicate points which are not
C$ adjacent cause no such problems.
C$
C$ N..............Number of original data points. N .LE. 1
C$ will cause immediate return, but FITPC2 will
C$ still perform correctly, always returning
C$ the first data point.
C$
C$ NSIGMA.........Number of tension parameters supplied in the
C$ array SIGMA(*). NSIGMA should be in the
C$ range 1..N.
C$
C$ SIGMA(*).......Tension parameters in range 0.0..infinity.
C$ The sign of each parameter is ignored. For
C$ generality, each original data point
C$ interval may have associated with it a
C$ separate tension parameter, although for
C$ most applications, a single one will
C$ suffice. If NSIGMA .LT. N, then
C$ SIGMA(NSIGMA) will be used for intervals
C$ NSIGMA..N.
C$
C$ The output arguments are:
C$
C$ ARCLEN.........Arc length of the curve (must be preserved
C$ for FITPC2).
C$ XPP(*),YPP(*)..Arrays of length N containing parameter
C$ values required for the interpolation by
C$ FITPC2.
C$
C$ The scratch argument is:
C$
C$ TEMP(*)........Working storage of at least 2*N locations;
C$ it is not subsequently required by FITPC2.
C$
C$ The effect of the tension parameter SIGMA is independent of
C$ the scale of the coordinates. For small SIGMA values (e.g.
C$ from 0 to about 5), the fitting function is a cubic in the
C$ coordinate u (= x or y) divided by a term linear in u; for
C$ SIGMA = 0, the function is exactly a cubic spline. As
C$ SIGMA becomes large (25 or greater), the fit approaches a
C$ straight-line interpolation. For SIGMA larger than
C$ approximately 10**t with a t-digit floating-point
C$ arithmetic, the fit becomes exact linear interpolation, so
C$ there is no advantage to making SIGMA extremely large. In
C$ fact, the returned parameters XPP(*) and YPP(*) for large
C$ SIGMA are of the order of SIGMA, so values of SIGMA near
C$ the machine overflow limit can lead to overflows in
C$ subsequent function evaluations. Large SIGMA values can
C$ also lead to underflow in this routine, but such underflows
C$ are harmless provided the host system quietly sets them to
C$ zero.
C$
C$ The code is developed according to the discussion in
C$
C$ S. Pruess, "Alternatives to the Exponential Spline in
C$ Tension", Math. Comp. 33, 1273-1281 (1979).
C$
C$ For t in the interval t(k)..t(k+1), Pruess' interpolant
C$ takes the form
C$
C$ 2
C$ s(t) = h (s'' F((t-t )/h ) + s'' F((t -t)/h )) + L (t)
C$ k k+1 k k k k+1 k k
C$
C$ with the linear interpolant given by
C$
C$ L (t) = u (t-t )/h + u (t -t)/h
C$ k k+1 k k k k+1 k
C$
C$ For large tension parameters, the first part becomes
C$ negligible, and the interpolant reduces to the linear term.
C$
C$ Pruess gives 4 functions F(t,p) dependent upon a parameter
C$ p which can be used to interpolate data with the tension of
C$ the fit proportional to p. Each of these is concave
C$ downwards and satisfies a number of common properties, some
C$ of which are:
C$
C$ F(t), F'(t), F''(t) are continuous for t in 0..1
C$ F''(t) >= 0 for t in 0..1
C$ F(0) = F(1) = 0
C$ F''(0) = 0
C$ F''(1) = 1
C$ For p set to its minimum value, F(t) = (t**3 - t)/6.
C$
C$ The first of Preuss' functions is that introduced by D.
C$ Schweikert and programmed by A.K. Cline, "Scalar and Planar
C$ Valued Curve Fitting Using Splines Under Tension", Comm.
C$ A.C.M. 17, 218-225 (1974). (Algorithm 476). Cline's
C$ algorithm contains sinh terms which are slow to evaluate
C$ and susceptible to premature underflow/overflow unless
C$ carefully programmed (which it was not; see Pruess'
C$ article, and also P. Rentrop, "An Algorithm for the
C$ Computation of the Exponential Spline", Numer. Math. 35,
C$ 81-93 (1980)). The function chosen here is Pruess's second
C$ one, a rational function due to H. Spaeth, which is defined
C$ for p in 0..infinity:
C$
C$ F(t,p) = (t**3/(1 + p*(1-t)) - t) / (2p**2 + 6p + 6)
C$
C$ With a simple rearrangement for p > 1, this function can be
C$ evaluated rapidly for any machine-representable p without
C$ danger of overflow, and any underflow is harmless.
C$
C$ For an arbitrary curve in the XY plane, instead of
C$ computing y = x(t), we compute x = x(t), y = y(t), where t
C$ is measured along the arc length of the curve. With k
C$ numbering the intervals 1 .. N-1 between adjacent pairs of
C$ N points, we define
C$
C$ t = arclength from point 1 to point k,
C$ k
C$
C$ t = arclength from point 1 to point (x,y), t <= t <= t.
C$ k k+1
C$
C$ h = t - t = arclength from point k to k+1
C$ k k+1 k
C$
C$ The interpolating polynomial for one coordinate, u(t), is
C$ then
C$
C$ 2
C$ s(t) = h (s'' F((t-t )/h ) + s'' F((t -t)/h )) + L (t)
C$ k k+1 k k k k+1 k k
C$
C$ with the linear interpolant given by
C$
C$ L (t) = u (t-t )/h + u (t -t)/h
C$ k k+1 k k k k+1 k
C$
C$ The second derivatives s''(1..N) are constants to be
C$ determined by boundary matching of derivatives and function
C$ values from adjacent interpolants at each point. This
C$ gives a symmetric tridiagonal system, As'' = c, with
C$ diagonal matrix elements
C$
C$ a = (h + h )F'(1),
C$ k k k-1
C$
C$ sub- and super-diagonal elements
C$
C$ b = -h F'(0),
C$ k k
C$
C$ and right-hand side
C$
C$ c = (u - u )/h - (u - u )/h
C$ k k+1 k k k k-1 k-1
C$
C$ Note that the tridiagonal matrix is the independent of
C$ u(t); in particular, it is the same for u = x, y, and z.
C$ Note also that we have N equations in N+2 unknowns, since
C$ a(1), a(N), c(1), and c(N) are undefined because they
C$ involve the unspecified values h(0) and h(N). For a
C$ non-closed curve, two additional conditions are required to
C$ define h(0) and h(N), or alternatively, two equations can
C$ be eliminated by defining values for s''(1) and s''(N).
C$ For the case of a closed curve treated in this routine, we
C$ can do the former, namely, since u(N+k) = u(k), set h(0) =
C$ h(N) = (t(1) - t(N)). This implies that we must in turn
C$ set s''(0) = s''(N), so that the first equation becomes
C$
C$ a s'' + b s'' + 0 + ... + 0 + b s'' = c
C$ 1 1 1 2 0 N 1
C$
C$ and the last
C$
C$ b s'' + 0 + ... + 0 + b s'' + a s'' = c
C$ 0 1 N-1 N-1 N N N
C$
C$ The spline equations for a closed curve thus require the
C$ solution of linear equations with an almost tridiagonal
C$ symmetric matrix, differing from a pure tridiagonal matrix
C$ in having two extra elements in the (N,1) and (1,N)
C$ positions. For diagonally dominant matrices, which we have
C$ in our case, the system is stable and can be solved without
C$ pivoting. Unfortunately, the extra elements destroy the
C$ extreme simplicity of the tridiagonal case for which the
C$ analytic solution can be written down by inspection. A
C$ full matrix LU decomposition is required, but it can be
C$ shown that the L and U matrices have a simple form. In
C$ fact, the decomposition A = LU has the general structure
C$ illustrated here for a 6 x 6 case
C$
C$ A = L U
C$
C$ | a1 b1 b0 | | 1 | | f1 b1 g1 |
C$ | b1 a2 b2 | | d1 1 | | f2 b2 g2 |
C$ | b2 a3 b3 | | d2 1 | | f3 b3 g3 |
C$ | b3 a4 b4 |=| d3 1 | | f4 b4 g4 |
C$ | b4 a5 b5 | | d4 1 | | f5 g5 |
C$ | b0 b5 a6 | | e1 e2 e3 e4 e5 1 | | f6 |
C$
C$ That is, given the almost tridiagonal matrix defined by 2
C$ N-vectors, a(1..N) and b(1..N-1), we can represent its LU
C$ decomposition by 5 N-vectors, the original b(1..N-2), plus
C$ d(1..N-2), e(1..N-1), f(1..N), and g(1..N-1). The solution
C$ of As = c can be obtained by a forward solution of Lw = c
C$ for w, then a backward solution of Us = w for s. The
C$ general formulae are
C$
C$ LU decomposition:
C$
C$ f(1) = a(1); g(1) = b0; e(1) = b0/f(1);
C$
C$ For k = 1 to N-3 by 1:
C$ d(k) = b(k)/f(k)
C$ f(k+1) = a(k+1) - d(k)*b(k)
C$ g(k+1) = -d(k)*g(k)
C$ e(k+1) = -e(k)*b(k)/f(k+1)
C$
C$ d(N-2) = b(N-2)/f(N-2)
C$ f(N-1) = a(N-1) - d(N-2)*b(N-2)
C$ g(N-1) = b(N-1) - d(N-2)*g(N-2)
C$ e(N-1) = (b(N-1) - e(N-2)*b(N-2))/f(N-1)
C$
C$ N-1
C$ f(N) = a(N) - SUM e(k)*g(k)
C$ k=1
C$
C$ Forward substitution of Lw = c:
C$
C$ w(1) = c(1)
C$
C$ For k = 1 to N-2 by 1:
C$ w(k+1) = c(k+1) - d(k)*w(k)
C$
C$ N-1
C$ w(N) = c(N) - SUM e(k)*w(k)
C$ k=1
C$
C$ Back substitution of Us = w:
C$
C$ s(N) = w(N)/f(N)
C$
C$ s(N-1) = (w(N-1) - g(N-1)*s(N))/f(N-1)
C$
C$ For k = N-2 to 1 by -1:
C$ s(k) = (w(k) - b(k)*s(k+1) - g(k)*s(N))/f(k)
C$
C$ Provided terms with negative subscripts are taken as 0,
C$ these formulas are valid for both N = 2 and N = 1. The
C$ work and storage required are
C$
C$ Step Adds Multiplies Divides Storage
C$ ---- ---- ---------- ------- -------
C$ form L,U 2N-1 4N-1 2N 6N
C$ solve Lw = c 2N-2 2N-2 0 N
C$ solve Us = w 2N-1 2N-1 N 0
C$
C$ Thus, we need about 6N adds, 8N multiplies, 3N divides, and
C$ 7N storage locations, or if p right-hand sides are present,
C$ so that the LU step need only be done once, (2+4p)N adds,
C$ (4+4p)N multiplies, (2+p)N divides, and (6+p)N storage.
C$
C$ For the pure tridiagonal case, b0 = 0, e(1..N-2) = 0,
C$ g(1..N-2) = 0, the figures are
C$
C$ Step Adds Multiplies Divides Storage
C$ ---- ---- ---------- ------- -------
C$ form L,U N-1 N-1 N-1 3N
C$ solve Lw = c N-1 N-1 0 N
C$ solve Us = w N-1 N-1 N 0
C$
C$ or about 3N adds, 3N multiplies, 2N divides, and 4N storage
C$ locations. With p right-hand sides, the requirements are
C$ (2p+1)N adds, (2p+1)N multiplies, (p+1)N divides, and
C$ (3+p)N storage locations. A direct solution in this case
C$ without computation of the LU decomposition takes the same
C$ work and storage, but destroys a(1..N), making it
C$ unsuitable for multiple right hand sides.
C$
C$ The presence of the two corner elements in the closed curve
C$ spline thus increases the work by factors of 2 for adds,
C$ (8/3) for multiplies, (3/2) for divides, and (7/4) for
C$ storage.
C$
C$ In the case of the parametric spline, we have one such
C$ equation As = c for each coordinate, with the same matrix
C$ A. This would then require 7N working locations for a 2-D
C$ spline, and 8N for a 3-D spline. However, if we solve the
C$ two or three equations simultaneously, the LU decomposition
C$ need not be preserved, and we can reduce the storage to 4N
C$ for 2-D and 5N for 3-D.
C$
C$ To see this, observe that merging the forward substitution
C$ with the LU decomposition step removes the need to save
C$ d(1..N-2) and e(1..N-1), since their elements can be used
C$ as they are produced and then discarded. We need only save
C$ wm(1..N)/f(1..N) (meaning element-by-element division)
C$ (overwriting cm(1..N), m = 1 or 2 for 2-D or 3-D spline),
C$ g(1..N-1)/f(1..N-1), and b(1..N-2)/f(1..N-2), and the back
C$ substitution step can produce sm(1..N) (overwriting on
C$ wm(1..N)), after which the storage locations for
C$ g(1..N-1)/f(1..N-1) and b(1..N-2)/f(1..N-2) can be freed.
C$ This is apparently what A.K. Cline did in ACM Algorithm
C$ 476, but he made no mention of the algorithm used for the
C$ linear system solution. This merging is not without cost,
C$ for an operation count in the code below shows that the
C$ total process takes about 10N adds, 12M multiplies, 6N
C$ divides, and 4N storage, compared to 10N adds, 10M
C$ multiplies, 4N divides, and 8N storage with the unmerged
C$ algorithm.
C$
C$ (29-JAN-83)