DIRECTORY CubicSplines, LocalSplines; LocalSplinesImpl: CEDAR PROGRAM EXPORTS LocalSplines = { 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: PROC [r: REAL] RETURNS [REAL] = INLINE { RETURN [ABS[r]] }; FourKnots: TYPE=ARRAY [0..3] OF FPCoords; MakeSpline: PUBLIC PROC [knots: KnotSequence, splineType: SplineType] RETURNS [CoeffsSequence] = TRUSTED { allCoeffs: CoeffsSequence; numknots: INTEGER _ knots.length; IF numknots <= 0 THEN SIGNAL TooFewKnots[numknots]; IF numknots < 3 THEN { 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 { allCoeffs[0].t1[1] _ knots[1][1]-knots[0][1]; allCoeffs[0].t1[2] _ knots[1][2]-knots[0][2]; }; RETURN[allCoeffs]; }; SELECT splineType FROM bsplineInterp => { interpKnots: KnotSequence _ KnotsToCPoints[knots]; allCoeffs _ MakeLocalSpline[interpKnots,bspline]; }; bezier => { 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 => { fourKnots[1] _ fourKnots[0]; fourKnots[2] _ knots[kindex]; kindex _ kindex+1; fourKnots[3] _ fourKnots[2]; }; =2 => { FOR i IN [1..2] DO fourKnots[i] _ knots[kindex]; kindex _ kindex+1; ENDLOOP; fourKnots[3] _ fourKnots[2]; }; >=3 => FOR i IN [1..3] DO fourKnots[i] _ knots[kindex]; kindex _ kindex+1; ENDLOOP; ENDCASE; ENDLOOP; }; IN [bspline..crspline] => allCoeffs _ MakeLocalSpline[knots,splineType]; ENDCASE => SIGNAL UnknownSpline[splineType]; RETURN[allCoeffs]; }; MakeLocalSpline: PROC [knots: KnotSequence, splineType: SplineType] RETURNS [CoeffsSequence] = TRUSTED { 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]; }; RealSequence: TYPE = REF RealSequenceRec; RealSequenceRec: TYPE = RECORD[element: SEQUENCE length:NAT OF REAL]; KnotsToCPoints: PROC [fpKnots: KnotSequence] RETURNS [fpCPoints: KnotSequence] = { 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; }; BSpline: PROC [fpKnots4: POINTER TO FourKnots] RETURNS [coeffs: Coeffs] = TRUSTED { FOR xyzw: INTEGER 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; }; CRSpline: PROC [fpKnots4: POINTER TO FourKnots] RETURNS [coeffs: Coeffs] = TRUSTED { a: REAL _5/E1; FOR xyzw: INTEGER 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; }; Bezier: PROC [fpKnots4: POINTER TO FourKnots] RETURNS [coeffs: Coeffs] = TRUSTED { FOR xyzw: INTEGER 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; }; ZeroSpline: PROC [fpKnots4: POINTER TO FourKnots] RETURNS [coeffs: Coeffs] = TRUSTED { FOR xyzw: INTEGER IN [1..NDIM] DO coeffs.t3[xyzw] _ 0; coeffs.t2[xyzw] _ 0; coeffs.t1[xyzw] _ 0; coeffs.t0[xyzw] _ 0; ENDLOOP; }; }. ΰLocalSplinesImpl.mesa Copyright Σ 1987 by Xerox Corporation. All rights reserved. m.stone September 26, 1980 5:09 PM Last Edited by: Stone, March 15, 1983 3:27 pm Russ Atkinson (RRA) February 2, 1987 10:32:54 pm PST 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 Perform 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 Κχ˜codešœ™K™