Today, I want to write about a pretty simple problem: finding the line where two planes intersect. At the end of the day, it boils down to a simple formula that's easy to find elsewhere online. However, I find the derivation interesting, and it serves as a nice introduction to some powerful techniques that can be used on harder problems like intersecting conic sections.

But first, in case you just want the formula: the planes \(\langle n_1, x\rangle + d_1 = 0\) and \(\langle n_2, x \rangle + d_2 = 0\) intersect along the line \(r_o + t r_d\), where \[\begin{aligned} r_d &:= n_1 \times n_2,\\ r_o &:= \frac {r_d \times(d_2n_1 - d_1n_2)} {\|r_d\|^2}. \end{aligned}\]

Prelude: intersecting two lines in the plane

When solving these sorts of problems, it is helpful to work in homogeneous coordinates. So we represent a point \((x, y) \in \mathbb{R}^2\) by the vector \((x, y, 1) \in \mathbb{R}^3\) (or any scalar multiple of this vector). And similarly, we represent the line \(ax + by + c = 0\) using the vector \((a, b, c) \in \mathbb{R}^3\) (or any scalar multiple of this vector). It's no coincidence that both points and lines share the same representation here, but that's a discussion for another day.

Suppose now that we want to find the intersection between two lines \(\ell_1\) and \(\ell_2\) (which are represented as vectors in \(\mathbb{R}^3\). A point \(p\) lies on lines \(\ell_1\) if \(\langle p, \ell_1\rangle = 0\), and similarly for \(\ell_2\), so their intersection must be a vector \(p \in \mathbb{R}^3\) which is simultaneously orthogonal to both \(\ell_1\) and \(\ell_2\). We can construct such a vector easily by taking the cross product \(\ell_1 \times \ell_2\).

Visually, each of our lines in \(\mathbb{R}^2\) becomes a plane passing through the origin when represented in homogeneous coordinates in \(\mathbb{R}^3\) and the point in which the lines intersect in \(\mathbb{R}^2\) becomes the line passing through the origin in \(\mathbb{R}^3\). This line can easily be found by taking the cross product of the planes' normal vectors, and yields homogeneous coordinates for the intersection point in \(\mathbb{R}^2\). Computing the intersection of these planes is particularly simple since both planes pass through the origin; computing the intersection between general planes will take us a bit more work later on.

Before we move on, though, I'll mention that this exact same construction also works to compute the line that connects two points. If we have points \(p_1, p_2\), then the line passing through these two points must be given by a vector \(\ell \in \mathbb{R}^3\) which is simultaneously orthogonal to \(p_2\) and \(p_2\), i.e. \(\ell = p_1 \times p_2\).

Lines in 3D

Now let's move to 3D. Since we've moved up a dimensions, points and planes are represented by vectors in \(\mathbb{R}^4\). But what about lines? One natural representation of lines is given by Plücker coordinates. If we consider the line traced out by \(r(t) = r_o + t r_d \in \mathbb{R^3}\), its Plücker coordinates are given by the vector \[(r_o \times r_d, r_d) \in \mathbb{R}^6.\] This formula looks a little strange at first, but it has many nice properties: for instance, if we had picked some other point on the line as the base point, say \(r_o + \lambda r_d\), then we would still get the same Plücker coordinates since \((r_o + \lambda r_d) \times r_d = r_o \times r_d\). Similarly, if we scale the line direction \(r_d\) by a constant \(\lambda\), then our Plücker coordinates also get scaled by \(\lambda\), but still represent the same line. And moreover, we can recover an equation from the line by its Plücker coordinates. We can recover the direction from the last 3 coordinates, and we can find a point on the line by applying the identity that \(r_d \times (r_o \times r_d) = \|r_d\|^2 r_o - \langle r_o, r_d\rangle r_d\), yielding the point on our line which is closest to the origin.

What if we're given two points \(p_1, p_2\) and we want the line \(\ell\) connecting them? To start with let's represent out points using ordinary coordinate vectors in \(x_1, x_2 \in \mathbb{R}^3\). We can think of our line as starting at \(x_1\), and proceeding in direction \(x_2 - x_1\), so we obtain Plücker coordinates \[\ell = ((x_2 - x_1) \times x_1, x_2-x_1) = (x_2 \times x_1, x_2 - x_1).\] But if you stare at this formula for a minute, and are familiar with wedge products, you may notice that this formula can be simplified to \(\ell = (x_2, 1) \wedge (x_1, 1).\) Indeed, if we write each \(p_i\) using homogeneous coordinates in \(\mathbb{R}^4\), the equation becomes \[\ell = p_2 \wedge p_1.\]

Planes in 3D

What about the plane \(P\) passing through three points \(p_1, p_2, p_3\)? We can construct its homogeneous coordinates as \(P = *(p_1 \wedge p_2 \wedge p_3)\), noting that \[\langle P, p_i \rangle = p_1 \wedge p_2 \wedge p_3 \wedge p_i = 0,\] for each of our points \(p_i\). Similarly, the plane passing through a point \(p\) and a line \(\ell\) is given by \(P = *(\ell \wedge p)\).

This formula will help us determine when a line \(\ell\) is contained in a plane \(P\): the line \(\ell\) lies in \(P\) if there is some point \(p\) such that \(P = *(\ell \wedge p)\). Using another wedge product identity, we can write this equation as \(P = -\iota_{p^\flat} *\ell\), which in turn is true if and only if \(P \wedge *\ell = 0\).

Intersecting planes in 3D

Now we can finally find the formula to intersect two planes in 3D. If we are given two planes \(P_1\) and \(P_2\), the line \(\ell\) contained in their intersection must satisfy the equations \[\begin{aligned}P_1 \wedge *\ell &= 0,\\P_2 \wedge *\ell &= 0.\end{aligned}\] It turns out that this is precisely the system of equations that we solved above to find the intersection of lines in 2D, just written in exterior algebra! And, it has the same solution (written in exterior algebra): \(\ell = *(P_1 \wedge P_2)\).

All that remains is to unpack all of the exterior algebra a formula in terms of more familiar vector algebra operations. Let's denote our planes as \(P_i = (n_i, d_i)\). Then, as we saw above, their wedge product is given by \[P_1 \wedge P_2 = (n_1 \times n_2, d_2 n_1 - d_1 n_2).\] Then the Hodge star simply swaps the components of this vectors, yielding \[\ell = *(P_1 \wedge P_2) = (d_2 n_1 - d_1 n_2, n_1 \times n_2).\] Finally, we can use the formula from earlier to convert these Plücker coordinates into a point on our line and a direction. Putting everything together, the planes \(\langle n_1, x\rangle + d_1 = 0\) and \(\langle n_2, x \rangle + d_2 = 0\) intersect along the line \(r_o + t r_d\), where \[\begin{aligned} r_d &:= n_1 \times n_2,\\ r_o &:= \frac {r_d \times(d_2n_1 - d_1n_2)} {\|r_d\|^2}. \end{aligned}\]

A more symmetric form

I'll quickly note at the end that everything I covered here is arguably cleaner if you use exterior algebra from the beginning. Rather than representing both points and hyperplanes using vectors in \(\mathbb{R}^{n+1}\), we can represent points as vectors in \(p \in \mathbb{R}^{n+1}\) and hyperplanes as vectors in the dual space \(P \in \Lambda^n \mathbb{R}^{n+1}\). And, rather than looking at the inner product between \(p\) and \(P\) (which is no longer well-defined), we can consider the primal-dual pairing \(\langle P, p\rangle = *(p \wedge P).\) One nice aspect of version of the theory is that then our formula for the line \(\ell\) between two points \(p_1, p_2\) is always given by \(\ell = p_1 \wedge p_2\), regardless of whether we are considering \(\mathbb{R}^2\), \(\mathbb{R}^2\), or \(\mathbb{R}^n\).

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)\)