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.

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.

Definition

Note that \[\begin{aligned} \text{Л}'(x) &= - \log |2 \sin x| \\ \text{Л}''(x) &= -\cot x \end{aligned}\]

The definition looks very strange. It's probably easiest to think of Л as being defined by the differential equation \[ \text{Л}''(x) = -\cot x \] subject to the initial conditions \[\begin{aligned} \text{Л}(0) &= 0\\ \text{Л}'\left(\frac \pi 2\right) &= 0\\ \end{aligned}\]

Useful Facts

A Circumradius Forumla

We begin with a seemingly-unrelated formula about the circumcircle of a triangle.

Let $t_{ijk}$ be a triangle with edge lengths $\ell_{ij}, \ell_{jk}, \ell_{ki}$ and angles $\alpha_{ij}^k, \alpha_{jk}^i, \alpha_{ki}^j$.

Energy Derivative

Futhermore, let $\lambda_{mn} = 2 \log \ell_{mn}$. We define a function \[ f( \lambda_{ij}, \lambda_{jk}, \lambda_{ki}) = \frac 12 \big( \alpha_{jk}^i \lambda_{jk} + \alpha_{ki}^j \lambda_{ki} + \alpha_{ij}^k \lambda_{ij}\big) + \text{Л}( \alpha_{jk}^i) + \text{Л}( \alpha_{ki}^j) + \text{Л}( \alpha_{ij}^k) \]

Möbius Transformations

Möbius transformations, also called fractional linear transformations, are complex functions of the form \[\varphi:z \mapsto \frac{az + b}{cz + d}\] where $a,b,c,d$ are complex numbers and $\det \mmat a b c d \neq 0$.

If we use projective coordinates on $\C$ (i.e. think of $z \in \C$ as the set of all $[u,v] \in \C^2$ with $z = \frac uv$), then the Möbius transformations become matrices \[\phi: \vvec z 1 \mapsto \mmat a b c d \vvec z 1 = \vvec {az + b}{cz + d} \sim \frac{az + b}{cz + d}\] You can check that multiplying together these matices gives the composition of the corresponding Möbius transformations. Using this perspective, we see that the nonzero determinant condition tells us that Möbius transformations are invertible.

If we multiply $a,b,c,d$ by the same nonzero complex number, then the Möbius transformation does not change. So Möbius transformations can be identified with the projective general linear group $PGL(2,\C)$. Since we are allowed to rescale the matrix entries, we can use the projective special linear group $PSL(2,\C)$ instead. $PSL(2,\C)$ is simple, so it must precisely capture the Möbius transformations.

Now, the projective coordinates don't just describe the complex plane. They add one point to the plane, $[1, 0]$, which turns the complex plane $\C$ into the Riemann sphere $\hat \C$. So a Möbius transformation maps the complex plane onto the sphere, applies some sort of transformation to the sphere, and maps the sphere back onto the plane. Two natural questions are: "What is the map between the plane and the sphere?" and "What maps do we apply to the sphere?".

What is the map between the plane and the sphere?

We can break this map from the plane to the sphere into multiple parts. First, we have a map \[ \begin{aligned} f & : \C \to \C^2 \cong \R^4\\ f & : z \mapsto \vvec z 1 \end{aligned} \] Next, we project from $\R^4$ to $\S^3$ by identifying positive scalar multiples of each other (since the image of $f$ does not include $0$, this is okay). This corresponds to identifying vectors in $\C^2$ which are positive real scalar multiples of each other.

Finally, we project from $\S^3$ to $\S^2$ with the Hopf fibration. This corresponds to identifying vectors in $\C^2$ which differ in phase (i.e. a complex number of unit norm).

Alternatively, we can think of this map as stereographic projection. If we first quotient $\C^2$ by phase, we can identify the image of the complex plane under $f$ with the plane of height-1 in $\C \times \R \cong \R^3$. Then, quotienting out by positive scalar multiplication is precisely stereographic projection onto the unit sphere.

What maps do we apply to the sphere?

We apply a linear map of determinant 1 to $\C^2$. TODO

Conformal Maps

