Interpolation cubique (courbe bézier...)

Vous débutez et vous avez besoin d'aide ? N'hésitez pas à poser vos questions
comtois
Messages : 5186
Inscription : mer. 21/janv./2004 17:48
Contact :

Re: Interpolation cubique (courbe bézier...)

Message par comtois »

Une petite modif pour définir les points avec la souris.

Code : Tout sélectionner

; This source file is part of OGRE
;     (Object-oriented Graphics Rendering Engine)
; For the latest info, see http://www.ogre3d.org/
;
; Copyright (c) 2000-2009 Torus Knot Software Ltd
InitSprite()
InitKeyboard()
InitMouse()
ExamineDesktops()
Scx = DesktopWidth(0)
Scy = DesktopHeight(0)
OpenScreen(Scx, Scy, DesktopDepth(0), "spline")

Structure Vector3
  x.f
  y.f
  z.f
EndStructure

Macro CopyVector3(a, b)
  a\x = b\x
  a\y = b\y
  a\z = b\z
EndMacro

Global Dim mCoeffs.f(3,3)
Global mAutoCalc = #True
Global NewList lPoints.Vector3()
Global Dim mPoints.Vector3(0)
Global Dim mTangents.Vector3(0)
Declare recalcTangents()
Declare interpolate2(*pos.Vector3, fromIndex, t.f)

; Set up matrix
; Hermite polynomial
mCoeffs(0,0) = 2
mCoeffs(0,1) = -2
mCoeffs(0,2) = 1
mCoeffs(0,3) = 1

mCoeffs(1,0) = -3
mCoeffs(1,1) = 3
mCoeffs(1,2) = -2
mCoeffs(1,3) = -1

mCoeffs(2,0) = 0
mCoeffs(2,1) = 0
mCoeffs(2,2) = 1
mCoeffs(2,3) = 0

mCoeffs(3,0) = 1
mCoeffs(3,1) = 0
mCoeffs(3,2) = 0
mCoeffs(3,3) = 0

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(*p.Vector3)
  AddElement(lPoints())
  lPoints()\x = *p\x
  lPoints()\y = *p\y
  lPoints()\z = *p\z
  
  Dim mPoints(ListSize(lPoints()))
  i = 0
  ForEach lPoints()
    mPoints(i) = lPoints()
    i + 1
    
  Next
  If (mAutoCalc)
    recalcTangents()
  EndIf
  
EndProcedure

