Fred Harthoorn
[email protected]
1
Tweedimensionale interpolatie
De inpterpolatie vindt plaats over 4 punten pi =
xi yi
!
(i = 0..3)
We trachten een vloeiende kromme aan te brengen door 4 punten in het platte vlak. We gaan uit van [p0 , p1 , p2 , p3 ] = [(x0 , y0 ), (x1 , y1 ), (x2 , y2 ), (x3 , y3 )] We cre¨eren twee kubische functies:
x = f (t) = a1 t3 + b1 t2 + c1 t + d1 y = g(t) = a2 t3 + b2 t2 + c2 t + d2
(1) (2)
We stellen dat xi = f (ti ) (i = 0, 1, 2, 3) yi = g(ti )
(3) (4)
De derde afgeleide van beide functies zijn constanten tussen t1 en t2 : x000 (t) = −x0 + 3x1 − 3x2 + x3 y 000 (t) = −y0 + 3y1 − 3y2 + y3
(5) (6)
We gaan beide constanten 3 keer integreren en dan krijgen we de oorspronkelijke 3e -graads kromme terug: Stel dat we 8 interpolatiepunten willen hebben. We delen dan x000 (t) door 8 en de volgende handeling voeren we 8 keer uit: x00j+1 (t) = x00j (t) + x000 (t) x0j+1 (t) = x0j (t) + x00j+1 (t) xj+1 (t) = xj (t) + x0j+1 (t)
(7) (8) (9)
00 yj+1 (t) = yj00 (t) + y 000 (t) 00 0 yj+1 (t) = yj0 (t) + yj+1 (t) 0 yj+1 (t) = yj (t) + yj+1 (t)
(10) (11) (12) 1
De re¨ele wereld beelden we af op een discrete ruimte zoals het computerscherm een vlak van bij voorbeeld 1024 bij 1024 pixels). Er moet dus een conversie plaats vinden van re¨ele getallen naar posities op het scherm, gehele getallen dus. Als we gehele getallen door 8 willen delen doen we 3 keer een shift rechts. Dit maakt de berekening stukken sneller. Hieronder volgt het algorithme opgeschreven in Delphi: TSpline = class(TGrafiek) private SplineLength : integer ; SplineArray : Array[0..SplineArrayLength-1] of Integer ; PlotDelay : integer ; SplinePts, Spline0, Spline1, Spline2, Spline3, DDDS, DDS, DS, S : integer ; tmpDDDS, TmpDDS, tmpDS, TmpS : integer ; protected public procedure SetSplineLength(n : integer) ; procedure InitSpline(Startw : Integer) ; procedure Spline(v : Integer) ; procedure extrapolation ; function SplineValue : integer ; procedure remember ; procedure restore ; function ddtSplineValue : integer ; function d2dtSplineValue : integer ; function d3dtSplineValue : integer ; procedure PlotSpline(x, y, wdth : integer ; Color : TColor) ; procedure PlotExtrapolation(j, size : integer ; clr : TColor) ; constructor Create(AOwner: TComponent); published property SplinePoints: Integer read SplineLength write SetSplineLength; property GetS : Integer Read S ; property GetdSdt : Integer Read DS ; property GetddSdt2 : Integer Read DDS ; property GetdddSdt3 : Integer Read DDDS ; end ;
constructor TSPline.Create(AOwner: TComponent); 2
begin inherited ; splinepts := 24 ; end; procedure TSpline.SetSplineLength(n : integer) ; begin SplineLength := n ; If SplineLength < 4 then SplineLength := 4 ; If SplineLength > 128 then SplineLength := 128 ; end ; procedure TSpline.InitSpline(Startw : Integer) ; var i : integer ; begin PlotDelay := SplineLength div 2 -1 ; for i := 0 to SplineLength-1 do SplineArray[i] := Startw ; Spline0 := Startw * SplineLength; Spline1 := Spline0; Spline2 := Spline0 ; Spline3 := Spline0 ; DDDS := 0 ; DDS := 0 ; DS := 0 ; S := Spline0 end ;
Procedure TSpline.Spline(v : Integer) ; var i : integer ; begin for i := 1 to SplineLength -1 do SplineArray[i-1] := SplineArray[i] ; SplineArray[SplineLength-1] := v ; Spline0 := Spline1 ; Spline1 := Spline2 ; Spline2 := Spline3 ; Spline3 := 0 ; for i := 0 to SplineLength-1 do Spline3 := Spline3 + SplineArray[i] ; DDDS := -Spline0 + 3 * Spline1 - 3 * Spline2 + Spline3 ; DDS := DDS + DDDS ; 3
DS := DS + DDS ; S := S + DS ; end ; procedure TSpline.extrapolation ; begin // DDS := DDS + DDDS ; DS := DS + DDS ; S := (S + DS) ; end ; procedure TSpline.remember ; begin TmpS := S ; TmpDS := DS ; TmpDDS := DDS ; TmpDDDS := DDDS ; end ; procedure begin DDDS := DDS := DS := S := end ;
TSpline.restore ; TmpDDDS ; TmpDDS ; TmpDS; TmpS ;
function TSpline.SplineValue : integer ; begin SplineValue := S div SplineLength end ; function TSpline.ddtSplineValue : integer ; begin ddtSplineValue := DS Div SplineLength * 4 end ; function TSpline.d2dtSplineValue : integer ; begin d2dtSplineValue := DDS 4
end ; function TSpline.d3dtSplineValue : integer ; begin d3dtSplineValue := DDDS end ; procedure TSPline.PlotSpline(x, y, wdth : integer ; Color : TColor) ; begin Spline(y) ; if x >= SplineLength then Plot(x - plotdelay, SplineValue, wdth, Color) ; end ; procedure TSPline.PlotExtrapolation(j, size : integer ; var i : integer ; begin begin Remember; for i := j to j + SplineLength div 2 + 10 do begin extrapolation; Plot(i-PlotDelay, SplineValue, Size, clr) ; end ; Restore end end ;
clr : TColor) ;
Hierbij geldt DDDS is de derde afgeleide, DDS de tweede, DS de eerste en S is de functie zelf.
5
