CS 15-458/858: Discrete Differential Geometry

Note: For the final assignment, you can do either this assignment or A5. Overall, you just need to be sure you completed A0 and A1, as well as 3 of the assignments A2–A6 (your choice which ones. See this Piazza post for more details.) Also, you cannot use late days on this assignment since it’s the last one.

Please implement the following routines in:

1. projects/vector-field-decomposition/tree-cotree.[js|cpp]:
   1. buildPrimalSpanningTree
   2. inPrimalSpanningTree
   3. buildDualSpanningCotree
   4. inDualSpanningCotree
   5. buildGenerators
2. projects/vector-field-decomposition/harmonic-bases.[js|cpp]:
   1. buildClosedPrimalOneForm
   2. compute

In addition, please implement either Hodge Decomposition, or Trivial Connections (on surfaces of genus 0):

Hodge Decomposition: Please implement the following routines in

1. projects/vector-field-decomposition/hodge-decomposition.js:
   1. constructor
   2. computeExactComponent
   3. computeCoExactComponent
   4. computeHarmonicComponent

Trivial Connections on Surfaces of Genus 0: Please implement computeConnections in projects/direction-field-design/trivial-connections.js. This file has function signatures for trivial connections on arbitrary surfaces, but for this assignment we are only requiring you to compute the connection form on simply-connected surfaces where you don’t have to worry about the period matrix. You are, of course, welcome to implement the general algorithm if you would like to!

Notes:

• Recall that homology generators are a set of loops which somehow describe all of the interesting loops on a surface. For example, here are the two homology generators for the torus.
• Your buildGenerators function should implement the Tree-Cotree algorithm, which is Algorithm 5 of section 8.2 of the notes.
• In tree-cotree.js the trees are stored using the dictionaries vertexParent and faceParent. You should store the primal spanning tree by storing the parent of vertex v in vertexParent[v]. The root should be its own parent. Similarly, you should store the dual spanning cotree by storing the parent of face f in faceParent[f].
• The Tree Cotree algorithm finds homology generators on the dual mesh. That is to say, you should take each dual edge which is not in either spanning tree, and form a loop by following its endpoint back up the dual cotree until they meet at the root of the dual cotree. Even though we’re thinking of these homology generators as loops on the dual mesh, we still store them as lists of ordinary halfedges, since edges in the dual mesh are in 1-1 correspondence with edges in the primal graph.
• buildClosedPrimalOneForm should use a homology generator to construct a closed primal 1-form as described in section 8.22 of the notes.
• The compute function in HarmonicBases should compute a harmonic basis using Algorithm 6 in section 8.2 of the notes. If you choose to implement Hodge decomposition for the assignment, you are welcome to use your Hodge decomposition code to solve for $d\alpha$. Otherwise, you can just ignore the hodgeDecomposition parameter and directly solve the linear system $\Delta \alpha = \delta \omega$ to find $\alpha$. (Recall that $\delta$ is the codifferential, which we talked a lot about in the Discrete Exterior Calculus assignment).

Notes (Hodge Decomposition):

