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

This post lists some useful formulas for evaluating the distortion of maps between triangle meshes. The distortions considered here are intrinsic, meaning they only depend on the change in triangle edge lengths. Conveniently, they also happen to have fairly simple formulas in terms of the triangle edge lengths.

Background

There are many ways to measure the distortion of a map \(\phi : M \to N\) between surfaces. The general strategy is to compute the Jacobian \(J_\phi\), which is a \(2\times 2\) matrix, and consider different functions of its singular values \(\sigma_1, \sigma_2\). For instance, the product \(\sigma_1\sigma_2\) measures the area distortion of \(\phi\) (as it is simply the determinant \(\det J_\phi\)), while the ratio \(\sigma_1 / \sigma_2\) measures the amount of anisotropic stretching induced by \(\phi\). (See e.g. section 2 of Khodakovsky et al. [ 2003; free version] for more details.)

A piecewise-linear map between triangle meshes, deforms each triangle by a linear map, so its distortion is constant on each triangle. Below, I give some formulas for computing such per-triangle distortions intrinsically, that is, computing the distortions using only the lengths of the initial and deformed triangles. In each formula, I refer to the triangle's vertices as \(i, j, k\). The initial edge lenths are denoted \(\ell_{ij}, \ell_{jk}, \ell_{ki}\), and the corner angles are denoted \(\alpha_i, \alpha_j, \alpha_k\). Quantities measured after deformation are denoted \(\tilde \ell_{ij}\), etc.

Area Distortion

The area distortion \(\sigma_1\sigma_2\) is simply given by the ratio of the deformed triangle's area to the original triangle's area. Using Heron's formula, one can show that the area of a triangle is given by \[\text{area}_{ijk} := \tfrac{1}{2\sqrt2}\sqrt{\left(\ell^2\right)^T \!\!\! A \ell^2},\] where \(\ell^2\) denotes the vector of squared edge lengths $(\ell_{ij}^2, \ell_{jk}^2, \ell_{ki}^2)^T$, and \(A\) is the matrix \[ A = \frac 12 \begin{pmatrix}-1 & 1 & 1 \\ 1 & -1 & 1\\ 1 & 1 & -1\end{pmatrix}.\] Hence, the area distortion is given by \[\sigma_1\sigma_2 = \sqrt{\frac{(\tilde \ell^2)^T A \tilde \ell^2}{(\ell^2)^T A \ell^2}}.\]

Symmetrized Anisotropic Distortion

Before considering the anisotropic distortion \(\sigma_1/\sigma_2\) itself, we begin with a symmetrized version \(\frac{\sigma_1}{\sigma_2} + \frac{\sigma_2}{\sigma_1}\), which is given by a similar formula: \[ \frac{\sigma_1}{\sigma_2} + \frac{\sigma_2}{\sigma_1} = \frac{\left(\ell^2\right)^T \!\!\! A \tilde \ell^2}{\sqrt{\left(\ell^2\right)^T \!\!\! A\, \ell^2}\sqrt{(\tilde \ell^2)^T \! A \tilde \ell^2}}. \] This formula can also be written directly in terms of the angles \(\alpha_i\) as \[ \frac{\sigma_1}{\sigma_2} + \frac{\sigma_2}{\sigma_1} = \begin{pmatrix} \cot \alpha_i\\\cot\alpha_j\\\cot\alpha_k\end{pmatrix}^T A^{-1} \begin{pmatrix} \cot \tilde\alpha_i\\\cot\tilde\alpha_j\\\cot\tilde\alpha_k\end{pmatrix}, \] where \(A^{-1}\) is given by \[ A^{-1} = \begin{pmatrix} 0 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \end{pmatrix}.\] The equivalence of these two expressions follows from the identity \[ \begin{pmatrix} \cot \alpha_i \\ \cot \alpha_j \\ \cot\alpha_k \end{pmatrix} = \frac {A^{-1}\ell^2} {\sqrt{\left(\tilde \ell^2\right)^T \!\!\! A \tilde \ell^2}}, \] (which can itself be derived using the law of cosines, the sine formula for area, and the definition that $\cot \alpha = \tfrac{\cos\alpha}{\sin\alpha}$.)

A derivation of the cotan formula for distortion, courtesy of Boris Springborn, can be found here, and connections to hyperbolic geometry are discussed, e.g., in p.11 of Joshua Bowman's PhD thesis with John H. Hubbard.

Anisotropic Distortion

The anisotropic distortion can easily be computed from the symmetrized anisotropic distortion by solving the quadratic equation \(y = x + \frac 1x\). Concretely, if the symmetrized distortion is given by \(d_s\), then the ordinary anisotropic distortion is given by \[ \frac{\sigma_1}{\sigma_2} = \frac 12 \left(d_s + \sqrt{d_s^2-4}\right) \]