Conformal maps have many nice properties. Directly from the definition, we can see that conformal maps preserve angles. I don't currently understand the more complicated nice properties.

Holomorphic Functions

The Cauchy-Riemnann equations give a necessary and sufficient condition for a function to be holomorphic. Let $f(x + iy) = u(x,y) + i v(x,y)$ where $u, v : \R^2 \to \R$. Then $f$ is holomorphic if and only if \[\begin{aligned} \pd u x &= \pd v y\\ \pd u y &= -\pd v x \end{aligned}\]

We can express the Cauchy-Riemnann equations nicely using the Wirtinger derivative. Let \[\pdo{\bar z} := \frac 12 \left(\pdo x + i \pdo y\right)\] Then the Cauchy-Riemann equations are simply the statment \[\pd f {\bar z} = 0\]

Now, let's look at $\C$ as a 2-dimensional real manifold with the standard metric using the identification $\C \sim \R^2$. The differential of $f = u+iv: \C \to \C$ is given by the jacobian \[f_* = \mmat {\pd ux} {\pd vx} {\pd uy} {\pd vy}\] By the Cauchy-Riemann equations, this must have the structure \[f_* = \mmat a {-b} b a\] where $a = \pd u x, b = \pd u y$ are real.

Now, we can compute the pullback of the Euclidean metric on $\C$. \[\begin{aligned} f^*g(v_1, v_2) &= g(f_*v_1, f_*v_2)\\ &= (f_*v_1)^T (f_*v_2)\\ &= v_1^T f_*^T f_* v_2 \end{aligned}\] So the pullback of the standard metric on $\C$ is given by $f_*^T f_*$. Using our expression for $f_*$, we see that \[\begin{aligned} f_*^Tf_* &= \mmat a b {-b} a \mmat a {-b} b a\\ &= \mmat{a^2 + b^2} 0 0 {a^2 + b^2}\\ &= (a^2 + b^2)\mathbb{I} \end{aligned}\] So the pullback of the metric is a scalar multiple of the metric. Thus, the Cauchy-Rimann equations tell us that holomorphic maps are confomal!

Conversely, suppose $f : \C \to \C$ is conformal. Let \[f_* = \mmat a b c d\] Again, the pullback of the metric is \[\begin{aligned} f_*^Tf_* &= \mmat a c b d \mmat a b c d\\ &= \mmat {a^2 + c^2} {ab+cd} {ab+cd} {b^2+d^2} \end{aligned}\] Since $f$ is conformal, we know that this is a scalar multiple of the identity. So $ab+cd = 0$, and $a^2 + c^2 = b^2 + d^2$.

You can solve this to find two solutions: $a=d, b=-c$ (in which case $f$ is holomorphic), or $a=-d, b=c$ (in which case $f$ is antiholomorphic).

So the only conformal maps $\C \to \C$ are holomorphic and antiholomorphic functions. In particular, the only orientation-preserving conformal maps are holomorphic functions.

Automorphisms of the Disk

This machinery of holomorphic functions allows us to nicely characterize the conformal automorphisms of the open unit disk (i.e. invertible conformal maps $D \to D$). As we saw above, it suffices to characterize holomorphic automorphisms of the disk.

Finally, we can classify holomorphic automorphisms of the disk. We can conclude that all orientation-preserving conformal maps of the unit disk to itself are given by these transformations.

Today's post is a little bit different than usual, but it's somewhat math-related, so I figured I'd post it here anyway. Recently, I folded Eric Joisel's origami dwarf

My folded dwarf

