Previous: fitpc2 Up: ../plot79_f.html Next: fitpc4


FITPC3

       SUBROUTINE  FITPC3 (N, X, Y, Z, XPP, YPP, ZPP, TEMP, ARCLEN,
      X     NSIGMA, SIGMA)
 C$    (Fit Tensioned P-Spline - 3-D Closed - Parameter Evaluation)
 C$    Fit a  set of  points forming  a closed  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.  The  curve is closed  in the  sense
 C$    that (X(N),Y(N),Z(N)) is connected to (X(1),Y(1),Z(1));  it
 C$    is not necessary to duplicate the first point at the end of
 C$    the arrays.
 C$
 C$    This  routine   computes  the   necessary  parameters   for
 C$    subsequent interpolation  at individual  points by  routine
 C$    FITPC4.
 C$
 C$    The input arguments are:
 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$    N..............Number of  original  data  points.  N .LE. 1
 C$                   will cause immediate return, but FITPC4 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 FITPC4).
 C$    XPP(*),
 C$    YPP(*),
 C$    ZPP(*).........Arrays  of  length  N  containing  parameter
 C$                   values required  for  the  interpolation  by
 C$                   FITPC4.
 C$
 C$    The scratch argument is:
 C$
 C$    TEMP(*)........Working  storage of at least  2*N locations;
 C$                   it is not subsequently required by FITPC4.
 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(*),  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(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 space  curve, instead  of computing  y  =
 C$    x(t), we compute x = x(t), y  = y(t), z = z(t), where t  is
 C$    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,z), 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.
 C$
 C$    (02-FEB-83)