-- CubicImpl.mesa
-- Last changed by Maureen Stone,  9-Sep-81 16:21:17
-- Written by Doug Wyatt, September 22, 1980  5:30 PM

DIRECTORY
	Cubic,
	Vector USING [Vec, Add, Sub, Mul, Div, Min, Max, In];

CubicImpl: CEDAR PROGRAM
IMPORTS Vector
EXPORTS Cubic = {
OPEN Cubic;


CoeffsToBezier: PUBLIC PROCEDURE[c: Coeffs] RETURNS[b: 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];
	};

BezierToCoeffs: PUBLIC PROCEDURE[b: Bezier] RETURNS[c: 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]];
	};

BezierPolygon: PUBLIC PROCEDURE[b: Bezier, epsilon: REAL,Proc: PROCEDURE[Vector.Vec]] = {
	lvlmax: CARDINAL=10;
	SubDivide: PROCEDURE[b: Bezier, lvl: CARDINAL] = {
		IF lvl>=lvlmax OR IsFlat[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];
	};

IsFlat: PROCEDURE[b: 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 < 1 THEN RETURN[TRUE];
	IF dy < dx THEN RETURN[
		ABS[d2.y-((d.y*d2.x)/d.x)]<eps AND ABS[d1.y-((d.y*d1.x)/d.x)]<eps]
	ELSE RETURN[
		ABS[d2.x-((d.x*d2.y)/d.y)]<eps AND ABS[d1.x-((d.x*d1.y)/d.y)]<eps];
	};

}.