-- CGSplineImpl.mesa
-- Last edit by Doug Wyatt, August 23, 1982 4:16 pm

-- Derived from RegularALSpline.mesa, by
-- P.Baudelaire, M.Stone, September 26, 1980 5:45 PM

DIRECTORY
CGSpline USING [Coeffs, Dim, ErrorType, Knots, SplineType],
CGStorage USING [pZone, qZone];

CGSplineImpl: CEDAR PROGRAM
IMPORTS CGStorage
EXPORTS CGSpline =
BEGIN OPEN CGSpline;

repZone: ZONE = CGStorage.qZone;
arrayZone: ZONE = CGStorage.pZone;

Error: PUBLIC ERROR[type: ErrorType] = CODE;

Array: TYPE = REF ArrayRep;
ArrayRep: TYPE = RECORD[SEQUENCE space: NAT OF REAL];

CoeffsArray: TYPE = REF CoeffsArrayRep;
CoeffsArrayRep: TYPE = RECORD[SEQUENCE space: NAT OF Coeffs];

Ref: TYPE = REF Rep;
Rep: PUBLIC TYPE = RECORD [
space: NAT,
h,a,b,c,d,r,s: Array ← NIL,
coeffs: CoeffsArray ← NIL
];

New: PUBLIC PROC[size: NAT] RETURNS[Ref] = {
self: Ref ← repZone.NEW[Rep ← [space: 0]];
Grow[self,size]; RETURN[self];
};

Grow: PROC[self: Ref, size: NAT] = {
IF size<=self.space THEN RETURN;
self.h ← arrayZone.NEW[ArrayRep[size]];
self.a ← arrayZone.NEW[ArrayRep[size]];
self.b ← arrayZone.NEW[ArrayRep[size]];
self.c ← arrayZone.NEW[ArrayRep[size]];
self.d ← arrayZone.NEW[ArrayRep[size]];
self.r ← arrayZone.NEW[ArrayRep[size]];
self.s ← arrayZone.NEW[ArrayRep[size]];
self.coeffs ← arrayZone.NEW[CoeffsArrayRep[size]];
self.space ← size;
};

--this procedure implements the natural spline algorithm
--with polygonal line length parametrization
--it expects splineType to be: naturalAL or cyclicAL

MakeSpline: PUBLIC PROCEDURE [self: Ref,
knots: Knots, splineType: SplineType, Proc: PROC[Coeffs]] =
BEGIN

xyzw: Dim;
numKnots: CARDINAL ← knots.length; --number of knots
lastKnot: CARDINAL ← numKnots-1; --max for knot index
lastK: CARDINAL ← lastKnot-1; --max for interval index (for coeffs and arrays h, a, b, c, d, r, & s)
size: CARDINAL ← (lastK+1);
k: CARDINAL;
h: Array ← NIL;
a: Array ← NIL;
b: Array ← NIL;
c: Array ← NIL;
d: Array ← NIL;
r: Array ← NIL;
s: Array ← NIL;
coeffs: CoeffsArray ← NIL;
deltaX, deltaY: REAL;

--first some error checks
IF numKnots <= 0 THEN ERROR Error[tooFewKnots];
-- IF splineType # naturalAL AND splineType # cyclicAL THEN
--  SIGNAL UnknownSpline[splineType];
IF splineType = cyclicAL THEN FOR xyzw IN Dim DO
 IF knots[0][xyzw] # knots[lastKnot][xyzw] THEN GOTO notCyclic;
 REPEAT
  notCyclic => ERROR Error[unmatchedEnds];
 ENDLOOP;

--special case for 1 and 2 points
IF numKnots < 3 THEN BEGIN
 coeffs: ARRAY [0..1) OF Coeffs;
 coeffs[0].t0 ← knots[0];
 coeffs[0].t1 ← coeffs[0].t2 ← coeffs[0].t3 ← [0,0];
 IF numKnots=2 THEN BEGIN
  coeffs[0].t1[X] ← knots[1][X]-knots[0][X];
  coeffs[0].t1[Y] ← knots[1][Y]-knots[0][Y];
  END;
 Proc[coeffs[0]];
 RETURN;
 END;

--allocate 4 or 7 arrays of reals for spline computation
--and 1 array of cubic coefficients to return
IF size>self.space THEN Grow[self,size];
coeffs ← self.coeffs;
h ← self.h;
a ← self.a;
b ← self.b;
d ← self.d;
IF splineType = cyclicAL THEN
 BEGIN
 c ← self.c;
 r ← self.r;
 s ← self.s;
 END;

--first get constant cubic coefficients
FOR k IN [0..lastK] DO
 coeffs[k].t0 ← knots[k];
 ENDLOOP;

--compute parametrization (polygonal line approximation)
--this assumes NDIM = 2
FOR k IN [0..lastK] DO
 deltaX ← knots[k+1][X] - knots[k][X];
 IF deltaX < 0 THEN deltaX ← -deltaX;
 deltaY ← knots[k+1][Y] - knots[k][Y];
 IF deltaY < 0 THEN deltaY ← -deltaY;
 h[k] ← IF deltaX > deltaY
   THEN (deltaX + deltaY/2)
   ELSE (deltaY + deltaX/2);
 --check overlapping knots!
 IF h[k] = 0 THEN h[k] ← 1;
 ENDLOOP;

SELECT splineType FROM
naturalAL => BEGIN
 --natural end conditions: precompute coefficients a
 a[0] ← 2 * (h[0] + h[1]);
 FOR k IN [1..lastK-1] DO
  a[k] ← 2 * (h[k] + h[k+1]) - h[k] * h[k] / a[k-1];
  ENDLOOP;
 END;
