<> <> <> DIRECTORY CubicSplines, LocalSplines; LocalSplinesImpl: CEDAR PROGRAM EXPORTS LocalSplines = BEGIN OPEN CubicSplines; <> ONE: REAL _ 1; TWO: REAL _ 2; E1: REAL _ 10; E2: REAL _ 100; E3: REAL _ 1000; E4: REAL _ 10000; TooFewKnots: PUBLIC SIGNAL[numknots: INTEGER] = CODE; UnknownSpline: PUBLIC SIGNAL[splineType: SplineType] = CODE; RABS: PROCEDURE[r: REAL] RETURNS[REAL] = INLINE BEGIN RETURN[IF r<0 THEN -r ELSE r] END; FourKnots: TYPE=ARRAY [0..3] OF FPCoords; MakeSpline: PUBLIC PROCEDURE [knots: KnotSequence, splineType: SplineType] RETURNS [CoeffsSequence] = TRUSTED BEGIN allCoeffs: CoeffsSequence; numknots: INTEGER _ knots.length; IF numknots <= 0 THEN SIGNAL TooFewKnots[numknots]; <> IF numknots < 3 THEN BEGIN allCoeffs _ NEW[CoeffsSequenceRec[1]]; allCoeffs[0].t0 _ knots[0]; allCoeffs[0].t1 _ [0,0]; allCoeffs[0].t2 _ [0,0]; allCoeffs[0].t3 _ [0,0]; IF numknots=2 THEN BEGIN allCoeffs[0].t1[1] _ knots[1][1]-knots[0][1]; allCoeffs[0].t1[2] _ knots[1][2]-knots[0][2]; END; RETURN[allCoeffs]; END; SELECT splineType FROM bsplineInterp => BEGIN interpKnots: KnotSequence _ KnotsToCPoints[knots]; allCoeffs _ MakeLocalSpline[interpKnots,bspline]; END; bezier => BEGIN kindex,cindex,i: INTEGER; fourKnots: FourKnots; <> IF numknots>=4 THEN i _ (numknots-1)/3 ELSE i _0; IF (numknots-1) MOD 3#0 THEN i _ i+1; --corresponds to knots hanging over allCoeffs _ NEW[CoeffsSequenceRec[i]]; cindex _ 0; <> FOR kindex IN [0..3] DO IF kindex >numknots-1 THEN EXIT; fourKnots[kindex] _ knots[kindex]; ENDLOOP; <> IF numknots=3 THEN fourKnots[3] _ fourKnots[2]; kindex _ 4; <> DO allCoeffs[cindex] _ Bezier[@fourKnots]; cindex _ cindex + 1; <> fourKnots[0] _ fourKnots[3]; SELECT numknots-kindex FROM <=0 => EXIT; =1 => BEGIN fourKnots[1] _ fourKnots[0]; fourKnots[2] _ knots[kindex]; kindex _ kindex+1; fourKnots[3] _ fourKnots[2]; END; =2 => BEGIN FOR i IN [1..2] DO fourKnots[i] _ knots[kindex]; kindex _ kindex+1; ENDLOOP; fourKnots[3] _ fourKnots[2]; END; >=3 => FOR i IN [1..3] DO fourKnots[i] _ knots[kindex]; kindex _ kindex+1; ENDLOOP; ENDCASE; ENDLOOP; END; IN [bspline..crspline] => allCoeffs _ MakeLocalSpline[knots,splineType]; ENDCASE => SIGNAL UnknownSpline[splineType]; RETURN[allCoeffs]; END; <> MakeLocalSpline: PROCEDURE [knots: KnotSequence, splineType: SplineType] RETURNS [CoeffsSequence] = TRUSTED BEGIN fourKnots: FourKnots; kindex,cindex,i: INTEGER; numknots: INTEGER _ knots.length; lastknot: INTEGER _ numknots-1; allCoeffs: CoeffsSequence; IF numknots >= 4 THEN i _ numknots-3 ELSE i _ 1; allCoeffs _ NEW[CoeffsSequenceRec[i]]; cindex _ 0; <> FOR kindex IN [0..3] DO IF kindex >lastknot THEN EXIT; fourKnots[kindex] _ knots[kindex]; ENDLOOP; IF numknots=3 THEN fourKnots[3] _ fourKnots[2]; kindex _ 4; <> DO allCoeffs[cindex] _ SELECT splineType FROM bspline => BSpline[@fourKnots], crspline => CRSpline[@fourKnots], ENDCASE => ZeroSpline[@fourKnots]; cindex _ cindex+1; IF kindex > lastknot THEN EXIT; <> fourKnots[0] _ fourKnots[1]; fourKnots[1] _ fourKnots[2]; fourKnots[2] _ fourKnots[3]; fourKnots[3] _ knots[kindex]; kindex _ kindex+1; ENDLOOP; RETURN[allCoeffs]; END; RealSequence: TYPE = REF RealSequenceRec; RealSequenceRec: TYPE = RECORD[element: SEQUENCE length:NAT OF REAL]; <> <> KnotsToCPoints: PROCEDURE[fpKnots: KnotSequence] RETURNS[fpCPoints: KnotSequence] = BEGIN delta: RealSequence; deltaMax,tolerence: REAL; xyzw,i,loopCount: INTEGER; numknots: INTEGER _ fpKnots.length; fpCPoints _ NEW[KnotSequenceRec[numknots+2]]; delta _ NEW[RealSequenceRec[numknots]]; <> FOR i IN [1..numknots] DO fpCPoints[i] _ fpKnots[i-1]; ENDLOOP; <> fpCPoints[0] _ fpCPoints[1]; fpCPoints[numknots+1] _ fpCPoints[numknots]; tolerence _ ONE/3; <

