The bezier drawing routines come in two parts.
The basic routines draw bezier curves following Knuth's description in
Metafont. The main routine is
draw_bezier(z0, z1, z2, z3)
where the zi are the control points, expressed in scaled integers,
with a fixed scale of 2^16.
Thus the gfxwin command is
bezier x0 y0 x1 y1 x2 y2 x3 y3
where the xi and yi are 32 bit integers, signed.
The first step is to subdivide the curve into octants. The main point
of this is that the group generated by x-y swaps and rotation by 90
degrees leaves pixels fixed, so that if S lies in this group, then the
Bezier curve through the S zi should be S applied to the Bezier curve
through the zi.
The subdivisions will have the property that x^', y^', x' - y' will all
be greater than or equal to 0, and will be flagged by an integer
telling which octant they belong to.
Thus, one curve becomes two in the upper plane (x' >= 0); then each of
these divides into two in the first quadrant (y' >= 0); then each of
these into octants (0 <= y' <= x').
A subroutine: find first t where B(t) goes from positive to negative
in [0, 1] if B is a quadratic Bezier function. (B = derivative of cubic)
Another subroutine: subdivision at t.
For each part, move the curve into the positive octant, and then
flag in three bits where it came from: mode bits NegY, NegX, SwapXY.
So in the end we are facing a chain of positive octant curves plus flags.
The flags affect the subroutine draw_pixel(x, y, mode) --- how pixels are
actually plotted.
For finding the crossing point of B(z0, z1, z2, t): first reduce to
case where z0 >= 0, min(z1, z2) < 0. Then continually bisect, choosing each
time the first subinterval where this remains true, or returning -1 if none.
This is \S 392 in Metafont.
In fact, following Knuth, I describe the structure of bezier control
as a set of
(*) vectors x[0..4], y[0..4]
(*) level ell 0..15
(*) octant 0..7
The point x[0], y[0] is just the origin of the curve, while the adjusted
z[1..3] = (x[1..3], y[1..3])
are the scaled differences 2^{ell}*(z[i] - z[i-1]), and z[4] is the
last point, the sum of z[0] + z[1] + z[2] + z[3] (invariably).
The octant is a flag made up of NegX, NegY, SwapXY.
====================================================================
The crossing routine (sections 391, 392 of Knuth)
We want to return the first value of t in [0, 1] where B(x0, x1, x2, t)
crosses from non-negative to negative, if it exists, or 1 if it never
does. The basic routine is this:
If x0 < 0, return(0).
If x0 = 0: if x1 < 0, return(0). Otherwise x1 >= 0. If x2 >= 0,
return(1).
If x0 > 0, x1 >= 0, x2 >= 0, return(1).
Otherwise x0 > 0 and min(x1, x2) < 0, or x0 = 0, x1 > 0, and x2 < 0.
Let t be the return value. Then t must be in (0, 1].
(It is still possible that there is no crossing, say if x1 < 0, x2 > 0,
and the minmium value lies above the x-axis.)
We are now reduced to case x0 >= 0, min(x1, x2) < 0.
Let t be the return value.
We keep bisecting the interval, choosing at each point the half where
this continues to hold, or exiting with t = 1 if neither does.
At start t = 0, ell = 0. Thus at the start of each loop, t represents
approximation at level ell to crossing, so if there is a crossing it
occurs in (t, t + 1/2^{ell}).
Loop: let x01 = half(x0 + x1), x02 = half(x1 + x2), x11 = half(x01 + x02).
then x0, x01, x11, x02, x2 represents the de Casteljau bisection.
If (x01 < 0 or x11 < 0) we replace x0 x1 x2 by x0 x01 x11, leave t
fixed, increment ell.
Otherwise x01 >= 0 and x11 >= 0; there is no crossing in the
first half. If x02 < 0 or x2 < 0 replace x0 x1 x2 by x11 x02 x2,
t by t + 1/2^{ell}, increment level ell. Else (there is no crossing)
return 1.
If ell < 16 loop, else exit loop. At end, return(t).
Knuth speeds this up a bit. He lets ell be implicit. He replaces
t by d = (1 + t)*2^{ell}, loops while d < 2^{16}, returns (d/2^{ell} - 1).
If we choose the first half interval, set d = 2*d, else d = 2*d + 1.
He works with y1 = x0 - x1, y2 = x1 - x2, and x0. The conditions
x0 >= 0, min(x1, x2) < 0 <=> max(y1, y1 + y2) > x0 >= 0.