--P.Baudelaire, M. Stone, September 26, 1980 5:20 PM
DIRECTORY
SplineDefs: FROM "SplineDefs",
SplineMemoryDefs: FROM "SplineMemoryDefs";

RegularUMSpline: PROGRAM IMPORTS SplineMemoryDefs EXPORTS SplineDefs =
BEGIN OPEN SplineDefs;

--this procedure implements the natural spline algorithm with unit knot mesh
--it expects splineType to be: naturalUM or cyclicUM
TooFewKnots: PUBLIC SIGNAL[numknots: INTEGER] = CODE;
UnknownSpline: PUBLIC SIGNAL[splineType: SplineType] = CODE;
UnmatchedEnds: PUBLIC SIGNAL = CODE;

MakeSpline: PUBLIC PROCEDURE
[knots: DESCRIPTOR FOR ARRAY OF FPCoords, splineType: SplineType]
RETURNS [DESCRIPTOR FOR ARRAY OF Coeffs] =
BEGIN

xyzw: CARDINAL[1..NDIM];
X: CARDINAL = 1;
Y: CARDINAL = 2;
numKnots: CARDINAL ← LENGTH[knots]; --number of knots
lastKnot: CARDINAL ← numKnots-1; --max for knot index
lastK: CARDINAL ← lastKnot-1; --max for interval index (for coeffs and arrays a, b, c, r, & s)
arraySize: CARDINAL ← (lastK+1)*SIZE[REAL];
coeffsSize: CARDINAL ← (lastK+1)*SIZE[Coeffs];
k: CARDINAL;
a: DESCRIPTOR FOR ARRAY OF REAL;
b: DESCRIPTOR FOR ARRAY OF REAL;
c: DESCRIPTOR FOR ARRAY OF REAL;
r: DESCRIPTOR FOR ARRAY OF REAL;
s: DESCRIPTOR FOR ARRAY OF REAL;
coeffs: DESCRIPTOR FOR ARRAY OF Coeffs;

--first some error checks
IF numKnots <= 0 THEN SIGNAL TooFewKnots[numKnots];
IF splineType # naturalUM AND splineType # cyclicUM THEN
SIGNAL UnknownSpline[splineType];
IF splineType = cyclicUM THEN FOR xyzw IN [1..NDIM] DO
IF knots[0][xyzw] # knots[lastKnot][xyzw] THEN GOTO notCyclic;
REPEAT
notCyclic => SIGNAL UnmatchedEnds[];
ENDLOOP;

--special case for 1 and 2 points
IF numKnots < 3 THEN BEGIN
coeffs ← DESCRIPTOR[SplineMemoryDefs.Allocate[SIZE[Coeffs]], 1];
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;
RETURN[coeffs];
END;

--from here on down can assume at least 3 points
--allocate 1 array of cubic coefficients to return
coeffs ← DESCRIPTOR[SplineMemoryDefs.Allocate[coeffsSize], lastK+1];
--allocate 2 or 5 arrays of reals for spline computation
a ← DESCRIPTOR[SplineMemoryDefs.Allocate[arraySize], lastK+1];
b ← DESCRIPTOR[SplineMemoryDefs.Allocate[arraySize], lastK+1];
IF splineType = cyclicUM THEN
BEGIN
c ← DESCRIPTOR[SplineMemoryDefs.Allocate[arraySize], lastK+1];
r ← DESCRIPTOR[SplineMemoryDefs.Allocate[arraySize], lastK+1];
s ← DESCRIPTOR[SplineMemoryDefs.Allocate[arraySize], lastK+1];
END;

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

--precompute coefficients a & c
a[0] ← 4;
FOR k IN [1..lastK-1] DO
a[k] ← 4 - 1/a[k-1];
ENDLOOP;
IF splineType = cyclicUM THEN
BEGIN
c[0] ← 1;
FOR k IN [1..lastK-1] DO
c[k] ← -c[k-1]/a[k-1];
ENDLOOP;
END;

FOR xyzw IN [1..NDIM] DO
--compute each dimension separately

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

--compute coefficient t2 (=second derivative/2)
SELECT splineType FROM
naturalUM => 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) - coeffs[k+1].t2[xyzw]) / a[k-1];
ENDLOOP;
END;
cyclicUM => BEGIN
p2: REAL;
p2 ← (6*(knots[lastKnot][xyzw] - 2*knots[lastKnot-1][xyzw] + knots[lastKnot-2][xyzw])
- s[0] - s[lastK-1]) / (r[0] + r[lastK-1] + 4);
FOR k IN [0..lastK-1] DO
coeffs[k].t2[xyzw] ← (r[k] * p2 + s[k]) / 2;
ENDLOOP;
coeffs[lastK].t2[xyzw] ← p2/2;
--note: coeffs[lastK+1].t2 = coeffs[0].t2
END;
ENDCASE;

--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] ← knots[k+1][xyzw] - knots[k][xyzw]
- (2*coeffs[k].t2[xyzw] + coeffs[k+1].t2[xyzw]) / 3;
ENDLOOP;
--note again: coeffs[lastK+1].t2 = coeffs[0].t2 (=0 for natural end conditions)
coeffs[lastK].t3[xyzw] ← (coeffs[0].t2[xyzw] - coeffs[lastK].t2[xyzw] ) / 3;
coeffs[lastK].t1[xyzw] ← knots[lastK+1][xyzw] - knots[lastK][xyzw]
- (2 * coeffs[lastK].t2[xyzw] + coeffs[0].t2[xyzw]) / 3;

ENDLOOP;

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

--return cubic coefficients
RETURN[coeffs];
END;

END.