CubicPathsImpl.mesa
Copyright Ó 1985, 1987 by Xerox Corporation. All rights reserved.
Written by: Maureen Stone August 5, 1985 8:49:21 pm PDT
Kurlander August 15, 1986 10:08:11 am PDT
Russ Atkinson (RRA) February 2, 1987 10:44:17 pm PST
Pier, August 5, 1987 3:50:54 pm PDT
Bier, February 6, 1989 4:44:55 pm PST
DIRECTORY
Cubic2, CubicPaths, CubicPathsExtras, CubicSplines, Imager, ImagerPath, ImagerTransformation, Polynomial, Real, RealFns, Vector2;
CubicPathsImpl:
CEDAR
PROGRAM
IMPORTS Cubic2, ImagerTransformation, Polynomial, RealFns, Vector2
EXPORTS CubicPaths, CubicPathsExtras
~ BEGIN OPEN CubicPaths;
Bezier: TYPE = Cubic2.Bezier;
Coeffs: TYPE = Cubic2.Coeffs;
ShortRealRootRec: TYPE = Polynomial.ShortRealRootRec;
VEC: TYPE = Imager.VEC;
X: NAT = CubicSplines.X;
Y: NAT = CubicSplines.Y;
PathFromCubic:
PUBLIC
PROC [coeffs: CubicSplines.CoeffsSequence, newBounds:
BOOL ←
TRUE]
RETURNS [Path] = {
path: Path ← NEW[PathRec];
path.cubics ← NEW[PathSequence[coeffs.length]];
FOR i:
NAT
IN [0..coeffs.length)
DO
cfs: Cubic2.Coeffs ← [
c3: [coeffs[i].t3[X], coeffs[i].t3[Y]],
c2: [coeffs[i].t2[X], coeffs[i].t2[Y]],
c1: [coeffs[i].t1[X], coeffs[i].t1[Y]],
c0: [coeffs[i].t0[X], coeffs[i].t0[Y]]];
path.cubics[i] ← Cubic2.CoeffsToBezier[cfs];
ENDLOOP;
IF newBounds THEN UpdateBounds[path];
RETURN [path];
};
EnumeratePath:
PUBLIC PROC [path: Path, moveTo: ImagerPath.MoveToProc, curveTo: ImagerPath.CurveToProc] = {
moveTo[path.cubics[0].b0];
FOR i:
NAT
IN [0..path.cubics.length)
DO
curveTo[path.cubics[i].b1, path.cubics[i].b2, path.cubics[i].b3];
ENDLOOP;
};
WalkPath:
PUBLIC PROC [path: Path, tol:
REAL, proc: PointProc] = {
subdivide:
PROC [bezier: Cubic2.Bezier]
RETURNS [quit:
BOOLEAN]= {
quit ← FALSE;
IF Cubic2.Flat[bezier, tol]
THEN {
len: REAL ← Vector2.Length[Vector2.Sub[bezier.b0, bezier.b3]];
IF len > 0 THEN {
step: REAL ← tol/len;
alpha: REAL ← step;
pt: VEC;
UNTIL alpha > 1
DO
pt ← Vector2.Add[Vector2.Mul[bezier.b3,alpha],Vector2.Mul[bezier.b0,1-alpha]];
IF proc[pt] THEN {quit ← TRUE; EXIT};
alpha ← alpha+step;
ENDLOOP;
}
ELSE quit ← proc[bezier.b3];
RETURN [quit]}
ELSE {
b1, b2: Cubic2.Bezier;
[b1,b2] ← Cubic2.Split[bezier];
IF NOT subdivide[b1] THEN [] ← subdivide[b2];
};
};
IF
NOT proc[path.cubics[0].b0]
THEN FOR i:
NAT
IN [0..path.cubics.length)
DO
[] ← subdivide[path.cubics[i]];
ENDLOOP;
};
PointOnPath:
PUBLIC PROC [pt:
VEC, path: Path, tol:
REAL]
RETURNS [
BOOLEAN] = {
tolSq: REAL ← tol*tol;
found: BOOLEAN ← FALSE;
hit: PointProc = {
IF Vector2.Square[Vector2.Sub[p, pt]] <= tolSq THEN found ← TRUE;
RETURN [found];
};
IF pt.x < path.bounds.x
OR pt.y < path.bounds.y
OR pt.x>(path.bounds.x+path.bounds.w)
OR pt.y > (path.bounds.y+path.bounds.h) THEN RETURN [FALSE]
ELSE WalkPath[path, tol, hit];
RETURN [found];
};
TranslatePath:
PUBLIC PROC [path: Path, offset:
VEC] = {
FOR i:
NAT
IN [0..path.cubics.length)
DO
path.cubics[i].b0 ← Vector2.Add[offset, path.cubics[i].b0];
path.cubics[i].b1 ← Vector2.Add[offset, path.cubics[i].b1];
path.cubics[i].b2 ← Vector2.Add[offset, path.cubics[i].b2];
path.cubics[i].b3 ← Vector2.Add[offset, path.cubics[i].b3];
ENDLOOP;
path.bounds.x ← path.bounds.x+offset.x;
path.bounds.y ← path.bounds.y+offset.y;
};
TransformPath:
PUBLIC PROC [path: Path, tranformation: ImagerTransformation.Transformation, newBounds:
BOOL ←
TRUE] = {
FOR i:
NAT
IN [0..path.cubics.length)
DO
path.cubics[i].b0 ← ImagerTransformation.Transform[tranformation, path.cubics[i].b0];
path.cubics[i].b1 ← ImagerTransformation.Transform[tranformation, path.cubics[i].b1];
path.cubics[i].b2 ← ImagerTransformation.Transform[tranformation, path.cubics[i].b2];
path.cubics[i].b3 ← ImagerTransformation.Transform[tranformation, path.cubics[i].b3];
ENDLOOP;
IF newBounds THEN UpdateBounds[path];
};
<<UpdateBounds:
PUBLIC PROC [path: Path] = {
minX, minY: REAL ← 10E6;
maxX, maxY: REAL ← -10E6;
pointProc: PointProc = {
IF p.x > maxX THEN maxX ← p.x;
IF p.y > maxY THEN maxY ← p.y;
IF p.x < minX THEN minX ← p.x;
IF p.y < minY THEN minY ← p.y;
};
WalkPath[path, 1.0, pointProc];
path.bounds ← [x: minX, y: minY, w: maxX-minX, h: maxY-minY];
};
>>
Evaluate:
PROC [a
3, a
2, a
1, a
0, b
3, b
2, b
1, b
0:
REAL, u:
REAL]
RETURNS [pt:
VEC] = {
pt.x ← (((a3)*u + a2)*u + a1)*u + a0;
pt.y ← (((b3)*u + b2)*u + b1)*u + b0;
};
EvaluateX:
PROC [coeffs: Cubic2.Coeffs, t:
REAL]
RETURNS [x:
REAL] = {
x ← (((coeffs.c3.x)*t + coeffs.c2.x)*t + coeffs.c1.x)*t + coeffs.c0.x;
};
EvaluateY:
PROC [coeffs: Cubic2.Coeffs, t:
REAL]
RETURNS [y:
REAL] = {
y ← (((coeffs.c3.y)*t + coeffs.c2.y)*t + coeffs.c1.y)*t + coeffs.c0.y;
};
UpdateBounds:
PUBLIC
PROC [path: Path] = {
path.cubics[0].b0;
minX, minY: REAL ← 10E6;
maxX, maxY: REAL ← -10E6;
loEndX, hiEndX, loEndY, hiEndY: REAL;
coeffs: Cubic2.Coeffs;
tArray: ARRAY [1..2] OF REAL;
ptX, ptY: REAL;
rootCount: NAT;
FOR i:
NAT
IN [0..path.cubics.length)
DO
coeffs ← Cubic2.BezierToCoeffs[path.cubics[i]];
Find vertical points (find dx/dt and solve for t)
[tArray, rootCount] ← QuadraticFormula[3.0*coeffs.c3.x, 2.0*coeffs.c2.x, coeffs.c1.x];
FOR j:
NAT
IN [1..rootCount]
DO
IF tArray[j] >= 0.0
AND tArray[j] <= 1.0
THEN {
ptX ← EvaluateX[coeffs, tArray[j]];
minX ← MIN[minX, ptX];
maxX ← MAX[maxX, ptX];
};
ENDLOOP;
Find horizontal points (find dy/dt and solve for t)
[tArray, rootCount] ← QuadraticFormula[3.0*coeffs.c3.y, 2.0*coeffs.c2.y, coeffs.c1.y];
FOR j:
NAT
IN [1..rootCount]
DO
IF tArray[j] >= 0.0
AND tArray[j] <= 1.0
THEN {
ptY ← EvaluateY[coeffs, tArray[j]];
minY ← MIN[minY, ptY];
maxY ← MAX[maxY, ptY];
};
ENDLOOP;
Consider Endpoints
loEndX ← path.cubics[i].b0.x;
loEndY ← path.cubics[i].b0.y;
hiEndX ← path.cubics[i].b3.x;
hiEndY ← path.cubics[i].b3.y;
minX ← MIN[minX, loEndX, hiEndX];
minY ← MIN[minY, loEndY, hiEndY];
maxX ← MAX[maxX, loEndX, hiEndX];
maxY ← MAX[maxY, loEndY, hiEndY];
ENDLOOP;
path.bounds ← [x: minX, y: minY, w: maxX-minX, h: maxY-minY];
};
CubicPathsExtras all below plus made UpdateBounds Public
ClosestPoint:
PUBLIC PROC [pt:
VEC, path: Path]
RETURNS [
VEC] = {
closest: VEC ← path.cubics[0].b0;
minD: REAL ← 10E6;
test: PointProc = {
thisD: REAL ← Vector2.Square[Vector2.Sub[p, pt]];
IF thisD < minD THEN {closest ← p; minD ← thisD};
};
WalkPath[path, 1.0, test];
RETURN [closest];
};
ClosestPointAnalytic:
PUBLIC
PROC [pt:
VEC, path: Path, tolerance:
REAL ← 9999.0]
RETURNS [closest:
VEC, success:
BOOL] = {
Find the fifth degree equation in the parameter u and solve using an iterative root finder.
The spline equation is first represented as:
x(u) = a3*u3 + a2*u2 + a1*u + a0
y(u) = b3*u3 + b2*u2 + b1*u + b0
for 0 d u d 1.
bezier: Bezier;
c: Coeffs;
a3, a2, a1, a0, b3, b2, b1, b0: REAL;
p0, p1, p2, p3: VEC;
bestDist2, thisDist2: REAL;
thisPt: VEC;
thisSuccess: BOOL;
bestDist2 ← Real.LargestNumber;
success ← FALSE;
FOR i:
NAT
IN [0..path.cubics.length)
DO
bezier ← path.cubics[i];
p0 ← bezier.b0; p1 ← bezier.b1; p2 ← bezier.b2; p3 ← bezier.b3;
a3 ← -p0.x + 3*p1.x - 3*p2.x + p3.x;
a2 ← 3*p0.x - 6*p1.x + 3*p2.x;
a1 ← -3*p0.x + 3*p1.x;
a0 ← p0.x;
b3 ← -p0.y + 3*p1.y - 3*p2.y + p3.y;
b2 ← 3*p0.y - 6*p1.y + 3*p2.y;
b1 ← -3*p0.y + 3*p1.y;
b0 ← p0.y;
c ← Cubic2.BezierToCoeffs[bezier];
[thisPt, thisDist2, thisSuccess] ← ClosestPointCubic[c.c3.x, c.c2.x, c.c1.x, c.c0.x, c.c3.y, c.c2.y, c.c1.y, c.c0.y, pt];
IF NOT thisSuccess THEN ERROR;
IF thisDist2 < bestDist2
THEN {
bestDist2 ← thisDist2;
closest ← thisPt;
success ← TRUE;
};
ENDLOOP;
};
ClosestPointCubic:
PROC [a3, a2, a1, a0, b3, b2, b1, b0:
REAL, q:
VEC]
RETURNS [closest:
VEC, bestD2:
REAL, success:
BOOL] = {
Given a cubic parametric curve
x(u) = a3*u3 + a2*u2 + a1*u + a0
y(u) = b3*u3 + b2*u2 + b1*u + b0,
for 0 d u d 1.
we wish to find the closest point on this curve to q.
First, we find the extrema of (x - q.x)2 + (y - q.y)2 by setting
x/u(x - q.x) + y/u(y - q.y) = 0.
This gives us up to 5 real roots. We eliminate those for which u < 0 or u > 1. We evaluate the spline at the remain values of u and compute (x - q.x)2 + (y - q.y)2. We also compute the distance from q to the spline endpoints at u = 0 and u = 1. We choose the best point.
c5, c4, c3, c2, c1, c0: REAL;
poly: Polynomial.Ref;
roots: ShortRealRootRec;
thisD2: REAL;
thisPt: VEC;
c5 ← 3*(a3*a3 + b3*b3);
c4 ← 5*(a2*a3 + b2*b3);
c3 ← 4*(a1*a3 + b1*b3) + 2*(a2*a2 + b2*b2);
c2 ← 3*(a1*a2 + b1*b2) + 3*(a3*(a0 - q.x) + b3*(b0 - q.y));
c1 ← 2*(a2*(a0 - q.x) + b2*(b0 - q.y)) + (a1*a1 + b1*b1);
c0 ← (a1*(a0 - q.x) + b1*(b0 - q.y));
poly ← Polynomial.Quintic[[c0, c1, c2, c3, c4, c5]];
roots ← CheapRealRootsInInterval[poly, 0.0, 1.0];
bestD2 ← Real.LargestNumber;
success ← FALSE;
FOR i:
NAT
IN [0..roots.nRoots)
DO
IF roots.realRoot[i] > 0.0
AND roots.realRoot[i] <= 1.0
THEN {
thisPt ← Evaluate[a3, a2, a1, a0, b3, b2, b1, b0, roots.realRoot[i]];
thisD2 ← Vector2.Square[Vector2.Sub[thisPt, q]];
IF thisD2 < bestD2
THEN {
bestD2 ← thisD2;
closest ← thisPt;
success ← TRUE;
};
};
ENDLOOP;
FOR u:
REAL ← 0, u+1.0
UNTIL u > 1.0
DO
thisPt ← Evaluate[a3, a2, a1, a0, b3, b2, b1, b0, u];
thisD2 ← Vector2.Square[Vector2.Sub[thisPt, q]];
IF thisD2 < bestD2
THEN {
bestD2 ← thisD2;
closest ← thisPt;
success ← TRUE;
};
ENDLOOP;
};
CopyPath:
PUBLIC PROC [path: Path]
RETURNS [Path] = {
new: Path ← NEW[PathRec];
new.cubics ← NEW[PathSequence[path.cubics.length]];
FOR i:
NAT
IN [0..path.cubics.length)
DO
new.cubics[i] ← path.cubics[i];
ENDLOOP;
new.bounds ← path.bounds;
RETURN [new];
};
ReversePath:
PUBLIC PROC [path: Path] = {
reverses it in place
last: NAT ← path.cubics.length-1;
FOR i:
NAT
IN [0..path.cubics.length/2)
DO
temp: Cubic2.Bezier ← path.cubics[i];
path.cubics[i] ← path.cubics[last-i];
path.cubics[last-i] ← temp;
ENDLOOP;
FOR i:
NAT
IN [0..path.cubics.length)
DO
temp: VEC ← path.cubics[i].b0;
path.cubics[i].b0 ← path.cubics[i].b3;
path.cubics[i].b3 ← temp;
temp ← path.cubics[i].b1;
path.cubics[i].b1 ← path.cubics[i].b2;
path.cubics[i].b2 ← temp;
ENDLOOP;
};
QuadraticFormula:
PUBLIC
PROC [a, b, c:
REAL]
RETURNS [roots:
ARRAY [1..2]
OF
REAL, rootCount:
NAT] = {
The solution to the equation "ax2+bx+c=0". If a=0 this is just a linear equation. If
c = 0, one root is zero and we solve a linear equation. Otherwise, we use the
quadratic formula in either the form [-b+-(b2-4ac)1/2]/2a or (2c)/[-b-+(b2-4ac)1/2]
depending on the sign of b. We will arrange the roots so that roots[1] < roots[2].
discriminant, temp: REAL;
IF a = 0
THEN {
[roots[1], rootCount] ← LinearFormula [b, c];
RETURN;
};
IF c = 0
THEN {
roots[1] ← 0.0;
[roots[2], rootCount] ← LinearFormula[a, b];
rootCount ← rootCount + 1;
IF roots[1] > roots[2]
THEN {
-- swap them
temp ← roots[1];
roots[1] ← roots[2];
roots[2] ← temp;
};
RETURN};
discriminant ← b*b - 4*a*c; -- 3 mult, 1 add
SELECT discriminant FROM
=0 => {rootCount ← 1;
roots[1] ← -b/(2*a);
RETURN}; -- 1 mult, 1 div (4 mult, 1 div, 1 add total)
<0 => {rootCount ← 0;
RETURN}; -- (3 mult, 1 add total)
>0 => {
sqrtRadical: REAL ← RealFns.SqRt[discriminant];
term: REAL;
rootCount ← 2;
term ← IF b < 0 THEN sqrtRadical - b ELSE -sqrtRadical - b;
roots[1] ← term/(2.0*a);
roots[2] ← (2.0*c)/term;
IF roots[1] > roots[2]
THEN {
-- swap them
temp ← roots[1];
roots[1] ← roots[2];
roots[2] ← temp;
};
RETURN}; -- 1 SqRt, 2 mult, 2 div, 2 add (1 SqRt, 5 mult, 2 div, 3 add total)
ENDCASE => ERROR;
}; -- end of QuadraticFormula
LinearFormula:
PROC [a, b:
REAL]
RETURNS [root:
REAL, rootCount:
NAT] = {
The solution to the equation "ax + b = 0".
IF a = 0
THEN {
IF b = 0 THEN {root ← 0.0; rootCount ← 0; RETURN}
ELSE {root ← 0.0; rootCount ← 0; RETURN} -- inconsistent case
}
ELSE {
IF b = 0 THEN {root ← 0.0; rootCount ← 1; RETURN}
ELSE {root ← -b/a; rootCount ← 1; RETURN};
};
};
CheapRealRootsInInterval:
PUBLIC
PROC [poly: Polynomial.Ref, lo, hi:
REAL]
RETURNS [roots: ShortRealRootRec] = {
Uses a modified Newton-Raphelson iteration to find the real roots of a poly. It finds only those roots in the interval [lo..hi]. If the interval is relatively small (e.g., [0..1]) this routine is significantly faster than Polynomial.CheapRootsInInterval.
d: NAT ← Polynomial.Degree[poly];
roots.nRoots ← 0;
IF d<=1
THEN {
IF d=1 THEN {roots.nRoots ← 1; roots.realRoot[0] ← -poly[0]/poly[1]}
}
ELSE {
savedCoeff: ARRAY [0..5] OF REAL;
FOR i: NAT IN [0..d] DO savedCoeff[i] ← poly[i] ENDLOOP;
Polynomial.Differentiate[poly];
BEGIN
extrema: ShortRealRootRec ← CheapRealRootsInInterval[poly, lo, hi];
x: REAL;
FOR i: NAT IN [0..d] DO poly[i] ← savedCoeff[i] ENDLOOP;
IF extrema.nRoots>0
THEN {
x ← RootBetween[poly, lo, extrema.realRoot[0]];
IF x<= extrema.realRoot[0] THEN {roots.nRoots ← 1; roots.realRoot[0] ← x};
FOR i:
NAT
IN [0..extrema.nRoots-1)
DO
x ← RootBetween[poly, extrema.realRoot[i], extrema.realRoot[i+1]];
IF x <= extrema.realRoot[i+1] THEN {roots.realRoot[roots.nRoots] ← x; roots.nRoots ← roots.nRoots + 1};
ENDLOOP;
x ← RootBetween[poly, extrema.realRoot[extrema.nRoots-1], hi];
IF x <= hi THEN {roots.realRoot[roots.nRoots] ← x; roots.nRoots ← roots.nRoots + 1}
}
ELSE {
x ← RootBetween[poly, lo, hi];
IF x <= hi THEN {roots.realRoot[0] ← x; roots.nRoots ← 1}
};
END
};
};
RootBetween:
PROC [poly: Polynomial.Ref, x0, x1:
REAL]
RETURNS [x: REAL] =
BEGIN
debug xi: ARRAY [0..50) OF REAL;
debug i: NAT ← 0;
OPEN Polynomial;
y0: REAL ← Polynomial.Eval[poly,x0];
y1: REAL ← Polynomial.Eval[poly,x1];
xx: REAL;
xx0: REAL ← IF y0<y1 THEN x0 ELSE x1;
xx1: REAL ← IF y0<y1 THEN x1 ELSE x0;
yy0: REAL ← MIN[y0,y1];
yy1: REAL ← MAX[y0,y1];
NewtonStep:
PROC [x:
REAL]
RETURNS [newx:
REAL] =
BEGIN
y: REAL ← Polynomial.Eval[poly, x];
deriv: REAL ← Polynomial.EvalDeriv[poly, x];
IF y<0 THEN {IF y>yy0 THEN {xx0 ← x; yy0 ← y}};
IF y>0 THEN {IF y<yy1 THEN {xx1 ← x; yy1 ← y}};
newx ← IF deriv=0 THEN x+y ELSE x-y/deriv;
END;
IF y0=0 THEN RETURN[x0];
IF y1=0 THEN RETURN[x1];
IF (y0 > 0 AND y1 > 0) OR (y0 < 0 AND y1 < 0) THEN RETURN [x1+100];
xx ← x ← (x0+x1)/2;
first try Newton's method
WHILE x
IN [x0..x1]
DO
newx:REAL ← NewtonStep[x];
IF x=newx THEN RETURN;
x ← NewtonStep[newx];
xx ← NewtonStep[xx];
debug xi[i] ← xx; i ← i+1;
IF xx=x THEN EXIT;
ENDLOOP;
If it oscillated or went out of range, try bisection
debug SIGNAL Bisection;
IF xx0
IN [x0..x1]
AND xx1
IN [x0..x1]
THEN
BEGIN
x0 ← MIN[xx0,xx1];
x1 ← MAX[xx0,xx1];
END;
y0 ← Polynomial.Eval[poly,x0];
y1 ← Polynomial.Eval[poly,x1];
THROUGH [0..500)
DO
y: REAL ← Polynomial.Eval[poly, x ← (x0+x1)/2];
IF x=x0
OR x=x1
THEN
{IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]};
IF (y > 0 AND y0 < 0) OR (y < 0 AND y0 > 0) THEN {x1 ← x; y1 ← y}
ELSE {x0 ← x; y0 ← y};
IF (y0 > 0 AND y1 > 0) OR (y0 < 0 AND y1 < 0) OR (y0=0 OR y1=0) THEN {IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]}
ENDLOOP;
ERROR;
END;
GetParam:
PUBLIC
PROC [bezier: Cubic2.Bezier, pt:
VEC]
RETURNS [u:
REAL] = {
Assumes that pt is a point on bezier. Finds a parameter, u, corresponding to that point.
p0, p1, p2, p3: VEC;
a,b,c,d: REAL;
ignoreX, ignoreY: BOOL ← FALSE;
polyX, polyY: Polynomial.Ref;
epsilon: REAL = 0.01;
rootsX, rootsY: ShortRealRootRec;
ShortRealRootRec: TYPE = RECORD [nRoots: NAT, realRoot: ARRAY [0..5) OF REAL];
p0 ← bezier.b0;
p1 ← bezier.b1;
p2 ← bezier.b2;
p3 ← bezier.b3;
d ← p0.x - pt.x;
c ← 3*(p1.x - p0.x);
b ← 3*p0.x - 6*p1.x + 3*p2.x;
a ← -p0.x + 3*(p1.x - p2.x) + p3.x;
IF ABS[d] < epsilon AND ABS[c] < epsilon AND ABS[b] < epsilon AND ABS[a] < epsilon THEN ignoreX ← TRUE
ELSE {
polyX ← Polynomial.Cubic[[d,c,b,a]];
rootsX ← CheapRealRootsInInterval[polyX, 0.0, 1.0];
};
d ← p0.y - pt.y;
c ← 3*(p1.y - p0.y);
b ← 3*p0.y - 6*p1.y + 3*p2.y;
a ← -p0.y + 3*(p1.y - p2.y) + p3.y;
IF ABS[d] < epsilon AND ABS[c] < epsilon AND ABS[b] < epsilon AND ABS[a] < epsilon THEN ignoreY ← TRUE
ELSE {
polyY ← Polynomial.Cubic[[d,c,b,a]];
rootsY ← CheapRealRootsInInterval[polyY, 0.0, 1.0];
};
IF ignoreX
THEN {
IF ignoreY THEN ERROR;
RETURN [rootsY.realRoot[0]];
}
ELSE
IF ignoreY
THEN {
RETURN [rootsX.realRoot[0]];
}
ELSE {
Take advantage of the roots of both polynomials to improve accuracy.
FOR i:
NAT
IN [0..rootsX.nRoots)
DO
FOR j:
NAT
IN [0..rootsY.nRoots)
DO
IF
ABS[rootsX.realRoot[i] - rootsY.realRoot[j]] < epsilon
THEN
RETURN [(rootsX.realRoot[i] + rootsY.realRoot[j]) / 2.0];
ENDLOOP;
ENDLOOP;
ERROR;
};
};
CubicMeetsLine:
PUBLIC
PROC [bezier: Bezier, a, b, c:
REAL]
RETURNS [points:
ARRAY [0..2]
OF
VEC, hitCount: [0..3], tangency:
ARRAY [0..2]
OF
BOOL] = {
The spline equation is represented as:
x(u) = a3*t3 + a2*t2 + a1*t + a0
y(u) = b3*t3 + b2*t2 + b1*t + b0
for 0 d u d 1.
The line is represented as:
ax + by + c = 0.
Thus, the points of intersection are at the t values that satisfy:
a*(a3*t3 + a2*t2 + a1*t + a0) + b*(b3*t3 + b2*t2 + b1*t + b0) + c = 0.
or equivalently:
(a*a3 + b*b3)*t3 + (a*a2 + b*b2)*t2 + (a*a1 + b*b1)*t + (a*a0 + b*b0 + c) = 0.
f: Coeffs;
poly: Polynomial.Ref;
roots: ShortRealRootRec;
lastT: REAL;
success: BOOL;
epsilon: REAL = 0.0001;
hitCount ← 0;
tangency[0] ← tangency[1] ← tangency[2] ← FALSE;
f ← Cubic2.BezierToCoeffs[bezier];
poly ← Polynomial.Cubic[[a*f.c0.x + b*f.c0.y + c, a*f.c1.x + b*f.c1.y, a*f.c2.x + b*f.c2.y, a*f.c3.x + b*f.c3.y]];
roots ← CheapRealRootsInInterval[poly, 0.0, 1.0];
lastT ← 2.0; -- a value outside the range [0..1];
success ← FALSE;
FOR i:
NAT
IN [0..roots.nRoots)
DO
IF roots.realRoot[i] > 0.0
AND roots.realRoot[i] <= 1.0
THEN {
IF
ABS[roots.realRoot[i]-lastT] > epsilon
THEN {
lastT ← roots.realRoot[i];
points[hitCount] ← Evaluate[f.c3.x, f.c2.x, f.c1.x, f.c0.x, f.c3.y, f.c2.y, f.c1.y, f.c0.y, lastT];
hitCount ← hitCount + 1;
}
ELSE {
IF lastT # 2.0 THEN tangency[hitCount-1] ← TRUE;
};
};
ENDLOOP;
};
END.