Version 26 of Cubic Splines

Updated 2005-05-08 16:39:44

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] ...

Greg Blair extended the code here in Tension Splines by adding several new splines with more controls.


 ##+##########################################################################
 #
 # 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

Lars H: A very nice little application, but not at all what my entry (which is now TIP #168) on the Tk 9.0 WishList was about. The problem is not to compute the control points of the splines, but to make the canvas understand something like a Postscript path (list of knot and control point coordinates for cubic Bezier polynomials). In want of support for this in the canvas widget, one could make do with some Tcl procedures for translating such a cubic path to something that the canvas can understand. I suppose this is the kind of thing that should go into tklib.

In the above code, the part of [Cubic::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.

AK: Just wondering, could MetaFont contains such algorithms ? It is a companion application to TeX, and Knuth used it generate his Computer Modern Fonts.

Lars H: No, it doesn't. MetaFont goes directly from Bezier curve to discrete pixel edges; there is no intermediate polygon representation (in that case a good thing, because it would introduce additional errors).


I found Keith's Cubis splines very useful, but missed the ability to work with a free sequence of X values (not monotonic increasing). Using his code, this was easy to implement by first collecting the length of the line between the points and then splining the X and Y values separately over the length. When the resulting curves are combined, you get the ability to spline just about every sequence of values. I added some code from Jonas Beskow who implemented catmull-Rom splines in Tcl. The code below displays both, just pick what suits you best:

# Adapted from http://www.cse.unsw.edu.au/~lambert/splines/source.html # (c) 2002 Jonas Beskow

# # Synopsis: # # spline::getSpline <vertex list> step # # Args: vertex list is of the form {x0 y0 x1 y1 ... xn yn} # step is the desired x-sampling step of the spline function # (catmull-rom splines, passes through the vertices) # # Result: spline points formatted as a list of the form {x0 y0 x1 y1 x2 y2} # the spline will start at x1 and end at x{n-1} # # spline::getPoint vertexlistvar x # # Args: vertexlistvar is a variable containing a list of vertices. # x is an x-coordinate lying between x1 and x{n-1} of the vertex list # The vertexlist will be pruned so vertices earlier in the list that # those needed to calculate the point at "x" are removed. # This will significantly speed up subsequent calls to getPoint for # long vertex lists. # # Result: returns the y-value of the spline function at point x. # # Anpassung an reellen Bedarf: Nur noch Catmull-Splines, und diese werden # nicht entlang der X-Achse, sondern der L‰ngenkoordinate des Splines auf- # gezogen. Dadurch entf‰llt die unpraktische Einschr‰nkung stetig steigen- # der X-Koordinaten! Allerdings geht dann alles nur noch halb so schnell. # # Peter K‰mpf 2005 # namespace eval spline {

    proc getPoint { vertexvar x } {
        upvar $vertexvar vertices

# # puts "vertexlist length: llength $vertices" # Scan through vertex list to find four consecutive points # where the input x lies between the middle two # After the scan the list is pruned to shorten the search next time #

        if {$x < [lindex $vertices 2] || $x > [lindex $vertices end-3]} {
            error "x-value out of range ([lindex $vertices 2] <= x <= [lindex $vertices end-3])" 
        }

#

        set i 0
        foreach \
          {vx(-2) vy(-2)} [lrange $vertices 0 end-6] \
          {vx(-1) vy(-1)} [lrange $vertices 2 end-4] \
          {vx(0)  vy(0)}  [lrange $vertices 4 end-2] \
          {vx(1)  vy(1)}  [lrange $vertices 6 end] {
            if {$vx(-1) <= $x && $vx(0) > $x} break
            incr i 2
        }

# # Pruning #

        if {$i > 0} {set vertices [lrange $vertices $i end]}

# # blend function for catmull-rom #

        set t [expr {($x-$vx(-1)) / ($vx(0)-$vx(-1))}]
        set y 0
        for {set j -2} {$j <= 1} {incr j} {
            switch -- $j {
                -2 {
                    set b [expr {((-$t+2)*$t-1)*$t/2.0}]
                }
                -1 {
                    set b [expr {(((3*$t-5)*$t)*$t+2)/2.0}]
                }
                0  {
                    set b [expr {((-3*$t+4)*$t+1)*$t/2.0}]
                }
                1  {
                    set b [expr {($t-1)*$t*$t/2.0}]
                }
            }
            set y [expr {$y + $b*$vy($j)}]
        }
        return $y
    }

# # Aufbau der Kurve aus der Vertex-Liste. Trick: Um die Anfangs- und Endsteigung # zu erhalten, wird die Liste am Anfang und Ende um einen extrapolierten Punkt # bereichert. #

    proc getSpline { vertices step } {
        set x0 [expr {1.1 * [lindex $vertices 0] - 0.1 * [lindex $vertices 2]}]
        set y0 [expr {1.1 * [lindex $vertices 1] - 0.1 * [lindex $vertices 3]}]
        set x1 [expr {1.1 * [lindex $vertices end-1] - 0.1 * [lindex $vertices end-3]}]
        set y1 [expr {1.1 * [lindex $vertices end]   - 0.1 * [lindex $vertices end-2]}]
        set vertices [linsert $vertices 0 $x0 $y0]
        lappend vertices $x1 $y1
        set result [list ]

# # Aufspalten in zwei Listen, die sich entlang der L‰ngskoordinate orientieren #

        set Laenge 0.0
        set xListe [list 0.0 $x0]
        set yListe [list 0.0 $y0]
        foreach {x1 y1} [lrange $vertices 2 end] {
            set x [expr {$x1 - $x0}]
            set y [expr {$y1 - $y0}]
            set Laenge [expr {$Laenge + sqrt($x*$x + $y*$y)}]
            lappend xListe $Laenge $x1
            lappend yListe $Laenge $y1
            set x0 $x1; set y0 $y1
        }

#

        set x0 [lindex $xListe 2]
        set x1 [lindex $xListe end-3]
        for {set x $x0} {$x < $x1} {set x [expr {$x+$step}]} {
            lappend result [getPoint xListe $x] [getPoint yListe $x]
        }
        return $result
    }

# # Aufspalten in zwei Listen, die sich entlang der L‰ngskoordinate orientieren # (c) Peter K‰mpf 2005 # # This is what you call to get the cubic spline coordinates with free # sequence of X values. # # Input: # vertices A list of x y pairs of the points to be connected # Precision Number of steps between two consecutive points # # Output: # Returns a list of all x y pairs which describe the spline. The # number of value pairs = (number of points - 1) x Precision + 1 # # Needs: CubicSpline, SolveTridiag #

    proc Cubic { vertices Precision } {
        set result [list ]
        set Laenge 0.0
        set x0     [lindex $vertices 0]
        set y0     [lindex $vertices 1]
        set xListe [list 0.0 $x0]
        set yListe [list 0.0 $y0]
        foreach {x1 y1} [lrange $vertices 2 end] {
            set x [expr {$x1 - $x0}]
            set y [expr {$y1 - $y0}]
            set Laenge [expr {$Laenge + sqrt($x*$x + $y*$y)}]
            lappend xListe $Laenge $x1
            lappend yListe $Laenge $y1
            set x0 $x1; set y0 $y1
        }

#

        set X [CubicSpline $xListe $Precision]
        set Y [CubicSpline $yListe $Precision]
        foreach {xLaenge xp} $X {yLaenge yp} $Y {
            lappend result $xp $yp
        }
        return $result
    }

##+########################################################################## # # CubicSpline - returns the x,y coordinates of the cubic spline using # XY control points. # vertices => {{x0 y0} {x1 y1} .... {xn yn}} # # XY points MUST BE SORTED by increasing x # # by Keith Vetter, March 2003 #

    proc CubicSpline { vertices {PRECISION 10} } {

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

        set idx 0
        foreach {x y} $vertices {
            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 {$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)}]
            }
            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 b1:n #

    proc 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)}]
        }
    }

}

