RegularALSplineImpl.mesa
Copyright Ó 1987, 1992 by Xerox Corporation. All rights reserved.
P.Baudelaire, M.Stone, September 26, 1980 5:45 PM
Last Edited by: Stone, March 15, 1983 3:27 pm
Russ Atkinson (RRA) February 2, 1987 10:27:48 pm PST
MakeSpline:
PUBLIC
PROC
[knots: KnotSequence, splineType: SplineType ]
RETURNS [CoeffsSequence] = {
RealSequence: TYPE = REF RealSequenceRec;
RealSequenceRec: TYPE = RECORD[element: SEQUENCE length:NAT OF REAL];
xyzw: CARDINAL[1..NDIM];
X: CARDINAL = 1;
Y: CARDINAL = 2;
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)
arraySize: CARDINAL ¬ (lastK+1)*SIZE[REAL];
coeffsSize: CARDINAL ¬ (lastK+1)*SIZE[Coeffs];
k: CARDINAL;
h: RealSequence;
a: RealSequence;
b: RealSequence;
c: RealSequence;
d: RealSequence;
r: RealSequence;
s: RealSequence;
coeffs: CoeffsSequence;
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 {
coeffs ¬ NEW[CoeffsSequenceRec[1]];
coeffs[0].t0 ¬ knots[0];
coeffs[0].t1 ¬ coeffs[0].t2 ¬ coeffs[0].t3 ¬ [0,0];
IF numKnots=2
THEN {
coeffs[0].t1[X] ¬ knots[1][X]-knots[0][X];
coeffs[0].t1[Y] ¬ knots[1][Y]-knots[0][Y];
};
RETURN[coeffs];
};
allocate 4 or 7 arrays of reals for spline computation
and 1 array of cubic coefficients to return
coeffs ¬ NEW[CoeffsSequenceRec[lastK+1]];
h ¬ NEW[RealSequenceRec[lastK+1]];
a ¬ NEW[RealSequenceRec[lastK+1]];
b ¬ NEW[RealSequenceRec[lastK+1]];
d ¬ NEW[RealSequenceRec[lastK+1]];
IF splineType = cyclicAL
THEN {
c ¬ NEW[RealSequenceRec[lastK+1]];
r ¬ NEW[RealSequenceRec[lastK+1]];
s ¬ NEW[RealSequenceRec[lastK+1]];
};
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 => {
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;
};
cyclicAL => {
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;
};
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 => {
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;
};
cyclicAL => {
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;
};
ENDCASE;
compute coefficient t2 (=second derivative/2)
SELECT splineType FROM
naturalAL => {
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;
};
cyclicAL => {
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;
};
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;
return cubic coefficients
RETURN[coeffs];
};