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

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

--this procedure implements the natural spline algorithm
--with polygonal line length parametrization
--it expects splineType to be: naturalAL or cyclicAL
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 h, a, b, c, d, r, & s)
arraySize: CARDINAL _ (lastK+1)*SIZE[REAL];
coeffsSize: CARDINAL _ (lastK+1)*SIZE[Coeffs];
k: CARDINAL;
h: DESCRIPTOR FOR ARRAY OF REAL;
a: DESCRIPTOR FOR ARRAY OF REAL;
b: DESCRIPTOR FOR ARRAY OF REAL;
c: DESCRIPTOR FOR ARRAY OF REAL;
d: DESCRIPTOR FOR ARRAY OF REAL;
r: DESCRIPTOR FOR ARRAY OF REAL;
s: DESCRIPTOR FOR ARRAY OF REAL;
coeffs: DESCRIPTOR FOR ARRAY OF Coeffs;
deltaX, deltaY: REAL;

--first some error checks
IF numKnots <= 0 THEN SIGNAL TooFewKnots[numKnots];
IF splineType # naturalAL AND splineType # cyclicAL THEN
	SIGNAL UnknownSpline[splineType];
IF splineType = cyclicAL 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;

--allocate 4 or 7 arrays of reals for spline computation
--and 1 array of cubic coefficients to return
coeffs _ DESCRIPTOR[SplineMemoryDefs.Allocate[coeffsSize], lastK+1];
h _ DESCRIPTOR[SplineMemoryDefs.Allocate[arraySize], lastK+1];
a _ DESCRIPTOR[SplineMemoryDefs.Allocate[arraySize], lastK+1];
b _ DESCRIPTOR[SplineMemoryDefs.Allocate[arraySize], lastK+1];
d _ DESCRIPTOR[SplineMemoryDefs.Allocate[arraySize], lastK+1];
IF splineType = cyclicAL 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;

--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;

FOR xyzw IN [1..NDIM] 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;

	--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
RETURN[coeffs];
END;

END.
(635)