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 M. Stone May 25, 1985 3:58:09 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]]]; }; }; Êü˜Jšœ™Jšœ%™%šÏk ˜ J˜J˜J˜J˜Jšœœ˜!J˜—Jšœœ˜ Jšœ˜Jšœœ˜,J˜Jšœœ˜šÏnœ œ(œœ˜nJšœœ"˜8šœœœœ˜CJ˜J˜J˜%J˜$Jšœ˜J˜J˜—J˜šœœœ˜&J˜_J˜Jšœ˜—šœœ˜Jšœœ˜J˜^J˜˜2J˜:—J˜J˜—šœ˜J˜.JšœL˜LJ˜—Jšœ˜ ˜J˜——šžœ œœ˜LJšœ œ!˜-Jš œœœœ$œ˜Qšœ œ˜J˜&J˜0J˜—Jšœ$˜(J˜—J˜J˜Jšœ-œb™’Jšœ œÏc;˜\š žœ œœœ œ˜KJšœœœœ˜-š˜Jšœœ œœ ˜&šœ˜Jšœœ œ ˜Jšœœœœ˜,Jšœ˜——Jšœ˜—š ž œ œœœœ˜Lš œœœœœ˜(Jšœ5˜;—J˜J˜Jš œœ œ œœ ˜3J˜šœ œ˜Jš œœœœœ˜:J˜ J˜—Jš œœœ œœ˜3J˜J˜—J™EJ™OJšŸœ<™QJš œœœœœ™=Jšœ œœ˜š ž œ œœœ œŸ˜mJšœœ ˜š˜Jšœ+œœ˜EJšœ˜JšœœœŸ˜CJšœ˜—Jšœ!™!Jšœ ˜ š˜Jš œœ+œœ˜BJšœ˜JšœœœŸ˜LJšœ˜—JšœŸ™ŸJšœ ˜ J˜—J˜Jš ž œœœœœ™:Jšœ(˜(šžœ œ˜.Jšœ˜J˜—š žœœœœœ˜FJšœA˜GJ˜—šÏb œœ˜'Jšœœ&˜.Jšœœ"˜)Jšœ˜ J˜—š  œœ˜'Jšœœ&˜.JšœœB˜IJšœ˜ J˜—š  œœ˜&J˜ Jšœœ˜J˜4Jšœœœœ˜0J˜—šžœœœœ˜VJ˜ Jšœ œ˜Jšœœ˜Jšœœ˜J˜#J˜#J˜J˜J˜J˜J˜J˜3Jšœœ ˜Jšœœ œ˜JšœD˜DJšœœ˜J˜J˜—J˜šž œœ œœ™QJ™!J™!Jšœœ™J™ JšœŸ'™FJ™\J™MJ™%J™&J™&J™J™J™šœœŸ™?J™"J™—šœ™J™J™J™J™J™5J™—J™J™—J˜—…—Œ