DIRECTORY CubicSplines, RegularALSpline; RegularALSplineImpl: CEDAR PROGRAM EXPORTS RegularALSpline = { OPEN CubicSplines; TooFewKnots: PUBLIC SIGNAL[numknots: INTEGER] = CODE; UnknownSpline: PUBLIC SIGNAL[splineType: SplineType] = CODE; UnmatchedEnds: PUBLIC SIGNAL = CODE; 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; 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; 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; 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]; }; 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]]; }; FOR k IN [0..lastK] DO coeffs[k].t0 _ knots[k]; ENDLOOP; 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); IF h[k] = 0 THEN h[k] _ 1; ENDLOOP; SELECT splineType FROM naturalAL => { 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 => { 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 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 => { 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 => { 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; SELECT splineType FROM naturalAL => { coeffs[0].t2[xyzw] _ 0; 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])); 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; 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; 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; 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; RETURN[coeffs]; }; }. ’RegularALSplineImpl.mesa Copyright Ó 1987 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 this procedure implements the natural spline algorithm with polygonal line length parametrization it expects splineType to be: naturalAL or cyclicAL max for interval index (for coeffs and arrays h, a, b, c, d, r, & s) first some error checks special case for 1 and 2 points allocate 4 or 7 arrays of reals for spline computation and 1 array of cubic coefficients to return first get constant cubic coefficients compute parametrization (polygonal line approximation) this assumes NDIM = 2 check overlapping knots! natural end conditions: precompute coefficients a cyclic end conditions: precompute coefficients a & c compute each dimension separately compute coefficient d natural end conditions: compute coefficients b cyclic end conditions: compute coefficients b, r, & s compute coefficient t2 (=second derivative/2) note: coeffs[lastK+1].t2.[xyzw] is 0 too but it is not stored, hence: note: coeffs[lastK+1].t2.[xyzw] = coeffs[0].t2[xyzw] compute coefficients t3 (=third derivative/6) & and t1 (=first derivative) note again: coeffs[lastK+1].t2[xyzw] = coeffs[0].t2[xyzw] (=0 for natural end conditions) now reparametrize coefficients to interval [0..1] return cubic coefficients ÊÃ˜codešœ™K™