• The procedure for Hodge Decomposition is given as Algorithm 3 in section 8.1 of the Notes.
• Note that the SparseMatrix class has an invertDiagonal function that you can use to invert diagonal matrices.
• For computeCoExactComponent, you should compute the coexact component $\delta \beta$ of a differential form $\omega$ by solving the equation $d\delta \beta = d\omega$. As stated in the notes, the most convenient way of doing this is to define a dual 0-form $\tilde \beta := *\beta$, and then to solve $d*d \tilde \beta = d \omega$. Then you can compute $\delta \beta$ using the fact that $\delta \beta = *d*\beta = *d\tilde\beta$. When doing these computations, you should keep tract of whether your forms are primal forms or dual forms, recalling that hodge stars take you from primal forms to dual forms, and hodge star inverses take you from dual forms to primal forms. (THis is discussed in detail in section 8.1.3 of the notes.
• Note that the 2-form Laplacian B is not necessarily positive definite, so you should use an LU solver instead of a Cholesky solver when solving systems involving B.

Notes (Trivial Connections):

• The trivial-connection-js file is structured for implementing the full trivial connections algorithm on arbitrary surfaces. Since that’s tricky, we’re only requiring that you compute connections on topological spheres. You should implement this in the computeConnections() function, and it should pass the tests labeled “computeConnections on a sphere”.
• Even though we’re working with spheres, it is still helpful to use the formulation of Trivial connections given in section 8.4.1 of the notes. In particular, you should solve the optimization problem given in Exercise 8.21.
• On a sphere, there are no harmonic 1-forms, so the $\gamma$ part of your Hodge decomposition will always be 0. Furthermore, the sphere has no homology generators, so the problem simplifies to \[\min_{\delta \beta} \;\|\delta\beta\|^2\;\text{s.t.}\; d\delta\beta = u\]
• Following the notes, we observe that the constraint $d\delta\beta = u$ determines $\beta$ up to a constant, and that constant is annihilated by $\delta$. So you can find the connection 1-form $\delta \beta$ with a single linear solve.

Submission Instructions
Please submit all source files to Gradescope by 5:59pm ET on May 21, per the usual submission instructions.

Warning: You cannot use late days on this assignment since it’s the last one.

Note: For the final assignment, you can do either this assignment or A5. Overall, you just need to be sure you completed A0 and A1, as well as 3 of the assignments A2–A6 (your choice which ones. See this Piazza post for more details.) Also, you cannot use late days on this assignment since it’s the last one.

In this assignment, you will investigate tools for working with tangent vector fields on surfaces. Tangent vector fields are central to classical differential geometry, and have many interesting applications. For this homework, we’ll look at one algorithm for designing vector fields, and along the way we’ll cover a lot of deep facts about surfaces.

There’s no PDF this week since the exercises are all from the notes. The notes will also give you all the background you’ll need to complete this assignment. It builds on a lot of the stuff we’ve already done in the class, especially discrete exterior calculus and the Laplacian.

Submit your written responses for Exercise 8.1 – Exercise 8.21 in the notes to Gradescope (per the usual submission instructions), except for Exercise 8.13.

Warning: You cannot use late days on this assignment since it’s the last one.

Note: For the final assignment, you can do either this assignment or A6. Overall, you just need to be sure you completed A0 and A1, as well as 3 of the assignments A2–A6 (your choice which ones. See this Piazza post for more details.) Also, you cannot use late days on this assignment since it’s the last one.

For the coding portion of your final assignment, you will implement the heat method, which is an algorithm for computing geodesic distance on curved surfaces. All of the details you need for implementation are described in Section 3 of the paper, up through and including Section 3.2. Note that you need only be concerned with the case of triangle meshes (not polygon meshes or point clouds); pay close attention to the paragraph labeled “Choice of Timestep.”

Please implement the following routines in:

1. projects/geodesic-distances/heat-method.[js|cpp]:
   1. constructor
   2. computeVectorField
   3. computeDivergence
   4. compute

Notes

• Refer to sections 3.2 of the paper for discretizations of Algorithm 1 (page 3).
• Recall that our Laplace matrix is positive semidefinite, which might differ from the sign convention the authors use.
• The tests for computeVectorField and computeDivergence depend on the A and F matrices you define in your constructor. So if you fail the tests but your functions look correct, check whether you have defined the flow and laplace matrices properly.
• Your solution should implement zero neumann boundary conditions (which are the “default behavior” of the cotan Laplacian) but feel free to tryout other Dirichlet and Neumann boundary conditions on your own.

Submission Instructions

Please submit your source files to Gradescope by 5:59pm ET on May 21, 2021.

Here’s the writeup for your last assignment, which plays the role of the final in the course. (Note that there is no final exam!)

Note: For the final assignment, you can do either this assignment or A6. Overall, you just need to be sure you completed A0 and A1, as well as 3 of the assignments A2–A6 (your choice which ones. See this Piazza post for more details.) Also, you cannot use late days on this assignment since it’s the last one.

This time, we’re taking off the “training wheels” and having you read a real paper, rather than course notes. Why? Because you’re ready for it! At this point you have all the fundamental knowledge you need to go out into the broader literature and start implementing all sorts of algorithms that are built on top of ideas from differential geometry. In fact, this particular algorithm is not much of a departure from things you’ve done already: solving simple equations involving the Laplacian on triangle meshes. As discussed in our lecture on the Laplacian, you’ll find many algorithms in digital geometry processing that have this flavor: compute some basic data (e.g., using a local formula at each vertex), solve a Laplace-like equation, compute some more basic data, and so on.

Your main references for this assignment will be:

• this video, which gives a brief (18-minute) overview of the algorithm, and
• this paper, which explains the algorithm in detail.

Written exercises for this assignment are found below. Submit your responses to Gradescope by 5:59pm ET, May 21, 2021.

The written part of your next assignment, on conformal surface flattening, is now available below. Conformal flattening is important for (among other things) making the connection between processing of 3D surfaces, and existing fast algorithms for 2D image processing. You’ll have the opportunity to implement one of these algorithms in the coding part of the assignment.

For the coding portion of this assignment, you will build the so-called “cotan-Laplace” matrix and start to see how it can be used for a broad range of surface processing tasks, including the Poisson equation and two kinds of curvature flow.

Getting Started

Please implement the following routines in

1. core/geometry.[cpp|js]:
   • laplaceMatrix
   • massMatrix
2. projects/poisson-problem/scalar-poisson-problem.[cpp|js]:
   • constructor
   • solve
3. projects/geometric-flow/mean-curvature-flow.[cpp|js]:
   • buildFlowOperator
   • integrate
4. projects/geometric-flow/modified-mean-curvature-flow.[cpp|js]:
   • constructor
   • buildFlowOperator

Notes

• Sections 6.4-6 of the course notes describe how to build the cotan-Laplace matrix and mass matrices, and outline how they can be used to solve equations on a mesh. In these applications you will be required to setup and solve a linear system of equations $Ax=b$ where $A$ is the Laplace matrix, or some slight modification thereof. Highly efficient numerical methods such as Cholesky Factorization can be used to solve such systems, but require A to be symmetric positive definite. Notice that the cotan-Laplace matrix described in the notes is negative semi-definite. To make it compatible for use with numerical methods like the Cholesky Factorization, your implementation of laplaceMatrix should instead produce a positive definite matrix, i.e., it should represent the expression
  $$(\Delta u)_i=\frac12 \sum_{ij}(\cot \alpha_{ij}+\cot \beta_{ij})(u_i–u_j).$$(Note that $u_i−u_j$ is reversed relative to the course notes.) To make this matrix strictly positive definite (rather than semidefinite), you should also add a small offset such as $10^{-8}$ to the diagonal of the matrix (which can be expressed in code as a floating point value 1e-8). This technique is known as Tikhonov regularization.
• The mass matrix is a diagonal matrix containing the barycentric dual area of each vertex.
• In the scalar Poisson problem, your goal is to discretize and solve the equation $\Delta \phi = \rho$ where $rho$ is a scalar density on vertices and $\Delta$ is the Laplace operator. Be careful about appropriately incorporating dual areas into the discretization of this equation (i.e., think about where and how the mass matrix should appear); also remember that the right-hand side cannot have a constant component (since then there is no solution).
• In your implementation of the implicit mean curvature flow algorithm, you can encode the surface $f:M \to \mathbb R^3$ as a single DenseMatrix of size $|V| \times 3$, and solve with the same (scalar) cotan-Laplace matrix used for the previous part. When constructing the flow operator, you should follow section 6.6 of the notes. But be careful – when we discretize equations of the form $\Delta f = \rho$, we create systems of the form $A f = M \rho$. So you’ll need to add in the mass matrix somewhere. Also, recall that our discrete Laplace matrix is the negative of the actual Laplacian.
• The modified mean curvature flow is nearly identical to standard mean curvature flow. The one and only difference is that you should not update the cotan-Laplace matrix each time step, i.e., you should always be using the cotans from the original input mesh. The mass matrix, however, must change on each iteration.

Submission Instructions

Please only turn in the source files geometry.[cpp|js], scalar-poisson-problem.[cpp|js], mean-curvature-flow.[cpp|js] and modified-mean-curvature-flow.[cpp|js] . These files should be submitted via the usual mechanism, as described on the assignments page.

These exercises will lead you through two different derivations for the cotan-Laplace operator. As we’ll discuss in class, this operator is basically the “Swiss army knife” of discrete differential geometry and digital geometry processing, opening the door to a huge number of interesting algorithms and applications. Note that this time the exercises all come from the course notes—you will need to read the accompanying notes in order to familiarize yourself with the necessary material (we’ll cover some of this stuff in greater depth later on, but you don’t need it for the assignment!)

For the coding portion of this assignment, you will implement various expressions for discrete curvatures and surfaces normals that you will derive in the written assignment. (However, the final expressions are given below in case you want to do the coding first.) Once implemented, you will be able to visualize these geometric quantities on a mesh. For simplicity, you may assume that the mesh has no boundary.

Getting Started
Please implement the following routines in core/geometry.[js|cpp]:

1. angle
2. dihedralAngle
3. vertexNormalAngleWeighted
4. vertexNormalSphereInscribed
5. vertexNormalAreaWeighted
6. vertexNormalGaussianCurvature
7. vertexNormalMeanCurvature
8. angleDefect
9. totalAngleDefect
10. scalarMeanCurvature
11. circumcentricDualArea
12. principalCurvatures

Notes 

1. The dihedral angle between the normals $N_{ijk}$ and $N_{ijl}$ of two adjacent faces $ijk$ and $ijl$ (respectively) is given by
$$ \theta_{ij} := \text{atan2}\left(\frac{e_{ij}}{\|e_{ij}\|} \cdot \left(N_{ijk} \times N_{jil}\right), N_{ijk} \cdot N_{jil}\right)$$

where $e_{ij}$ is the vector from vertex $i$ to vertex $j$.

2. The formulas for the angle weighted normal, sphere inscribed normal, area weighted normal, discrete Gaussian curvature normal and discrete mean curvature normal at vertex $i$ are
\[\begin{aligned}
N_i^\phi &:= \sum_{ijk \in F} \phi_i^{jk}N_{ijk}\\
N_i^S &:= \sum_{ijk \in F} \frac{e_{ij} \times e_{ik}} {\|e_{ij}\|^2\|e_{ik}\|^2}\\
N_i^A &:= \sum_{ijk \in F} A_{ijk}N_{ijk}\\
KN_i &= \frac 12 \sum_{ij \in E} \frac{\theta_{ij}}{\|e_{ij}\|}e_{ij}\\
HN_i &= \frac 12 \sum_{ij \in E}\left(\cot\left(\alpha_k^{ij}\right) + \cot\left(\beta_l^{ij}\right)\right)e_{ij}
\end{aligned}
\]

where $\phi_i^{jk}$ is the interior angle between edges $e_{ij}$ and $e_{ik}$, and $A_{ijk}$ is the area of face $ijk$. Note that sums are taken only over elements (edges or faces) containing vertex $i$. Normalize the final value of all your normal vectors before returning them.

3. The circumcentric dual area at vertex $i$ is given by
\[A_i := \frac 18 \sum_{ijk \in F} \|e_{ik}\|^2\cot\left(\alpha_j^{ki}\right) + \|e_{ij}\|^2\cot\left(\beta_k^{ij}\right)\]

4. The discrete scalar Gaussian curvature (also known as angle defect) and discrete scalar mean curvature at vertex $i$ are given by
\[\begin{aligned}
K_i &:= 2\pi – \sum_{ijk \in F} \phi_i^{jk}\\
H_i &:= \frac 12 \sum_{ij \in E} \theta_{ij}\|e_{ij}\|
\end{aligned}
\]

Note that these quantities are discrete 2-forms, i.e., they represent the total Gaussian and mean curvature integrated over a neighborhood of a vertex. They can be converted to pointwise quantities (i.e., discrete 0-forms at vertices) by dividing them by the  circumcentric dual area of the vertex (i.e., by applying the discrete Hodge star).

5. You are required to derive expressions for the principal curvatures $\kappa_1$ and $\kappa_2$ in exercise 4 of the written assignment. Your implementation of principalCurvatures should return the (pointwise) minimum and maximum principal curvature values at a vertex (in that order).

Submission Instructions

Please submit only the source code file geometry.js or geometry.cpp (depending on whether you’re using JavaScript or C++). This means your code should not depend on edits made to any other files. Then follow the usual hand-in instructions.

The written portion of Assignment 2 can be found below. It takes a look at the curvature of smooth and discrete surfaces, which we have been talking about in lecture.

If you’re typesetting this assignment, feel free to use this template.