Version 10 of Tension Splines
Updated 2004-04-17 00:12:302003-12-08 - new page
2003-12-17 - added bezier curves
---
http://mini.net/files/CubicSplinesWithTension.png (How do I upload an image to http://mini.net/files ? ftp+gftp+sftp+rcp all fail! Wiki entry http://wiki.tcl.tk/861 discusses the use of images in wiki pages, but nowhere documents how to upload an image to the wiki site.
GPS: You can't upload images to the files directory unless you are Richard Suchenwirth or have other connections. I think this was done to prevent inappropriate content from showing up. :) )
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!
# 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 http://mini.net/tcl/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\nby Greg Blair September 2003\[email protected]\nGUI derived from\nCubic Splines\nby 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.
}