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. ncompiler localspline m.stone September 26, 1980 5:09 PM Last Edited by: Stone, March 15, 1983 3:27 pm here are some REAL constants used for assigning real constants. special case for 1 and 2 points allocate for allCoeffs set up first 4 knots, duplicating if necessary only case where duplication is necessary since check above for numknots<3 need to do this at least once Take 4th knot + next 3 knots. Special case if less than 3 left take knots and Spline procedure name and draw the curve set up first 4 knots need to do this at least once shift the segment over Preform Inversion routine to get control points from drawn points From Yamaguci's paper inits set the first and last control points, and the tolerence Main Loop iterate and collect maximum delta set the first and last control points End Check take 4 knots and make the BSpline Coeficient matrix take 4 knots and make the Catmul-Rom Coeficient matrix tan = a(P2 - P1) where a is chosen to be .5 take 4 knots and make the Bezier Coeficient matrix ZeroSpline fills the coeff matrix with zeros. Mostly for a pesky SELECT expression, but might be generally useful Ê(˜J˜Jšœ™Jšœ#™#Jšœ-™-šÏk ˜ J˜ Jšœ ˜ —J˜Jšœ œœ˜7Jšœœ˜J˜Jšœ?™?Jšœœœœ˜Jš œœœœ œ ˜AJš œ œœœœ˜5Jšœœœœ˜