Other distortions

Once you have computed the distortions \(\sigma_1\sigma_2\) and \(\sigma_1 / \sigma_2\), then you can find the values of the two singular values \(\sigma_1\) and \(\sigma_2\) by multiplying and dividing the distortion values respectively. These singular values can then be used to evaluate many other measurements of distortion.

For discussion of intrinsic calculations for distortion in the context of elasticity, see Appendix B of Sassen et al. [2020; arxiv version]

Everyone who discusses Möbius transformations mentions that they map circles to circles, but it can be hard to find concrete equations describing exactly how a given circle is changed by a given Möbius transformation.

Image of a Möbius transformation from Möbius transformations revealed by Arnold & Rogness.

I recently had to derive a formula for the curvature of a circle after applying a certain Möbius transformation, and the resulting formula was surprisingly simple. Suppose we start with a circle centered at \(x \in \mathbb{R}^2\) with radius \(r > 0\). If we apply a Möbius transformation which fixes the unit circle and sends some point \(z\) inside the unit disk to the origin, then the circle's curvature becomes \[ \tilde r = \frac{1-\|z\|^2}{1 + 2\langle x, z \rangle + (\|x\|^2 - r^2)\|z\|^2}r. \]

In particular, I was interested in Möbius transformations fixing the unit circle. If we ignore rotations, which don't affect the curvature of circles anyway, such a Möbius transformations is determined entirely by the point \(z\) which is sent to zero: \[f_z(p) := \frac{p-z}{1-\bar zp}.\] (See e.g. the wikipedia article.) Note that the inverse of \(f_z\) is given by \(f_{-z}\), since \(f_z(f_{-z}(p)) = p\).

Suppose we start with a circle of radius \(r\) centered at a point \(x \in \mathbb{R}^2\). We can encode our circle as the zero level set of a quadratic form \[Q(p) := a \|p\|^2 - \langle b, p \rangle + c = 0,\] with coefficients \(a = 1\), \(b = 2x\), and \(c = \|x\|^2 - r^2\). We can recover the center from \(Q\) as \(x = \tfrac c {2a}\), and the radius as \(r = \tfrac 1 {2a} \sqrt{\|b\|^2 - 4 ac}\).

Now we can determine how the Möbius transformation \(f_{-z}\) affects our circle by finding the zero set of the transformed map \(Q(f_{z}(p))\). The calculation is easier to do if we express \(Q(p)\) using complex numbers as \(Q(p) := a \bar p p - \Re(\bar b p) + c\). If we substitute in the definition of \(f_z(p)\) and do some nasty algebra, we find that \[\begin{aligned}\|1-\bar zp\|^2Q(f_z(p)) &= \left((a + \Re(\bar b z) + c \|z\|^2\right) \|p\|^2\\ &\quad- \Re\left((2 a \bar z + \bar b + b \bar z^2 + 2 c z)p\right)\\ &\quad+ (a \|z\|^2 + \Re(\bar b z) + c) \end{aligned}.\]

That is, up to a scalar multiple, \(Q(f_z(p))\) is itself a quadratic form with coefficients \[\begin{aligned} \tilde a &:= a + \Re(\bar b z) + c \|z\|^2\\ \tilde b &:= 2az + b + \bar b z^2 + 2 c \bar z\\ \tilde c &:= a \|z\|^2 + \Re(\bar b z) + c. \end{aligned}\]

We can now extract the transformed radius from \(Q(f_z(p))\). Another gnarly calculation shows that \[ \|\tilde b\|^2 - 4 \tilde a \tilde c = (\|b\|^2 - 4 a c)\left(1-\|z\|^2\right)^2 = 4a^2 r^2 \left(1-\|a\|^2\right)^2. \] Therefore, the radius of the transformed circle is given by \[ \tilde r = \frac{1-\|z\|^2}{1 + \tfrac 1a\langle b, z \rangle + \tfrac ca \|z\|^2}r, \] as promised above. (The sign on the \(\langle b z\rangle\) term is flipped in the earlier expression above since there we applied \(f_z\) to the circle, rather than \(f_{-z}\).)

The \(1-\|z\|^2\) term suggests that there may be a slick proof of this formula using hyperbolic geometry: Möbius transformations which fix the unit circle are precisely the isometries of the Poincaré disk model of the hyperbolic plane, whose metric tensor is given by: \[ds^2 = \frac{4 dz^2}{\left(1-\|z\|^2\right)^2}.\] But I'm happy with this concrete calculation for now.