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.