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