CurvatureImpl.mesa
M. Stone February 13, 1985 4:41:55 pm PST
DIRECTORY
Filters,
Seq,
Curvature;
CurvatureImpl: CEDAR PROGRAM = {
DValues:
PROC = {
i: INTEGER ← 0;
n: INTEGER ← values.l;
IF n<2 THEN RETURN;
maxDValue ← Real.SmallestNormalizedNumber;
minDValue ← Real.LargestNumber;
dValues ← NEW[ValueSeq[n]];
FOR i IN [0..n-1) DO
dvalue.v is difference between adjacent values. dvalue.s is half way between
the s of adjacent values
dValues[i].v ← values[i+1].v-values[i].v;
dValues[i].s ← (values[i+1].s+values[i].s)/2;
IF dValues[i].v >maxDValue THEN maxDValue ← dValues[i].v;
IF dValues[i].v <minDValue THEN minDValue ← dValues[i].v;
ENDLOOP;
dValues[n-1].v ← values[0].v-values[n-1].v;
don't have the information to do this calculation exactly right
dValues[n-1].s ← values[n-1].s;
IF dValues[n-1].v >maxDValue THEN maxDValue ← dValues[n-1].v;
IF dValues[n-1].v <minDValue THEN minDValue ← dValues[n-1].v;
PushReal[maxDValue];
PushReal[minDValue];
};
this is angle/(m0+m1)/2 in degrees
CheckCurvature: SampleProc = {
v0: Vec ← Sub[sa.xy,Prev[sa].xy];
v1: Vec ← Sub[Next[sa].xy,sa.xy];
m0: REAL ← ABS[sa.s-Prev[sa].s];
m1: REAL ← ABS[Next[sa].s-sa.s];
x: REAL ← Dot[v0,v1];
y: REAL ← Cross[v0,v1];
sa.value ← RealFns.ArcTanDeg[y,x];
sa.value ← .5*sa.value/(m0+m1);
};
this is angle/(m0+m1)/2 in radians
CCurveAngle: SampleProc = {
v0: Vec ← Sub[sa.xy,Prev[sa].xy];
v1: Vec ← Sub[Next[sa].xy,sa.xy];
m0: REAL ← ABS[sa.s-Prev[sa].s];
m1: REAL ← ABS[Next[sa].s-sa.s];
x: REAL ← Dot[v0,v1];
y: REAL ← Cross[v0,v1];
sa.value ← RealFns.ArcTan[y,x];
sa.value ← .5*sa.value/(m0+m1);
};
FindCircleCenter:
PROCEDURE[p0,p1,p2:Vector.Vec]
RETURNS [o:Vector.Vec,valid:
BOOLEAN]=
{d1,d2:Vec;
d20,d21,det:REAL;
x20,y20,x21,y21,x22,y22:REAL;
eps:REAL𡤀.001;
d1←Mul[Sub[p0,p1],2];
d2←Mul[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�t[Matrix[d1.x,d1.y,d2.x,d2.y]];
IF ABS [det] < eps THEN
{RETURN[Vec[0,0],FALSE]} ELSE
{o.x←(d20*d2.y-d21*d1.y)/det; o.y←(d1.x*d21-d2.x*d20)/det;
RETURN[o,TRUE];};
};
ArcLength: SampleProc = {
sa.value ← (Next[sa].s+sa.s)/2;
};
}.