Module Spline

Comtois · Post by **Comtois** » Sat May 16, 2015 5:39 pm

ModuleSpline.pb

DeclareModule Spline
  Declare AddPoint(Spline, x.f, y.f, z.f)
  Declare Clear(Spline)
  Declare Compute(Spline, t.f)
  Declare CountPoints(Spline)
  Declare Create()
  Declare Free(Spline)
  Declare.f PointX(Spline, index)
  Declare.f PointY(Spline, index)
  Declare.f PointZ(Spline, index)
  Declare.f X(Spline)
  Declare.f Y(Spline)
  Declare.f Z(Spline)
EndDeclareModule

Module Spline
  #Epsilon = 0.01  

  Structure Vector3
    x.f
    y.f
    z.f
  EndStructure
  
  Structure s_Spline
    AutoCalc.i
    Ret.Vector3
    Array Points.Vector3(0)
    Array Tangents.Vector3(0)
    Array Coeffs.f(3,3)
  EndStructure
  
  Declare recalcTangents(*This.s_Spline)
  Declare interpolate2(*This.s_Spline, fromIndex, t.f)
  
  Macro CopyVector3(a, b)
    a\x = b\x
    a\y = b\y
    a\z = b\z
  EndMacro

  Macro EqualVector3(a, b)
    a\x<=b\x+#Epsilon And a\x>=b\x-#Epsilon And a\y<=b\y+#Epsilon And a\y>=b\y-#Epsilon And a\z<=b\z+#Epsilon And a\z>=b\z-#Epsilon
  EndMacro  
  
  Procedure Create()
    Protected *Spline.s_Spline = AllocateMemory(SizeOf(s_Spline))
    If *Spline
      InitializeStructure(*Spline, s_Spline)
      
      ; Set up matrix Hermite polynomial
      With *Spline
        \Coeffs(0,0) = 2
        \Coeffs(0,1) = -2
        \Coeffs(0,2) = 1
        \Coeffs(0,3) = 1
        
        \Coeffs(1,0) = -3
        \Coeffs(1,1) = 3
        \Coeffs(1,2) = -2
        \Coeffs(1,3) = -1
        
        \Coeffs(2,0) = 0
        \Coeffs(2,1) = 0
        \Coeffs(2,2) = 1
        \Coeffs(2,3) = 0
        
        \Coeffs(3,0) = 1
        \Coeffs(3,1) = 0
        \Coeffs(3,2) = 0
        \Coeffs(3,3) = 0
        
        \AutoCalc = #True
      EndWith
      ProcedureReturn *Spline
    Else
      ProcedureReturn #False
    EndIf  
  EndProcedure
  
  Procedure Free(*This.s_Spline)
    If *This
      ClearStructure(*This, s_Spline)
      FreeMemory(*This)
    EndIf 
  EndProcedure
  
  Procedure Clear(*This.s_Spline)
    If *This
      Dim *This\Points.Vector3(0)
      Dim *This\Tangents.Vector3(0)
    EndIf
  EndProcedure
  
  Procedure concatenate(Array r.f(2), Array m.f(2), Array m2.f(2))
    r(0,0) = m(0,0) * m2(0,0) + m(0,1) * m2(1,0) + m(0,2) * m2(2,0) + m(0,3) * m2(3,0);
    r(0,1) = m(0,0) * m2(0,1) + m(0,1) * m2(1,1) + m(0,2) * m2(2,1) + m(0,3) * m2(3,1);
    r(0,2) = m(0,0) * m2(0,2) + m(0,1) * m2(1,2) + m(0,2) * m2(2,2) + m(0,3) * m2(3,2);
    r(0,3) = m(0,0) * m2(0,3) + m(0,1) * m2(1,3) + m(0,2) * m2(2,3) + m(0,3) * m2(3,3);
    
    r(1,0) = m(1,0) * m2(0,0) + m(1,1) * m2(1,0) + m(1,2) * m2(2,0) + m(1,3) * m2(3,0);
    r(1,1) = m(1,0) * m2(0,1) + m(1,1) * m2(1,1) + m(1,2) * m2(2,1) + m(1,3) * m2(3,1);
    r(1,2) = m(1,0) * m2(0,2) + m(1,1) * m2(1,2) + m(1,2) * m2(2,2) + m(1,3) * m2(3,2);
    r(1,3) = m(1,0) * m2(0,3) + m(1,1) * m2(1,3) + m(1,2) * m2(2,3) + m(1,3) * m2(3,3);
    
    r(2,0) = m(2,0) * m2(0,0) + m(2,1) * m2(1,0) + m(2,2) * m2(2,0) + m(2,3) * m2(3,0);
    r(2,1) = m(2,0) * m2(0,1) + m(2,1) * m2(1,1) + m(2,2) * m2(2,1) + m(2,3) * m2(3,1);
    r(2,2) = m(2,0) * m2(0,2) + m(2,1) * m2(1,2) + m(2,2) * m2(2,2) + m(2,3) * m2(3,2);
    r(2,3) = m(2,0) * m2(0,3) + m(2,1) * m2(1,3) + m(2,2) * m2(2,3) + m(2,3) * m2(3,3);
    
    r(3,0) = m(3,0) * m2(0,0) + m(3,1) * m2(1,0) + m(3,2) * m2(2,0) + m(3,3) * m2(3,0);
    r(3,1) = m(3,0) * m2(0,1) + m(3,1) * m2(1,1) + m(3,2) * m2(2,1) + m(3,3) * m2(3,1);
    r(3,2) = m(3,0) * m2(0,2) + m(3,1) * m2(1,2) + m(3,2) * m2(2,2) + m(3,3) * m2(3,2);
    r(3,3) = m(3,0) * m2(0,3) + m(3,1) * m2(1,3) + m(3,2) * m2(2,3) + m(3,3) * m2(3,3);
  EndProcedure
  
  Procedure Multiplication(Array r.f(1), Array v.f(1), Array m.f(2))
    r(0) = m(0,0) * v(0) + m(1,0) * v(1) + m(2,0) * v(2) + m(3,0) * v(3)
    r(1) = m(0,1) * v(0) + m(1,1) * v(1) + m(2,1) * v(2) + m(3,1) * v(3)
    r(2) = m(0,2) * v(0) + m(1,2) * v(1) + m(2,2) * v(2) + m(3,2) * v(3)
    r(3) = m(0,3) * v(0) + m(1,3) * v(1) + m(2,3) * v(2) + m(3,3) * v(3)
  EndProcedure
  
  Procedure addPoint(*This.s_Spline, x.f, y.f, z.f)
    Size = ArraySize(*This\Points())
    ReDim *This\Points(Size+1)
    *This\Points(Size)\x = x
    *This\Points(Size)\y = y
    *This\Points(Size)\z = z
    
    If (*This\AutoCalc)
      recalcTangents(*This)
    EndIf
    
  EndProcedure
  
  Procedure Compute(*This.s_Spline, t.f)
    
    ;Currently assumes points are evenly spaced, will cause velocity
    ;change where this is Not the Case
    ;TODO: base on arclength?
    
    
    ; Work out which segment this is in
    fSeg.f = t * (ArraySize(*This\Points()) - 1)
    segIdx = Int(fSeg)
    ; Apportion t
    t = fSeg - segIdx
    interpolate2(*This, segIdx, t)
    
  EndProcedure
  
  Procedure interpolate2(*This.s_Spline, fromIndex, t.f)
    
    ;// Bounds check
    ;assert (fromIndex < mPoints.size() && "fromIndex out of bounds")
    If ArraySize(*This\Points()) = 0 Or fromIndex >=ArraySize(*This\Points())
      ProcedureReturn
    EndIf
    
    If (fromIndex + 1) = ArraySize(*This\Points())
      
      ;// Duff request, cannot blend To nothing
      ;// Just Return source
      CopyVector3(*This\Ret, *This\Points(fromIndex))
      ProcedureReturn
      
    EndIf
    
    ;// Fast special cases
    If t = 0.0
      
      CopyVector3(*This\Ret, *This\Points(fromIndex))
      ProcedureReturn
      
    ElseIf t = 1.0
      CopyVector3(*This\Ret, *This\Points(fromIndex + 1))
      ProcedureReturn
    EndIf
    
    ;// Real interpolation
    ;// Form a vector of powers of t
    Define.f t2, t3
    t2 = t * t
    t3 = t2 * t
    Dim powers.f(3)
    powers(0) = t3
    powers(1) = t2
    powers(2) = t
    powers(3) = 1
    
    ;// Algorithm is ret = powers * Coeffs * Matrix4(point1, point2, tangent1, tangent2)
    Define.Vector3 point1, point2, tan1, tan2
    CopyVector3(point1, *This\Points(fromIndex))
    CopyVector3(point2, *This\Points(fromIndex+1))
    CopyVector3(tan1  , *This\Tangents(fromIndex))
    CopyVector3(tan2  , *This\Tangents(fromIndex+1))
    
    Dim pt.f(3,3)
    
    pt(0,0) = point1\x
    pt(0,1) = point1\y
    pt(0,2) = point1\z
    pt(0,3) = 1.0
    
    pt(1,0) = point2\x
    pt(1,1) = point2\y
    pt(1,2) = point2\z
    pt(1,3) = 1.0
    
    pt(2,0) = tan1\x
    pt(2,1) = tan1\y
    pt(2,2) = tan1\z
    pt(2,3) = 1.0
    
    pt(3,0) = tan2\x
    pt(3,1) = tan2\y
    pt(3,2) = tan2\z
    pt(3,3) = 1.0
    
    Dim ret.f(3)
    Dim r.f(3,3)
    ;ret = powers * Coeffs * pt;
    concatenate(r(), *This\Coeffs(), pt())
    Multiplication(ret(), powers(), r())
    
    *This\Ret\x = ret(0)   
    *This\Ret\y = ret(1)
    *This\Ret\z = ret(2)
    
  EndProcedure
  
  Procedure recalcTangents(*This.s_Spline)
    
    ;         // Catmull-Rom approach
    ;         //
    ;         // tangent[i] = 0.5 * (point[i+1] - point[i-1])
    ;         //
    ;         // Assume endpoint tangents are parallel With line With neighbour
    
    Define.i i, numPoints, isClosed
    
    numPoints = ArraySize(*This\Points())
    If numPoints < 2
      
      ; Can't do anything yet
      ProcedureReturn
    EndIf
    
    ; Closed Or open?
    
    If EqualVector3(*This\Points(0), *This\Points(numPoints-1))  
      
      isClosed = #True
      
    Else
      
      isClosed = #False
    EndIf
    
    Dim *This\Tangents.Vector3(numPoints)
        
    For i = 0 To numPoints-1
      
      If i =0
        
        ; Special Case start
        If isClosed
          ; Use numPoints-2 since numPoints-1 is the last point And == [0]
          *This\Tangents(i)\x = 0.5 * (*This\Points(1)\x - *This\Points(numPoints-2)\x)
          *This\Tangents(i)\y = 0.5 * (*This\Points(1)\y - *This\Points(numPoints-2)\y)
          *This\Tangents(i)\z = 0.5 * (*This\Points(1)\z - *This\Points(numPoints-2)\z)
        Else
          
          *This\Tangents(i)\x = 0.5 * (*This\Points(1)\x - *This\Points(0)\x)
          *This\Tangents(i)\y = 0.5 * (*This\Points(1)\y - *This\Points(0)\y)
          *This\Tangents(i)\z = 0.5 * (*This\Points(1)\z - *This\Points(0)\z)
        EndIf
        
      ElseIf i = numPoints-1
        
        ; Special Case End
        If isClosed
          
          ; Use same tangent As already calculated For [0]
          *This\Tangents(i)\x = *This\Tangents(0)\x
          *This\Tangents(i)\y = *This\Tangents(0)\y
          *This\Tangents(i)\z = *This\Tangents(0)\z
          
        Else
          
          *This\Tangents(i)\x = 0.5 * (*This\Points(i)\x - *This\Points(i-1)\x)
          *This\Tangents(i)\y = 0.5 * (*This\Points(i)\y - *This\Points(i-1)\y)
          *This\Tangents(i)\z = 0.5 * (*This\Points(i)\z - *This\Points(i-1)\z)
        EndIf
        
      Else
        *This\Tangents(i)\x = 0.5 * (*This\Points(i+1)\x - *This\Points(i-1)\x)
        *This\Tangents(i)\y = 0.5 * (*This\Points(i+1)\y - *This\Points(i-1)\y)
        *This\Tangents(i)\z = 0.5 * (*This\Points(i+1)\z - *This\Points(i-1)\z)
      EndIf
      
    Next
    
  EndProcedure
  
  Procedure CountPoints(*This.s_Spline)
    If *This
      ProcedureReturn ArraySize(*This\Points())
    EndIf  
  EndProcedure
  
  Procedure.f PointX(*This.s_Spline, index)
    If *This
      If index >= 0 And index < ArraySize(*This\Points())
        ProcedureReturn *This\Points(index)\x
      EndIf  
    EndIf  
  EndProcedure
  
  Procedure.f PointY(*This.s_Spline, index)
    If *This
      If index >= 0 And index < ArraySize(*This\Points())
        ProcedureReturn *This\Points(index)\y
      EndIf  
    EndIf  
  EndProcedure
  
  Procedure.f PointZ(*This.s_Spline, index)
    If *This
      If index >= 0 And index < ArraySize(*This\Points())
        ProcedureReturn *This\Points(index)\z
      EndIf  
    EndIf  
  EndProcedure
  
  Procedure.f X(*This.s_Spline)
    If *This
      ProcedureReturn *This\Ret\x
    EndIf  
  EndProcedure
  
  Procedure.f Y(*This.s_Spline)
    If *This
      ProcedureReturn *This\Ret\y
    EndIf  
  EndProcedure
  
  Procedure.f Z(*This.s_Spline)
    If *This
      ProcedureReturn *This\Ret\z
    EndIf  
  EndProcedure
