-- CubicImpl.mesa
-- Last changed by Doug Wyatt, September 7, 1980 12:56 PM

DIRECTORY
Cubic,
RealFns USING [SqRt],
Vector USING [Vec, Add, Sub, Mul, Div, Min, Max, In, Dot],
Memory USING [zone];

CubicImpl: PROGRAM
IMPORTS Memory,Vector,RealFns
EXPORTS Cubic = {
OPEN Cubic,Vector;

zone: UNCOUNTED ZONE = Memory.zone;

CCoeffsToBezier: PUBLIC PROCEDURE[c: POINTER TO Coeffs,
b: POINTER TO Bezier] = {
OPEN Vector,c,b;
-- Compute the bezier control points of the cubic
-- These points form a convex hull of the curve
-- b0 is the starting point, b3 is the ending point
-- b0 ← c0;
-- b1 ← c0 + c1/3;
-- b2 ← c0 + 2*c1/3 + c2/3; [ = b1 + (c1+c2)/3 ]
-- b3 ← c0 + c1 + c2 + c3;
b0←c0;
b1←Add[c0,Div[c1,3]];
b2←Add[b1,Div[Add[c1,c2],3]];
b3←Add[Add[Add[c0,c1],c2],c3];
};

CBezierToCoeffs: PUBLIC PROCEDURE[b: POINTER TO Bezier,
c: POINTER TO Coeffs] = {
OPEN Vector,b,c;
temp: Vec;
-- Compute the cubic coefficients from the Bezier points.
-- c0 ← b0;
-- c1 ← 3(b1-b0);
-- c2 ← 3(b0-2b1+b2); [ = 3(b2-b1) - c1 ]
-- c3 ← -b0 +3(b1-b2) + b3; [ = b3 - (b0 + 3(b2-b1)) ]
c0 ← b0;
c1 ← Mul[Sub[b1,b0],3];
c2 ← Sub[(temp←Mul[Sub[b2,b1],3]),c1];
c3 ← Sub[b3,Add[b0,temp]];
};

CBezierPolygon: PUBLIC PROCEDURE[b: POINTER TO Bezier, epsilon: REAL,
Proc: PROCEDURE[Vector.Vec]] = {
lvlmax: CARDINAL=10;
SubDivide: PROCEDURE[b: POINTER TO Bezier, lvl: CARDINAL] = {
IF lvl>=lvlmax OR FlatTest[b,epsilon] THEN Proc[b.b3]
ELSE {
OPEN Vector;
q1,q2,q3,qp1,qp2,q: Vec;
bb: Bezier;
q1←Mul[Add[b.b0,b.b1],0.5];
q2←Mul[Add[b.b1,b.b2],0.5];
q3←Mul[Add[b.b2,b.b3],0.5];
qp1←Mul[Add[q1,q2],0.5];
qp2←Mul[Add[q2,q3],0.5];
q←Mul[Add[qp1,qp2],0.5];
bb←[b0: b.b0, b1: q1, b2: qp1, b3: q];
SubDivide[@bb,lvl+1];
bb←[b0: q, b1: qp2, b2: q3, b3: b.b3];
SubDivide[@bb,lvl+1];
};
};
SubDivide[b,1];
};

CBezierRelax: PUBLIC PROCEDURE[b: POINTER TO Bezier, costheta: REAL,
Proc: PROCEDURE[Bezier]] = {
lvlmax: CARDINAL=10;
SubDivide: PROCEDURE[b: POINTER TO Bezier, lvl: CARDINAL] = {
IF lvl>=lvlmax OR ObliqueTest[b,costheta] THEN Proc[b↑]
ELSE {
OPEN Vector;
q1,q2,q3,qp1,qp2,q: Vec;
bb: Bezier;
q1←Mul[Add[b.b0,b.b1],0.5];
q2←Mul[Add[b.b1,b.b2],0.5];
q3←Mul[Add[b.b2,b.b3],0.5];
qp1←Mul[Add[q1,q2],0.5];
qp2←Mul[Add[q2,q3],0.5];
q←Mul[Add[qp1,qp2],0.5];
bb←[b0: b.b0, b1: q1, b2: qp1, b3: q];
SubDivide[@bb,lvl+1];
bb←[b0: q, b1: qp2, b2: q3, b3: b.b3];
SubDivide[@bb,lvl+1];
};
};
SubDivide[b,1];
};

FlatTest: PROCEDURE[b: POINTER TO Bezier, eps: REAL]
RETURNS[BOOLEAN] = INLINE {
OPEN Vector;
dx,dy: REAL;
d1,d2,d,bl,bh: Vec;
oh: Vec=[0.5,0.5];
bh←Add[Max[b.b0,b.b3],oh];
bl←Sub[Min[b.b0,b.b3],oh];
IF ~In[b.b1,bl,bh] OR ~In[b.b2,bl,bh] THEN RETURN[FALSE];
d1←Sub[b.b1,b.b0];
d2←Sub[b.b2,b.b0];
d←Sub[b.b3,b.b0];
dx←ABS[d.x]; dy←ABS[d.y];
IF dx+dy < ceps THEN RETURN[TRUE];
IF dy < dx THEN RETURN[
ABS[d2.y-((d.y*d2.x)/d.x)]<ceps AND ABS[d1.y-((d.y*d1.x)/d.x)]<ceps]
ELSE RETURN[
ABS[d2.x-((d.x*d2.y)/d.y)]<ceps AND ABS[d1.x-((d.x*d1.y)/d.y)]<ceps];
};

ObliqueTest: PROCEDURE[b: POINTER TO Bezier, cang: REAL]
RETURNS[BOOLEAN] = INLINE {
OPEN Vector;
d1,d2,d3:Vec;
d1←Norm[Sub[b.b1,b.b0]];
d2←Norm[Sub[b.b2,b.b1]];
d3←Norm[Sub[b.b3,b.b2]];
RETURN [Dot[d1,d2] > cang AND Dot[d2,d3] > cang];
};

SetEpsilon:PUBLIC PROCEDURE[e:REAL]=
{ceps←e;};

Mag:PROCEDURE[v:Vec] RETURNS [m:REAL]=
BEGIN
m←RealFns.SqRt[v.x*v.x+v.y*v.y];
END;

Norm:PROCEDURE[v:Vec] RETURNS [nv:Vec]=
BEGIN
m:REAL;
IF (m←Mag[v]) = 0 THEN RETURN[Vec[0,0]];
nv.x←v.x/m; nv.y←v.y/m;
END;



ceps:REAL←1.5;

}.