CubicImpl.mesa
Copyright Ó 1985, 1992 by Xerox Corporation. All rights reserved.
Last changed by Maureen Stone, 9-Sep-81 16:21:17
Doug Wyatt, September 5, 1985 1:31:43 pm PDT
DIRECTORY
Cubic USING [Bezier, Coeffs],
Vector USING [Add, Div, In, Max, Min, Mul, Sub, VEC];
CubicImpl: CEDAR PROGRAM
IMPORTS Vector
EXPORTS Cubic = {
OPEN Cubic;
VEC: TYPE ~ Vector.VEC;
CoeffsToBezier: PUBLIC PROCEDURE[c: Coeffs] RETURNS [b: Bezier] = {
OPEN Vector, c;
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;
b.b0¬c0;
b.b1¬Add[c0,Div[c1,3]];
b.b2¬Add[b.b1,Div[Add[c1,c2],3]];
b.b3¬Add[Add[Add[c0,c1],c2],c3];
};
BezierToCoeffs: PUBLIC PROCEDURE[b: Bezier] RETURNS [c: Coeffs] = {
OPEN Vector,b;
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)) ]
c.c0 ¬ b0;
c.c1 ¬ Mul[Sub[b1,b0],3];
c.c2 ¬ Sub[(temp¬Mul[Sub[b2,b1],3]),c.c1];
c.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: PUBLIC PROCEDURE[b: Bezier, eps: REAL]
RETURNS[BOOLEAN] = {
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];
};
}.