Cubic2Impl.mesa
Copyright Ó 1982, 1985, 1987 by Xerox Corporation. All rights reserved.
Last changed by Doug Wyatt, August 23, 1982 4:12 pm
Eric Nickell, June 27, 1985 10:16:03 pm PDT
Stone, June 29, 1985 10:40:13 am PDT
Russ Atkinson (RRA) February 2, 1987 10:18:12 pm PST
Bier, May 11, 1989 10:27:02 pm PDT
DIRECTORY
Cubic2, Cubic2Extras, Vector2, Imager;
Cubic2Impl:
CEDAR
PROGRAM
IMPORTS Vector2
EXPORTS Cubic2, Cubic2Extras = {
OPEN Vector2, Cubic2;
VEC: TYPE = Imager.VEC;
CoeffsToBezier:
PUBLIC
PROC [coeffs: Coeffs]
RETURNS [bezier: Bezier] = {
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] = {
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];
};
AlphaSplit:
PUBLIC
PROC [bezier: Bezier, alpha:
REAL]
RETURNS[Bezier,Bezier] = {
Splits one bezier cubic into two bezier cubics at the point u = alpha, where u is the standard bezier parameter.
Fraction:
PROC[q,r:
VEC, alpha:
REAL, beta:
REAL]
RETURNS [
VEC] = {
RETURN[[q.x*beta+r.x*alpha, q.y*beta+r.y*alpha]] };
a,b,c,ab,bc,p: VEC;
oneMinusAlpha: REAL;
oneMinusAlpha ← 1.0-alpha;
a ← Fraction[bezier.b0, bezier.b1, alpha, oneMinusAlpha];
b ← Fraction[bezier.b1, bezier.b2, alpha, oneMinusAlpha];
c ← Fraction[bezier.b2, bezier.b3, alpha, oneMinusAlpha];
ab ← Fraction[a, b, alpha, oneMinusAlpha];
bc ← Fraction[b, c, alpha, oneMinusAlpha];
p ← Fraction[ab, bc, alpha, oneMinusAlpha];
RETURN[[bezier.b0, a, ab, p],[p, bc, c, bezier.b3]];
};
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] = {
Replaced by Bier, February 6, 1989. The new version using epsilon instead of 0.5 for the first rejection test, and has more illuminating comments.
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];
};
};
[Artwork node; type 'ArtworkInterpress on' to command tool]
Flat:
PUBLIC
PROC [bezier: Bezier, epsilon:
REAL]
RETURNS [
BOOL] = {
A Bezier segment is considered flat if it can be approximated by a straight line from its first control point to its last with a guaranteed maximum error less than epsilon.
Trivial Rejection 1: If b1 or b2 are more than epsilon from the bounding box of b0 and b3, then we assume the curve is not flat.
Trivial Rejection 2: If b0 is within epsilon of b3 then assume the curve is flat. This relies on the assumption that the curve is not a very high frequency one in this region.
dx,dy: REAL;
d1,d2,d,bl,bh: VEC;
oh:
VEC=[epsilon, epsilon];
Trivial Rejection 1
bh ← Vector2.Add[UpperRight[bezier.b0,bezier.b3],oh];
bl ← Vector2.Sub[LowerLeft[bezier.b0,bezier.b3],oh];
IF ~In[bezier.b1,bl,bh]
OR ~In[bezier.b2,bl,bh]
THEN
RETURN[
FALSE];
Trivial Rejection 2
d ← Vector2.Sub[bezier.b3,bezier.b0];
dx ← ABS[d.x]; dy ← ABS[d.y];
IF dx+dy < epsilon
THEN
RETURN[
TRUE];
Compute Pseudo-Altitudes, A and B
d1 ← Vector2.Sub[bezier.b1,bezier.b0];
d2 ← Vector2.Sub[bezier.b2,bezier.b0];
IF dy < dx
THEN {
Compute Vertical Altitudes as shown in the figure.
dydx: REAL ← d.y/d.x;
RETURN[ABS[d2.y-d2.x*dydx]<epsilon AND ABS[d1.y-d1.x*dydx]<epsilon] }
ELSE {
Compute Horizontal Altitudes (not shown).
dxdy: REAL ← d.x/d.y;
RETURN[ABS[d2.x-d2.y*dxdy]<epsilon AND ABS[d1.x-d1.y*dxdy]<epsilon] };
};
LowerLeft:
PROC[a:
VEC, b:
VEC]
RETURNS[
VEC] =
INLINE {
RETURN[[
MIN[a.x,b.x],
MIN[a.y,b.y]]]};
UpperRight:
PROC[a:
VEC, b:
VEC]
RETURNS[
VEC] =
INLINE {
RETURN[[
MAX[a.x,b.x],
MAX[a.y,b.y]]] };
Min:
PROC [a:
VEC, b:
VEC]
RETURNS [
VEC] =
INLINE {
RETURN[[MIN[a.x,b.x],MIN[a.y,b.y]]];
};
Max:
PROC [a:
VEC, b:
VEC]
RETURNS [
VEC] =
INLINE {
RETURN[[MAX[a.x,b.x],MAX[a.y,b.y]]];
};
In:
PROC [a:
VEC, b:
VEC, c:
VEC]
RETURNS [
BOOLEAN] =
INLINE {
RETURN[a.x IN[b.x..c.x] AND a.y IN[b.y..c.y]];
};
}.