# Assignment 3 (Coding): The Laplacian and Curvature Flow

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.js:
• laplaceMatrix
• massMatrix
2. projects/poisson-problem/scalar-poisson-problem.js:
• constructor
• solve
3. projects/geometric-flow/mean-curvature-flow.js:
• buildFlowOperator
• integrate
4. projects/geometric-flow/modified-mean-curvature-flow.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 = \tfrac{1}{2}\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 \rightarrow \mathbb{R}^3$ as a single DenseMatrix of size $|V| x 3$, and solve with the same (scalar) cotan-Laplace matrix used for the previous part.
• 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 rename your geometry.js, scalar-poisson-problem.js, mean-curvature-flow.js and modified-mean-curvature-flow.js files to geometry.txt, scalar-poisson-problem.txt, mean-curvature-flow.txt and modified-mean-curvature-flow.txt (respectively) and put them in a single zip file called solution.zip. These files and your solution to the written exercises should be submitted together in a single email to Geometry.Collective@gmail.com with the subject line DDG17A3.

## 31 thoughts on “Assignment 3 (Coding): The Laplacian and Curvature Flow”

1. yangxio says:

Coding 2. projects/poisson-problem/scalar-poisson-problem.js:
I am doing the coding assignment and realize that in addition to “solve(rho)”, we also need to implement “constructor(geometry)” as well in order for the code to run. Would you please confirm this point as well? And similarly, whenever we are asked to implement a routine, are we supposed to figure out on our own all necessary dependencies and fill them out as necessary?
Thanks for the clarification.

1. rohansawhney says:

Yes, you are supposed to make very minor edits to the constructor in projects/poisson-problem/scalar-poisson-problem.js – this should be as simple as replacing the placeholders with calls to functions you’ve already implemented. Notice that the constructor builds important matrices such as the cotan Laplace ahead of time, before any calls to solve. This is an important optimization for solving a linear system of equations Ax = b where the right hand side b changes frequently but the matrix A remains the same.

In general, you won’t have to figure out the dependencies between functions, except with the constructor in the same file. I’ll make sure to mention this in the assignment writeups.

2. yangxio says:

Question out of curiosity:
I want to see how badly/well (depending on perspectives) our discrete flows approximate the theoretically continuous ones when applied to the classical setup of soap bubble on two loops, especially for the critical case, the local minima, and the point at which the bubble breaks. So could you please let us know how to create our own meshes (.obj and other necessary files) as inputs for “index.html”? Or even better, could you please make such files for us to play? I believe that it is better for us to see a simple visual example whose “solution” is well-known (intuitively) to compare with.
(P.S. to the bunny: No offense, you are cute! But soap bubble is more classical and we can make it in the bathroom.)

1. rohansawhney says:

If you want to simulate something similar to what’s happening in this video (https://www.youtube.com/watch?v=mziis4pbBOw), then you can try out your discrete flows on these cylindrical meshes (https://drive.google.com/a/andrew.cmu.edu/file/d/0B_HXoJovrElnb3JVZzJhSzdFUjQ/view?usp=sharing, https://drive.google.com/a/andrew.cmu.edu/file/d/0B_HXoJovrElnYlJHLThkbFdtQWc/view?usp=sharing). However, you will have to modify your flow operator so that the boundary vertices remain fixed. You can model your own meshes with tools such as blender (https://www.blender.org). Also, the “input” folder in the code repo has a bunch of meshes you can load and try your flows on. Let us know how it goes!

3. anobani says:

i have coded the laplace Matrix according the specification, however i am not passing the test, for the entries of adjacency, should the all be +1 or -1, as specified we need positive definite matrix by multiplying the entries by -1, can’t we just set the entries to 1 ?

4. plav says:

I am stuck in the beginning of the exercise.
1) Is the cotan-laplace matrix the same as the laplacian matrix? The textbook cites this sparse matrix as the representation of the cotan-Laplace operator… but the Laplacian matrix itself is only the matrix of degrees – adjacency. I can imagine how an algorithm could use laplacian matrix to calculate the cotan-laplace, but I do not know if this is what I should code.

1. plav says:

I have just read again… so it has nothing to do with the degree – adjacency matrix.

1. rohansawhney says:

The adjacency matrix in the notes is not the cotan-Laplace matrix. Refer to the formula in the assignment writeup above to implement the cotan-Laplace matrix.

5. ahagerty says:

My code for mean-curvature-flow is working, but I must be confused on some aspect of modified-mean-curvature-flow.