EndModule

And an example (use LeftButton to create a path, RightButton to clear Spline)

Code: Select all

XIncludeFile "ModuleSpline.pb"

InitSprite()
InitKeyboard()
InitMouse()
ExamineDesktops()
Scx = DesktopWidth(0)
Scy = DesktopHeight(0)
OpenWindow(0,0,0,Scx, Scy, "Test Spline", #PB_Window_BorderLess)
OpenWindowedScreen(WindowID(0), 0,0,Scx,Scy)

MySpline = Spline::Create()

xd.f = 0
yd.f = 0
pas.f = 1
time.f = 0

CreateSprite(0, 64,64)
StartDrawing(SpriteOutput(0))
Line(0,0,1,64, RGB(255,255,0))
LineXY(0,16,64,32, RGB(255,255,0))
LineXY(0,64-16,64,32, RGB(255,255,0))
FillArea(32,32,RGB(255,255,0),RGB(85,85,0))
StopDrawing()

Repeat
  While WindowEvent()
  Wend
  
  ClearScreen(0)
  ExamineMouse()
  If MouseButton(1)
    If flag = 0
      Flag = 1
      MouseX.f = MouseX()
      MouseY.f = MouseY()
      Spline::AddPoint(MySpline, MouseX, 0, MouseY)
      If xd = 0 And yd = 0
        xd = MouseX
        yd = MouseY
      EndIf
    EndIf
  Else
    flag = 0
  EndIf
  If MouseButton(2)
    Spline::Clear(MySpline)
  EndIf
  StartDrawing(ScreenOutput())
  For t=0 To 300
    
    Spline::Compute(MySpline, t/300.0)
    x.f = Spline::X(MySpline)
    z.f = Spline::Z(MySpline)
    If t > 0
      LineXY(xd, yd, x, z, $AABBCC)
    EndIf
    xd = x
    yd = z   
  Next t
  
  For i=0 To Spline::CountPoints(MySpline)-1
    Circle(Spline::PointX(MySpline, i), Spline::PointZ(MySpline, i), 2, $FF)
  Next   
  
  LineXY(0, MouseY(), Scx, MouseY(), $AAAAAA)
  LineXY(MouseX(), 0, MouseX(), Scy, $AAAAAA)
  StopDrawing()
  
  ;-Sprite
  Spline::Compute(MySpline, time)
  x.f = Spline::X(MySpline)
  z.f = Spline::Z(MySpline)
  Spline::Compute(MySpline, time-0.001)
  x1.f = Spline::X(MySpline)
  z1.f = Spline::Z(MySpline)
  DisplayTransparentSprite(0,x-32, z-32)  
  
  DirectionX.f = x-x1
  DirectionZ.f = z-z1 
  Angle.f = Degree(ATan2(DirectionX, DirectionZ))
  
  RotateSprite(0,Angle,#PB_Absolute)
  time + pas * 0.02 / Spline::CountPoints(MySpline)
  
  If time > 1
    Time = 1
    pas = - pas
  ElseIf time < 0
    Time = 0
    pas = - pas
  EndIf
  
  ExamineKeyboard()
  FlipBuffers()
  
Until KeyboardPushed(#PB_Key_Escape)

Spline::Free(MySpline)

[EDIT]
Fixed Procedure Free()
New example (rotation sprite added)

[EDIT2]
Fixed Procedure recalcTangents() --> added EqualVector3()

Post by **Fred** » Mon May 18, 2015 1:15 pm

Nice code !

IdeasVacuum · Post by **IdeasVacuum** » Mon May 18, 2015 3:42 pm

...should be added to PB as a function

Comtois · Post by **Comtois** » Mon May 18, 2015 3:46 pm

IdeasVacuum wrote:...should be added to PB as a function

Done...

juror · Post by **juror** » Mon May 18, 2015 6:29 pm

Most impressive.

davido · Post by **davido** » Mon May 18, 2015 8:42 pm

Great. Thank you very much.

Little John · Post by **Little John** » Mon May 18, 2015 8:44 pm

Hi,

very cool!

However, ATM I've got two questions.

There are different kinds of splines. What kind of splines exactly does your code produce?
I'm not able to detect it from the code. But when I'd use the code in a program, I'd like to know what exactly my program does.

Can your module also be used with simple 2d graphics on an image gadget?

Comtois · Post by **Comtois** » Mon May 18, 2015 10:04 pm

Little John wrote: There are different kinds of splines. What kind of splines exactly does your code produce?

i have adapted Ogre3D's code for PureBasic :

Ogre3D wrote:Detailed Description

A very simple spline class which implements the Catmull-Rom class of splines.

Remarks
Splines are bendy lines. You define a series of points, and the spline forms a smoother line between the points to eliminate the sharp angles.

Catmull-Rom splines are a specialisation of the general Hermite spline. With a Hermite spline, you define the start and end point of the line, and 2 tangents, one at the start of the line and one at the end. The Catmull-Rom spline simplifies this by just asking you to define a series of points, and the tangents are created for you.

Little John wrote:Can your module also be used with simple 2d graphics on an image gadget?

You get a Vector3 use it like you want. In my example i use it to draw lines and circles on screen, and moving a sprite. You can use it to move an Entity 3D etc...

applePi · Post by **applePi** » Tue May 19, 2015 8:59 am

thanks Comtois,
i have tested it with y=sin(x) for period 0 to 2Pi with steps 0.1 , added inside the mouse click capture , and it has great visual appeal
click the mouse once to fill the array with y=sin(x) data

Code: Select all

XIncludeFile "ModuleSpline.pb"

InitSprite()
InitKeyboard()
InitMouse()
ExamineDesktops()
Scx = DesktopWidth(0)
Scy = DesktopHeight(0)
OpenWindow(0,0,0,Scx, Scy, "Test Spline", #PB_Window_BorderLess)
OpenWindowedScreen(WindowID(0), 0,0,Scx,Scy)

MySpline = Spline::Create()

xd.f = 0
yd.f = 0
pas.f = 1
time.f = 0

CreateSprite(0, 64,64)
StartDrawing(SpriteOutput(0))
Line(0,0,1,64, RGB(255,255,0))
LineXY(0,16,64,32, RGB(255,255,0))
LineXY(0,64-16,64,32, RGB(255,255,0))
FillArea(32,32,RGB(255,255,0),RGB(85,85,0))
StopDrawing()


Repeat
  While WindowEvent()
  Wend
  
  ClearScreen(0)
  ExamineMouse()
  If MouseButton(1)
    If flag = 0;;;;;;
      Flag = 1
      While xx.f<=2*#PI
        xx.f+0.1
        ;MouseX.f = MouseX()
        MouseX.f = xx*100 + 50
        ;MouseY.f = MouseY()
        MouseY.f = Sin(xx)* 200 + 300
      
      ;Debug MouseX: Debug MouseY:Debug "----------------------"
      Spline::AddPoint(MySpline, MouseX, 0, MouseY)
      If xd = 0 And yd = 0
        xd = MouseX
        yd = MouseY
      EndIf
      Wend
    EndIf;;;;;
  Else
    flag = 0
  EndIf
  If MouseButton(2)
    Spline::Clear(MySpline)
  EndIf
  StartDrawing(ScreenOutput())
  For t=0 To 300
    
    Spline::Compute(MySpline, t/300.0)
    x.f = Spline::X(MySpline)
    z.f = Spline::Z(MySpline)
    If t > 0
      LineXY(xd, yd, x, z, $AABBCC)
    EndIf
    xd = x
    yd = z   
  Next t
  
  For i=0 To Spline::CountPoints(MySpline)-1
    Circle(Spline::PointX(MySpline, i), Spline::PointZ(MySpline, i), 2, $FF)
  Next   
  
  LineXY(0, MouseY(), Scx, MouseY(), $AAAAAA)
  LineXY(MouseX(), 0, MouseX(), Scy, $AAAAAA)
  StopDrawing()
  
  ;-Sprite
  Spline::Compute(MySpline, time)
  x.f = Spline::X(MySpline)
  z.f = Spline::Z(MySpline)
  Spline::Compute(MySpline, time-0.001)
  x1.f = Spline::X(MySpline)
  z1.f = Spline::Z(MySpline)
  DisplayTransparentSprite(0,x-32, z-32)  
  
  DirectionX.f = x-x1
  DirectionZ.f = z-z1 
  Angle.f = Degree(ATan2(DirectionX, DirectionZ))
  
  RotateSprite(0,Angle,#PB_Absolute)
  time + pas * 0.2 / Spline::CountPoints(MySpline)
  
  If time > 1
    Time = 1
    pas = - pas
  ElseIf time < 0
    Time = 0
    pas = - pas
  EndIf
  
  ExamineKeyboard()
  FlipBuffers()
  
Until KeyboardPushed(#PB_Key_Escape)

Spline::Free(MySpline)

Little John · Post by **Little John** » Tue May 19, 2015 10:31 am

@Comtois:
Thank you for the explanations. And for the code, of course!

PureBasic Forums - English

Module Spline

Module Spline

Re: Module Spline

Re: Module Spline

Re: Module Spline

Re: Module Spline

Re: Module Spline

Re: Module Spline

Re: Module Spline

Re: Module Spline

Re: Module Spline