-- CGSplineImpl.mesa -- Last edit by Doug Wyatt, August 23, 1982 4:16 pm -- Derived from RegularALSpline.mesa, by -- P.Baudelaire, M.Stone, September 26, 1980 5:45 PM DIRECTORY CGSpline USING [Coeffs, Dim, ErrorType, Knots, SplineType], CGStorage USING [pZone, qZone]; CGSplineImpl: CEDAR PROGRAM IMPORTS CGStorage EXPORTS CGSpline = BEGIN OPEN CGSpline; repZone: ZONE = CGStorage.qZone; arrayZone: ZONE = CGStorage.pZone; Error: PUBLIC ERROR[type: ErrorType] = CODE; Array: TYPE = REF ArrayRep; ArrayRep: TYPE = RECORD[SEQUENCE space: NAT OF REAL]; CoeffsArray: TYPE = REF CoeffsArrayRep; CoeffsArrayRep: TYPE = RECORD[SEQUENCE space: NAT OF Coeffs]; Ref: TYPE = REF Rep; Rep: PUBLIC TYPE = RECORD [ space: NAT, h,a,b,c,d,r,s: Array _ NIL, coeffs: CoeffsArray _ NIL ]; New: PUBLIC PROC[size: NAT] RETURNS[Ref] = { self: Ref _ repZone.NEW[Rep _ [space: 0]]; Grow[self,size]; RETURN[self]; }; Grow: PROC[self: Ref, size: NAT] = { IF size<=self.space THEN RETURN; self.h _ arrayZone.NEW[ArrayRep[size]]; self.a _ arrayZone.NEW[ArrayRep[size]]; self.b _ arrayZone.NEW[ArrayRep[size]]; self.c _ arrayZone.NEW[ArrayRep[size]]; self.d _ arrayZone.NEW[ArrayRep[size]]; self.r _ arrayZone.NEW[ArrayRep[size]]; self.s _ arrayZone.NEW[ArrayRep[size]]; self.coeffs _ arrayZone.NEW[CoeffsArrayRep[size]]; self.space _ size; }; --this procedure implements the natural spline algorithm --with polygonal line length parametrization --it expects splineType to be: naturalAL or cyclicAL MakeSpline: PUBLIC PROCEDURE [self: Ref, knots: Knots, splineType: SplineType, Proc: PROC[Coeffs]] = BEGIN xyzw: Dim; 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 h, a, b, c, d, r, & s) size: CARDINAL _ (lastK+1); k: CARDINAL; h: Array _ NIL; a: Array _ NIL; b: Array _ NIL; c: Array _ NIL; d: Array _ NIL; r: Array _ NIL; s: Array _ NIL; coeffs: CoeffsArray _ NIL; deltaX, deltaY: REAL; --first some error checks IF numKnots <= 0 THEN ERROR Error[tooFewKnots]; -- IF splineType # naturalAL AND splineType # cyclicAL THEN -- SIGNAL UnknownSpline[splineType]; IF splineType = cyclicAL THEN FOR xyzw IN Dim DO IF knots[0][xyzw] # knots[lastKnot][xyzw] THEN GOTO notCyclic; REPEAT notCyclic => ERROR Error[unmatchedEnds]; ENDLOOP; --special case for 1 and 2 points IF numKnots < 3 THEN BEGIN coeffs: ARRAY [0..1) OF Coeffs; 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; Proc[coeffs[0]]; RETURN; END; --allocate 4 or 7 arrays of reals for spline computation --and 1 array of cubic coefficients to return IF size>self.space THEN Grow[self,size]; coeffs _ self.coeffs; h _ self.h; a _ self.a; b _ self.b; d _ self.d; IF splineType = cyclicAL THEN BEGIN c _ self.c; r _ self.r; s _ self.s; END; --first get constant cubic coefficients FOR k IN [0..lastK] DO coeffs[k].t0 _ knots[k]; ENDLOOP; --compute parametrization (polygonal line approximation) --this assumes NDIM = 2 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); --check overlapping knots! IF h[k] = 0 THEN h[k] _ 1; ENDLOOP; SELECT splineType FROM naturalAL => BEGIN --natural end conditions: precompute coefficients a 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; END; cyclicAL => BEGIN --cyclic end conditions: precompute coefficients a & c 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; END; ENDCASE => ERROR Error[bug]; FOR xyzw IN Dim DO --compute each dimension separately --compute coefficient d 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 => BEGIN --natural end conditions: compute coefficients b 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; END; cyclicAL => BEGIN --cyclic end conditions: compute coefficients b, r, & s 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; END; ENDCASE; --compute coefficient t2 (=second derivative/2) SELECT splineType FROM naturalAL => BEGIN coeffs[0].t2[xyzw] _ 0; --note: coeffs[lastK+1].t2.[xyzw] is 0 too but it is not stored, hence: 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; END; cyclicAL => BEGIN 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])); --note: coeffs[lastK+1].t2.[xyzw] = coeffs[0].t2[xyzw] coeffs[lastK].t2[xyzw] _ p2/2; FOR k IN [0..lastK-1] DO coeffs[k].t2[xyzw] _ (r[k] * p2 + s[k]) / 2; ENDLOOP; END; ENDCASE => ERROR Error[bug]; --compute coefficients t3 (=third derivative/6) & and t1 (=first derivative) 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; --note again: coeffs[lastK+1].t2[xyzw] = coeffs[0].t2[xyzw] (=0 for natural end conditions) 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; --now reparametrize coefficients to interval [0..1] 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; --release computation arrays -- SplineMemoryDefs.Free[BASE[h]]; -- SplineMemoryDefs.Free[BASE[a]]; -- SplineMemoryDefs.Free[BASE[b]]; -- SplineMemoryDefs.Free[BASE[d]]; -- IF splineType = cyclicAL THEN -- BEGIN -- SplineMemoryDefs.Free[BASE[c]]; -- SplineMemoryDefs.Free[BASE[r]]; -- SplineMemoryDefs.Free[BASE[s]]; -- END; --return cubic coefficients FOR i: CARDINAL IN[0..size) DO Proc[coeffs[i]] ENDLOOP; -- RETURN[coeffs]; END; END. Κ)– "Mesa" style˜IprocšφΟcIœ_œΟk œ žœ:žœ žœžœžœ žœ žœžœžœžœžœžœžœ žœžœžœžœžœžœžœžœžœžœ!žœžœžœžœžœžœžœ žœžœžœ žœžœžœΟnœžœžœžœžœ žœ'žœŸœžœžœžœžœžœžœ'žœ'žœ'žœ'žœ'žœ'žœ,žœ3›œŸ œžœž œ5Ÿœžœ žœžœœ žœœžœGœžœžœ žœ žœ žœ žœ žœ žœ žœžœžœžœžœžœbžœžœžœžœžœžœ)žœžœ žœžœžœ"žœžœžœ žœžœYžœ žœžœžœ žœ žœžœ žœ žœžœžœžœgžœžœXžœžœžœ)žœ(žœžœ žœžœQžœžœ žœžœ žœžœ žœ(žœ žœžœ žœžœžœžœœžœ žœ žœžœ žœžœ4œžœžœžœ8žœžœžœ7œ žœžœžœ<žœžœžœžœ(žœžœžœžœžœžœžœ$œœžœžœ žœ8žœžœžœžœ žœžœ1œžœžœžœ#žœžœžœžœ-žœžœžœ8œ"žœžœžœ#žœžœžœžœ/žœ$žœž œžœžœXžœžœžœ0œžœ žœžœHœ<žœž œžœžœSžœžœžœžœ©7œ#žœžœžœ4žœžœžœžœMœžœžœžœ£žœ\œΆ4œžœžœ žœ§žœžœΙœžœžœžœ žœžœžœžœ˜φ6—…—x§