> loopCount _ 0; DO deltaMax _ 0; <> FOR xyzw IN [1..NDIM] DO FOR i IN [1..numknots] DO delta[i-1] _ fpKnots[i-1][xyzw] -fpCPoints[i][xyzw] +fpKnots[i-1][xyzw]/2-(fpCPoints[i-1][xyzw]+fpCPoints[i+1][xyzw])/4; IF RABS[delta[i-1]] >= deltaMax THEN deltaMax _ RABS[delta[i-1]]; ENDLOOP; FOR i IN [1..numknots] DO fpCPoints[i][xyzw] _ delta[i-1]+fpCPoints[i][xyzw]; ENDLOOP; ENDLOOP; <> fpCPoints[0] _ fpCPoints[1]; fpCPoints[numknots+1] _ fpCPoints[numknots]; <> IF RABS[deltaMax] <= tolerence THEN EXIT; ENDLOOP; END; <> BSpline: PROCEDURE[fpKnots4: POINTER TO FourKnots] RETURNS[coeffs: Coeffs] = TRUSTED BEGIN xyzw: INTEGER; FOR xyzw IN [1..NDIM] DO coeffs.t3[xyzw] _ (-fpKnots4[0][xyzw] +3*fpKnots4[1][xyzw] -3*fpKnots4[2][xyzw] +fpKnots4[3][xyzw])/6; coeffs.t2[xyzw] _ (3*fpKnots4[0][xyzw]-6*fpKnots4[1][xyzw]+3*fpKnots4[2][xyzw])/6; coeffs.t1[xyzw] _ (-3*fpKnots4[0][xyzw]+3*fpKnots4[2][xyzw])/6; coeffs.t0[xyzw] _ (fpKnots4[0][xyzw]+4*fpKnots4[1][xyzw]+fpKnots4[2][xyzw])/6; ENDLOOP; END; <> <> CRSpline: PROCEDURE[fpKnots4: POINTER TO FourKnots] RETURNS[coeffs: Coeffs] = TRUSTED BEGIN xyzw: INTEGER; a: REAL _5/E1; FOR xyzw IN [1..NDIM] DO coeffs.t3[xyzw] _ (-a*fpKnots4[0][xyzw] +(2-a)*fpKnots4[1][xyzw] +(a-2)*fpKnots4[2][xyzw] +a*fpKnots4[3][xyzw]); coeffs.t2[xyzw] _ (2*a*fpKnots4[0][xyzw]+(a-3)*fpKnots4[1][xyzw]+(3-2*a)*fpKnots4[2][xyzw]-a*fpKnots4[3][xyzw]); coeffs.t1[xyzw] _ (-a*fpKnots4[0][xyzw]+a*fpKnots4[2][xyzw]); coeffs.t0[xyzw] _ fpKnots4[1][xyzw]; ENDLOOP; END; <> Bezier: PROCEDURE[fpKnots4: POINTER TO FourKnots] RETURNS[coeffs: Coeffs] = TRUSTED BEGIN xyzw: INTEGER; FOR xyzw IN [1..NDIM] DO coeffs.t3[xyzw] _ -fpKnots4[0][xyzw] +3*fpKnots4[1][xyzw] -3*fpKnots4[2][xyzw] +fpKnots4[3][xyzw]; coeffs.t2[xyzw] _ 3*fpKnots4[0][xyzw]-6*fpKnots4[1][xyzw]+3*fpKnots4[2][xyzw]; coeffs.t1[xyzw] _ -3*fpKnots4[0][xyzw]+3*fpKnots4[1][xyzw]; coeffs.t0[xyzw] _ fpKnots4[0][xyzw]; ENDLOOP; END; <> <> ZeroSpline: PROCEDURE[fpKnots4: POINTER TO FourKnots] RETURNS[coeffs: Coeffs] = TRUSTED BEGIN xyzw: INTEGER; FOR xyzw IN [1..NDIM] DO coeffs.t3[xyzw] _ 0; coeffs.t2[xyzw] _ 0; coeffs.t1[xyzw] _ 0; coeffs.t0[xyzw] _ 0; ENDLOOP; END; END.