Bézier approximation of conics I want to draw conics in some software. Sounds simple enough, but the combination of special cases and noisy numerics cause me tons of grief, so I turn to you looking for help.
The setup is as follows: I've got a conic section, described by a quadratic form $Q$ as the set of points $\{(x,y)\mid Q(x,y)=0\}$. Internally I'm using homogeneous coordinates and a symmetric matrix for the conic, but a quadratic equation with $6$ coefficients works just as well. The portion of this conic that is visible to the user is bounded by a rectangle, e.g. $\{(x,y)\mid 0\le x\le w,0\le y\le h\}$. Now I want to find a set of (probably quadratic) Bézier curves which approximate this conic reasonably well. I'm not talking about rational Bézier curves here. I want my algorithm to be able to deal with various special cases (see below) and imprecise floating-point computations.
Special cases
There are a number of possible cases:


*

*Ellipse

*Hyperbola

*Parabola (center at infinity!)

*Circle (no well-defined principal axes!)

*Single point (e.g. circle of radius zero)

*Pair of lines (matrix has rank two)

*Pair of lines, one of them the line at infinity

*Double line (matrix has rank one)


It is easy to characterize all of these in terms of the coefficients of the quadratic form. But numerics make this tricky, as the conic would usually not fall exactly into one of these classes. What's worse, due to rounding errors it may behave differently in different situations. So for example a pair of lines might look like a hyperbola at first, but when trying to find it's vertical tangents or some such, then the two coinciding solutions might become slightly imaginary instead of real and distinct.
One may be tempted to allow for numeric errors by considering values smaller than some $\varepsilon$ as exactly zero. But getting those $\varepsilon$ right is really hard: too large a value will let a hyperbola look like a pair of lines when it should not, while too small a value will not be enough to exclude the problems associated with these special cases.
It would be best if all the special cases could be included with one of the not so special cases. A pair of lines should be like a hyperbola. A circle or a single point should be like an ellipse. A parabola could be considered either a hyperbola or an ellipse, it should not matter. If the algorithm were robust enough to correctly draw all the relevant almost-special cases, then it should be doing the right thing for the really special cases as well. The double line is really hard, so even a solution which doesn't handle that correctly would be a huge step forward.
It would be best to have some form of parsimonious algorithm: one which computes information from the quadratic form sparingly in order to avoid inconsistencies. Anything that can be derived from known information should be concluded from that to stay consistent.
Outline of one approach
One approach I keep coming back to looks as follows, at a high level:


*

*Compute some special points: intersections with the boundary, points where the tangents are horizontal or vertical, things like this.

*Connect these special points so that the combinatorics match those of the conic arcs to be drawn.

*Draw each conic arc, starting at the points of intersection. Recursively add more points until the approximation is good enough.


One way to do the first step in a consistent way would be by evaluating the sign of the quadratic form for each of the rectangle corners. If the endpoints of a rectangle edge have different signs, there has to be one point of intersection in between them, so pick the more likely one. If the endpoints have equal sign, there are zero or two points of intersection, and both of these choices should be consistent with the rest so simply rely on the numeric result.
But how can you tell which special points to connect? Each point on the cyclic boundary has to connect to either the point before or the one after, so there are only two possible solutions, but so far I haven't come up with an easy way to distinguish between these. Adding further special points like those with special tangents makes this harder.
Once the points to connect have been identified, there might be a question of how to connect them, i.e. which direction to go. But if I managed to include all the horizontal and vertical tangencies, then the arc should cover at most a quarter turn, and the direction should be obvious from this.
During recursive refinement, one will often want to intersect some line with the conic. How can we tell which of the two points of intersection actually belongs to the arc that's currently being processed?
One also has to determine control points between these points on the conic. After some experiments with cubic splines matching curvature, I decided to go for quadratic. Finding tangent directions from the gradient is easy, and intersecting these gives the control point. Except if the tangents are almost identical, so the point of intersection is poorly defined…
Literature
I doubt I'm the first one to consider these questions.
But so far my research into literature hasn't come up with anything promising. Some people do pixel-based scan conversion. Some use rational Bézier splines. Some parametrize the curve like a circle or parabola, which works very poorly for degenerate cases.  So if you know of a useful reference for this, please share it.
 A: Here is how I would approach this.
First perform the centering/reduction of the conic to be able to use simple parametric equations. [The degenerate cases can be processed specifically using points and line segments.] Working with the parametric equation reduces parts of the problem to 1D, which makes it much more tractable.
At the same time, the viewing window becomes a rotated rectangle. You can see it as the intersection of four clipping half planes.
Solve the clipping problem for a single half plane. The intersections with the line will tell you a range of the parameter which gets excluded/kept. (One or two intervals.)
Then you can combine the four clipping planes and consider the intersections of the allowed ranges of the parameter. This will give you up to four intervals, hence up to four arcs.
Now you have to approximate an arc of a conic with one of more Bezier, say quadratic. By constructing the tangents at the endpoints and intersecting them, you get three control points that define the Bezier arc.
You must find a way to assess the closeness of the Bezier to the true arc. An expedite solution can be to sample a few points along the Bezier and estimate their distance to the conic. The exact distance is difficult to compute but you have enough with an approximation (the so-called algebraic distance).
When the total distance exceeds a threshold, subdivide the conic arc and redo the approximation. (Do this recursively.)
A: This is too long for a comment, and it may not be very helpful, but I see a few fundamental problems with what you're doing.
The first is the use of algebraic coefficients to represent the curves. The relationships between these coefficients and geometric parameters are very complex. That being the case, it's very difficult to find $\epsilon$ values that have some geometric meaning. And you always need these epsilon values if you're using floating-point arithmetic.
The second problem is that you seem to be trying to distinguish between cases that are fundamentally indistinguishable in floating point arithmetic. There is really no difference between a circle and an almost-circular ellipse, or between a very pointy hyperbola and a pair of lines. I'm not sure what you're going to do with the curves, but in most applications the calculation will give you the same answer regardless of whether you use the hyperbola or the pair of lines, so distinguishing between them is pointless (and futile). If you're willing to use Bezier curve approximations, this underscores the idea that the exact forms of the curves don't matter very much to you.
Lastly, it sounds like one of the main things you want to do with the curves is draw them. Implicit equations like $Q(x,y) = 0$ are not well suited to vector-based drawing. Parametric equations are much better. There are ways to obtain parametric equations for conic section curves represented in implicit form, but they have lots of $\epsilon$ problems. You'd be much better off if you used parametric equations from the outset. The disadvantage is that you will need special cases -- there is no unified way to represent all conic section curves with a single type of parametric equation. Rational quadratics come close, but still won't give you both branches of a hyperbola, or a pair of lines. But identifying the special cases from the beginning is much better than trying to recover them from an implicit equation later.
If you just want to draw the curves on a graphics display, and you're forced to use implicit equations, then you can generate pixels directly, without first generating vectors. 
A good reference on this topic is a paper written by E. T. Y Lee in 1987, entitled "The rational Bézier representation for conics". There have been a few others since then, but none of them are significantly better, AFAICS. Once you have rational quadratics, approximation by polynomial quadratics or polynomial cubics is straightforward using something like the test-and-subdivide approach suggested by Yves.
