
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 := 0] = {
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<from THEN SIGNAL OpenCurve;
d ¬ 0;
IF to<from THEN {
FOR i: NAT IN (from..sc.length) DO d ¬ d+dist[i]; ENDLOOP;
from ¬ 0;
};
FOR i: NAT IN (from..to] DO d ¬ d+dist[i]; ENDLOOP;
};

None: PUBLIC SIGNAL = CODE;
FindSequence: PUBLIC PROC[sc: SampledCurve, start: NAT, condition: ConditionProc] RETURNS[ft, lt: NAT] = {	--
index: NAT ¬ start;
DO
IF condition[sc[Wrap[sc,index ! OpenCurve => 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 Ó 1985, 1992 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]]];
};
};

Êi  •NewlineDelimiter
–(cedarcode) style™ code™Kšœ	Ïeœ6™BK™%K™,—K˜ šÏk	˜	K˜K˜K˜K˜Kšœžœ˜!—K˜ KšÐblœžœž˜ Kšžœ˜Kšžœžœ˜,K˜ Kšžœžœ
žœ˜šÏnœžœ(žœžœ˜nKšœžœ"˜8šœžœžœžœžœ	žœ˜CKšœžœ˜Kšœžœ˜K˜%K˜$K˜K˜K˜ —K˜šžœžœžœž˜&K˜_K˜Kšžœ˜—šžœžœ˜Kšœžœ˜K˜^K˜˜2K˜:—K˜K˜—šžœ˜K˜.K˜LK˜—Kšžœ˜˜K˜ ——š œžœžœ˜LKšœ	žœ!˜-Kš
žœžœžœžœ$žœ˜Qšžœžœ˜K˜&K˜0K˜—Kšžœ$˜(K˜—K˜ K˜ Kšœ-žœb™’KšœžœÏc;˜\š œžœžœžœ˜PKšžœžœžœžœ˜-šž˜Kšžœžœžœžœ
˜&šžœ˜Kšœžœ	žœ˜Kšžœžœžœžœ˜,K˜——K˜—š
 	œžœžœžœžœ˜Lš	œžœžœžœžœ˜(Kšžœ5˜;—K˜K˜Kš
žœžœžœ	žœžœ˜3K˜šžœ	žœ˜Kš
žœžœžœžœžœ˜:K˜	K˜—Kš
žœžœžœžœžœ˜3K˜K˜ —K™EK™OKš¡œ<™QKš	œžœžœžœžœ™=Kšœžœžœ˜š œžœžœžœ	žœ¡˜mKšœžœ	˜šž˜Kšžœ+žœžœž˜EKšžœ˜Kšžœžœžœ¡˜CKšžœ˜—K™!K˜šž˜Kš	žœžœ+žœžœž˜BKšžœ˜Kšžœžœžœ¡˜LKšžœ˜—K™ŸK˜K˜—K˜ Kš œžœžœžœžœžœ™:K˜(š œžœ˜.K˜K˜—š œžœžœžœžœžœ˜FKšžœA˜GK˜—šÏbœžœ˜'Kšœžœ&˜.Kšœžœ"˜)Kšžœ˜
K˜—š¢œžœ˜'Kšœžœ&˜.KšœžœB˜IKšžœ˜
K˜—š¢œžœ˜&Kšœžœ˜Kšœžœ˜K˜4Kšžœžœžœžœ˜0K˜—š œžœžœžœ
žœžœ˜VKšœžœ˜Kšœžœ˜Kšœžœ˜Kšœžœ˜K˜#K˜#K˜K˜K˜K˜K˜K˜3Kšžœžœž˜Kšœžœžœžœž˜K˜DKšžœžœ˜K˜K˜ —K˜š 	œžœžœ
žœžœžœ™QKšœžœ™!Kšœžœ™!Kšœžœ™Kšœžœ™ Kšœ¡'™FK™\K™MK™%K™&K™&K™K™K™šžœžœ¡™?K™"K™—šžœ™K™K™K™K™K™5K™—K™K™ —K˜ — …—    ’    