--compiler localspline --m.stone September 26, 1980 5:09 PM -- Last Edited by: Stone, March 15, 1983 3:27 pm DIRECTORY SplineDefs: FROM "SplineDefs"; LocalSpline: PROGRAM EXPORTS SplineDefs = BEGIN OPEN SplineDefs; --here are some REAL constants used for assigning real constants. 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] = BEGIN allCoeffs: CoeffsSequence; numknots: INTEGER _ knots.length; IF numknots <= 0 THEN SIGNAL TooFewKnots[numknots]; --special case for 1 and 2 points 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; --allocate for allCoeffs 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; --set up first 4 knots, duplicating if necessary FOR kindex IN [0..3] DO IF kindex >numknots-1 THEN EXIT; fourKnots[kindex] _ knots[kindex]; ENDLOOP; --only case where duplication is necessary since check above for numknots<3 IF numknots=3 THEN fourKnots[3] _ fourKnots[2]; kindex _ 4; --need to do this at least once DO allCoeffs[cindex] _ Bezier[@fourKnots]; cindex _ cindex + 1; --Take 4th knot + next 3 knots. Special case if less than 3 left 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; --take knots and Spline procedure name and draw the curve MakeLocalSpline: PROCEDURE [knots: KnotSequence, splineType: SplineType] RETURNS [CoeffsSequence] = 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; --set up first 4 knots 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; --need to do this at least once DO allCoeffs[cindex] _ SELECT splineType FROM bspline => BSpline[@fourKnots], crspline => CRSpline[@fourKnots], ENDCASE => ZeroSpline[@fourKnots]; cindex _ cindex+1; IF kindex > lastknot THEN EXIT; --shift the segment over 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]; --Preform Inversion routine to get control points from drawn points --From Yamaguci's paper 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]]; --inits FOR i IN [1..numknots] DO fpCPoints[i] _ fpKnots[i-1]; ENDLOOP; --set the first and last control points, and the tolerence fpCPoints[0] _ fpCPoints[1]; fpCPoints[numknots+1] _ fpCPoints[numknots]; tolerence _ ONE/3; --Main Loop loopCount _ 0; DO deltaMax _ 0; --iterate and collect maximum delta 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; --set the first and last control points fpCPoints[0] _ fpCPoints[1]; fpCPoints[numknots+1] _ fpCPoints[numknots]; --End Check IF RABS[deltaMax] <= tolerence THEN EXIT; ENDLOOP; END; --take 4 knots and make the BSpline Coeficient matrix BSpline: PROCEDURE[fpKnots4: POINTER TO FourKnots] RETURNS[coeffs: Coeffs] = 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; --take 4 knots and make the Catmul-Rom Coeficient matrix --tan = a(P2 - P1) where a is chosen to be .5 CRSpline: PROCEDURE[fpKnots4: POINTER TO FourKnots] RETURNS[coeffs: Coeffs] = 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; --take 4 knots and make the Bezier Coeficient matrix Bezier: PROCEDURE[fpKnots4: POINTER TO FourKnots] RETURNS[coeffs: Coeffs] = 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 fills the coeff matrix with zeros. Mostly for a pesky SELECT expression, but --might be generally useful ZeroSpline: PROCEDURE[fpKnots4: POINTER TO FourKnots] RETURNS[coeffs: Coeffs] = 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. ʘJšóÏcnÏk œžœžœžœžœžœBžœžœžœžœ žœ žœ žœ žœžœžœ žœžœžœžœžœžœž œžœžœžœžœžœžœžœžœžœžœžœžœžœ Ïn œžœž œ0žœžœ&žœžœžœžœ"žœžœžœžœ„žœ žœžœcžœžœžœžœ žœžœlžœžœžœœžœ žœžœ žœžœžœ $œžœ&1œžœžœžœžœžœžœ+žœLœžœ žœ, œžœDBœ#žœžœ žœžœ…žœžœžœžœžœDžœ)žœžœžœžœDžœžœžœžœžœIžœžœžœ žœ:Ÿœž œ1žœžœ(žœ žœžœ,žœžœžœžœ%žœžœžœžœžœžœ'žœžœ žœ+ žœžœ žœIžœ1žœžœžœœŽžœžœ žœžœžœ#žœžœ žœžœžœžœ\Ÿœž œžœžœ*žœžœ žœ žœ'žœœžœžœžœžœ;œYžœ œžœ$œžœžœžœžœžœžœžœ…žœžœžœ žœžœžœžœžœ8žœžœ(œL œžœžœžœžœžœžœ6Ÿœž œ žœžœ žœžœžœžœžœžœžœÏžœžœgŸœž œ žœžœ žœžœžœžœžœžœžœžœËžœžœ5Ÿœž œ žœžœ žœžœžœžœžœžœžœ™žœžœwŸ œž œ žœžœ žœžœžœžœžœžœžœZžœžœžœ˜‚8—…—