(you can see Joisel's rough instructions for the dwarf here). One interesting feature of the dwarf that Joisel points out in his instructions is that the dwarf is folded out of a 28 by 28 grid. As Joisel observes, usually origami models use grids whose dimensions are powers of 2 - it's simple to fold a piece of paper in half repeatedly to obtain an 8 by 8 grid, or a 32 by 32 grid. But 28 by 28 is trickier. In fact, Joisel advises you to use a ruler to form the grid instead of bothering to fold 28ths by hand. But it turns out that it's not so hard to fold 28ths after all. That's what I'm writing about today. But before jumping straight into folding 28ths, we'll start with a slightly easier topic.

Folding a Square in Thirds

Here is a nice little folding sequence to fold a piece of paper in thirds. First, take your square and fold it in half. Unfold, and you are left with a vertical crease cutting the square in half.

Next, fold and unfold the square in half diagonally. Now, fold and unfold a crease from the bottom right corner to the middle of the top edge. And now you're done! The two diagonal creases intersect each other at a point one third of the way across the paper!

Why Does This Work?

There's probably some sort of clever argument you can make using Euclidean geometry and similar triangles and the like to show that this algorithm really does find you one third of the paper. But I think it's easier to use coordinate geometry instead. Let's imagine our square of paper as living in the plane so that its right edge is the $y$ axis and the bottom edge is the $x$ axis.

Note that the two diagonal creases that we folded lie along the lines $y = x + 1$ and $y = -2x$. Now, we can solve for their intersection. \[\begin{aligned} x + 1 &= -2x\\ 3x &= -1\\ x &= -1/3 \end{aligned}\] we find that they intersect at $x = -1/3$. That's one third of the way across the paper!

Generalizing to Arbitrary Fractions

The simple fact that $3 = 2+1$ played a crucial role in our proof above. The $2x$ on the right hand side and the single $x$ on the left hand side combined to give us a factor of $3x$. And that $3$ became the denominator of $1/3$. So what would happen if our right hand side were $-4x$ instead of $-2x$?

\[\begin{aligned} x + 1 &= -4x\\ 5x &= -1\\ x &= -1/5 \end{aligned}\]

Then, instead of finding a point $1/3$ of the way across the page, we would find a point $1/5$ of the way across the page! In general, if we can fold a line with a slope of $-n$, we can then fold the paper into segments of width $1/(n+1)$.

And how do we fold a line of slope $-n$? Earlier we folded a line with slope $-2$ by first folding the paper in half, and then folding a diagonal cutting one of the halves in half. This creates a line of slope $-2$ because half of a square is a $2:1$ rectangle, and its diagonal has slope $-2$. Similarly, we can use an $n:1$ rectangle to fold a diagonal of slope $-n$.

So given a fold $1/n$ of the way across the paper, we can find $1/(n+1)$ as follows: Suppose we start with a square that has a crease $1/n$ of the way across.

Next, fold and unfold the square in half diagonally. Now, fold and unfold a crease from the bottom right corner to the top of our starting crease. And now you're done! The two diagonal creases intersect each other at a point $1/(n+1)$ of the way across the paper!

Folding 28ths

This procedure gives us a straightforward, if tedious method of folding a square into 28ths: First fold it in half, then find $1/3$, then use $1/3$ to find $1/4$, then use $1/4$ to find $1/5$, and so on, until we finally use $1/27$ to find $1/28$. Of course, this is a terrible idea for several reasons. It would take a long time to fold, and would leave countless extra creases on your square. With a little bit of thought, we can fold 28ths with far less effort, and making minimal extra creases.

$28 = 4 \cdot 7$. Folding things in quarters is easy: just fold in half twice. So the only difficult part of folding 28ths is folding 7ths. $7 = 6 + 1$, so we can obtain $1/7$ by first folding $1/6$. And $1/6$ is just half of $1/3$, which we already know how to fold. Here is whole folding sequence:

First, take your square and fold it in half. Unfold, and you are left with a vertical crease cutting the square in half.

Next, fold and unfold the square in half diagonally. You only need to make a strong crease in the top right. Now, fold the diagonal from the bottom right corner to the middle of the top edge. Make a pinch where this crease intersects the other diagonal crease. As we saw earlier, this intersection is $1/3$ across the paper. Now, fold the right edge to the intersection you just made, pinching at the top of the paper. This creates a pinch $1/6$ across the paper. Now, fold the a diagonal from the bottom right corner to the top of the $1/6$ pinch you just made. Pinch where this diagonal intersects your original diagonal. This intersection is $1/7$ across the paper. Finally, you can fold the right edge to your $1/7$ intersection to create $1/14$, and you can fold the right edge to the $1/14$ crease to create $1/28$. Then you're done!