cyclicAL => BEGIN
 --cyclic end conditions: precompute coefficients a & c
 a[0] ← 2 * (h[0] + h[lastK]);
 FOR k IN [1..lastK-1] DO
  a[k] ← 2 * (h[k-1] + h[k]) - h[k-1] * h[k-1] / a[k-1];
  ENDLOOP;
 c[0] ← h[lastK];
 FOR k IN [1..lastK-1] DO
  c[k] ← - h[k-1] * c[k-1] / a[k-1];
  ENDLOOP;
 END;
 ENDCASE => ERROR Error[bug];

FOR xyzw IN Dim DO
--compute each dimension separately

 --compute coefficient d
 FOR k IN [0..lastK] DO
  d[k] ← (knots[k+1][xyzw] - knots[k][xyzw]) / h[k];
  ENDLOOP;

 IF numKnots > 2 THEN SELECT splineType FROM
 naturalAL => BEGIN
  --natural end conditions: compute coefficients b
  FOR k IN [0..lastK-1] DO
   b[k] ← 6 * (d[k+1] - d[k]);
   ENDLOOP;
  FOR k IN [1..lastK-1] DO
   b[k] ← b[k] - h[k] * b[k-1] / a[k-1];
   ENDLOOP;
  END;
 cyclicAL => BEGIN
  --cyclic end conditions: compute coefficients b, r, & s
  b[0] ← 6 * (d[0] - d[lastK]);
  FOR k IN [1..lastK-1] DO
   b[k] ← 6 * (d[k] - d[k-1]);
   ENDLOOP;
  FOR k IN [1..lastK-1] DO
   b[k] ← b[k] - h[k-1] * b[k-1] / a[k-1];
   ENDLOOP;
  r[lastK] ← 1;
  s[lastK] ← 0;
  FOR k DECREASING IN [0..lastK-1] DO
   r[k] ← - (h[k] * r[k+1] + c[k]) / a[k];
   s[k] ← (b[k] - h[k] * s[k+1]) / a[k];
   ENDLOOP;
  END;
 ENDCASE;

 --compute coefficient t2 (=second derivative/2)
 SELECT splineType FROM
 naturalAL => BEGIN
  coeffs[0].t2[xyzw] ← 0;
  --note: coeffs[lastK+1].t2.[xyzw] is 0 too but it is not stored, hence:
  coeffs[lastK].t2[xyzw] ← b[lastK-1] / (2 * a[lastK-1]);
  FOR k DECREASING IN [1..lastK-1] DO
   coeffs[k].t2[xyzw] ← ((b[k-1] / 2) - h[k] * coeffs[k+1].t2[xyzw]) / a[k-1];
   ENDLOOP;
  END;
 cyclicAL => BEGIN
  p2: REAL;
  p2 ← 6*(d[lastK] - d[lastK-1]) - h[lastK] * s[0] - h[lastK-1] * s[lastK-1];
  p2 ← p2 / (h[lastK] * r[0] + h[lastK-1] * r[lastK-1] + 2 * (h[lastK] + h[lastK-1]));
  --note: coeffs[lastK+1].t2.[xyzw] = coeffs[0].t2[xyzw]
  coeffs[lastK].t2[xyzw] ← p2/2;
  FOR k IN [0..lastK-1] DO
   coeffs[k].t2[xyzw] ← (r[k] * p2 + s[k]) / 2;
   ENDLOOP;
  END;
 ENDCASE => ERROR Error[bug];

 --compute coefficients t3 (=third derivative/6) & and t1 (=first derivative)
 FOR k IN [0..lastK-1] DO
  coeffs[k].t3[xyzw] ← (coeffs[k+1].t2[xyzw] - coeffs[k].t2[xyzw]) / 3;
  coeffs[k].t1[xyzw] ← d[k] - h[k] * (2*coeffs[k].t2[xyzw] + coeffs[k+1].t2[xyzw]) / 3;
  ENDLOOP;
 --note again: coeffs[lastK+1].t2[xyzw] = coeffs[0].t2[xyzw] (=0 for natural end conditions)
 coeffs[lastK].t3[xyzw] ← (coeffs[0].t2[xyzw] - coeffs[lastK].t2[xyzw]) / 3;
 coeffs[lastK].t1[xyzw] ← d[lastK] - h[lastK] * (2 * coeffs[lastK].t2[xyzw] + coeffs[0].t2[xyzw]) / 3;

 --now reparametrize coefficients to interval [0..1]
 FOR k IN [0..lastK] DO
  coeffs[k].t3[xyzw] ← h[k] * h[k] * coeffs[k].t3[xyzw];
  coeffs[k].t2[xyzw] ← h[k] * h[k] * coeffs[k].t2[xyzw];
  coeffs[k].t1[xyzw] ← h[k] * coeffs[k].t1[xyzw];
  ENDLOOP;

 ENDLOOP;

--release computation arrays
-- SplineMemoryDefs.Free[BASE[h]];
-- SplineMemoryDefs.Free[BASE[a]];
-- SplineMemoryDefs.Free[BASE[b]];
-- SplineMemoryDefs.Free[BASE[d]];
-- IF splineType = cyclicAL THEN
--  BEGIN
--  SplineMemoryDefs.Free[BASE[c]];
--  SplineMemoryDefs.Free[BASE[r]];
--  SplineMemoryDefs.Free[BASE[s]];
--  END;

--return cubic coefficients
FOR i: CARDINAL IN[0..size) DO Proc[coeffs[i]] ENDLOOP;
-- RETURN[coeffs];
END;

END.