4

I am looking for a detailed description of an algorithm for the classical problem of computing the intersection of two conic curves. The curves are given by two equations of the form:

$$ a x^2 + b y^2 + c xy + d x+e y + f = 0 $$

I need the algorithm to extend my graphics library MetaPict. If possible I prefer references with pseudo-code (or real code), but other references are also welcome.

mvw
  • 34,562
soegaard
  • 181

4 Answers4

8

You can work by considering the pencil containing the two conics ($\lambda f(x,y)+(1-\lambda) g(x,y)=0$) and find its degenerate element, which is formed of two straight lines. Then intersecting the straight lines with one of the conics is straightforward.

Another approach is by eliminating $y$ to obtain a fourth degree polynomial equation.

Combine the two equations to eliminate $y^2$. Then you can express $y$ in terms of a rational function of $x$. $$ a x^2 + b y^2 + c xy + d x+e y + f = 0\\ a' x^2 + b' y^2 + c' xy + d' x+e' y + f' = 0\\ $$ $$ (ab'-a'b) x^2 + (cb'-c'b) xy + (db'-d'b) x+(eb'-e'b) y + (fb'-f'b) = 0,\\ y=-\frac{(ab'-a'b) x^2 +(db'-d'b) x+ (fb'-f'b)}{(cb'-c'b) x +(eb'-e'b)}=\frac{a''x^2+d''x+f''}{c''x+e''}. $$ Substituting in one of the original equations, you get a univariate polynomial.

$$(a x^2+ d x+ f)(c''x+e'')^2+ b (a''x^2+d''x+f'')^2 - (c x +e) (c''x+e'')(a''x^2+d''x+f'')= 0 .$$

This can be solved by closed formulas.

  • 1
    The book by Richter-Gebert use the pencil method. Finding $lambda$ requires one to solve a third degree equation (I wrote a solver following Kahan's description of a Newton-Raphson algorithm tailored to the cubic equation). Kahan shows that the classical closed formula can be numerically unstable, thus a iterative method is used.

    Then the degenerate form must be split. This may or may not give two real lines (a double line and a point is also possible). In the book the work is done with homogenous coordinates (possible complex) so there are no explicit special cases in the calculation.

     – soegaard May 08 '16 at 18:25
  • 1
    @soegaard: if I am right, the two approaches are essentially one, though this is not immediately obvious. The resolution of the quartic works by factorization in two quadratics, and this factorization requires the solution of a cubic. You can parallel some steps of the pencil approach. Resolution of a quadratic is indeed not seamless. –  May 08 '16 at 19:51
  • 1
    Did Kahan say why the closed formula was unstable? Specifically, what makes him believe the instability is inherent to the formula and not cough his implementation? The solution to the cubic equation is just arithmetic, and a couple of simple roots. For conics, the relations among the coefficients are not arbitrary, and we know what they are. I fail to see how I can be prevented from using this information to construct a stable implementation. But I could be mistaken. – Slumberland May 18 '16 at 01:14
  • @Slumberland Bit late, but I didn't see your comment before. The Kahan paper is here: https://people.eecs.berkeley.edu/~wkahan/Math128/Cubic.pdf – soegaard Dec 02 '18 at 19:39
4

The book "Perspectives on Projective Geometry" by Jürgen Richter-Gebert contains a detailed description of both theory and algorithms concerning calculations with conics. The book is light of examples though.

In the algorithm for splitting a degenerate conic into two lines, there is missing a minus sign. Page 190 step 3 should be: beta = sqrt( - B_{i,i} ).

For an example see the answers Intersection of conics using matrix representation and Decomposition of a degenerate conic which use a slightly different algorithm, but is nevertheless very helpful.

Community
  • 1
soegaard
  • 181
2

I have a complete solution which works fine for me but should undergo a bit more testing. I made it because Geogebra is terrible at finding roots.

My approach will benefit only modestly from parallelization. I am not convinced matrices are an improvement, but for a graphics library it is surely worth checking. My experience is this: symbolically, they make the overall problem easier to write and organize, but numerically they make the problem sequence longer, and more prone to error. Specifically, the steps are sequential, there is no way to avoid solving a cubic equation, and the matrices are not direct transformations, but represent ordinary linear systems which have to be solved. Moreover, matrix inverse cannot in general be assumed.

Instead, I simply solve the fourth order equation directly. I focus on a solid, efficient calculation method for all roots (real and imaginary) of a third order polynomial. From there, the fourth order solution is basically free. This delimits the problem, allowing me to focus on identifying and mitigating the sources of numerical instability.

For example, consider the case that the two sections are touching (tangent). The author of Intersection of conics using matrix representation flags this as a potential source of instability. But represented as an ordinary polynomial, this problem is well defined: we have simply a double root: f(x) = 0 and f'(x) = 0.

I believe my method is stable, and I have spent time insuring that difficult cases are handled correctly, but I have not been exhaustive. We should test it more carefully for implementation in a library. If you are interested, let me know. I am happy to post it here.

Community
  • 1
Slumberland
  • 181
1

WLOG, one of the conics is $$x^2\pm y^2=1$$

because you can center/reduce one of the equations and let the other conic undergo the same transformation. (There is also the limit parabolic case, which I am not considering.)

Then you can solve the quadratic (hyperbolic-)trigonometric equation

$$a\cos^2\theta+2b\cos\theta\sin\theta+c\sin^2\theta+2d\cos\theta+2e\sin\theta+f=0$$ (or with $\cosh,\sinh$).

This can be rationalized by the Weierstrass transformation, giving

$$a\left(\frac{1-t^2}{1+t^2}\right)^2+2b\frac{1-t^2}{1+t^2}\frac{2t}{1+t^2}+c\left(\frac{2t}{1+t^2}\right)^2+2d\frac{1-t^2}{1+t^2}+2e\frac{2t}{1+t^2}+f=0$$

which is quartic. (Similar formula in the hyperbolic case.)