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
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;
};
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.
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;
};
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];
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;
now index is the first TRUE value
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;
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);
lt ¬ index-1;
};
CurvatureProc: TYPE = PROC[[dIn,dOut: VEC] RETURNS [REAL];
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];};
};
}.
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]]];
};
};