DIRECTORY CubicSplines, RegularUMSpline; RegularUMSplineImpl: CEDAR PROGRAM EXPORTS RegularUMSpline = BEGIN OPEN CubicSplines; TooFewKnots: PUBLIC SIGNAL[numknots: INTEGER] = CODE; UnknownSpline: PUBLIC SIGNAL[splineType: SplineType] = CODE; UnmatchedEnds: PUBLIC SIGNAL = CODE; MakeSpline: PUBLIC PROCEDURE [knots: KnotSequence, splineType: SplineType] RETURNS [CoeffsSequence] = BEGIN 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 BEGIN coeffs _ NEW[CoeffsSequenceRec[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; coeffs _ NEW[CoeffsSequenceRec[lastK+1]]; a _ NEW[RealSequenceRec[lastK+1]]; b _ NEW[RealSequenceRec[lastK+1]]; IF splineType = cyclicUM THEN BEGIN c _ NEW[RealSequenceRec[lastK+1]]; r _ NEW[RealSequenceRec[lastK+1]]; s _ NEW[RealSequenceRec[lastK+1]]; END; 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 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 IF numKnots > 2 THEN SELECT splineType FROM naturalUM => BEGIN 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 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; SELECT splineType FROM naturalUM => BEGIN 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; 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; END; 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]; END; END. 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 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 Êp˜J˜Jšœ3™3Jšœ-™-Jšœ6™6šÏk ˜ J˜ J˜—J˜Jšœ œœ˜šœœ˜˜AJ˜—Jšœ˜—Jšœ˜—šœ˜Jšœ5™5J˜Bšœœ˜˜AJ˜—Jšœ˜—J˜ J˜ šœ œœ˜#J˜ J˜Jšœ˜—Jšœ˜—Jšœ˜—˜Jšœ-™-Jšœ˜šœ ˜J˜JšœE™EJ˜7šœ œœ˜#J˜EJšœ˜—Jšœ˜—šœ˜Jšœœ˜ ˜UJ˜/—šœœ˜J˜,Jšœ˜—J˜Jšœ(™(Jšœ˜—Jšœ˜—˜JšœJ™Jšœœ˜J˜E˜6J˜4—Jšœ˜—JšœN™NJ˜M˜BJ˜8——˜Jšœ˜—J˜Jšœ™Jšœ ˜Jšœ˜J˜Jšœ˜J˜J˜—…—`Ø