Tension Splines

Greg Blair

2003-12-08 - new page

2003-12-17 - added bezier curves


https://web.archive.org/web/20070208151609/www3.sympatico.ca/gblair/CubicSplinesWithTension.png

I needed some spline with tension tcl code. I hacked Keith Vetter's Cubic Splines GUI and added some tension algorithms. The reference where each algorithm came from is listed in the code.

Enjoy!

BAJ You can run this code via Jacl/Swank/Java Web Start at http://www.onemoonscientific.com/swank/tensionspline.jnlp


 # the next line restarts using wish \
 exec wish "$0" "$@" -visual truecolor

 ## Cubic Tension Splines: TensionSpline.tcl

 ###############################################################################################
 ## Greg Blair 2003-09-29                                                                     ##
 ## I took Keith Vetters Cubic Spline code and added TC and TCB tension, cardinal,            ##
 ## Catmull Rom splines as well as linear and cosine interpolation                            ##
 ## using C code I found at:                                                                  ##
 ##                                                                                           ##
 ## http://astronomy.swin.edu.au/~pbourke/analysis/interpolation/ and                         ##
 ## http://www.cubic.org/~submissive/sourcerer/hermite.htm                                    ##
 ##                                                                                           ##
 ## I hand converted the c code to TCL, and extended Keith's test GUI for the new curve types ##
 ##                                                                                           ##
 ## I have include the comments from the websites as documentation to the methods implemented ##
 ###############################################################################################

 ##Keith Vetter 2003-03-07 - one feature often noted as missing from tk the ability of the canvas to do cubic splines.

 ##Here is a routine PolyLine::CubicSpline that takes a list of (x,y) values of control points and returns a
 ## list points on the cubic spline curve that's suitable to be used by: canvas create line [PolyLine::CubicSplint $xy] ...

 ##+##########################################################################
 #
 # CubicSpline -- routines to generate the coordinates of a cubic spline
 # given a list of control points. Also included is demo/test harness
 # by Keith Vetter
 #
 # Revisions:
 # KPV Mar 07, 2003 - initial revision
 #

 namespace eval Cubic {}
 namespace eval PolyLine {}

 ##+##########################################################################
 #
 # CubicSpline - returns the x,y coordinates of the cubic spline using
 # xy control points.
 # xy => {{x0 y0} {x1 y1} .... {xn yn}}
 #
 # XY points MUST BE SORTED by increasing x
 #
 proc PolyLine::CubicSpline {xy {PRECISION 10}} {

    set np [expr {[llength $xy] / 2}]
    if {$np <= 1} return

    set idx 0
    foreach {x y} $xy {
        set X($idx) $x
        set Y($idx) $y
        incr idx
    }

    for {set i 1; set last $X(0)} {$i < $np} {set last $X($i); incr i} {
        set h($i) [expr {double($X($i) - $last)}]
        if {$h($i) == 0} return
        if {$h($i) < 0} return ;# ERROR not sorted
    }

    if {$np > 2} {
        for {set i 1} {$i < $np-1} {incr i} {
            set i2 [expr {$i + 1}]
            set i0 [expr {$i - 1}]
            set diag($i) [expr {($h($i) + $h($i2))/3.0}]
            set sup($i) [expr {$h($i2) / 6.0}]
            set sub($i) [expr {$h($i) / 6.0}]
            set a($i) [expr {($Y($i2) - $Y($i))/$h($i2) - 
                                 ($Y($i) - $Y($i0)) / $h($i)}]
        }
        PolyLine::SolveTridiag sub diag sup a [expr {$np - 2}]
    }
    set a(0) [set a([expr {$np - 1}]) 0]

    # Now generate the point list
    set xy [list $X(0) $Y(0)]
    for {set i 1} {$i < $np} {incr i} {
        set i0 [expr {$i - 1}]
        for {set j 1} {$j <= $PRECISION} {incr j} {
            set t1 [expr {($h($i) * $j) / $PRECISION}]
            set t2 [expr {$h($i) - $t1}]
            set y [expr {((-$a($i0)/6 * ($t2 + $h($i)) * $t1 + 
                               $Y($i0))* $t2 + (-$a($i)/6 * 
                               ($t1+$h($i)) * $t2 + $Y($i)) * $t1)/$h($i)}]
            set t [expr {$X($i0) + $t1}]
            lappend xy $t $y
        }
    }
    return $xy
 }
 ##+##########################################################################
 # SolveTriDiag -- solves the linear system for tridiagoal NxN matrix A
 # using Gaussian elimination (no pivoting). Since A is sparse, we pass
 # in three diagonals:
 #     sub(i)  => a(i,i-1)    diag(i) => a(i,i)    sup(i)  => a(i,i+1)
 #
 # Result is returned in b[1:n]
 #
 proc PolyLine::SolveTridiag {N_sub N_diag N_sup N_b n} {
    upvar 1 $N_sub sub
    upvar 1 $N_diag diag
    upvar 1 $N_sup sup
    upvar 1 $N_b b

    # Factorization and forward substitution
    for {set i 2} {$i <= $n} {incr i} {
        set i0 [expr {$i - 1}]
        set sub($i) [expr {$sub($i) / $diag($i0)}]
        set diag($i) [expr {$diag($i) - $sub($i) * $sup($i0)}]
        set b($i) [expr {$b($i) - $sub($i) * $b($i0)}]
    }
    set b($n) [expr {$b($n) / $diag($n)}]
    for {set i [expr {$n - 1}]} {$i >= 1} {incr i -1} {
        set i2 [expr {$i + 1}]
        set b($i) [expr {($b($i) - $sup($i) * $b($i2)) / $diag($i)}]
    }
 }

 #########################################################################
 ## from: http://astronomy.swin.edu.au/~pbourke/analysis/interpolation/ ##
 #########################################################################

 ## Discussed here are a number of interpolation methods, this is by no
 ## means an exhaustive list but the methods shown tend to be those in
 ## common use in computer graphics. The main attributes is that they
 ## are easy to compute and are stable. Interpolation as used here is
 ## different to "smoothing", the techniques discussed here have the
 ## characteristic that the estimated curve passes through all the
 ## given points. The idea is that the points are in some sense correct
 ## and lie on an underlying but unknown curve, the problem is to be
 ## able to estimate the values of the curve at any position between
 ## the known points.

 ## Linear interpolation is the simplest method of getting values at
 ## positions in between the data points. The points are simply joined
 ## by straight line segments. Each segment (bounded by two data points)
 ## can be interpolated independently. The parameter mu defines where
 ## to estimate the value on the interpolated line, it is 0 at the
 ## first point and 1 and the second point. For interpolated values
 ## between the two points mu ranges between 0 and 1. Values of mu
 ## outside this range result in extrapolation. This convention is
 ## followed for all the subsequent methods below. As with subsequent
 ## more interesting methods, a snippet of plain C code will server
 ## to describe the mathematics.

 proc LinearInterpolate {y1 y2 mu} {
 # http://astronomy.swin.edu.au/~pbourke/analysis/interpolation/
   ## return [expr {$y1*(1-$mu)+$y2*$mu}] ;# 2 multiplies 1 add, 1 subtract
   return [expr {$y1+$mu*($y2-$y1)}]      ;# 1 multiply   1 add, 1 subtract
 }

 ## Linear interpolation results in discontinuities at each point.
 ## Often a smoother interpolating function is desirable, perhaps
 ## the simplest is cosine interpolation. A suitable orientated
 ## piece of a cosine function serves to provide a smooth transition
 ## between adjacent segments.

 proc CosineInterpolate {y1 y2 mu} {
 # http://astronomy.swin.edu.au/~pbourke/analysis/interpolation/

 ## Paul Burke, on his website does mentions:
 ## By a cute trick the cosine interpolation reverts to linear
 ## if applied independently to each coordinate.

 ## must linear interpolate in x coord, cosine interp in y coord

   set mu2 [expr {(1.-cos(3.14159265358979323846*$mu))/2.}]
   ## return [expr {$y1*(1.-$mu2)+$y2*$mu2}]
   return [expr {$y1+$mu2*($y2-$y1)}]
 }

 ## Cubic interpolation is the simplest method that offers true
 ## continuity between the segments. As such it requires more
 ## than just the two endpoints of the segment but also the two
 ## points on either side of them. So the function requires 4
 ## points in all labeled y0, y1, y2, and y3, in the code below.
 ## mu still behaves the same way for interpolating between
 ## the segment y1 to y2. This does raise issues for how to
 ## interpolate between the first and last segments. In the
 ## examples here I just haven't bothered. A common solution
 ## is the dream up two extra points at the start and end of
 ## the sequence, the new points are created so that they
 ## have a slope equal to the slope of the start or end segment.

 proc CubicInterpolate {y0 y1 y2 y3 mu} {
 # http://astronomy.swin.edu.au/~pbourke/analysis/interpolation/
   set mu2 [expr {$mu*$mu}]
   set a0 [expr {$y3 - $y2 - $y0 + $y1}]
   set a1 [expr {$y0 - $y1 - $a0}]
   set a2 [expr {$y2 - $y0}]
   set a3 $y1

   return [expr {$a0*$mu*$mu2+$a1*$mu2+$a2*$mu+$a3}]
 }

 ## Hermite interpolation like cubic requires 4 points so that
 ## it can achieve a higher degree of continuity. In addition
 ## it has nice tension and biasing controls. Tension can be
 ## used to tighten up the curvature at the known points. The
 ## bias is used to twist the curve about the known points.
 ## The examples shown here have the default tension and
 ## bias values of 0, it will be left as an exercise for
 ## the reader to explore different tension and bias values.

 ##
 ##   Tension: 1 is high, 0 normal, -1 is low
 ##   Bias: 0 is even,
 ##         positive is towards first segment,
 ##         negative towards the other

 proc HermiteInterpolate {y0 y1 y2 y3 mu tension bias} {
 # http://astronomy.swin.edu.au/~pbourke/analysis/interpolation/

 ## GB this is a Kochanek-Bartels Spline
 ## (also called TCB-Splines, for tension,continuity,bias)
 ## GB with continuity c set to zero

 ## GB note variables m0,m1 are the slopes at y1,y2

   set mu2 [expr {$mu * $mu}]
   set mu3 [expr {$mu2 * $mu}]
   set m0  [expr {($y1-$y0)*(1+$bias)*(1-$tension)/2 + 
                  ($y2-$y1)*(1-$bias)*(1-$tension)/2}]
   set m1  [expr {($y2-$y1)*(1+$bias)*(1-$tension)/2 + 
                  ($y3-$y2)*(1-$bias)*(1-$tension)/2}]
   set a0  [expr { 2*$mu3 - 3*$mu2 + 1}]
   set a1  [expr {   $mu3 - 2*$mu2 + $mu}]
   set a2  [expr {   $mu3 -   $mu2}]
   set a3  [expr {-2*$mu3 + 3*$mu2}]

   return [expr {$a0*$y1+$a1*$m0+$a2*$m1+$a3*$y2}]
 }

 ## While you may think the above cases were 2 dimensional,
 ## they are just 1 dimensional interpolation (the horizontal
 ## axis is linear). In most cases the interpolation can be
 ## extended into higher dimensions simply by applying it to
 ## each of the x,y,z coordinates independently. This is
 ## shown on the right for 3 dimensions for all but the
 ## cosine interpolation. By a cute trick the cosine
 ## interpolation reverts to linear if applied independently
 ## to each coordinate.

 proc HermiteSpline {P1 P2 T1 T2 s} {
 # http://www.cubic.org/~submissive/sourcerer/hermite.htm

 #   Vector S: The interpolation-point and it's powers up to 3:
 #   Vector C: The parameters of our hermite curve:
 #   Matrix h: The matrix form of the 4 hermite polynomials:

 #   S = | s^3 s^2 s 1 |        |  2  -2   1   1 |      | P1 |
 #                          h = | -3   3  -2  -1 |  C = | P2 |
 #                              |  0   0   1   0 |      | T1 |
 #                              |  1   0   0   0 |      | T2 |

 #  To calculate a point on the curve you build the Vector S,
 #  multiply it with the matrix h and then multiply with C.
 #    P = S * h * C

    if {[info exist HermiteCache($s,h1)]} {
        set h1 $HermiteCache($s,h1)
        set h2 $HermiteCache($s,h2)
        set h3 $HermiteCache($s,h3)
        set h4 $HermiteCache($s,h4)
    } else {
        set s2 [expr {$s*$s}]
        set s3 [expr {$s*$s2}]
        set h1 [expr { 2.*$s3 - 3.*$s2 + 1.}]
        set h2 [expr {-2.*$s3 + 3.*$s2}]
        set h3 [expr {    $s3 - 2.*$s2+$s}]
        set h4 [expr {    $s3 -    $s2}]
        set HermiteCache($s,h1) $h1
        set HermiteCache($s,h2) $h2
        set HermiteCache($s,h3) $h3
        set HermiteCache($s,h4) $h4
    }
    return [expr {$h1*$P1 + $h2 *$P2 + $h3* $T1 + $h4*$T2}]
 }

 proc CardinalSpline {P0 P1 P2 P3 a s} {
 # http://www.cubic.org/~submissive/sourcerer/hermite.htm

 # Cardinal splines are just a subset of the hermite curves.
 # They don't need the tangent points because they will be calculated from
 # the control points. We'll lose some of the flexibility of the hermite
 # curves, but as a tradeoff the curves will be much easier to use.
 # The formula for the tangents for cardinal splines is:

 #          T   =  a * ( P     -  P      )
 #            i            i+1      i-1

 # a is a constant which affects the tightness of the curve
 # (a should be between 0 and 1, but this is not a must)

    set T1 [expr {$a*($P2-$P0)}]
    set T2 [expr {$a*($P3-$P1)}]
    return [HermiteSpline $P1 $P2 $T1 $T2 $s]
 }

 proc CatmullRomSpline {P0 P1 P2 P3 s} {
 # http://www.cubic.org/~submissive/sourcerer/hermite.htm

 # The Catmull-Rom spline is again just a subset of the cardinal splines.
 # You only have to define a as 0.5, and you can draw and interpolate
 # Catmull-Rom splines.
 #
 #           T   =  0.5 * ( P     -  P      )
 #            i              i+1      i-1

    return [CardinalSpline $P0 $P1 $P2 $P3 0.5 $s]
 }

 proc Lerp {a b mu} {
 # -------------------------------------------------
 # http://www.cubic.org/~submissive/sourcerer/bezier.htm
 # simple linear interpolation between two points
 # -------------------------------------------------
    set ax [lindex $a 0]
    set ay [lindex $a 1]
    set bx [lindex $b 0]
    set by [lindex $b 1]
    return [list [expr {$ax + ($bx-$ax)*$mu}] [expr {$ay + ($by-$ay)*$mu}] ]
 }

 proc BezierSpline {a b c d mu} {
 # --------------------------------------------------------
 # http://www.cubic.org/~submissive/sourcerer/bezier.htm
 # evaluate a point on a bezier-curve. mu goes from 0 to 1.0
 # --------------------------------------------------------

    set  ab   [Lerp $a    $b    $mu]
    set  bc   [Lerp $b    $c    $mu]
    set  cd   [Lerp $c    $d    $mu]
    set  abbc [Lerp $ab   $bc   $mu]
    set  bccd [Lerp $bc   $cd   $mu]
    return    [Lerp $abbc $bccd $mu]
 }

 proc HornerBezier {degree xc yc t} {
    upvar 1 $xc xcoeff
    upvar 1 $yc ycoeff

 # -------------------------------------------------------------------------
 # Greg Blair, algorithm derived from:
 # hornbez(degree, coeff, t)
 # Curves and Surface for Computer Aided Geometric Design, A Practical Guide
 # Farin, Gerald, Academic Press, 2nd Edition, 1990, page 48
 # -------------------------------------------------------------------------

 ## uses  Horner's scheme to compute one coordinate
 ## value of a  Bezier curve. Has to be called
 ## for each coordinate  (x,y, and/or z) of a control polygon.
 ## Input:   degree: degree of curve.
 ##          coeff:  array with coefficients of curve.
 ##          t:      parameter value.
 ## Output: coordinate value.

        set t1 [expr {1.0 - $t}]
    set fact 1.0
        set n_choose_i 1

        set x [expr {$xcoeff(0)*$t1}]
        set y [expr {$ycoeff(0)*$t1}]
        for {set i 1} {$i < $degree} {incr i} {
                set fact [expr {$fact*$t}]
                set n_choose_i [expr {$n_choose_i*($degree-$i+1)/$i}] ;#  always int!
                set x [expr {$x + $fact*$n_choose_i*$xcoeff($i)*$t1}]
                set y [expr {$y + $fact*$n_choose_i*$ycoeff($i)*$t1}]
        }
        set x [expr {$x + $fact*$t*$xcoeff($degree)}]
        set y [expr {$y + $fact*$t*$ycoeff($degree)}]

        return [list $x $y]
 }

 proc TCBSpline {P0 P1 P2 P3 s t c b} {
 # http://www.cubic.org/~submissive/sourcerer/hermite.htm
 ##################################################################
 ## note by Greg Blair:                                          ##
 ##                                                              ##
 ## TCB splines due to Doris Kochanek, a UofWaterloo MSc student ##
 ##                                                              ##
 ## These TCB spline equations can also be found on page 433 of  ##
 ## Splines for use in Computer Graphics & Geometric Modelling   ##
 ## Bartels, Beatty, Barsky, Morgan Kaufman, 1987                ##
 ##################################################################

 # The Kochanek-Bartels Splines
 # (also called TCB-Splines, for tension,continuity,bias)
 # Now we're going down to the guts of curve interpolation:
 # The kb-splines (mostly known from Autodesk's 3d-Studio Max and Newtek's Lightwave)
 # are nothing more than hermite curves and a handfull of formulas to calculate the tangents.

 # These curves have been introduced by D. Kochanek and R. Bartels in 1984 to
 # give animators more control over keyframe animation.
 # They introduced three control-values for each keyframe point:

 # Tension:    How sharply does the curve bend?
 # Continuity: How rapid is the change in speed and direction?
 # Bias:       What is the direction of the curve as it passes through the keypoint?

 # I won't try to derive the tangent-formulas here. I think just giving you something you can use
 # is a better idea. (if you're interested you can ask me. I can write it down and send it to you
 # via email.)

 #    The "incoming" Tangent equation:

 #              (1-t)*(1-c)*(1+b)
 #    TS    =   -----------------  * ( P   -  P    )
 #      i              2                i      i-1

 #              (1-t)*(1+c)*(1-b)
 #          +   -----------------  * ( P   -  P    )
 #                     2                i+1    i

 #    The "outgoing" Tangent equation:

 #              (1-t)*(1+c)*(1+b)
 #    TD    =   -----------------  * ( P   -  P    )
 #      i              2                i      i-1

 #              (1-t)*(1-c)*(1-b)
 #          +   -----------------  * ( P   -  P    )
 #                     2                i+1    i

 #    When you want to interpolate the curve you should use this vector:

 #        |  P(i)    |
 #    C = |  P(i+1)  |
 #        |  TD(i)   |
 #        |  TS(i+1) |

 # You might notice that you always need the previous and next point
 # if you want to calculate the curve. This might be a problem when you try
 # to calculate keyframe data from Lightwave or 3D-Studio. I don't
 # know exactly how these programs handle the cases of the first and last point, but
 # there are enough sources available on the internet. Just search around
 # a little bit. (Newtek has a good developer section. You can download
 # the origignal Lightwave motion code on their web-site).

    set TD1 [expr {(1.-$t)*(1.-$c)*(1.+$b)/2*($P1 - $P0) 
                 + (1.-$t)*(1.+$c)*(1.-$b)/2*($P2 - $P1)}]
    set TS2 [expr {(1.-$t)*(1.+$c)*(1.+$b)/2*($P2 - $P1) 
                 + (1.-$t)*(1.-$c)*(1.-$b)/2*($P3 - $P2)}]

    return [HermiteSpline $P1 $P2 $TD1 $TS2 $s]

 }

 proc KeyFrameTCBSpline {P0 P1 P2 P3 s t c b N0 N1 N2} {
 # http://www.cubic.org/~submissive/sourcerer/hermite.htm

 # Speed Control
 # If you write yourself keyframe-interpolation code and put it into a program
 # you'll notice one problem:

 # Unless you have your keyframes in fixed intervals you will have a sudden change
 # of speed and direction whenever you pass a keyframe-point.
 # This can be avoided if you take the number of key-positions (frames) between two keyframes into account:

 # N is the number of frames (seconds, whatever) between two keypoints.

 #                       2 * N
 #                            i-1
 #   TD  =  TD *     ---------------       adjustment of outgoing tangent
 #     i      i          N    + N
 #                        i-1    i

 #                       2 * N
 #                            i
 #   TS  =  TS *     ---------------       adjustment of incomming tangent
 #     i      i          N    + N
 #                        i-1    i

    set TD1 [expr {(1.-$t)*(1.-$c)*(1.+$b)/2*($P1 - $P0) 
                 + (1.-$t)*(1.+$c)*(1.-$b)/2*($P2 - $P1)}]
    set TS2 [expr {(1.-$t)*(1.+$c)*(1.+$b)/2*($P2 - $P1) 
                 + (1.-$t)*(1.-$c)*(1.-$b)/2*($P3 - $P2)}]

    set D [expr {$N0 + $N1}]
    if {$D} {
        set TD1 [expr {$TD1 * 2.*$N0/$D}]
    }
    set D [expr {$N1+$N2}]
    if {$D} {
        set TS2 [expr {$TS2 * 2.*$N2/$D}]
    }
    return [HermiteSpline $P1 $P2 $TD1 $TS2 $s]
 }

 proc PolyLine::Linear {xy {PRECISION 10}} {

    set np [expr {[llength $xy] / 2}]
    if {$np <= 1} return

    set idx 0
    foreach {x y} $xy {
        set X($idx) $x
        set Y($idx) $y
        incr idx
    }

    # Now generate the point list
    set xy [list $X(0) $Y(0)]
    for {set i 1} {$i < $np} {incr i} {
        set i0 [expr {$i - 1}]
        for {set j 0} {$j <= $PRECISION} {incr j} {
            set mu [expr {double($j) / double($PRECISION)}]
            set x [LinearInterpolate $X($i0) $X($i) $mu]
            set y [LinearInterpolate $Y($i0) $Y($i) $mu]
            lappend xy $x $y
        }
    }
    return $xy
 }

 proc PolyLine::Cosine {xy {PRECISION 10}} {

    set np [expr {[llength $xy] / 2}]
    if {$np <= 1} return

    set idx 0
    foreach {x y} $xy {
        set X($idx) $x
        set Y($idx) $y
        incr idx
    }

    # Now generate the point list
    set xy [list $X(0) $Y(0)]
    for {set i 1} {$i < $np} {incr i} {
        set i0 [expr {$i - 1}]
        for {set j 0} {$j <= $PRECISION} {incr j} {
            set mu [expr {double($j) / double($PRECISION)}]
            ## Paul Burke, on the website does mention:
            ## By a cute trick the cosine interpolation reverts to linear
            ## if applied independently to each coordinate.

            ## Now
            set x [CosineInterpolate $X($i0) $X($i) $mu]
            ## does generate linear interpolation

            ## linear interpolation on x and cosine interpolation on y
            set x [LinearInterpolate $X($i0) $X($i) $mu]
            set y [CosineInterpolate $Y($i0) $Y($i) $mu]
            lappend xy $x $y
        }
    }
    return $xy
 }

 proc PolyLine::CatmullRom {xy {PRECISION 10}} {

    set np [expr {[llength $xy] / 2}]
    if {$np <= 1} return

    set idx 0
    foreach {x y} $xy {
        set X($idx) $x
        set Y($idx) $y
        incr idx
        if {$idx == 1} {
            set X($idx) $x
            set Y($idx) $y
            incr idx
        }
    }
    ## duplicate last point
    set X($idx) $x
    set Y($idx) $y
    incr idx

    # Now generate the point list
    set xy [list $X(0) $Y(0)]
    for {set i 1} {$i < $np} {incr i} {
        set i0 [expr {$i - 1}]
        set i2 [expr {$i + 1}]
        set i3 [expr {$i + 2}]
        for {set j 0} {$j <= $PRECISION} {incr j} {
            set mu [expr {double($j) / double($PRECISION)}]
            set x [CatmullRomSpline $X($i0) $X($i) $X($i2) $X($i3) $mu]
            set y [CatmullRomSpline $Y($i0) $Y($i) $Y($i2) $Y($i3) $mu]
            lappend xy $x $y
        }
    }
    return $xy
 }

 proc PolyLine::Bezier {xy {PRECISION 10}} {

    set np [expr {[llength $xy] / 2}]
    if {$np < 4} return

    set idx 0
    foreach {x y} $xy {
        set X($idx) $x
        set Y($idx) $y
        incr idx
    }

    set xy {}

    set idx 0
    while {[expr {$idx+4}] <= $np} {
        set a [list $X($idx) $Y($idx)]; incr idx
        set b [list $X($idx) $Y($idx)]; incr idx
        set c [list $X($idx) $Y($idx)]; incr idx
        set d [list $X($idx) $Y($idx)];# incr idx   ;# use last pt as 1st pt of next segment
        for {set j 0} {$j <= $PRECISION} {incr j} {
            set mu [expr {double($j) / double($PRECISION)}]
            set pt [BezierSpline $a $b $c $d $mu]
            lappend xy [lindex $pt 0] [lindex $pt 1]
        }
    }
    return $xy
 }

 proc PolyLine::FDBezier {xy {PRECISION 10}} { ;# forward differencing Bezier calc.
 ## http://www.niksula.cs.hut.fi/~hkankaan/Homepages/bezierfast.html
 ## Calculating bezier curves is quite processor intensive task, but
 ## this algorithm will make calculating at least 10 times faster.
 ## Especially with bezier surfaces of bezier spaces, you will need the speed.

 ## All polynomial functions can be calculated via forward differencing.
 ## The key idea is to start at some point of the function, move forwards
 ## at constant step and use Taylor series to calculate the next value.

 ## Further reading

 ## The algo I presented here may not be the best one, but it works fine
 ## for most 3d applications. If you need more, read these:

 ## Adaptive forward differencing for rendering curves and surfaces;
 ## Sheue-Ling Lien, Michael Shantz and Vaughan Pratt;
 ## Proceedings of the 14th annual conference on Computer graphics, 1987, Pages 111 - 118

 ## Rendering trimmed NURBS with adaptive forward differencing;
 ## M. Shantz and Sheue-Ling Chang;
 ## Proceedings of the 15th annual conference on Computer graphics, 1988, Pages 189 - 198

 ## Rendering cubic curves and surfaces with integer adaptive forward differencing;
 ## S.-L. Chang and M. S. R. Rocchetti;
 ## Conference proceedings on Computer graphics, 1989, Pages 157 - 166

    set np [expr {[llength $xy] / 2}]
    if {$np < 4} return

    set idx 0
    foreach {x y} $xy {
        set X($idx) $x
        set Y($idx) $y
        incr idx
    }

    set t [expr {1.0 / double($PRECISION)}]
    set temp [expr { $t * $t}]

    set idx 0
    set xy {}
    while {[expr {$idx+4}] <= $np} {    ;# segment loop
        set X0 $X($idx)
        set Y0 $Y($idx); incr idx
        set X1 $X($idx)
        set Y1 $Y($idx); incr idx
        set X2 $X($idx)
        set Y2 $Y($idx); incr idx
        set X3 $X($idx)
        set Y3 $Y($idx);# re-use last pt for 1st pt of next segment

        set xf $X0
        set yf $Y0
        set xfd         [expr { 3 * ($X1 - $X0) * $t}]
        set yfd         [expr { 3 * ($Y1 - $Y0) * $t}]
        set xfdd_per_2  [expr { 3 * ($X0 - 2 * $X1 + $X2) * $temp}]
        set yfdd_per_2  [expr { 3 * ($Y0 - 2 * $Y1 + $Y2) * $temp}]
        set xfddd_per_2 [expr { 3 * (3 * ($X1 - $X2) + $X3 - $X0) * $temp * $t}]
        set yfddd_per_2 [expr { 3 * (3 * ($Y1 - $Y2) + $Y3 - $Y0) * $temp * $t}]

        set xfddd       [expr { $xfddd_per_2 + $xfddd_per_2}]
        set yfddd       [expr { $yfddd_per_2 + $yfddd_per_2}]
        set xfdd        [expr { $xfdd_per_2 + $xfdd_per_2}]
        set yfdd        [expr { $yfdd_per_2 + $yfdd_per_2}]
        set xfddd_per_6 [expr { $xfddd_per_2 * (1.0 / 3)}]
        set yfddd_per_6 [expr { $yfddd_per_2 * (1.0 / 3)}]

        for {set j $PRECISION} {$j} {incr j -1} {   ;# pt in segment loop
            lappend xy $xf $yf

            set xf         [expr { $xf + $xfd + $xfdd_per_2 + $xfddd_per_6}]
            set yf         [expr { $yf + $yfd + $yfdd_per_2 + $yfddd_per_6}]
            set xfd        [expr { $xfd + $xfdd + $xfddd_per_2}]
            set yfd        [expr { $yfd + $yfdd + $yfddd_per_2}]
            set xfdd       [expr { $xfdd + $xfddd}]
            set yfdd       [expr { $yfdd + $yfddd}]
            set xfdd_per_2 [expr { $xfdd_per_2 + $xfddd_per_2}]
            set yfdd_per_2 [expr { $yfdd_per_2 + $yfddd_per_2}]
        }
    }
    lappend xy $X3 $Y3

    return $xy
 }

 proc PolyLine::HornBezier {xy {PRECISION 10}} {

    set np [expr {[llength $xy] / 2}]
    if {$np < 3} return

    set idx 0
    foreach {x y} $xy {
        set X($idx) $x
        set Y($idx) $y
        incr idx
    }

    set degree [expr { $np -1 }]

    # Now generate the point list
    # 1st knot
    set xy [list $X(0) $Y(0)]
    # interpolate for intermediate points
    for {set j 1} {$j < $PRECISION} {incr j} {
        set t [expr {double($j) / double($PRECISION)}]
        set pt [HornerBezier $degree X Y $t]
        lappend xy [lindex $pt 0] [lindex $pt 1]
    }
    # 2nd knot
    lappend xy $X($degree) $Y($degree)

    return $xy
 }

 proc PolyLine::Cardinal {xy a {PRECISION 10}} {

    set np [expr {[llength $xy] / 2}]
    if {$np <= 1} return

    set idx 0
    foreach {x y} $xy {
        set X($idx) $x
        set Y($idx) $y
        incr idx
        if {$idx == 1} {
            set X($idx) $x
            set Y($idx) $y
            incr idx
        }
    }
    ## duplicate last point
    set X($idx) $x
    set Y($idx) $y
    incr idx

    # Now generate the point list
    set xy [list $X(0) $Y(0)]
    for {set i 1} {$i < $np} {incr i} {
        set i0 [expr {$i - 1}]
        set i2 [expr {$i + 1}]
        set i3 [expr {$i + 2}]
        for {set j 0} {$j <= $PRECISION} {incr j} {
            set mu [expr {double($j) / double($PRECISION)}]
            set x [CardinalSpline $X($i0) $X($i) $X($i2) $X($i3) $a $mu]
            set y [CardinalSpline $Y($i0) $Y($i) $Y($i2) $Y($i3) $a $mu]
            lappend xy $x $y
        }
    }
    return $xy
 }

 proc PolyLine::Cubic {xy {PRECISION 10}} {

    set np [expr {[llength $xy] / 2}]
    if {$np <= 1} return

    set idx 0
    foreach {x y} $xy {
        set X($idx) $x
        set Y($idx) $y
        incr idx
        if {$idx == 1} {
            set X($idx) $x
            set Y($idx) $y
            incr idx
        }
    }
    ## duplicate last point
    set X($idx) $x
    set Y($idx) $y
    incr idx

    # Now generate the point list
    set xy [list $X(0) $Y(0)]
    for {set i 1} {$i < $np} {incr i} {
        set i0 [expr {$i - 1}]
        set i2 [expr {$i + 1}]
        set i3 [expr {$i + 2}]
        for {set j 0} {$j <= $PRECISION} {incr j} {
            set mu [expr {double($j) / double($PRECISION)}]
            set x [CubicInterpolate $X($i0) $X($i) $X($i2) $X($i3) $mu]
            set y [CubicInterpolate $Y($i0) $Y($i) $Y($i2) $Y($i3) $mu]
            lappend xy $x $y
        }
    }
    return $xy
 }

 proc PolyLine::Hermite {xy tension bias {PRECISION 10}} {

    set np [expr {[llength $xy] / 2}]
    if {$np <= 1} return

    set idx 0
    foreach {x y} $xy {
        set X($idx) $x
        set Y($idx) $y
        incr idx
        if {$idx == 1} {
            set X($idx) $x
            set Y($idx) $y
            incr idx
        }
    }
    ## duplicate last point
    set X($idx) $x
    set Y($idx) $y
    incr idx

    # Now generate the point list
    set xy [list $X(0) $Y(0)]
    for {set i 1} {$i < $np} {incr i} {
        set i0 [expr {$i - 1}]
        set i2 [expr {$i + 1}]
        set i3 [expr {$i + 2}]
        for {set j 0} {$j <= $PRECISION} {incr j} {
            set mu [expr {double($j) / double($PRECISION)}]
            set x [HermiteInterpolate $X($i0) $X($i) $X($i2) $X($i3) $mu $tension $bias]
            set y [HermiteInterpolate $Y($i0) $Y($i) $Y($i2) $Y($i3) $mu $tension $bias]
            lappend xy $x $y
        }
    }
    return $xy
 }
 proc PolyLine::TCB {xy tension continuity bias {PRECISION 10}} {

    set np [expr {[llength $xy] / 2}]
    if {$np <= 1} return

    set idx 0
    foreach {x y} $xy {
        set X($idx) $x
        set Y($idx) $y
        incr idx
        if {$idx == 1} {
            set X($idx) $x
            set Y($idx) $y
            incr idx
        }
    }
    ## duplicate last point
    set X($idx) $x
    set Y($idx) $y
    incr idx

    # Now generate the point list
    set xy [list $X(0) $Y(0)]
    for {set i 1} {$i < $np} {incr i} {
        set i0 [expr {$i - 1}]
        set i2 [expr {$i + 1}]
        set i3 [expr {$i + 2}]
        for {set j 0} {$j <= $PRECISION} {incr j} {
            set mu [expr {double($j) / double($PRECISION)}]
            set x [TCBSpline $X($i0) $X($i) $X($i2) $X($i3) $mu $tension $continuity $bias]
            set y [TCBSpline $Y($i0) $Y($i) $Y($i2) $Y($i3) $mu $tension $continuity $bias]
            lappend xy $x $y
        }
    }
    return $xy
 }
 proc PolyLine::KeyFrameTCB {xy tension continuity bias {PRECISION 10}} {

    set np [expr {[llength $xy] / 2}]
    if {$np <= 1} return

    set idx 0
    foreach {x y} $xy {
        set X($idx) $x
        set Y($idx) $y
        incr idx
        if {$idx == 1} { ;## duplicate 1st pt
            set X($idx) $x
            set Y($idx) $y
            incr idx
        }
    }
    ## duplicate last point
    set X($idx) $x
    set Y($idx) $y
    incr idx

    # Now generate the point list
    set xy [list $X(0) $Y(0)]
    for {set i 1} {$i < $np} {incr i} {
        set i0 [expr {$i - 1}]
        set i2 [expr {$i + 1}]
        set i3 [expr {$i + 2}]
        set X0 $X($i0)
        set X1 $X($i)
        set X2 $X($i2)
        set X3 $X($i3)
        set Y0 $Y($i0)
        set Y1 $Y($i)
        set Y2 $Y($i2)
        set Y3 $Y($i3)
        set N0 [expr {$X1 - $X0}]
        set N1 [expr {$X2 - $X1}]
        set N2 [expr {$X3 - $X3}]
        for {set j 0} {$j <= $PRECISION} {incr j} {
            set mu [expr {double($j) / double($PRECISION)}]
            set x [KeyFrameTCBSpline $X0 $X1 $X2 $X3 $mu $tension $continuity $bias $N0 $N1 $N2]
            set y [KeyFrameTCBSpline $Y0 $Y1 $Y2 $Y3 $mu $tension $continuity $bias $N0 $N1 $N2]
            lappend xy $x $y
        }
    }
    return $xy
 }
 if {[file tail [info script]] == [file tail $argv0]} {
 ##Now to demonstrate and test the code you can use the following:

 ################################################################
 ################################################################
 #
 # Test harness and demo code
 #

 package require Tk

 set S(title) "Cubic Splines with Tension"
 set S(r) 5                                      ;# Control point size
 set S(w) 5                                      ;# Line width
 set S(precision) 10
 set INITPOINTS {-200 0 -80 -125 30 100 200 0}

 proc Reset {} {
    global S

    set S(tension) 0
    set S(continuity) 0
    set S(bias) 0
    set S(a) 0.5
    DrawCurve
 }

 proc BrushedMetal {c width height} {
    ## from wiki https://wiki.tcl-lang.org/9776

    return
    ## for {set row 0} {$row < $height} {incr row 1} {
        ## set line_color [expr {45000+int(1000000*rand())%3000}]
        ## $c create line 0 $row $width $row -width 1 
                ## -fill [format "#%04x%04x%04x" 
                ## $line_color $line_color $line_color]
    ## }
 }

 proc DoDisplay {} {
    global S INITPOINTS

    wm title . $S(title)
    pack [frame .ctrl -relief ridge -bd 2 -padx 5 -pady 5] \
        -side right -fill both -ipady 5
    pack [frame .screen -bd 2 -relief raised] -side top -fill both -expand 1

    canvas .c -relief raised -borderwidth 0 -height 500 -width 500
    BrushedMetal .c 500 500
    pack .c   -in .screen -side top    -fill both -expand 1
    bind all <Alt-c> {catch {console show}}
    bind .c <Configure> {ReCenter %W %h %w}
    .c bind p <B1-Motion> [list MouseMove %x %y]

    bind all <Escape> {exit}

    DoCtrlFrame
    AddCtrlPoint $INITPOINTS
    update
 }
 proc DoCtrlFrame {} {
    frame .buttons
    pack [button .buttons.add   -text "Add Pt"     -command AddCtrlPoint] -side left
    pack [button .buttons.dele  -text "Del Pt"  -command DeleteCtrlPoint] -side left
    pack [button .buttons.clear -text "Clear Pts"  -command ClearCtrlPoint] -side left

    set font {Helvetica 14 bold}
    scale .prec -width 10 -orient h -variable S(precision) -font $font \
        -label Precision: -relief ridge -from 1 -to 40 -command DrawCurve

    scale .bias -width 10 -orient h -variable S(bias) -font $font \
        -label "B - Bias:" -relief ridge -from -1. -to 1. -resolution 0.01 -command DrawCurve

    scale .tens -width 10 -orient h -variable S(tension) -font $font \
        -label "T - Tension:" -relief ridge -from -1. -to 1. -resolution 0.01 -command DrawCurve

    scale .cont -width 10 -orient h -variable S(continuity) -font $font \
        -label "C - Continuity:" -relief ridge -from -1. -to 1. -resolution 0.01 -command DrawCurve

    scale .a    -width 10 -orient h -variable S(a) -font $font \
        -label "Cardinal Spline - A:" -relief ridge -from 0. -to 1. -resolution 0.01 -command DrawCurve

    button .reset -text Reset -command {Reset}

    labelframe .rb -text "INTERPOLATOR" -relief ridge -bd 2

    set rbList {} ;# radio button list
    lappend rbList {line "Linear"}
    lappend rbList {cosi "Cosine"}
    lappend rbList {cubp "Cubic Polynomial"}
    lappend rbList {cubs "Cubic Spline - Keith Vetter"}
    lappend rbList {catm "CatmullRom Spline"}
    lappend rbList {card "Cardinal Spline (A)"}
    lappend rbList {tbsp "TB Spline"}
    lappend rbList {tcbs "TCB Spline"}
    lappend rbList {bezi "Piecewise Bezier Spline"}
    lappend rbList {fbez "Fwd Dif Piecewise Bezier Spline"}
    lappend rbList {hbez "Horner Bezier Spline"}
    ## KeyFrame splines are drawing a linear line, not a spline for 1st segment
    ## I guess our assumption that we can test these splines by assuming
    ## the number of frames between key frames is the difference in X coordinates
    ## is not adequate.
    ## remove from GUI:
    ## lappend rbList {kfs  "Key Frame TCB Spline"}

    set i 1
    foreach entry $rbList {
        set key  [lindex $entry 0]
        set text [lindex $entry 1]
        pack [frame .rb.l$i]
        radiobutton .rb.l$i.$key -variable S(type) -value $key -font $font \
            -width 30 -text $text -command DrawCurve -anchor w
        pack .rb.l$i.$key -side top -pady 2 -anchor w -fill x
        incr i
    }

    pack [button .buttons.about -text About -command About] -side left
    pack [button .buttons.exit  -text Exit  -command exit] -side left

    set row 1
    grid .buttons   -in .ctrl -row $row -sticky ew ; incr row
    #grid .add   -in .ctrl -row $row -sticky ew
    #grid .dele  -in .ctrl -row $row -sticky ew
    #grid .clear -in .ctrl -row $row -sticky ew ; incr row
    grid .prec  -in .ctrl -row $row -sticky ew ; incr row
    grid .tens  -in .ctrl -row $row -sticky ew ; incr row
    grid .cont  -in .ctrl -row $row -sticky ew ; incr row
    grid .bias  -in .ctrl -row $row -sticky ew ; incr row
    grid .a     -in .ctrl -row $row -sticky ew ; incr row
    grid .reset -in .ctrl -row $row -sticky ew ; incr row
    grid .rb    -in .ctrl -row $row -sticky ew ; incr row
    #grid .about -in .ctrl -row $row -sticky ew ; incr row
    #grid .exit  -in .ctrl -row $row -sticky ew ; incr row
 }
 proc ReCenter {W h w} {                   ;# Called by configure event
    foreach h2 [expr {$h / 2}] w2 [expr {$w / 2}] break
    $W config -scrollregion [list -$w2 -$h2 $w2 $h2]
    BrushedMetal $W [expr {2 * $w2}] [expr {2 * $h2}]
 }
 proc MakeBox {x y r} {
    return [list [expr {$x-$r}] [expr {$y-$r}] [expr {$x+$r}] [expr {$y+$r}]]
 }
 proc MouseMove {X Y} {
    regexp {p([0-9]+)} [.c itemcget current -tag] => who
    set X [.c canvasx $X] ; set Y [.c canvasy $Y]
    foreach x $::P($who,x)   y $::P($who,y) break
    foreach ::P($who,x) $X   ::P($who,y) $Y break
    .c move p$who [expr {$X - $x}] [expr {$Y - $y}]
    DrawCurve
 }
 proc AddCtrlPoint {{xy {}}} {
    global P S

    set np [llength [array names P *x]]

    if {$xy == {}} {
        set w [expr {[winfo width .c] - 50}]
        set xy [list [expr {$w * rand() - $w/2}] [expr {50 * rand() - 25}]]
    }
    foreach {x y} $xy {
        set P($np,x) $x
        set P($np,y) $y
        .c create oval [MakeBox $x $y $S(r)] -tag [list p p$np] -fill yellow
        incr np
    }
    DrawCurve
 }
 proc DeleteCtrlPoint {} {
    global P

    set np [llength [array names P *x]]
    if {$np == 0} return
    incr np -1

    # Always delete the rightmost control point
    # swap RIGHTMOST and NP then delete NP
    set rightmost [lindex [lindex [SortPoints] end] end]
    .c delete p$rightmost
    .c itemconfig p$np -tag [list p p$rightmost]
    set P($rightmost,x) $P($np,x)
    set P($rightmost,y) $P($np,y)

    unset P($np,x)
    unset P($np,y)

    DrawCurve
 }
 proc ClearCtrlPoint {} {
    global P
    .c delete p
    catch {unset P}
    DrawCurve
 }
 proc SortPoints {} {
    global P

    set np [llength [array names P *x]]
    set xy {}
    for {set i 0} {$i < $np} {incr i} {
        lappend xy [list $P($i,x) $P($i,y) $i]
    }
    set xy2 [lsort -real -index 0 $xy]
    return $xy2
 }
 proc DrawCurve {args} {
    global S

    set xy {}
    foreach pt [SortPoints] {                   ;# Flatten point list
        foreach {x y} $pt break
        lappend xy $x $y
    }

    switch $::S(type) {
        "line" { set xy [PolyLine::Linear      $xy $::S(precision)] }
        "cosi" { set xy [PolyLine::Cosine      $xy $::S(precision)] }
        "cubp" { set xy [PolyLine::Cubic       $xy $::S(precision)] }
        "cubs" { set xy [PolyLine::CubicSpline $xy $::S(precision)] }
        "catm" { set xy [PolyLine::CatmullRom  $xy $::S(precision)] }
        "bezi" { set xy [PolyLine::Bezier      $xy $::S(precision)] }
        "fbez" { set xy [PolyLine::FDBezier    $xy $::S(precision)] }
        "hbez" { set xy [PolyLine::HornBezier  $xy $::S(precision)] }
        "tbsp" { set xy [PolyLine::Hermite     $xy $::S(tension) $::S(bias) $::S(precision)] }
        "tcbs" { set xy [PolyLine::TCB         $xy $::S(tension) $::S(continuity) $::S(bias) $::S(precision)] }
        "card" { set xy [PolyLine::Cardinal    $xy $::S(a) $::S(precision)] }
        "kfs"  { set xy [PolyLine::KeyFrameTCB $xy $::S(tension) $::S(continuity) $::S(bias) $::S(precision)] }
    }
    .c delete cubic
    if {$xy == {}} return
    .c create line $xy -tag cubic -width $::S(w)
    .c lower cubic
 }
 proc About {} {
    set msg "Tension Splines
 by Greg Blair September 2003
 [email protected]
 GUI derived from
 Cubic Splines
 by Keith Vetter, March 2003"
    tk_messageBox -title About -message $msg
 }
 ################################################################
 ################################################################
 DoDisplay
 set S(type) line
 Reset
 DrawCurve

 ##In the above code, the part of [PolyLine::CubicSpline] that comes after Now generate
 ##the point list does such a conversion, but it is rather crude. Approximation with
 ##polygons is not bad as such (indeed, it is what Postscript does internally when it
 ##flattens paths), but in general one would want the number of line segments to be
 ##adjusted to the actual curve, so that the "precision" really would give a bound on
 ##the approximation error. I don't know how one would achieve this effectively, but
 ##there are probably good algorithms in the literature.

 }

uniquename 2013aug18

Since images at 'external sites' (such as www3.sympatico.ca above) tend to go dead, here is a screenshot image 'locally stored' on the disk drives of the server of this wiki. This image was created from the code above, but to fit the GUI on a netbook screen (max of 600 pixels high), the size 14 font was changed to 7 --- and a '-pady 0' parameter was added to 'button' and 'radiobutton' definitions. Also the '-width' (height) of the 'scale' widgets was reduced to 6 pixels.

blair_TensionSplines_wiki10530_cosine_screenshot_838x546.jpg