<> <> <> <> <> 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]; }; }.