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
DIRECTORY
CubicSplines,
LocalSplines;
LocalSplinesImpl: CEDAR PROGRAM EXPORTS LocalSplines = { OPEN CubicSplines;
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: 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];
special case for 1 and 2 points
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;
allocate for allCoeffs
IF numknots>=4 THEN i ← (numknots-1)/3 ELSE i 𡤀 
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 =>
{
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 {
take knots and Spline procedure name and draw the curve
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];
};
RealSequence: TYPE = REF RealSequenceRec;
RealSequenceRec: TYPE = RECORD[element: SEQUENCE length:NAT OF REAL];
KnotsToCPoints: PROC [fpKnots: KnotSequence] RETURNS [fpCPoints: KnotSequence] = {
Perform Inversion routine to get control points from drawn points
From Yamaguci's paper
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;
};
BSpline: PROC [fpKnots4: POINTER TO FourKnots] RETURNS [coeffs: Coeffs] = TRUSTED {
take 4 knots and make the BSpline Coeficient matrix
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 {
take 4 knots and make the Catmul-Rom Coeficient matrix
tan = a(P2 - P1) where a is chosen to be .5
a: REAL 𡤅/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 {
take 4 knots and make the Bezier Coeficient matrix
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 {
ZeroSpline fills the coeff matrix with zeros. Mostly for a pesky SELECT expression, but might be generally useful
FOR xyzw: INTEGER IN [1..NDIM] DO
coeffs.t3[xyzw] ← 0;
coeffs.t2[xyzw] ← 0;
coeffs.t1[xyzw] ← 0;
coeffs.t0[xyzw] ← 0;
ENDLOOP;
};
}.