Previous: fitpo2 Up: ../plot79_f.html Next: fitpo4
SUBROUTINE FITPO3 (N,X,Y,Z, XPP,YPP,ZPP, TEMP, ARCLEN, NEPOPT,
X NSIGMA, SIGMA)
C$ (Fit Tensioned P-Spline - 3-D Open - Parameter Evaluation)
C$ Fit a set of points forming an open space curve to a
C$ parametric function dependent upon a set of continuous
C$ tension parameters, SIGMA(*), which can be used to adjust
C$ the shape of the fit.
C$
C$ The input arguments are:
C$
C$ N..............Number of original data points. N .LE. 1
C$ will cause immediate return, but FITPO4 will
C$ still perform correctly, always returning
C$ the first data point.
C$
C$ X(*),
C$ Y(*),
C$ Z(*)...........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$ 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$ NEPOPT......... .LE. 0 - provide internal estimates of
C$ endpoint derivative conditions
C$ .EQ. 1 - endpoint conditions from first
C$ derivatives supplied in XPP(1),
C$ YPP(1), ZPP(1), XPP(N), YPP(N),
C$ and ZPP(N) at entry.
C$ .GE. 2 - endpoint conditions from second
C$ derivatives supplied in XPP(1),
C$ YPP(1), ZPP(1), XPP(N), YPP(N),
C$ and ZPP(N) at entry.
C$
C$ The output arguments are:
C$
C$ ARCLEN.........Arc length of the curve (must be preserved
C$ for FITPO4).
C$ XPP(*),
C$ YPP(*),
C$ ZPP(*).........Arrays of length N containing parameter
C$ values required for the interpolation by
C$ FITPO4.
C$
C$ The scratch argument is:
C$
C$ TEMP(*)........Working storage of at least N locations; it
C$ is not subsequently required by FITPO4.
C$
C$ When derivative estimates for the endpoints are provided in
C$ XPP(1), YPP(1), ZPP(1), XPP(N), YPP(N), and ZPP(N), these
C$ derivatives are with respect to ARC LENGTH.
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, y, z) 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(*), YPP(*), and ZPP(*)
C$ for large SIGMA are of the order of SIGMA, so values of
C$ SIGMA near the machine overflow limit can lead to overflows
C$ in 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(i)..t(i+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$ i i+1 i i i+1 i i
C$
C$ with the linear interpolant given by
C$
C$ L (t) = u (t -t)/h + u (t -t)/h
C$ i i+1 i i i i+1 i
C$
C$ For large tension parameters, the first part becomes
C$ negligible, and the interpolant reduces to the linear term.
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). In
C$ this routine, we choose the latter course, allowing the
C$ user to provide either first or second derivative
C$ estimates, or to automatically supply values from a
C$ polynomial fitted to the endpoints. Since equations 1, 2,
C$ N-1, and N in the tridiagonal system are
C$
C$ a(1)s''(1) + b(1)s''(2) = c(1)
C$
C$ b(1)s''(1) + a(2)s''(2) + b(2)s''(3) = c(2)
C$
C$ b(N-2)s''(N-2) + a(N-1)s''(N-1) + b(N-1)s''(N) = c(N-1)
C$
C$ b(N-1)s''(N-1) + a(N )s''(N) = c(N)
C$
C$ we can discard the first and last and rewrite the other two
C$ as
C$
C$ a(2)s''(2) + b(2)s''(3) = c(2 ) - b(1 )s''(1)
C$
C$ b(N-2)s''(N-2) + a(N-1)s''(N-1) = c(N-1) - b(N-1)s''(N)
C$
C$ This gives us a tridiagonal system with N-2 equations in
C$ N-2 unknowns which can be solved directly.
C$
C$ It is often more convenient and intuitive to provide first
C$ derivative estimates. The first derivative of the
C$ interpolant on the interval k..k+1 is given by
C$
C$
C$ s'(t) = h (s'' F'((t-t )/h ) - s'' F'((t -t)/h ))
C$ k k k+1 k k k k+1 k
C$
C$ + (u - u )/h
C$ k+1 k k
C$
C$ At (x(1),y(1),z(1)), k = 1 and t = 0, and at
C$ (x(N),y(N),z(N)), k = N-1, t = 1, so we have
C$
C$ s'(0) = h (s'' F'(0) - s'' F'(1)) + (u - u )/h
C$ 1 1 2 1 2 1 1
C$
C$ s' (1) = h (s'' F'(1) - s'' F'(0)) + (u - u )/h
C$ N-1 N-1 N N-1 N N-1 N-1
C$
C$ These can be rearranged to define s''(1) and s''(N) which
C$ can be substituted into equations 2 and N-1, resulting in
C$ modified values for a(2), c(2), a(N-1), and c(N-1).
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 as ACM Algorithm
C$ 476; it contains sinh terms which are slow to evaluate and
C$ susceptible to premature underflow/overflow unless
C$ carefully programmed. The function chosen here is Pruess's
C$ second one, a rational function due to H. Spaeth, which is
C$ defined 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$ (20-JUL-89)