--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 𡤀 
  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 𡤅/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.