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


# Synopsis:
#
# spline::getSpline <vertex list> step ?method?
#
# Args:   vertex list is of the form {x0 y0 x1 y1 ... xn yn} 
#         step is the desired x-sampling step of the spline function
#         method can be cmr (catmull-rom splines, passes through the 
#           vertices) or bsp (b-splines, does not pass throug 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 ?-method method? ?-prune pruneflag?
#
# 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 
#       If pruneflag is set to true (default) the vertexlist will be pruned so
#         vertices earlier in the list than those needed to calculate the 
#         point at "x" are removed. This will significantly speed up 
#         subsequent calls to getPoint for long vertex lists (assuming that x 
#         increases for each subsequent call)
#       method is either cmr or bsp
#       
# Result: returns the y-value of the spline function at point x.
#         
# Gui: the last section of the file contains a gui for testing purposes. 
#      To run it, comment out the return statement.
#      In the GUI, right button adds points, left button moves points. 
#      note: new points must be added to the right of the last point!

namespace eval spline {
 variable a
 set a(HALFPI) [expr 2*atan(1)]
} 


# blend function for b-spline
proc spline::b-bsp {i t} {
 switch -- $i {
  -2 {expr {(((-$t+3)*$t-3)*$t+1)/6.}}
  -1 {expr {(((3*$t-6)*$t)*$t+4)/6.}}
  0  {expr {(((-3*$t+3)*$t+3)*$t+1)/6.}}
  1  {expr {$t*$t*$t/6.}}
 }
}

# blend function for catmull-rom
proc spline::b-cmr {i t} {
 switch -- $i {
  -2 {expr {((-$t+2)*$t-1)*$t/2.}}
  -1 {expr {(((3*$t-5)*$t)*$t+2)/2.}}
  0  {expr {((-3*$t+4)*$t+1)*$t/2.}}
  1  {expr {($t-1)*$t*$t/2.}}
 }
}

proc spline::getPoint {vertexvar x args} {
 upvar $vertexvar vertices
 # puts [info level 0]
 # 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

 array set a {
  -method cmr
  -prune true
 }
 array set a $args

 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
 }
 if [string is true $a(-prune)] {
  if {$i>0} {set vertices [lrange $vertices $i end]}
 }
 set t [expr {($x-$vx(-1))/($vx(0)-$vx(-1))}]
 set y 0
 for {set j -2} {$j<=1} {incr j} {
  set b [b-$a(-method) $j $t]
  set y [expr {$y+$b*$vy($j)}]
 }
 return $y
}

# b-spline or catmull-rom      
proc spline::getSpline {vertices step {method cmr}} {

 set x0 [expr {$step*ceil(1.0*[lindex $vertices 2]/$step)}]
 set x1 [expr {$step*floor(1.0*[lindex $vertices end-3]/$step)}]
 for {set x $x0} {$x<$x1} {set x [expr {$x+$step}]} {
  lappend result $x [getPoint vertices $x -method $method -prune 1]
 }
 return $result
}

# experimental home-cooked splines
proc spline::generate2 {vertices step} {
 puts [info level 0]
 variable a
 set np [expr [llength $vertices]/2]

 set i 0
 foreach {x y} $vertices {
  set v(x,$i) $x
  set v(y,$i) $y
  incr i
 }
 for {set i 1;set ip 0; set in 2} {$i<$np-1} {incr i; incr ip; incr in} {
  set v(k,$i) [expr {($v(y,$in)-$v(y,$ip))/($v(x,$in)-$v(x,$ip))}]
 }
 
 set x0 [expr {$step*ceil(1.0*$v(x,1)/$step)}]
 set x1 [expr {$step*floor(1.0*$v(x,[expr $np-2])/$step)}]
 #puts x0=$x0
 #puts x1=$x1
 set seg 0
 for {set x $x0} {$x<$x1} {set x [expr {$x+$step}]} {
  #puts x=$x
  while {$v(x,$seg)<=$x} {incr seg} ;#find next vertex above $x
  set pseg [expr {$seg-1}]
  set prevx $v(x,$pseg)
  set nextx $v(x,$seg)
  set prevy $v(y,$pseg)
  set nexty $v(y,$seg)
  set prevk $v(k,$pseg)
  set nextk $v(k,$seg)
  set t [expr {($x-$prevx)/($nextx-$prevx)}]
#  set w0 [expr {cos($t*$a(HALFPI))}]
#  set w1 [expr {sin($t*$a(HALFPI))}]
  set ww [expr {cos($t*$a(HALFPI))}]
  set w0 [expr {$ww*$ww}]
  set w1 [expr {1-$w0}]
#  if {$t<.5} {set w0 1;set w1 0} else {set w0 0; set w1 1}
  set w0 [expr 1-$t]
  set w1 $t

  puts t=$t:w0=$w0,w1=$w1
  set y [expr {$w0*($prevy+$prevk*($x-$prevx)) + $w1*($nexty+$nextk*($x-$nextx))}]
  lappend result $x $y
 }
 return $result
}

return

#------------------------------------------------------------------------------
# test code bewlow
namespace eval spline::test {}

proc spline::test::movePt {c id x y} {
 $c coords $id $x $y $x $y
 updateSpline $c
}
proc spline::test::addPt {c x y} {
 variable 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 spline::test::updateSpline {c} {
 variable points
 variable a
 foreach id $points {
  foreach {x y x y} [$c coords $id] break
  lappend vertices $x $y
 }
 if {[llength $vertices] < 8} return
# $c coords jb [generate2 $vertices 2]
 puts [time {
  $c coords cmr [spline::getSpline $vertices 2 cmr]
 }]
 $c coords bsp [spline::getSpline $vertices 2 bsp]
}

proc spline::test::gui {vertices step} {
 variable points
 catch {destroy .c}
 pack [canvas .c -width 800 -bg black] -expand 1 -fill both
 
 
 set points ""
 .c create line 0 0 0 0 -fill yellow -tags [list spline cmr]
 .c create line 0 0 0 0 -fill green -tags [list spline bsp]
 .c create line 0 0 0 0 -fill blue -tags [list spline jb]
 foreach {x y} $vertices {
  addPt .c $x $y
 }
 bind .c <ButtonRelease-3> [namespace code [list addPt .c %x %y]]
 updateSpline .c
}

spline::test::gui [list 0 100 50 200 100 50 200 70 300 100] 1]