Procedure interpolate(*pos.Vector3, 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 * (ListSize(lPoints()) - 1)
  segIdx = Int(fSeg)
  ; Apportion t
  t = fSeg - segIdx
  interpolate2(*pos, segIdx, t)
  
EndProcedure

Procedure interpolate2(*pos.Vector3, fromIndex, t.f)
  
  ;// Bounds check
  ;assert (fromIndex < mPoints.size() && "fromIndex out of bounds")
  If ListSize(lPoints()) = 0 Or fromIndex >= ListSize(lPoints())
    ProcedureReturn 
  EndIf
  
  If (fromIndex + 1) = ListSize(lPoints())
    
    ;// Duff request, cannot blend To nothing
    ;// Just Return source
    CopyVector3(*pos, mPoints(fromIndex))
    ProcedureReturn
    
  EndIf
  
  ;// Fast special cases
  If t = 0.0
    
    CopyVector3(*pos, mPoints(fromIndex))
    ProcedureReturn
    
  ElseIf t = 1.0
    CopyVector3(*pos, mPoints(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 * mCoeffs * Matrix4(point1, point2, tangent1, tangent2)
  Define.Vector3 point1, point2, tan1, tan2
  CopyVector3(point1, mPoints(fromIndex))
  CopyVector3(point2, mPoints(fromIndex+1))
  CopyVector3(tan1  , mTangents(fromIndex))
  CopyVector3(tan2  , mTangents(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 * mCoeffs * pt;
  concatenate(r(), mCoeffs(), pt())
  Multiplication(ret(), powers(), r())
  
  *pos\x = ret(0)   
  *pos\y = ret(1)
  *pos\z = ret(2)
  ;ProcedureReturn Vector3(ret.x, ret.y, ret.z);
  
EndProcedure

Procedure recalcTangents()
  
  ;         // 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 = ListSize(lPoints())
  If numPoints < 2
    
    ; Can't do anything yet
    ProcedureReturn
  EndIf
  
  ; Closed Or open?
  
  If mPoints(0) = mPoints(numPoints-1)
    
    isClosed = #True
    
  Else
    
    isClosed = #False
  EndIf
  
  Dim mTangents.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]
        mTangents(i)\x = 0.5 * (mPoints(1)\x - mPoints(numPoints-2)\x)
        mTangents(i)\y = 0.5 * (mPoints(1)\y - mPoints(numPoints-2)\y)
        mTangents(i)\z = 0.5 * (mPoints(1)\z - mPoints(numPoints-2)\z)
      Else
        
        mTangents(i)\x = 0.5 * (mPoints(1)\x - mPoints(0)\x)
        mTangents(i)\y = 0.5 * (mPoints(1)\y - mPoints(0)\y)
        mTangents(i)\z = 0.5 * (mPoints(1)\z - mPoints(0)\z)
      EndIf
      
    ElseIf i = numPoints-1
      
      ; Special Case End
      If isClosed
        
        ; Use same tangent As already calculated For [0]
        mTangents(i)\x = mTangents(0)\x
        mTangents(i)\y = mTangents(0)\y
        mTangents(i)\z = mTangents(0)\z
        
      Else
        
        mTangents(i)\x = 0.5 * (mPoints(i)\x - mPoints(i-1)\x)
        mTangents(i)\y = 0.5 * (mPoints(i)\y - mPoints(i-1)\y)
        mTangents(i)\z = 0.5 * (mPoints(i)\z - mPoints(i-1)\z)
      EndIf
      
    Else
      mTangents(i)\x = 0.5 * (mPoints(i+1)\x - mPoints(i-1)\x)
      mTangents(i)\y = 0.5 * (mPoints(i+1)\y - mPoints(i-1)\y)
      mTangents(i)\z = 0.5 * (mPoints(i+1)\z - mPoints(i-1)\z)
    EndIf
    
  Next
  
EndProcedure

;- Ajoute quelques points et ensuite on trace la spline
p.Vector3
pos.Vector3

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

Repeat
  ClearScreen(0)
  ExamineMouse()
  If MouseButton(1)
    If flag = 0
      Flag = 1
      p\x = MouseX()
      p\z = MouseY()
      addPoint(@p)
      If xd = 0 And yd = 0
        xd = p\x
        yd = p\z
      EndIf
    EndIf
  Else
    flag = 0
  EndIf 
  StartDrawing(ScreenOutput())
  For t=0 To 300
    
    interpolate(@pos, t/300.0)
    If t > 0
      LineXY(xd, yd, pos\x, pos\z, $AABBCC)
    EndIf
    xd = pos\x
    yd = pos\z   
  Next t
  interpolate(@pos, time)
  Circle(pos\x, pos\z, 10, $FFAAAA)
  
   time + pas * 0.02 / ListSize(lPoints())
  
  If time > 1
    Time = 1
    pas = - pas
  ElseIf time < 0
    Time = 0
    pas = - pas
  EndIf
  
  For i=0 To ListSize(lPoints())-1
    Circle(mPoints(i)\x, mPoints(i)\z, 2, $FF)
  Next   
  
  LineXY(0, MouseY(), Scx, MouseY(), $AAAAAA)
  LineXY(MouseX(), 0, MouseX(), Scy, $AAAAAA)
  StopDrawing()
  
  ExamineKeyboard()
  FlipBuffers()
Until KeyboardPushed(#PB_Key_Escape)
Dernière modification par comtois le jeu. 26/juil./2012 10:41, modifié 2 fois.
http://purebasic.developpez.com/
Je ne réponds à aucune question technique en PV, utilisez le forum, il est fait pour ça, et la réponse peut profiter à tous.
Mesa
Messages : 1126
Inscription : mer. 14/sept./2011 16:59

Re: Interpolation cubique (courbe bézier...)

Message par Mesa »

J'ai une erreur, index de tableau hors limite ligne 159.

Mesa.
comtois
Messages : 5186
Inscription : mer. 21/janv./2004 17:48
Contact :

Re: Interpolation cubique (courbe bézier...)

Message par comtois »

C'est exact, le problème apparait avec le débogueur activé. J'avais fait la correction chez moi et pas ici ! C'est corrigé, tu peux tester à nouveau :)
http://purebasic.developpez.com/
Je ne réponds à aucune question technique en PV, utilisez le forum, il est fait pour ça, et la réponse peut profiter à tous.
Répondre