Intersections of Conics

Suppose we have two conic sections and we want to find their intersections. If the conics are circles, then this is easy, and can be done in closed form. However, if we have, say, a pair of hyperbolas then things are harder. For instance, a pair of hyperbolas can intersect in four points, unlike a pair of circles which generally intersect in two points.

However, there is still a pretty slick algorithm for computing the intersection points between any two conic sections, which I'll discuss in this post. A detailed description can be found in Perspectives on Projective Geometry by Jürgen Richter-Gebert.

Matrix representation of conic sections

It turns out to be much easier to do these kind of calculations on conics if we represent them as symmetric matrices via homogeneous coordinates. Explicitly, a conic can be written as the set of points \(x, y\) satisfying an equation of the form \[Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0.\] This equation can be rewritten in matrix notation as \[\begin{pmatrix}x & y & 1\end{pmatrix} \begin{pmatrix}A & B/2 & D/2 \\ B/2 & C & E/2 \\ D/2 & E/2 & F\end{pmatrix} \begin{pmatrix}x\\y\\1\end{pmatrix}=0,\] and we can identify our conic with the matrix in the middle of this equation. Writing conics as matrices allows us to work with them algebraically: for instance, we can add together two conics by adding together their matrices. This kind of algebraic manifpulation of conics is the key to the algorithm for computing intersections that I describe below.

Computing the intersections

Suppose we have two conics given by symmetric matrices \(Q_1, Q_2 \in \mathbb{R}^{3 \times 3}\) with intersection points \(p_1, \ldots, p_4\). Explicitly, this means that for each \(p_i\), we have \[p_i^T Q_1 p_i = p_i^T Q_2 p_i = 0.\]

But by linearity, this means that for any coefficients \(\lambda, \mu \in \mathbb{R}\), we also have \[p_i^T (\lambda Q_1 + \mu Q_2) p_i = 0,\] defining a whole family of conic sections that also pass through these four points. If we plot some of the other conics generated from the two hyperbolas shown above, we see that we find some hyperbolas, some ellipses, and even a pair of diagonal lines all passing through our four intersection points.

Those two diagonal lines turn out to be the key to locating the intersection points \(p_i\). After all, it's pretty easy to solve for the intersection between a conic and a line – you just have to solve a quadratic equation.

Formally, the two diagonal lines make up a degenerate hyperbola. A matrix represents a degenerate conic if its determinant is zero, so we can find this degenerate hyperbola by solving \(\det(\lambda Q_1 + \mu Q_2) = 0\), which is a cubic equation (since \(Q_1\) and \(Q_2\) are 3 \(\times\) 3 matrices). Then we simply need to identify the two lines which compose this degenerate hyperbola and intersect them with one of the original hyperbolas. In summary, the algorithm is as follows:

1. Find \(\lambda, \mu \neq 0\) such that \(\det(\lambda Q_1 + \mu Q_2) = 0\).
2. Decompose the degenerate hyperbola \(\lambda Q_1 + \mu Q_2\) into a pair of lines \(g, h\).
3. Intersect \(g\) and \(h\) with \(Q_1\) to identify all four intersection points.

Below, I'll discuss each of the steps in some more detail.

Finding the degenerate hyperbola

Richter-Gebert gives a detailed algorithm for finding the roots of a cubic equation in homogeneous coordinates, but if you have access to a linear algebra library, you can find the roots more easily. First, if \(Q_2\) is itself degenerate (i.e. \(\det(Q_2) = 0\)), then taking \(\lambda = 0, \mu = 1\) gives us the solution that we desire. On the other hand, if \(Q_2\) is nondegenerate (i.e. invertible), then we can set \(\lambda = 1\) and multiply through by \(Q_2^{-1}\) and instead solve \[\det\left(\mu I - \left(-Q_1 Q_2^{-1}\right)\right)=0.\] Now, this problem may not obviously look easier, but there is a big advantage: the solutions \(\mu\) are precisely the eigenvalues of the matrix \(-Q_1Q_2^{-1}\), which can easily be computed by standard linear algebra libraries!

Decomposing the degenerate hyperbola into lines

For now, suppose we have some degenerate hyperbola represented by a symmetric matrix \(A \in \mathbb{R}^{3\times3}\). This degenerate hyperbola is really just a pair of lines, which can be represented by two vectors \(g, h \in \mathbb{R}^3\). Explicitly, a point \(x\) lies on the first line if \(g^Tx = 0\), and similarly for \(h\). Hence, \(x\) lies on the union of the two lines whenever \(x^T gh^T x = 0\). So if our symmetric matrix \(A\) represents this pair of lines, then (up to scale), \(A\) must equal the symmetrized matrix \(gh^T + hg^T\). Richter-Gebert gives a neat algorithm for computing \(g\) and \(h\) from \(A\). First, he notes that the cofactor matrix of \(A\) is given by \(A^\triangle = -(g \times h)(g\times h)^T\). And furthermore, a direct calculation shows that \(gh^T - hg^T\) is simply the cross product matrix \([g\times h]_\times\), and hence \(2gh^T = A + [g\times h]_\times\).

Concretely, then, \(g\) and \(h\) may be obtained form \(A\) as follows:

1. \(B \gets A^\triangle\)
2. \(i \gets\) the index of a nonzero diagonal entry of \(B\)
3. \(\beta \gets \sqrt{-B(i,i)}\)
4. \(p \gets B(:, i) / \beta\)
5. \(C \gets A + [p]_\times\)
6. \(i, j \gets\) the index of a nonzero entry of \(C\)
7. \(g = C(i, :), h = C(:, j)\)

Intersecting a line with a conic

Suppose we wish to intersect the conic \(A \in \mathbb{R}^{3\times 3}\) with the line \(g \in \mathbb{R}^3\). Perhaps unsurprisingly, Richter-Gebert also has a neat algorithm for intersecting a conic with a line in homogeneous coordinates. I find the motivation for this procedure to be the most mysterious, although the steps themselves are fairly simple. First, we assume that \(g(2) \neq 0\), reindexing the entries if necessary. Then, we do the following:

1. \(B \gets [g]_\times^T A [g]_\times\)
2. \(\alpha \gets \frac 1 {g(2)} \sqrt{B(0, 1)^2 - B(0, 0) B(1, 1)}\)
3. \(C \gets B + \alpha [g]_\times\)
4. \(i, j \gets\) the index of a nonzero entry of \(C\)
5. \(g = C(i, :), h = C(:, j)\)