#------------------------------------------------------------------------------ # test code #

proc main {} {

    global points
    global Step
    global Precision
    global cattime
    global contime
    set vertices [list 10 100  50 190 100  50 180 150 \
                      280 200 330  80 240 100 330 180 \
                      420 300 320 400 180 280]
    set Step       3
    set Precision 10
    font create big -family Helvetica -size 18 -weight bold

# # Open window with canvas on the left #

    wm title . "Splinetest"
    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 300 -width 500
    pack   .c -in .screen -side top -fill both -expand 1
    bind all <Alt-c> {console show}

# # and control widgets on the right side #

    label .cat -text "Catmull-Rom Splines" -fg blue -font big
    scale .cats -orient h -variable Step -command {updateSpline .c} \
        -label Stepsize: -relief ridge -from 1 -to 50
    label .catt -text "Time to draw"
    entry .cate -textvariable cattime -width 24 -bd 0 -justify right

#

    label .con -text "Cubic Splines" -fg green -font big
    scale .cons -orient h -variable Precision -command {updateSpline .c} \
        -label Precision: -relief ridge -from 1 -to 20
    label .cont -text "Time to draw"
    entry .cone -textvariable contime -width 24 -bd 0 -justify right

#

    message .m -width 200 \
        -text "Click and drag to move points;\
               Double-click to add a point;\
               Shift-Double-click to remove a point."
    button .about -text About -command About -width 12

# # Plot both lines #

    set points ""
    set contime [time { .c create line [spline::Cubic $vertices $Precision] \
        -width 2 -fill green -tags [list spline cube]}]
    set cattime [time { .c create line [spline::getSpline $vertices $Step] \
        -width 2 -fill blue -tags [list spline cmr]}]
    foreach {x y} $vertices {
        addPt .c $x $y
    }
    bind .c <Double-ButtonPress-1> [namespace code [list addPt .c %x %y]]
    bind .c <Shift-Double-ButtonPress-1> [namespace code [list rmPt .c %x %y]]

# # Position the widgets #

    grid .cat   -in .ctrl -row 1 -sticky ew -pady 10
    grid .cats  -in .ctrl -row 2 -sticky ew -pady 10
    grid .catt  -in .ctrl -row 3 -sticky ew
    grid .cate  -in .ctrl -row 4 -sticky ew -padx 10

#

    grid .con   -in .ctrl -row 6 -sticky ew -pady 10
    grid .cons  -in .ctrl -row 7 -sticky ew -pady 10
    grid .cont  -in .ctrl -row 8 -sticky ew
    grid .cone  -in .ctrl -row 9 -sticky ew -padx 10

#

    grid .m     -in .ctrl -row 11 -sticky ew -pady 30
    grid rowconfigure .ctrl {1 2 3 4 6 7 8 9 11 13} -weight 1
    grid rowconfigure .ctrl {0 5 10 12}  -weight 1 -minsize 10

#

    grid .about -in .ctrl -row 13 -pady 10
    update

}

