DIRECTORY Seq, SampledCurves, Complex, Vector, RealFns USING[ArcTanDeg, ArcTan]; SampledCurvesImpl: CEDAR PROGRAM IMPORTS Vector, RealFns EXPORTS SampledCurves = {OPEN SampledCurves; VEC: TYPE = Vector.VEC; SampledCurveFromSamples: PUBLIC PROC [samples: Seq.ComplexSequence, closed: BOOLEAN] RETURNS[SampledCurve] = { sc: SampledCurve _ NEW[SampledCurveRec[samples.length]]; computeValues: PROC[p0,p1,p2: VEC] RETURNS[t: VEC, ds, k: REAL] = { d1: VEC _ Vector.Sub[p1,p0]; d2: VEC _ Vector.Sub[p2,p1]; ds _ Vector.Mag[d1] + Vector.Mag[d2]; t _ Vector.Unit[Vector.Add[d1, d2]]; k _ ComputeK[d1,d2]; }; sc.closed _ closed; FOR i: NAT IN (0..samples.length-1) DO [t: sc[i].t, ds: sc[i].ds, k: sc[i].k] _ computeValues[samples[i-1], samples[i], samples[i+1]]; sc[i].p _ samples[i]; ENDLOOP; IF closed THEN { last: NAT _ samples.length-1; [t: sc[0].t, ds: sc[0].ds, k: sc[0].k] _ computeValues[samples[last], samples[0], samples[1]]; sc[0].p _ samples[0]; [t: sc[last].t, ds: sc[last].ds, k: sc[last].k] _ computeValues[samples[last-1], samples[last], samples[0]]; sc[last].p _ samples[last]; } ELSE { sc[0] _ [p: samples[0], t: [0,0], k: 0, ds:0]; sc[samples.length-1] _ [p: samples[samples.length-1], t: [0,0], k: 0, ds:0]; }; RETURN[sc]; }; DeltaK: PUBLIC PROC[sc: SampledCurve] RETURNS [deltas: Seq.RealSequence] = { deltas _ NEW[Seq.RealSequenceRec[sc.length]]; FOR i: NAT IN (0..deltas.length-1) DO deltas[i] _ sc[i+1].k - sc[i-1].k; ENDLOOP; IF sc.closed THEN { deltas[0] _ sc[1].k-sc[sc.length-1].k; deltas[sc.length-1] _ sc[0].k-sc[sc.length-2].k; } ELSE deltas[0] _deltas[sc.length-1] _ 0; }; OpenCurve: PUBLIC SIGNAL = CODE; --raised when you try an illegal operation on an open curve Wrap: PUBLIC PROC[sc: SampledCurve, index: INT] RETURNS [newIndex: NAT] = { IF index IN [0..sc.length) THEN RETURN[index] ELSE IF NOT sc.closed THEN SIGNAL OpenCurve ELSE { new: INT _ index MOD sc.length; RETURN[IF new<0 THEN sc.length+new ELSE new] }; }; ArcLength: PUBLIC PROC[sc: SampledCurve, from, to: NAT] RETURNS[d: REAL] = { dist: PROC[index: NAT] RETURNS[REAL] = { RETURN[Vector.Mag[Vector.Sub[sc[index].p,sc[index-1].p]]]}; from _ Wrap[sc, from]; to _ Wrap[sc, to]; IF NOT sc.closed AND to SIGNAL None]]] THEN EXIT ELSE index _ index+1; IF index-start > sc.length THEN SIGNAL None; --never found anything ENDLOOP; ft _ index; DO IF NOT condition[sc[Wrap[sc,index ! OpenCurve => EXIT]]] THEN EXIT ELSE index _ index+1; IF index-start > sc.length THEN {index _ start; EXIT}; --whole curve is TRUE ENDLOOP; lt _ index-1; }; ComputeK: CurvatureProc _ DThetaDegrees; SetKProc: PUBLIC PROC[proc: CurvatureProc] = { ComputeK _ proc; }; FindAngleDegrees: PUBLIC PROC[dIn,dOut: Vector.VEC] RETURNS [REAL] = { RETURN[RealFns.ArcTanDeg[Vector.Cross[dIn,dOut],Vector.Dot[dIn,dOut]]]; }; DThetaDegrees: PUBLIC CurvatureProc = { ds: REAL _ Vector.Mag[dIn] + Vector.Mag[dOut]; k: REAL _ FindAngleDegrees[dIn, dOut]/ds; RETURN[k]; }; DThetaRadians: PUBLIC CurvatureProc = { ds: REAL _ Vector.Mag[dIn] + Vector.Mag[dOut]; k: REAL _ RealFns.ArcTan[Vector.Cross[dIn,dOut],Vector.Dot[dIn,dOut]]/ds; RETURN[k]; }; CircleCenter: PUBLIC CurvatureProc = { center: VEC; valid: BOOLEAN; [center, valid] _ FindCircleCenter[[0,0],dIn, dOut]; RETURN[IF valid THEN Vector.Mag[center] ELSE 0]; }; FindCircleCenter: PUBLIC PROCEDURE[p0,p1,p2: VEC] RETURNS [center: VEC,valid:BOOLEAN]= {d1,d2:VEC; d20,d21,det:REAL; x20,y20,x21,y21,x22,y22:REAL; eps:REAL_0.001; d1_Vector.Mul[Vector.Sub[p0,p1],2]; d2_Vector.Mul[Vector.Sub[p1,p2],2]; x20_p0.x*p0.x; y20_p0.y*p0.y; x21_p1.x*p1.x; y21_p1.y*p1.y; x22_p2.x*p2.x; y22_p2.y*p2.y; d20_x20-x21+y20-y21; d21_x21-x22+y21-y22; det_Vector.Det[Vector.Matrix[d1.x,d1.y,d2.x,d2.y]]; IF ABS [det] < eps THEN {RETURN[VEC[0,0],FALSE]} ELSE {center.x_(d20*d2.y-d21*d1.y)/det; center.y_(d1.x*d21-d2.x*d20)/det; RETURN[center,TRUE];}; }; }. SampledCurvesImpl.mesa Copyright c 1985 by Xerox Corporation. All rights reserved. M. Stone May 25, 1985 3:58:09 pm PDT Doug Wyatt, September 5, 1985 2:39:48 pm PDT Wrap takes an index and returns the newIndex MOD curve length if the curve is closed, does the end test and raises OpenCurve if the curve is open. This finds a sequence of TRUE values as defined by the ConditionProc. ft is the index of the first TRUE value, lt is the index of the last TRUE value Uses unwrapped values ie. start, ft, lt can be outside the range [0..bool.length) ConditionProc: TYPE = PROC[point: PointRec] RETURNS[BOOLEAN]; now index is the first TRUE value now index is the first FALSE value. Three cases: index is truly first false value; index is sc.length+1 (open curve); index is start (whole curve is marked); CurvatureProc: TYPE = PROC[[dIn,dOut: VEC] RETURNS [REAL]; CircleTan: PROC [a, b, c: Complex.VEC, maxSin: REAL] RETURNS [t: Complex.VEC] = { d1: Vector.VEC _ Vector.Sub[b,a]; d2: Vector.VEC _ Vector.Sub[c,b]; a1,b1,c1,a2,b2,c2,det: REAL; midAB,midBC,center: Complex.VEC; [d1,d2] _ TestCorner[a, b, c]; --will signal TooSharp if it's a corner The center of the circle is the intersection of the two perpendicular bisectors of d1 and d2 ax+by = c is a line. Find two lines, use Cramer's rule to find intersection. Normal to b-center is tangent vector. midAB _ Vector.Div[Vector.Add[a,b],2]; midBC _ Vector.Div[Vector.Add[b,c],2]; a1 _ - d1.x; b1 _ - d1.y; a2 _ - d2.x; b2 _ - d2.y; det _ a1*b2-a2*b1; IF RealFns.AlmostZero[det, -100] THEN { --3 pts nearly colinear t _ Vector.Unit[Complex.Sub[c,a]]; } ELSE { c1 _ a1*midAB.x+b1*midAB.y; c2 _ a2*midBC.x+b2*midBC.y; center.x _ (c1*b2-c2*b1)/det; center.y _ (c2*a1-c1*a2)/det; t _ Vector.Unit[Vector.Normal[Vector.Sub[b,center]]]; }; }; Êz˜codešœ™Kšœ Ïmœ1™