LocalSplinesImpl.mesa
Copyright Ó 1987, 1992 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 ¬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 =>
{
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;
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;
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 ¬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 {
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;
};
}.