In constructor I define this.A using geometry.laplaceMatrix, and then in buildFlowOperator I use this.A to define the flow operator using the same formula as in mean-curvature-flow. However, the tests aren’t going through.

My understanding is that the only difference should be that when defining the flow operator, we don’t want to rebuild the Laplace matrix, that’s why we store the original Laplace matrix in this.A instead of calling geometry.laplaceMatrix inside of buildFlowOperator.

Am I missing something? There’s so little code that I don’t understand where I could have a bug.

1. rohansawhney says:

Seems like you are on the right track. It’s difficult for me to say what’s going wrong without looking through your code. I’d be happy to help you further in office hours 🙂

1. ahagerty says:

Yup, looks like my error was just something to do with javascript… I had to basically “clone” this.A instead of using it directly, once I did that everything seems to work.

2. yangxio says:

Maybe there is some bug like pass-by-reference vs pass-by-value issue that possibly messes the Laplace matrix up.

3. minkyuk says:

So did simply making this.A non-alias solve the issue?

6. anobani says:

i am trying to solve the buildFlowOperator, but i feel lost in the equations, since we need to define the curvature we can relate to the mass matrix since it encodes the barrycentricArea, in which can define the change of the curvature according to the change of that area- in which is pushing the vertex with or against the normal, so can i assume that $f_0$ in the equation $(f_h – f_0)/h = \Delta f$ is the initial mass matrix ( please correct me if this logic is wrong). i still don’t know the logic of the flow of the equations in the course notes page (94) , as the following:
1. $(f_h – f_0)/h = \Delta f$
2. $f_h = f_0+ h\Delta f_0$

in the second equation why we have $\Delta f_0$ rather than $\Delta f$. and i am still by reasoning over the equations and the use of Laplace matrix to solve this, will it help me find the value of $f_h$ ?

1. anobani says:

my apology, just mixed between the mass matrix and the Laplace, but still want to ask about the two equations and the relation of $\Delta f$ and $\Delta f_0$

1. rohansawhney says:

1) $f$ represents the positions of the vertices of the mesh – in the discrete case, it is a|V|x3 matrix.
2) $\Delta$ in Section 6.6 of the course notes equals the inverse of the mass matrix times the cotan-Laplace given above. This matrix is not necessarily symmetric, which means that it cannot be decomposed with Cholesky factorization (However, LU factorization will still work). If you read Section 6.3 carefully, you’ll notice that a common trick is to move the mass matrix to the right hand side. More specifically, $\Delta x = b$ can be rewritten as $Lx = Mb$ where $L$ is cotan Laplace matrix and $M$ is the mass matrix.
3) We have the freedom to choose the argument of $\Delta$ in the equation $(f_h – f_0)/h = \Delta f$. Setting $f$ to $f_0$ gives us the forward euler (https://en.wikipedia.org/wiki/Euler_method) integration scheme, while setting it to $f_h$ gives us the backward euler integration (https://en.wikipedia.org/wiki/Backward_Euler_method). It turns out that backward Euler is a lot more stable than forward Euler. We would like you to use backward Euler for this assignment.

7. plav says:

Hi Rohan,
For me, the few paragraphs in the book plus this description of the problem are not enough. Follow some topics:

1) I assume that $(I- h\Delta)f_h = f_0$ is the system we are trying to solve.
2) If (1) is correct, the buildFlowOperator should solve this: $(I- h\Delta)$. However, there should be no mass involved here…. and the function has the mass matrix as the argument.
3) if (1) is correct, I can solve this system as I did for the Poisson.
3.1 In this case, I should guarantee that $f$ has derivative 0, by subtracting it by $\bar{f}$.
3.2 Then, I should solve the system as $Ax = -M(f – \bar{f})$
3.3 In this case, the method integrate(h) is doing exactly (3.2)
4) When you say I should treat $f$ as a single DenseMatrix of size $|V|\times 3$, do you mean that I should create a $|V|\times 1$ matrix with an array of the respective array $[x, y, z]$ of the point at each matrix position*, or we should decompose this array and effectively create a $|V|\times 3$ DenseMatrix with numbers at each matrix position.

Can you check if this is the right idea?
If it is not, can you give me some hint?

Thanks

*assuming this is possible in javascript

1. rohansawhney says:

1) Yes, this is the system we are trying to solve.
2) buildFlowOperator should construct $I – h\Delta$, but look at my comments to anobani’s questions – I describe how $\Delta$ (as described in the notes) can be decomposed. Also, be careful about sign here: our cotan-Lapalce matrix is positive definite, rather than negative semi definite as in the notes.
3) No, we only need to remove the constant component from the right hand side for a Poisson problem with the $\Delta$ operator. This analysis does not apply to arbitrary systems of equations $Ax = b$. For mean curvature flow, $A$ does not equal $\Delta$. Simply solve $(I – h\Delta)f_h = f_0$ using Cholesky factorization.
4) The DenseMatrix class cannot store arbitrary Javascript objects such as arrays. As the writeup suggests, you should encode the x, y, z components of $f$ as a |V|x3 DenseMatrix.

1. plav says:

Ok. It makes sense.

However, now that I can sketch this operation, I am lost with the geometry and mesh datastructures.

1) Originally, when I talked with you about vertexIndex, you told me that it is a dictionary that was necessary because the ordering in the mesh.vertices could be different from some real index. Apparently, this real index could change, so this dictionary was necessary to preserve the original mesh.vertices array intact. However, I just checked the function that creates the dictionary and it is basically a more efficient way to recover the indices in mesh.vertices. So there is only one indexing system.

2) The description says that geometry.positions is an array containing the position of each vertex in a mesh. However, it seems to be a dictionary that uses a mesh.vertices[i] as the key and maps to its position (see lines 19-22 of geometry.js)

3) if (1) and (2) are true, the indexing of the rows of the matrices I use to solve the system of equations does not need to involve the dictionary vertexIndex as I am using a loop and I already have the index of the vertex.

All these details of the code + the details of the math make me stay for hours to solve a small bit of the assignment.

1. plav says:

correcting:
3) if (1) and (2) are true, as I am using a loop (let vertex of positions) the dict vertexIndex is necessary to convert the key of positions (which is a vertex) to its index. The resulting indexing system is exactly the same as the indexing for row and cols of the matrices we are creating (Mass, Laplacian, etc.).

1. plav says:

I am stuck here….geometry.positions is an array, but takes a vertex as an index? I checked the class vertex, it has an internal attribute index, but it is only returned as a string representation of the object….

2. rohansawhney says:

1) Yes, all of this is true. The framework has the ability to use different indexing when mapping mesh vertices, edges, faces etc to positions in a matrix, but right now it defaults to a single simple scheme. The reason for this is that none of the current projects in the codebase require any fancy indexing, but that doesn’t mean that the codebase shouldn’t have this ability for future projects!

2) Yes, geometry.positions is a dictionary that maps objects instantiated from the Vertex class to its position (a Vector object) in 3D.

3) Sure, but in general its not good idea to write code that “fights the system” or ignores its abstractions. For example, since vertexIndex is a visible property of the mean curvature flow class, you can’t be certain that its not used elsewhere (like in the rendering code).

8. apoorva says:

So ‘removing the constant component’ would mean ‘ensuring that the entries of the vector add up to 0’ in the discrete case, right? I’m trying to apply this transformation to M * rho.

1. pmichel1 says:

There is a more principled way of removing the constant component at the beginning of the lecture notes on the laplacian

1. apoorva says:

Are you referring to the expression $$\left(\int\limits_{M}\rho~dV\right)/|M|$$? That’s what I’m doing in the discrete case to compute M*rhoBar (my M*rhobar is a vector all of whose entries are [sum of the entries of M*rho divided by the total area]). I’m not sure if I got this right conceptually, because it doesn’t seem to work…

2. apoorva says:

When I do this, my solution is very close to the reference solution, but it’s not within the 1e-6 bound. Should I try harder to get it closer?

1. rohansawhney says:

A few people are having this problem. Make sure you are adding the positive 1e-8 shift only to the diagonal terms of the Laplace matrix.

1. apoorva says:

Yes, I’m only adding it to the diagonal. I’m not sure why there’s disagreement with the tests (the picture of the bunny looks okay when I add densities to it).
Here’s what I’m doing –
I compute M*rho, and I subtract it from a vector of size numVertices each of whose entries is sum(M*rho)/numVertices. Then I use this resultant vector as the right-hand side when I call this.A.chol().solvePositiveDefinite(…).

1. rohansawhney says:

You’re almost there – look at the smooth expression for rhobar you posted in the comments above and compare it with your discrete calculation, do they agree? What is $|M|$ in the smooth expression?

9. intrepidowl says:

1. For the mean curvature flow, what are the implications for the matrix in the writeup being opposite positive/negative semidefinite? Does that mean we need to multiply all the Laplacians we use by -1 because we reversed sign when constructing the laplaceMatrix?

2. For the mean curvature flow, does the normalize() call in the function matter? I haven’t gotten different results from moving it around.