Re: Interpolation cubique (courbe bézier...)
Publié : mar. 24/juil./2012 10:52
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)