proc movePt { c id x y } {

    $c coords $id $x $y $x $y
    updateSpline $c

}

proc addPt { c x y } {

    global points
    set id [.c create oval $x $y $x $y -width 6 -outline red]
    .c bind $id <B1-Motion> [namespace code [list movePt .c $id %x %y]]
    lappend points $id
    updateSpline $c

}

proc rmPt { c x y } {

    global points
    set Index 0

# puts "Punkt an $x $y entfernen"

    foreach id $points {
        foreach {xp yp xp yp} [$c coords $id] break
        if {[expr {abs($x - $xp)}] < 3 &&  [expr {abs($y - $yp)}] < 3} {
            set points [lreplace $points $Index $Index]
            $c delete $id
            break
        }
        incr Index
    }
    updateSpline $c

}

proc updateSpline { c {a 0} } {

    global points
    global Step
    global Precision
    global cattime
    global contime

#

    foreach id $points {
        foreach {x y x y} [$c coords $id] break
        lappend vertices $x $y
    }
    if {[llength $vertices] < 4} return
    set contime [time { $c coords cube [spline::Cubic $vertices $Precision] }]
    set cattime [time { $c coords cmr  [spline::getSpline $vertices $Step] }]

}

proc About {} {

    set msg "Cubic Splines by Keith Vetter, March 2003 \
        \nCatmull-Rom Splines by Jonas Beskow, 2002 \
        \nFree order of X points by Peter K‰mpf, 2005"
    tk_messageBox -title About -message $msg

}

# main


[ Category Graphics | Category Mathematics ]