-- CGCubicImpl.mesa
-- Last changed by Doug Wyatt, August 23, 1982 4:12 pm

DIRECTORY
CGCubic USING [Bezier, Coeffs],
CGVector USING [Add, Sub, Mul, Div, Min, Max, In],
GraphicsBasic USING [Vec];

CGCubicImpl: CEDAR PROGRAM
IMPORTS CGVector
EXPORTS CGCubic = {
OPEN Vector: CGVector, CGCubic, GraphicsBasic;

CoeffsToBezier: PUBLIC PROC[coeffs: Coeffs] RETURNS[bezier: Bezier] = {
OPEN Vector;
-- 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;
bezier.b0 ← coeffs.c0;
bezier.b1 ← Add[coeffs.c0,Div[coeffs.c1,3]];
bezier.b2 ← Add[bezier.b1,Div[Add[coeffs.c1,coeffs.c2],3]];
bezier.b3 ← Add[Add[Add[coeffs.c0,coeffs.c1],coeffs.c2],coeffs.c3];
RETURN[bezier];
};

BezierToCoeffs: PUBLIC PROC[bezier: Bezier] RETURNS[coeffs: Coeffs] = {
OPEN Vector;
-- 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)) ]
t: Vec ← Mul[Sub[bezier.b2,bezier.b1],3];
coeffs.c0 ← bezier.b0;
coeffs.c1 ← Mul[Sub[bezier.b1,bezier.b0],3];
coeffs.c2 ← Sub[t,coeffs.c1];
coeffs.c3 ← Sub[bezier.b3,Add[bezier.b0,t]];
RETURN[coeffs];
};

Split: PUBLIC PROC[bezier: Bezier] RETURNS[Bezier,Bezier] = {
a,b,c,ab,bc,p: Vec;
Mid: PROC[u,v: Vec] RETURNS[Vec] = INLINE { RETURN[[(u.x+v.x)/2,(u.y+v.y)/2]] };
a ← Mid[bezier.b0,bezier.b1];
b ← Mid[bezier.b1,bezier.b2];
c ← Mid[bezier.b2,bezier.b3];
ab ← Mid[a,b];
bc ← Mid[b,c];
p ← Mid[ab,bc];
RETURN[[bezier.b0, a, ab, p],[p, bc, c, bezier.b3]];
};

Flat: PUBLIC PROC[bezier: Bezier, epsilon: REAL] RETURNS[BOOLEAN] = {
OPEN Vector;
dx,dy: REAL;
d1,d2,d,bl,bh: Vec;
oh: Vec=[0.5,0.5];
bh ← Add[Max[bezier.b0,bezier.b3],oh];
bl ← Sub[Min[bezier.b0,bezier.b3],oh];
IF ~In[bezier.b1,bl,bh] OR ~In[bezier.b2,bl,bh] THEN RETURN[FALSE];
d1 ← Sub[bezier.b1,bezier.b0];
d2 ← Sub[bezier.b2,bezier.b0];
d ← Sub[bezier.b3,bezier.b0];
dx ← ABS[d.x]; dy ← ABS[d.y];
IF dx+dy < 1 THEN RETURN[TRUE];
IF dy < dx THEN {
dydx: REAL ← d.y/d.x;
RETURN[ABS[d2.y-d2.x*dydx]<epsilon AND ABS[d1.y-d1.x*dydx]<epsilon] }
ELSE {
dxdy: REAL ← d.x/d.y;
RETURN[ABS[d2.x-d2.y*dxdy]<epsilon AND ABS[d1.x-d1.y*dxdy]<epsilon] };
};

}.