DIRECTORY CubicSplines, RegularUMSpline; RegularUMSplineImpl: CEDAR PROGRAM EXPORTS RegularUMSpline = { 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; --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: RealSequence; b: RealSequence; c: RealSequence; r: RealSequence; s: RealSequence; coeffs: CoeffsSequence; 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; 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]]; a _ NEW[RealSequenceRec[lastK+1]]; b _ NEW[RealSequenceRec[lastK+1]]; IF splineType = cyclicUM 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; a[0] _ 4; FOR k IN [1..lastK-1] DO a[k] _ 4 - 1/a[k-1]; ENDLOOP; IF splineType = cyclicUM THEN { c[0] _ 1; FOR k IN [1..lastK-1] DO c[k] _ -c[k-1]/a[k-1]; ENDLOOP; }; FOR xyzw IN [1..NDIM] DO IF numKnots > 2 THEN SELECT splineType FROM naturalUM => { 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; }; cyclicUM => { 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; }; ENDCASE; SELECT splineType FROM naturalUM => { 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) - coeffs[k+1].t2[xyzw]) / a[k-1]; ENDLOOP; }; cyclicUM => { 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; }; 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] _ knots[k+1][xyzw] - knots[k][xyzw] - (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] _ knots[lastK+1][xyzw] - knots[lastK][xyzw] - (2 * coeffs[lastK].t2[xyzw] + coeffs[0].t2[xyzw]) / 3; ENDLOOP; RETURN[coeffs]; }; }. ”RegularUMSplineImpl.mesa Copyright Σ 1987 by Xerox Corporation. All rights reserved. P.Baudelaire, M. Stone, September 26, 1980 5:20 PM Last Edited by: Stone, March 15, 1983 3:27 pm Last Edited by: Pier, February 13, 1984 3:44:08 pm PST Russ Atkinson (RRA) February 2, 1987 10:25:46 pm PST this procedure implements the natural spline algorithm with unit knot mesh it expects splineType to be: naturalUM or cyclicUM first some error checks special case for 1 and 2 points from here on down can assume at least 3 points allocate 1 array of cubic coefficients to return allocate 2 or 5 arrays of reals for spline computation first get constant cubic coefficients precompute coefficients a & c compute each dimension separately 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 = coeffs[0].t2 compute coefficients t3 (=third derivative/6) & and t1 (=first derivative) note again: coeffs[lastK+1].t2 = coeffs[0].t2 (=0 for natural end conditions) return cubic coefficients Κ\˜code™K™š˜Kšœ œ˜$—Kšœ˜—K˜Kšœ™šœœ˜Kšœ œ˜#K˜K˜3šœ œ˜K˜*K˜*Kšœ˜—Kšœ ˜Kšœ˜—K˜Kšœ.™.Kšœ0™0Kšœ œ˜)Kšœ6™6Kšœœ˜"Kšœœ˜"šœœ˜Kšœœ˜"Kšœœ˜"Kšœœ˜"Kšœ˜—K˜Kšœ%™%šœœ ˜K˜Kšœ˜—K˜Kšœ™K˜ šœœ˜K˜Kšœ˜—šœœ˜K˜ šœœ˜K˜Kšœ˜—Kšœ˜—K˜Kšœœœ˜Kšœ!™!K˜šœœœ ˜+šœ˜Kšœ.™.K˜>šœœ˜˜AK˜—Kšœ˜—Kšœ˜—šœ ˜ Kšœ5™5K˜Bšœœ˜˜AK˜—Kšœ˜—K˜ K˜ šœ œœ˜#K˜ K˜Kšœ˜—Kšœ˜—Kšœ˜—˜Kšœ-™-Kšœ ˜šœ˜K˜KšœE™EK˜7šœ œœ˜#K˜EKšœ˜—Kšœ˜—šœ ˜ Kšœœ˜ ˜UK˜/—šœœ˜K˜,Kšœ˜—K˜Kšœ(™(Kšœ˜—Kšœ˜—˜KšœJ™Jšœœ˜K˜E˜6K˜4—Kšœ˜—KšœN™NK˜M˜BK˜8——˜Kšœ˜—K˜Kšœ™Kšœ ˜Kšœ˜—K˜Kšœ˜K˜K˜—…—