[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 '''Cubic::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 [[Cubic::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 {}
 
 ##+##########################################################################
 # 
 # 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 Cubic::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)}]
        }
        Cubic::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 Cubic::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)}]
    }
 }
 
----
Now to demonstrate and test the code you can use the following:

----
 ################################################################
 ################################################################
 #
 # Test harness and demo code
 #
 
 package require Tk
 
 set S(title) "Cubic Spline"
 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 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
    pack .c   -in .screen -side top    -fill both -expand 1
    bind all <Alt-c> {console show}
    bind .c <Configure> {ReCenter %W %h %w}
    .c bind p <B1-Motion> [list MouseMove %x %y]
 
    DoCtrlFrame
    AddCtrlPoint $INITPOINTS
    update
 }
 proc DoCtrlFrame {} {
    button .add  -text "Add Point"   -command AddCtrlPoint -bd 4
    button .dele -text "Delete Point" -command DeleteCtrlPoint -bd 4
    button .clear -text "Clear Points"  -command ClearCtrlPoint
    .add configure  -font "[font actual [.add cget -font]] -weight bold"
    .dele configure -font [.add cget -font]
    .clear configure  -font [.add cget -font]
 
    scale .prec -orient h -variable S(precision) -font [.add cget -font] \
        -label Precision: -relief ridge -from 1 -to 20 -command DrawCurve
    button .about -text About -command About
    
    grid .add   -in .ctrl -row 1 -sticky ew
    grid .dele  -in .ctrl -row 2 -sticky ew
    grid .clear -in .ctrl -row 3 -sticky ew -pady 10
    grid .prec  -in .ctrl -row 4 -sticky ew -pady 30
    grid rowconfigure .ctrl 50 -weight 1
 
    grid .about   -in .ctrl -row 100 -sticky ew
 }
 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]
 }
 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 ClearCtrlPoints {} {
    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
    }
    
    set xy [Cubic::CubicSpline $xy $::S(precision)]
    .c delete cubic
    if {$xy == {}} return
    .c create line $xy -tag cubic -width $::S(w)
    .c lower cubic
 }
 proc About {} {
    set msg "Cubic Splines\nby Keith Vetter, March 2003"
    tk_messageBox -title About -message $msg
 }
 ################################################################
 ################################################################
 DoDisplay
 DrawCurve

----
[Category Graphics] | [Category Mathematics]