(This note was originally written by Keenan Crane; at this point all responsibility for content and/or errors lies with Peter Schröder)
Parallel Transport
Suppose we have a tangent vector \(u = df(X)\) sitting on an immersed surface \(f(M)\). How do we move from one point of the surface to another while preserving the direction of \(u\)? If \(f(M)\) is completely flat (like the plane itself) then the most natural thing is to slide \(u\) from one point to the other along a straight path—keeping the angle with some reference direction fixed—to obtain a new vector \(u^\prime\). This process is called parallel transport, and the tangent vectors \(u\) and \(u^\prime\) are, as usual, said to be parallel. Parallel transport on a curved surface is a bit trickier. If we keep \(u\) pointing in the same direction, then it ceases to be tangent and now sticks out the side. On the other hand, if we instead keep \(u\) flat against the surface then we no longer have a consistent, global reference direction. Overall, the notion of “same direction” is not very well-defined!
Since there is no notion of “same direction” when comparing vectors on different tangent spaces of a surface, the notion of parallel becomes very flexible. On a surface, one explicitly specifies a parallel transport map along a curve \(\gamma\) connecting point \(p\) to point \(q\)
\[P_\gamma:T_pM\to T_qM\]
that “teleports” vectors from the tangent plane \(T_pM\) to the tangent plane \(T_qM\) isometrically (vectors can be arbitrarily rotated, but not stretched or sheared). We could say that by definition two vectors \(X\in T_pM\) and \(Y\in T_qM\) are parallel (along \(\gamma\)) if \(P_\gamma(X)=Y\). Note that the notion of parallelity is path dependent: the result of parallel transport depends on the path \(\gamma\) connecting \(p\) and \(q\).
The parallel transport map \(P_{\gamma}\) is required to satisfy the axiom that \(P_{\gamma_1\cdot\gamma_2} = P_{\gamma_1}\circ P_{\gamma_2}\) where $\gamma_1\cdot\gamma_2$ is the concatenation of paths, and \(P_{\gamma_1}\circ P_{\gamma_2}\) is the concatenation of linear maps. Other than that, \(P_{\gamma}\) can be designed arbitrarily. A parallel transport defines how to “connect” a tangent space to its neighboring tangent spaces. Hence a choice of parallel transport is also called a connection.
There are a few special types of connections. If the parallel transport is defined in such a way that is path independent, then we call it a flat connection or trivial connection. Another special connection which perhaps mimics parallelity for the flat space the most is the Levi-Civita connection. In the Levi-Civita connection, if \(\gamma\) is a geodesic connecting \(p\) and \(q\), then the angle made between \(X\) and \(\gamma’\) at \(p\) is the same as that made between \(P_{\gamma}X\) and \(\gamma’\) at \(q\).
Remark. An alternative way of describing connection is to describe what it means for vectors to be parallel infinitesimally. Given a connection, i.e. a notion of parallel transport, we can measure the rate of change of a vector field \(Y\) along a certain direction \(X\) (directional derivative of a vector field)
\[\nabla_X Y:={d\over dt}P_{\gamma(t)}^{-1}(Y(\gamma(t)))\Big|_{t=0},\quad \gamma(0)=p, \gamma'(0) = X.\]
A vector field \(Y\) is parallel along a curve \(\gamma\) if \(\nabla_{\gamma'(t)}Y = 0\) for all \(t\). The derivative operator \(\nabla\) is called the covariant derivative or the connection derivative with respect to a connection.
Now, any other connection is different from a connection by a 1-form. How? Suppose \(\tilde P_\gamma\) is the parallel transport of this new connection, then generically we may write \(\tilde P_\gamma = e^{\int_{\gamma} \omega J}\circ P_\gamma\); that is, the new parallel transport \(\tilde P\) is done by first parallel transport using \(P\), and then perform an additional rotation by angle \(\int_\gamma\omega\) for some angle-valued 1-form \(\omega\). Infinitesimally, two connection derivatives are differ by
\[\tilde\nabla = \nabla – \omega J.\]
Therefore, if we already have a canonical connection \(\nabla\), say Levi-Civita connection, then all other connections are parameterized by an angle-valued 1-form \(\omega\).
Discrete Connections
How should we specify a connection in the discrete setting? Well, for a given a pair of triangles \((f_{ijk},f_{jil})\), we can imagine rigidly unfolding them into the plane, translating a vector from one to the other, applying a rotation by some small angle \(\theta_{ij}\), and then rigidly “refolding” these triangles into their initial configuration, as illustrated above. In other words, we can describe a connection on a triangle mesh via a single angle \(\varphi_{ij} \in \mathbb R\) for each oriented dual edge in our mesh. We should also make sure that \(\varphi_{ji} = -\varphi_{ij}\). In other words, the motion we make going from face \(f_{jil}\) to face \(f_{ijk}\) should be the opposite of the motion from \(f_{ijk}\) to \(f_{jil}\). Enforcing symmetry ensures that our notion of “parallel” is consistent no matter which direction we travel. The whole collection of angles \(\varphi \in \mathbb R^{|E|}\) is called a discrete connection.
By the way, does this object sound familiar? It should! In particular, we have a single number per oriented dual edge, which changes sign when we change orientation. In other words, \(\varphi\) is a real-valued, dual, discrete \(1\)-form.
The Levi-Civita Connection
In terms of the picture above, we said that an angle \(\varphi_{ij}=0\) means “just translate; don’t rotate.” If we set all of our angles to zero, we get a very special connection called the Levi-Civita connection. (Those with some geometry background should note that a discrete connection really encodes the deviation from Levi-Civita; it should not be thought of as the connection itself.) The Levi-Civita connection effectively tries to “twist” a tangent vector as little as possible as it moves it from one point to the next. There are many ways to describe the Levi-Civita connection in the smooth setting, but a particularly nice geometric description is given by Kobayashi:
Theorem (Kobayashi): The Levi-Civita connection on a smooth surface is the pullback under the Gauss map of the Levi-Civita connection on the sphere.
What does this statement mean? First, recall that the Gauss map \(N: M \rightarrow S^2\) takes a point on the surface to its corresponding unit normal—this normal can also be thought of as a point on the unit sphere. And what’s the Levi-Civita connection on the sphere? Well, we said that Levi-Civita tries to “twist” vectors as little as possible. On a sphere, it’s not hard to see that the motion of least twist looks like a rotation of the tangent plane along a great arc in the direction \(Z\) of parallel transport. More explicitly, we want a rotation around the axis \(N \times Z\), where \(N\) is the normal of our initial tangent plane. Altogether, then, Kobayashi’s theorem says the following. If we start out with a tangent vector \(\tilde{X}\) on our surface and want to transport it in the direction \(\tilde{Z}\), we should first find the tangent plane with normal \(N\) on the sphere, and the two corresponding tangent vectors \(X\) and \(Z\). (Extrinsically, of course, these are just the same two vectors!) We can then apply an (infinitesimal) rotation along the great arc in the direction \(Z\), dragging \(X\) along with us.
Exercise: Use Kobayashi’s theorem to justify the “unfold, translate, refold” procedure that is used to define the discrete Levi-Civita connection.
Holonomy
At this point you may be wondering what all this stuff has to do with vector field design. Well, once we define a connection on our mesh, there’s an easy way to construct a vector field: start out with an initial vector, parallel transport it to its neighbors using the connection, and repeat until you’ve covered the surface (as depicted above). One thing to notice is that the vector field we end up with is completely determined by our choice of connection. In effect, we can design vector fields by instead designing connections.
However, something can go wrong here: depending on which connection we use, the procedure above may not provide a consistent description of any vector field. For instance, consider the planar mesh below, and a connection that applies a rotation of \(18^\circ\) as we cross each edge in counter-clockwise order. By the time we get back to the beginning, we’ve rotated our initial vector (1) by only \(5 \times 18^\circ = 90^\circ\). In other words, our connection would have us believe that (1) and (6) are parallel vectors!
This phenomenon is referred to as the holonomy of the connection. More generally, holonomy is the difference in angle between an initial and final vector that has been transported around a closed loop. (This definition applies in both the discrete and smooth setting.)
Trivial Connections
To construct a consistently-defined vector field, we must ensure that our connection has zero holonomy around every loop. Such a connection is called a trivial connection. In fact, the following exercise shows that this condition is sufficient to guarantee consistency everywhere:
Exercise: Show that parallel transport by a trivial connection is path-independent. Hint: consider two different paths from point \(a\) to point \(b\).
As a result we can forget about the particular paths along which vectors are transported, and can again imagine that we simply “teleport” them directly from one point to another. If we then reconstruct a vector field via a trivial connection, we get a parallel vector field, i.e., a field where (at least according to the connection) every vector is parallel to every other vector. In a sense, parallel vector fields on surfaces are a generalization of constant vector fields in the plane. But actually, the following exercise shows that any vector field can be considered parallel—as long as we choose the right connection:
Exercise: Show that every discrete vector field (i.e., a vector per face) is parallel with respect to some trivial discrete connection. Hint: think about the difference between vectors in adjacent triangles.
Curvature of a Connection
We can use a trivial connection to define a vector field, but how do we find a trivial connection? The first thing you might try is the Levi-Civita connection—after all, it has a simple, straightforward definition. Sadly, the Levi-Civita connection is not in general trivial:
Exercise: Show that the holonomy of the discrete Levi-Civita connection around the boundary of any dual cell equals the angle defect of the enclosed vertex.
Therefore, unless our mesh is developable, Levi-Civita will exhibit some non-zero amount of holonomy. Actually, you may recall that angle defect is used to define a discrete notion of Gaussian curvature. We can also use a connection to determine curvature—in particular, the curvature of a connection (smooth or discrete) over a topological disk \(D \subset M\) is given by the holonomy around the region boundary \(\partial D\).
Exercise: Show that these two notions of curvature are the same, i.e., show that the curvature of the discrete Levi-Civita connection over any disk \(D\) equals the total discrete Gaussian curvature over that region.
Curvature gives us one tool to test whether a connection is trivial. In particular, a trivial connection must have zero curvature everywhere. For this reason it’s reasonable to say that every trivial connection is “flat.” But is every flat connection also trivial? Well, remember that the curvature of a connection is defined in terms of the holonomy around region boundaries. Any such boundary is called a contractible loop because it can be continuously deformed to a point without “catching” on anything:
In general, there may also be noncontractible loops on a surface that cannot be described as the boundary of any disk. For instance, consider the loop \(\gamma\) pictured on the torus to the left:
In general a surface of genus \(g\) will have \(2g\) “basic types” of noncontractible loops called generators. More precisely, two loops are said to be homotopic if we can get from one to the other by simply sliding it along the surface without ever breaking it. No two distinct generators are homotopic to each other, and what’s more, we can connect multiple copies of the generators to “generate” any noncontractible loop on the surface. For instance, consider the loop \(\gamma^3\), which consists of three copies of \(\gamma\) joined end-to-end. (Formally, the space of loops together with the operation of concatenation describe the first homology group on the surface.)
If we want to check if a connection is trivial, we need to know that it has trivial holonomy around both contractible and noncontractible loops. Equivalently: it must have zero curvature and trivial holonomy around noncontractible loops. As you’re about to demonstrate, though, we don’t need to check all the loops—just a small collection of basis loops.
Exercise: Show that the holonomy around any discrete loop is determined by the curvature at each vertex and the holonomy around a collection of \(2g\) generators.
Singularities
There’s one more issue we run into when trying to find a trivial connection. You may remember the Gauss-Bonnet theorem, which says that \(\sum_{v \in V} d(v) = 2\pi\chi\), i.e., the total Gaussian curvature over a surface equals \(2\pi\) times the Euler characteristic \(\chi\). In fact, this theorem holds if we replace the Gaussian curvature with the curvature of any connection (not just Levi-Civita). But something’s odd here: didn’t we say that a trivial connection should have zero holonomy—hence zero curvature? So unless \(\chi=0\) (i.e., \(M\) is a torus) we have a problem!
Fortunately the solution is simple: we can permit our connection to exhibit nonzero holonomy (hence nonzero curvature) around some loops, as long as this holonomy is an integer multiple of \(2\pi\). Geometrically, a vector parallel transported around any closed loop will still end up back where it started, even if it “spins around” some whole number of times \(k\) along the way. Any vertex where \(k \ne 0\) is called a singularity (see below for some examples). As we’ll see in a moment, singularities actually make it easier to design vector fields with the desired appearance, since one can control the global appearance of the field using only a small number of degrees of freedom.
Vector Field Design
Now on to the fun part: designing vector fields. At this point, you’ve already written most of the code you’ll need! But let’s take a look at the details. To keep things simple we’re going to assume that \(M\) is a topological sphere, so you can forget about non-contractible loops for now.
Our goal is to find a connection \(1\)-form \(\varphi\) such that the holonomy around every loop is zero. If we let \(\varphi_{ij} = 0\) for every dual edge \(e^\star_{ij}\), then the holonomy around any dual cell will be equal to the integrated Gaussian curvature over that cell, which we’ll denote by \(K \in \mathbb R^{|V|}\). Therefore, we need to find angles \(\varphi_{ij}\) such that\[d_0^T\varphi = -K,\]i.e., the holonomy around the boundary of every dual cell should exactly cancel the Gaussian curvature. (Notice that \(d_0^T\) is the discrete exterior derivative on dual \(1\)-forms.) We also need to incorporate singularities. That’s easy enough: we can just ask that the angles \(\varphi_{ij}\) cancel the existing Gaussian curvature, but add curvature corresponding to singularities:\[d_0^T\varphi = -K+2\pi k.\]Here \(k \in \mathbb Z^{|V|}\) is a vector of integers encoding the type of singularity we want at each vertex. To “design” a vector field, then, we can just set \(k\) to a nonzero value wherever we want a singularity. Of course, we need to make sure that \(\sum_i k_i = \chi\) so that we do not violate Gauss-Bonnet.
We now have a nice linear system whose solution gives us a trivial connection with a prescribed set of singularities. One last question, though: is the solution unique? Well, our connection is determined by one angle per edge, and we have one equation to satisfy per dual cell—or equivalently, one per vertex. But since we have roughly three times as many edges as vertices (which you showed earlier on!), this system is underdetermined. In other words, there are many different trivial connections on our surface. Which one gives us the “nicest” vector field? While there’s no completely objective answer to this question, the most reasonable thing may be to look for the trivial connection closest to Levi-Civita. Why? Well, remember that Levi-Civita “twists” vectors as little as possible, so we’re effectively asking for the smoothest vector field. Computationally, then, we need to find the solution to the last equation with minimum norm (since the angles \(\varphi_{ij}\) already encode the deviation from Levi-Civita). As a result, we get the optimization problem\[
\begin{array}{rl}
\underset{\varphi \in \mathbb R^{|E|}}{\min} & ||\varphi||^2 \\
s.t. & d_0^T\varphi = -K+2\pi k.
\end{array}
\]One way to solve this problem would be to use some kind of steepest descent method. However, we can be a bit more clever here by recognizing that this optimization problem is equivalent to looking for the solution to \(d_0^T\varphi = -K+2\pi k\) that has no component in the null space of \(d_0^T\)—any other solution will have larger norm.
Exercise: Show that the column space of \(d_1^T\) is a subspace of the null space of \(d_0^T\). Hint: what happens when you apply \(d\) twice?
Hence, to get the smoothest trivial connection with the prescribed curvature we could (1) compute any solution \(\tilde{\varphi}\) to \(d_0^T\varphi = -K+2\pi k\), then (2) project out the null space component by computing \(\varphi = \tilde{\varphi} – d_1^T(d_1 d_1^T)^{-1} d_1 \tilde{\varphi}\). Overall, then, we get a trivial connection by solving two nice, sparse linear systems. Sounds pretty good, and in fact that’s how the algorithm was originally proposed way back in 2010. But it turns out there’s an even nicer, more elegant way to compute trivial connections, using Helmholtz-Hodge decomposition. (This formulation will also make it a little easier to work with surfaces of nontrivial topology.)
Trivial Connections++
So far, we’ve been working with a \(1\)-form \(\varphi\), which describes the deviation of our connection from Levi-Civita. Just for fun, let’s rewrite this problem in the smooth setting, on a closed surface \(M\) of genus \(g\). Our measure of smoothness is still the total deviation from Levi-Civita, which we can again write as \(||\varphi||^2\). We still have the same constraint on simply-connected cycles, namely \(d\varphi = u\) where \(u = -K + 2\pi k\) (in the smooth setting, we can think of \(k\) as a sum of Dirac deltas). This time around, we’ll also consider the holonomy around the \(2g\) generators \(\Gamma_i\). Again, we’ll use the \(1\)-form \(\varphi\) to “cancel” the holonomy we experience in the smooth setting, i.e., we’ll enforce the constraint\[
\int_{\Gamma_i} \varphi = v_i,
\]where \(-v_i\) is the holonomy we get by walking around the loop \(\Gamma_i\). (In the discrete setting we can measure this quantity as before: transport a vector around each loop by unfolding, sliding, and refolding without any extra in-plane rotation.) Overall, then, we get the optimization problem\[
\begin{array}{rl}
\displaystyle\min_\varphi & ||\varphi||^2 \\
\mathrm{s.t.} & d\varphi = u, \\
& \int_{\Gamma_i} \varphi = v_i,\ i=1,\ldots,2g.
\end{array}
\]Like any other \(1\)-form, \(\varphi\) has a Hodge decomposition, i.e., we can write it as\[
\varphi = d\alpha + \delta\beta + h
\]for some \(0\)-form \(\alpha\), \(2\)-form \(\beta\), and harmonic \(1\)-form \(h\). This expression can be used to simplify our optimization problem:
\[
\begin{array}{rl}
\displaystyle\min_\varphi & ||\delta\beta||^2 + ||h||^2 \\
\mathrm{s.t.} & d\delta\beta = u, \\
& \int_{\Gamma_i} \delta\beta + h = v_i,\ i=1,\ldots,2g.
\end{array}
\]
There are a couple nice things about this reformulation. First, we can find the coexact part by simply solving the linear equation \(d\delta\beta = u\), or equivalently \(d * d *\beta = u\). As with Hodge decomposition, we can make this system even nicer by making the substitution \(\tilde{\beta} := *\beta\), in which case we end up with a standard scalar Poisson problem\[
\Delta \tilde{\beta} = u.
\]We can then recover \(\beta\) itself via \(\beta = *\tilde{\beta}\), as before. Note that the solution \(\tilde{\beta}\) is unique up to a constant, since on a compact domain the only harmonic 0-forms are the constant functions (as you showed when we studied the Poisson equation). Of course, since constants are in the kernel of \(\delta\), we still get a uniquely-determined coexact part \(\delta\beta\).
The second nice thing about this formulation is that we can directly solve for the harmonic part \(\gamma\) by solving the system of linear equations\[
\int_{\Gamma_i} h = v_i – \int_{\Gamma_i} \delta\beta,
\]i.e., since \(\delta\beta\) is uniquely determined by Poisson problem for \(\beta\), we can just move it to the right-hand side.
A slightly nicer way to write this latter system is using the period matrix of our surface. Let \(\Gamma_1,\ldots,\Gamma_{2g}\) be a collection of homology generators, and let \(h_1,\ldots,h_{2g}\) be a basis for the harmonic 1-forms. The period matrix \(\mathsf{P} \in \mathbb R^{2g \times 2g}\) is then given by\[
\mathsf{P}_{ij} = \int_{\Gamma_i} h_j,
\]i.e., it measures how much each harmonic basis “lines up” with each generating cycle. Period matrices have an intimate relationship with the conformal structure of a surface, which we discussed when looking at parameterization. But we don’t have time to talk about that now—we have to compute vector fields! In the discrete setting, we can compute the entries of the period matrix by simply summing up \(1\)-form values over generating cycles. In other words, if \(\Gamma_i\) is a collection of dual edges forming a loop, and \(h_i\) is a dual discrete 1-form (i.e., a value per dual edge), then we have\[
\mathsf{P}_{ij} = \sum_{e^\star_k \in \Gamma_i} (h_j)_k.
\]
Once we have the period matrix, we can express our constraint on generator holonomy as follows. Let \(\mathsf{z} \in \mathbb R^{2g}\) be the coefficients of the harmonic component \(h\) with respect to our basis of harmonic \(1\)-forms, i.e.,\[
h = \sum_{i=1}^{2g} \mathsf{z}_i h_i.
\]Also, let \(\tilde{\mathsf{v}} \in \mathbb R^{2g}\) be the right-hand side of our constraint equation, encoding both the desired generator holonomy as well as the integrals of the coexact part along each generator:\[
\tilde{\mathsf{v}}_i = \mathsf{v}_i – \int_{\Gamma_i} \delta\beta.
\]Then the harmonic component can be found by solving the \(2g \times 2g\) linear system\[
\mathsf{P}\mathsf{z} = \tilde{\mathsf{v}},
\]where the period matrix \(\mathsf{P}\) is constant (i.e., it depends only on the mesh geometry and not the configuration of singularities or generator holonomy).
Overall, then, we have the following algorithm for computing the smoothest vector field on a simplicial surface with a prescribed collection of singularities:
Trivial-Connection++
Require: Vector \(k \in \mathbb{Z}^{|V|}\) of singularity indices adding up to \(2\pi\chi\)
- Solve \(\Delta \tilde{\beta} = u\)
- Solve \(\mathsf{P}\mathsf{z} = \tilde{\mathsf{v}}\)
- \(\varphi \leftarrow \delta\beta + h\)
The resulting \(1\)-form can be used to produce a unit vector field via the procedure described previously. Note that the most expensive part of the entire algorithm is prefactoring the cotan-Laplace matrix, which is subsequently used to compute both the harmonic \(1\)-form bases and to update the potential \(\beta\). In comparison, all other steps (finding generating loops, etc.) have a negligible cost, and moreover can be computed just once upon initialization (e.g., the period matrix \(\mathsf{P}\)). In short, finding the smoothest vector field with prescribed singularities costs about as much as solving a single scalar Poisson problem! If you’ve been paying attention, you’ll notice that this statement is kind of a theme in these notes: if treated correctly, many of the fundamental geometry processing tasks we’re interested in basically boil down to solving a Poisson equation. (This outcome is particularly nice, since in principle we can use the same prefactorization for many different applications!)
