Written by Deniz Ozbay, Tal Rastopchin, and Alexander Rougellis
Surface Deformation and Remeshing
Representing a curved surface with a mesh of polygons can prove to be difficult because, given that meshes are discrete surfaces, there are times when small changes to the triangulation of the surface can cause big changes to the overall geometry and measured quantities (as can be seen with the example of the Schwarz Lantern). Sometimes we want to deform surfaces, and the deformation of the surface can result in a “bad” triangulation. Representing these surfaces with such deformations is done by optimizing these deformations \(\delta S\) on the given surface by solving \(\min_{\delta S} f(S + \delta S)\) where \(S\) is that given surface. When the real valued function \(f: S \rightarrow \mathbb{R}\) gives us surface area, our optimization will minimize the surface area, which is also known as mean curvature flow. Minimizing surfaces can be done with gradient descent (given a small triangle mesh) and Newton’s Method (given a large mesh). After minimizing, if the surface is determined to be “bad”, then remeshing is needed using remeshing operations such as edge flip, edge split, and edge collapse (to name a few basic operations that could be used).
Mean Curvature Flow
Given that one of the goals of this project is to optimize functions on meshes, our first task was to put together a simple implementation of mean curvature flow in MATLAB. Professor Etienne Vouga introduced mean curvature flow as a deformation of a surface that minimizes surface area. He explained that if we have a function for the surface area of a mesh, as well as the gradient, we can use an optimization method like gradient descent in order to compute the deformation induced by the mean curvature flow. After our introduction to geometry processing during the tutorial week we knew that we could use the doublarea function to write a function that returned the surface area of a mesh. However, computing the gradient of this function is tricky—what would even be the domain of the surface area function? If the doublearea function relies on computing triangle areas using the cross product, and we are summing the result over all triangles in a mesh, is the domain some sort of “collection of triangles” that we are summing over?
To answer this question, Professor Etienne Vouga pointed us to the paper “Can Mean-Curvature Flow be Modified to be Non-singular?” and explained that we could express the gradient of the surface area function as the Laplacian applied to the \(x\), \(y\), and \(z\) vertex coordinate columns. In particular, the paper explained that “informally, mean-curvature flow can be thought of as a flow that pushes a point on a surface towards the average position of its neighbors.” The paper specifically explains that when \(M\) is a two dimensional manifold, \(\Phi_t : M \rightarrow \mathbb{R}^3\) is a smooth family of immersions of the manifold \(M\), and \(g_t( \cdot, \cdot)\) is a metric induced by the immersion at time \(t\), we have that \(\Phi_t\) is a solution to the mean-curvature flow if
\(\frac{\partial \Phi_t}{\partial t} = \Delta _t \Phi _t\)
where \(\Delta_t\) is the Laplace-Beltrami operator defined with respect to the metric \(g_t\). If we interpret the Laplacian as a local averaging operator, this equation exactly captures the idea that mean-curvature flow is a deformation that pushes a point on a surface towards the average position of its neighbors.
For our implementation, Professor Vouga explained that instead of worrying about the metric induced by the flow we could just compute the Laplacian matrix for the first step and use it throughout the entire simulation. One thing we learned was that sometimes for computation and derivation it can be easier to look at the same function, like surface area, from a bunch of different perspectives. For example, once we knew that the gradient of the surface area was the Laplacian of the \(V\) matrix, we could compute the Hessian by reshaping the \(V\) matrix into a column vector and reshaping the Laplacian matrix into a larger block matrix with 3 copies of the Laplacian matrix along the diagonal. Computing the Hessian in this way allowed us to also try to implement a Newton’s method approach to minimize surface area.
After a lot of discussion, programming, testing, and playing around with initial conditions, we finally got our MATLAB implementation of mean curvature flow to work.
A figure displaying the open cylinder mesh used as a seed surface for the mean curvature flow. 
A figure displaying the catenoid mesh resulting from the mean curvature flow.
The figure on the left is a cylinder that was our initial surface and the figure on the right is the resulting minimal surface produced by the mean curvature flow. It was really cool to see the end result—the resulting minimal surfaces can be very pretty. We’re excited to learn more about optimization on surfaces as well as how remeshing could potentially improve this simulation process.
Current and Future Work
Currently, we are working on implementing Newton’s method and the gradient descent method in the same program as remeshing. Our algorithm runs Newton’s method as long as it decreases the surface area. If a Newton step does not decrease the surface area, the algorithm then switches to the gradient descent method. Then it remeshes the triangles using edge flipping to obtain a Delaunay triangulation. To build this code, we made use of numerous different structures, like an edge face matrix, edge vertex matrix and a vertex face matrix. Then, we were able to implement the edge flipping condition. At this step we first find the neighboring two faces of the edge we check to flip. We check if the angles at the third vertices of the corresponding two faces is bigger than \(\pi\) radians, and if so, the edge is flipped to get the remeshing closer to a Delaunay triangulation. This allows us to optimize the surface and do the remeshing at the same time, which we think is very cool!
While implementing our algorithm, we did lots of debugging and got some funny shapes. Of course, one of the best parts about trying to make the code work was working together as a team. We had a lot of fun trying to decipher what went wrong when we got surfaces compressing to lines then to funny looking disks. Our next step is to work on coming up with algorithms for edge split and edge collapse, so we can fully implement remeshing on the shape. We are really excited to see what the final product will look like, to try lots of cool surfaces and see what the optimized remeshing for each will be!