Code archives/Algorithms/Natural Cubic Spline

; Natural Cubic Spline 
; Based on a java-applet by Tim Lambert
;
; Drag the red dots around with the mouse

Graphics 800,600,0,2
SetBuffer BackBuffer()

Type Cubic
 Field a#
 Field b#
 Field c#
 Field d#
End Type

Const Points = 8 ; Number of points
Const Steps = 24 ; Number of substeps between points (Higher number = Smoother curve = Slower drawing)

;Arrays for calculations
Const n=Points-1  
Dim xp#(n)
Dim yp#(n)
Dim xc.Cubic(n)
Dim yc.Cubic(n)
Dim gamma#(n+1)
Dim delta#(n+1)
Dim D#(n+1)

;Place the points
For i=0 To n
 xp#(i)=i*(GraphicsWidth()/Points+1)+((GraphicsWidth()/Points+1)/2)
 yp#(i)=GraphicsHeight()/2
Next

;Init the CubicsArrays
InitCubics()

CalcCubics()
Cls
DrawSpline()
Flip

Dragging=False

While Not KeyHit(1)

 If MouseDown(1) And Dragging=False Then
  mx=MouseX():my=MouseY()
  For i=0 To n
   px=xp(i):py=yp(i)
   If mx>px-3 And mx<px+3 And my>py-3 And my<py+3 Then
    DragPoint=i
	Dragging=True
    Exit
   EndIf
  Next
 EndIf 

 If MouseDown(1)=False And Dragging=True Then
  Dragging=False
 EndIf

 If Dragging Then
  xp#(DragPoint)=MouseX()
  yp#(DragPoint)=MouseY()
  CalcCubics()
  Cls
  DrawSpline()	
  Flip
 EndIf 

Wend
End


Function DrawSpline()
 x#=CubicEval(xc(0),0.0)
 y#=CubicEval(yc(0),0.0)

  Color 255,255,255

For i=0 To n-1
 For j=1 To Steps
  oldx#=x#
  oldy#=y#
  stp#=Float(j)/Float(Steps)
  x#=CubicEval(xc(i),stp#)
  y#=CubicEval(yc(i),stp#)
  Line oldx#,oldy#,x#,y#
 Next
Next

  Color 255,0,0 
  For i=0 To n
   Oval xp#(i)-3,yp#(i)-3,6,6,True
  Next
End Function


Function CalcCubics()

gamma(0) = 1.0/2.0
For i=1 To n-1
 gamma(i) = 1/(4-gamma(i-1))
Next
gamma(n) = 1/(2-gamma(n-1))

; First the X-points

delta(0) = 3*(xp(1)-xp(0))*gamma(0)
For i=1 To n-1
 delta(i) = (3*(xp(i+1)-xp(i-1))-delta(i-1))*gamma(i)
Next
delta(n)=(3*(xp(n)-xp(n-1))-delta(n-1))*gamma(n)

D(n)=delta(n)
For i=n-1 To 0 Step -1
 D(i) = delta(i) - gamma(i)*D(i+1)
Next

For i=0 To n-1
 xc(i)\a# = xp(i)
 xc(i)\b# = D(i)
 xc(i)\c# = 3 * (xp(i+1)-xp(i)) - 2*D(i) - D(i+1)
 xc(i)\d# = 2 * (xp(i)-xp(i+1)) + D(i) + D(i+1)
Next

; Then the Y-points

delta(0) = 3*(yp(1)-yp(0))*gamma(0)
For i=1 To n-1
 delta(i) = (3*(yp(i+1)-yp(i-1))-delta(i-1))*gamma(i)
Next
delta(n)=(3*(yp(n)-yp(n-1))-delta(n-1))*gamma(n)

D(n)=delta(n)
For i=n-1 To 0 Step -1
 D(i) = delta(i) - gamma(i)*D(i+1)
Next

For i=0 To n-1
 yc(i)\a# = yp(i)
 yc(i)\b# = D(i)
 yc(i)\c# = 3 * (yp(i+1)-yp(i)) - 2*D(i) - D(i+1)
 yc(i)\d# = 2 * (yp(i)-yp(i+1)) + D(i) + D(i+1)
Next

End Function

Function CubicEval#(Cubic.Cubic,u#)
 ta#=Cubic\a#
 tb#=Cubic\b#
 tc#=Cubic\c#
 td#=Cubic\d#	

 Return (((td#*u#)+tc#)*u#+tb#)*u#+ta#
End Function

Function InitCubics()
 For i=0 To n
  xc.Cubic(i) = New Cubic
  yc.Cubic(i) = New Cubic
 Next
End Function