Demo SGI research projects

Minimal Surfaces, But With Saddle Points

By Natasha Diederen, Alice Mehalek, Zhecheng Wang, and Olga Guțan

This week we worked on extending the results described here. We learned an array of new techniques and enhanced existing skills that we had developed the week(s) before. Here is some of the work we have accomplished since the last blog post. 


One of the improvements we made was to create a tiling function which created an \( n^3 \) grid of our triply periodic surfaces, so that we were better able to visualise them as a structure. We started off with a surface inside a \( [-1, 1]^3 \) cube, and then imagined an \(n^3\) grid (pictured below). To make a duplicate of the surface in one of these grid cubes, we considered how much the surface would need to be translated in the \(x\), \(y\) and \(z\) directions. For example to duplicate the original surface in the black cube into the green square, we would need to shift all the \(x\)-coordinates in the mesh by two (since the cube is length two in all directions) and leave the \(y\)- and \(z\)-coordinates unchanged. Similarly, to duplicate the original surface into the purple square, we would need to shift the all \(x\)-coordinates in the mesh by two, all the \(y\)-coordinates by two, and all the \(z\)-coordinates by \(2n\).

Figure 1. A visualization of the the surface tiling.

To copy the surface \(n^3\) times into the right grid cubes, we need find all the unique permutations of groups of three vectors chosen from \((0,2,4, \dots, 2n)\) and add them to the vertices matrix of the of the mesh. To update the face data, we add multiples of the number of vertices each time we duplicate into a new cube. 

With this function in place, we can now see our triply periodic surfaces morphing as a larger structure.

Figure 2. A 3x3x3 Tiling of the Surface.


A skill we continued developing and something we have grown to enjoy, is what we affectionately call “blendering.” To speak in technical terms, we use the open-source software Blender to generate triangle meshes that we, then, use as tests for our code. 

For context: Blender is a free and open-source 3D computer graphics software tool set used for a plethora of purposes, such as creating animated films, 3D printed models, motion graphics, virtual reality, and computer games. It includes many features and it, truly, has endless possibilities. We used one small compartment of it: mesh creation and mesh editing, but we look forward to perhaps experiencing more of its possibilities in the future.

We seek to create shapes that are non-manifold; mathematically, this means that there exist local regions in our surface that are not homeomorphic to a subset of the Euclidean space. In other words, non-manifold shapes contain junctions where more than two faces meet at an edge, or more than two faces meet at a vertex without sharing an edge.

Figure 3. Manifold versus nonmanifold edges and vertices. Source.

This is intellectually intriguing to consider, because most standard geometry processing methods and techniques do not consider this class of shapes. As such, most algorithms and ideas need to be redesigned for non-manifold surfaces. 

Our Blender work consisted of a lot of trial-and-error. None of us had used Blender before, so the learning curve was steep. Yet, despite the occasional frustration, we persevered. With each hour worked, we increased our understanding and expertise, and in the end we were able to generate surfaces we were quite proud of. Most importantly, these triangle meshes have been valuable input for our algorithm and have helped us explore in more detail the differences between manifold and non-manifold surfaces. 

Figure 4. Manifold and Nonmanifold Periodic Surfaces.


The new fellows joining this week came from a previous project on minimal surface led by Stephanie Wang, which used Houdini as a basis for generating minimal surfaces. Thus, we decided we could use Houdini to carry out some physical deformations, to see how non-manifold surfaces performed compared to similar manifold surfaces. 

We used a standard Houdini vellum solver with some modifications to simulate a section of our surface falling under gravity. Below are some of the simulations we created.

Figure 5. A Nonmanifold and a Manifold Surface Experiencing Gravity.

Newton’s Method

When we were running Pinkhall and Polthier’s algorithm on our surfaces, we noticed that that algorithm would not stop at a local saddle point such as the Schwarz P surface, but would instead run until there was only a thin strip of area left, which is indeed a minimum, but not a very useful one.

Therefore, we switched to Newton’s Method to solve our optimization problem.

We define the triangle surface area as an energy: let the three vertices of a triangle be \(\mathbf{q}_1\), \(\mathbf{q}_2\), \(\mathbf{q}_3\). Then \(E = \frac{1}{2} \|\mathbf{n}\| \), where \(\mathbf{n} = (\mathbf{q}_2-\mathbf{q}_1) \times (\mathbf{q}_3-\mathbf{q}_1)\) is the surface area normal. Then we derive its Jacobian and Hessian, and construct the Jacobian and Hessian for all faces in the mesh.

However, this optimization scheme still did not converge to the desired minimum, perhaps because our initialization is far from the solution. Additionally, one of our project mentors implemented the approach in C++ and, similarly, no results ensued. Later, we tried to add line search to Newton’s Method, but also no luck.

Although our algorithm still does not converge to some minimal surfaces which we know to exist, it has generated the following fun bugs.

Figure 6.
Figure 7.
Figure 8.

Future Work

In the previous blog post, we discussed studying the physical properties of nonmanifold TPMS. Over the past week, we used the Vellum Solver in Houdini and explored some of the differences in physical properties between manifold and nonmanifold surfaces. However, this is just the beginning — we can continue to expand our work in that direction.

Additional goals may include writing a script to generate many possible initial conditions, then converting the initial conditions into minimal surfaces, either by using the Pinkall and Polthier algorithm, or by implementing another minimal-surface-generating algorithm. 

More work can be done on enumerating all of the possible nonmanifold structures that the Pinkall and Polthier algorithm generates. The researchers can, then, categorize the structures based on their geometric or physical properties. As mentioned last week, this is still an open problem


We would like to thank our mentors Etienne Vouga, Nicholas Sharp, and Josh Vekhter for their patient guidance and the numerous hours they spent helping us debug our Matlab code, even when the answers were not immediately obvious to any of us. A special thanks goes to Stephanie Wang, who shared her Houdini expertise with us and, thus, made a lot of our visualizations possible. We would also like to thank our Teaching Assistant Erik Amezquita

SGI research projects

Robust computation of the Hausdorff distance between triangle meshes

Authors: Bryce Van Ross, Talant Talipov, Deniz Ozbay

The SGI project titled Robust computation of the Hausdorff distance between triangle meshes originally was planned for a 2 week duration, and due to common interest in continuing, was extended to 3 weeks. This research was led by mentor Dr. Leonardo Sacht of the Department of Mathematics of the Federal University of Santa Catarina (UFSC), Brazil. Accompanying support was TA Erik Amezquita, and the project team consisted of SGI fellows, math (under)graduate students, Deniz Ozbay, Talant Talipov, and Bryce Van Ross. Here is an introduction to the idea of our project. The following is a summary of our research and our experiences, regarding computation of the Hausdorff distance. 


Given two triangle meshes A, B in , the following are defined:

1-sided Hausdorff distance h: 

$$h(A, B) = \max\{d(x, B) : x\in A\} = \max\{\min\{\|x-y\|: x\in A\}: y\in B\}$$

where d is the Euclidean distance and |x-y| is the Euclidean norm. Note that h, in general, is not symmetric. In this sense, h differs from our intuitive sense of distance, being bidirectional. So, h(B, A) can potentially be a smaller (or larger) distance in contrast to h(A, B). This motivates the need for an absolute Hausdorff distance.

2-sided Hausdorff distance H: 

$$H(A,B) = \max\{h(A, B), h(B, A)\}$$

By definition, H is symmetric. Again, by definition, H depends on h. Thus, to yield accurate distance values for H, we must be confident in our implementation and precision of computing h.

For more mathematical explanation of the Hausdorff distance, please refer to this Wikipedia documentation and this YouTube video.


Objects are geometrically complex and so too can their measurements be. There are different ways to compare meshes to each other via a range of geometry processing techniques and geometric properties. Distance is often a common metric of mesh comparison, but the conventional notion of distance is at times limited in scope. See Figures 1 and 2 below.

Figure 1: This distance is limited to the red vertices, ignoring other points of the triangles.
Figure 2: This distance ignores the spatial positions of the triangles. So, the distance is skewed to the points of the triangles, and not the contribution of the space between the triangles.

Figures from the Hausdorff distance between convex polygons. 

Our research focuses on robustly (efficiently, for all possible cases) computing the Hausdorff distance h for pairs of triangular meshes of objects. The Hausdorff distance h is fundamentally a maximum distance among a set of desirable distances between 2 meshes. These desirable distances are minimum distances of all possible vectors resulting from points taken from the first mesh to the second mesh. 

Why is h significant? If h tends to 0, then this implies that our meshes, and the objects they represent, are very similar. Strong similarity indicates minimal change(s) from one mesh to the other, with possible dynamics being a slight deformation, rotation, translation, compression, stretch, etc. However, if h >> 0, then this implies that our meshes, and the objects they represent, are dissimilar. Weak similarity indicates significant change(s) from one mesh to the other, associated with the earlier dynamics. Intuitively, h depends on the strength of ideal correspondence from triangle to triangle, between the meshes. In summary, h serves as a means of calculating the similarity between triangle meshes by maximally separating the meshes according to the collection of all minimum pointwise-mesh distances. 


The Hausdorff distance can be used for a variety of geometry processing purposes. Primary domains include computer vision, computer graphics, digital fabrication, 3D-printing, and modeling, among others. 

Consider computer vision, an area vital to image processing, having a multitude of technological applications in our daily lives. It is often desirable to identify a best-candidate target relative to some initial mesh template. In reference to the set of points within the template, the Hausdorff distance can be computed for each potential target. The target with the minimum Hausdorff distance would qualify as being the best fit, ideally being a close approximation to the template object. The Hausdorff distance plays a similar role relative to other areas of computer science, engineering, animation, etc. See Figure 3 and 4, below, for a general representation of h.

Figure 3: Hausdorff distance h corresponds to the dotted lined distance of the left image. In the right image, h is found in the black shaded region of the green triangle via the Branch and Bound Method. This image is found in Figure 1 of the initial reading provided by Dr. Leonardo Sacht.
Figure 4: Hausdorff distance h corresponds to the solid lined distance of the left image. This distance is from the furthest “leftmost” point of the first mesh (armadillo) to the closest “leftmost” point of the second mesh. This image is found in Figure 5 of the initial reading provided by Dr. Leonardo Sacht.

Branch and Bound Method

Our goal was to implement the branch-and-bound method for calculating H. The main idea is to calculate the individual upper bounds of Hausdorff distance for each triangle meshes of mesh A and the common lower bound for the whole mesh A. If the upper bound of some triangle is greater than the general lower bound, then this face is discarded and the remaining ones are subdivided. So, we have these 3 main steps:

1. Calculating the lower bound

We define the lower bound as the minimum of the distances from all the vertices of mesh A to mesh B. Firstly, we choose the vertex P on mesh A. Secondly, we compute the distances from point P  to all the faces of mesh B. The actual distance from point P to mesh B is the minimum of the distances that were calculated the step before. For more theoretical details you should check this blog post:

The implementation of this part:

Calculating the distance from the point P to each triangle mesh T of mesh B is a very complicated process. So, I would like not to show the code and only describe it. The main features that should be considered during this computation is the position of point P relatively to the triangle T. For example, if the projection of point P on the triangle’s plane lies inside the triangle, then the distance from point P to triangle is just the length of the corresponding normal vector. In the other cases it could be the distance to the closest edge or vertex of triangle T.

During testing this code our team faced the problem: the calculating point-triangle distance takes too much time. So, we created the bounding-sphere idea. Instead of computing the point-triangle distance we decided to compute point-sphere distance, which is very simple. But what sphere should we pick? The first idea that sprung to our minds was the sphere that is generated by a circumscribed circle of the triangle T. But the computation of its center is also complicated. That is why we chose the center of mass M as the center of the sphere and maximal distance from M to each vertex of triangle T. So, if the distance from the point P to this sphere is greater than the actual minimum, then the distance to the corresponding triangle is exactly not the minimum. This trick made our code work approximately 30% faster. This is the realization:

2. Calculating the upper bounds

Overall, the upper bound is derived by the distances between the vertices and  triangle inequality. For more theoretical details you should check this blog post:

During testing the code from this page on big meshes our team faced the memory problem. On the grounds that we tried to store the lengths that we already computed, it took too much memory. That is why we decided just compute these length one more time, even though it takes a little bit more time (it is not crucial):

3. Discarding and subdividing

The faces are subdivided in the following way:  we add the midpoints and triangles that are generated by the previous vertices and these new points. In the end, we have 4 new faces instead of the old one. For more theoretical details you should check this blog post:


Below are our results for two simple meshes, first one being a sphere mesh and the second one being the simple mesh found in the blog post linked under the “Discarding and subdividing” section.

Figure 5: Results for a sphere mesh with different tolerance levels
Figure 6: Results for a smaller, simple mesh with different tolerance levels


The intuition behind how to determine the Hausdorff distance is relatively simple, however the implementation of computing this distance isn’t trivial. Among the 3 tasks of this Hausdorff distance algorithm (finding the lower bound, finding the upper bound, and finding a triangular subdivision routine), the latter two tasks were dependent on finding the lower bound. We at first thought that the subdivision routine would be the most complicated process, and the lower bound would be the easiest. We were very wrong: the lower bound was actually the most challenging aspect of our code. Finding vertex-vertex distances was the easiest aspect of this computation. Given that in MATLAB triangular meshes are represented as faces of vertex points, it is difficult to identify specific non-vertex points of some triangle. To account for this, we at first used computations dependent on finding a series of normals amongst points. This yielded a functional, albeit slow, algorithm. Upon comparison, the lower bounds computation was the cause of this latency and needed to be improved. At this point, it was a matter of finding a different implementation. It was fun brainstorming with each other about possible ways to do this. It was more fun to consider counterexamples to our initial ideas, and to proceed from there. At a certain point, we incorporated geometric ideas (using the centroid of triangles) and topological ideas (using the closed balls) to find a subset of relevant faces relative to some vertex of the first mesh, instead of having to consider all possible faces of the secondary mesh. Bryce’s part was having to mathematically prove his argument, for it to be correct, but only to find out later it would be not feasible to implement (oh well). Currently, we are trying to resolve bugs, but throughout the entire process we learned a lot, had fun working with each other, and gained a stronger understanding of techniques used within geometry processing.

In conclusion, it was really fascinating to see the connection between the theoretical ideas and the implementation of an algorithm, especially how a theoretically simple algorithm can be very difficult to implement. We were able to learn more about the topological and geometric ideas behind the problem as well as the coding part of the project. We look forward to finding more robust ways to code our algorithm, and learning more about the mathematics behind this seemingly simple measure of the geometric difference between two meshes.


Wrapping Up SGI 2021

After six weeks of intensive research and tutorials on applied geometry, we finally are ready to wrap up SGI 2021 and send our Fellows back to their home institutions. Directing this program has been one of the most rewarding experiences of my career, and it has been a pleasure seeing our students advance as scientists, mathematicians, and supportive community members.

SGI’s success is entirely thanks to a huge team of volunteers whose time and energy made the program possible.

Below, I acknowledge the many colleagues who have participated in the planning, leadership, and day-to-day activities of SGI 2021. The sheer length of this list is an inspiring reflection of the supportive community we enjoy in geometry research.

To start, SGI 2021 was only possible with the support of our generous sponsors, whose donations allowed us to offer a stipend to each SGI Fellow commensurate with a summer internship:

SGI was organized over the previous year with guidance from colleagues at a variety of institutions worldwide. A Steering Committee provided opinions and advice throughout the planning process:

Within MIT, several faculty and staff members contributed substantial effort to make the program a success. My team in the Geometric Data Processing (GDP) Group provided advice and volunteer support, from feedback on how to structure the program to hand-packing 72 boxes to ship to our Fellows and volunteers; GDP admin Mieke Moran organized payments and many other key aspects that made the program run smoothly. Our EECS Department Chair Prof. Asu Ozdaglar, AI+D Chair Prof. Antonio Torralba, and CSAIL Director Prof. Daniela Rus advocated for the program and provided support and encouragement as SGI obtained final approval within the MIT administration. CSAIL Assistant Director Carmen Finn provided critical last-minute help to make sure our Fellows were paid on time. Prof. Frédo Durand provided much-needed advice—and allowed me to vent—at several junctures.

SGI 2021 received far more applications than anticipated, and our final cadre of 34 Fellows and 29 additional tutorial week participants was painfully difficult to select. Our admissions committee carefully read all the applications:

The first week of SGI centered around five days of tutorials to get our Fellows up to speed on geometry processing research. Each tutorial day was organized by a different volunteer, who took charge of the content for the entire day and generated valuable course materials we can reuse in the future:

  • Day 1: Dr. Oded Stein (MIT), basic techniques in geometry processing
  • Day 2: Hsueh-Ti (Derek) Liu (University of Toronto) and Jiayi Eris Zhang (University of Toronto and Stanford), shape deformation
  • Day 3: Silvia Sellán (University of Toronto), shape representations
  • Day 4: Michal Edelstein (Technion) and Abhishek Sharma (École Polytechnique), shape correspondence
  • Day 5: Prof. Amir Vaxman (Utrecht University), directional fields

The remaining five weeks of SGI included a host of 1-3 week research projects, each involving experienced mentors working closely with SGI Fellow. Our full list of projects and project mentors is as follows:

  • Dr. Itzik Ben-Shabat: Self-supervised normal estimation using shape implicit neural representation (August 16-August 27)
  • Prof. Mikhail Bessmeltsev and Prof. Ed Chien: Singularity-Free Frame Field Design for Line Drawing Vectorization (July 26-August 6)
  • Dr. Tolga Birdal and Prof. Nina Miolane: Uncertainty Aware 3D Multi-Way Matching via Soft Functional Maps (July 26-August 6)
  • Prof. David Bommes and Dr. Pierre-Alexandre Beaufort: Quadrilateral and hexahedral mesh optimization with locally injective maps (July 26-August 6)
  • Prof. Marcel Campen: Improved 2D Higher-Order Meshing (July 26-July 30)
  • Prof. Keenan Crane: Surface Parameterization via Intrinsic Triangulations (August 9-August 27)
  • Dr. Matheus Gadelha: Learning Classifiers of Parametric Implicit Functions (August 16-August 27)
  • Prof. Lin Gao and Jie Yang: Unsupervised partial symmetry detection for 3D models with geometric deep learning (August 16-August 27)
  • Christian Hafner and Prof. Bernd Bickel: Joints for Elastic Strips (August 9-August 13)
  • Yijiang Huang and Prof. Caitlin Mueller: Design optimization via shape morphing (August 16-August 27)
  • Dr. Xiangru Huang: Moving Object Detection from consecutive LiDAR Point Clouds (August 23-August 27)
  • Prof. Michael Kazhdan: Multigrid on meshes (July 26-July 30)
  • Prof. Paul Kry and Alexander Mercier-Aubin: Incompressible flow on meshes (July 26-August 6)
  • Prof. Kathryn Leonard: 2D shape complexity (July 26-July 30)
  • Prof. David Levin: Optimal Interlocking Parts via Implicit Shape Optimizations (July 26-August 6)
  • Angran Li, Kuanren Qian, and Prof. Yongjie Jessica Zhang: Geometric Modeling for Isogeometric Analysis with Engineering Applications (August 2-August 6)
  • David Palmer: Bayesian Rotation Synchronization (August 2-August 13); Planar-faced and other specialized hexahedral meshes (August 16-August 27)
  • Prof. Jorg Peters (The beauty of) Semi-structured splines (August 9-August 13)
  • Alison Pouplin and Dimitris Kalatzis: Learning projection of hierarchical data with a variational auto-encoder onto the Klein disk (July 26-August 6)
  • Prof. Leonardo Sacht: Robust computation of the Hausdorff distance between triangle meshes (August 9-August 20)
  • Prof. Yusuf Sahillioglu: Cut optimization for parameterization (August 2-August 13)
  • Josua Sassen and Prof. Martin Rumpf: Mesh simplification driven by shell elasticity (August 9-August 20)
  • Dr. Nicholas Sharp, Prof. Etienne Vouga, Josh Vekhter: Nonmanifold Periodic Minimal Surfaces (August 9-August 27)
  • Dr. Tal Shnitzer-Dery: Self-similarity loss for shape descriptor learning in correspondence problems (August 9-August 13)
  • Dr. Ankita Shukla and Prof. Pavan Turaga: Role of orthogonality in deep representation learning for eco-conservation applications (August 9-August 13)
  • Prof. Noah Snavely: Reconstructing the Moon and Earth in 3D from the World’s Photos (August 9-August 13)
  • Prof. Justin Solomon: Anisotropic Schrödinger Bridges (August 16-August 27)
  • Prof. Marco Tarini: Better Volume-encoded parametrizations (August 2-August 13)
  • Prof. Amir Vaxman: High-order directional field design (July 26-August 6)
  • Prof. Etienne Vouga: Differentiable Remeshing (July 26-August 6)
  • Dr. Stephanie Wang: Discrete Laplacian, area functional, and minimal surfaces (August 16-August 20)
  • Paul Zhang: Classifying hexahedral mesh singular vertices (July 26-August 6); Subdivision Surface Mesh Fitting (August 16-August 27)

An intrepid team of TAs helped our participants learn new topics, refined the tutorial activities, and supported the research projects:

Each week of SGI, we had multiple guest speakers drop by to share their research and experiences, and to give advice to the SGI Fellows:

  • Prof. Katia Bertoldi, Harvard: Multistable inflatable origami from deployable structures to robots
  • Prof. Michael Bronstein, Twitter/Imperial College London: Geometric Deep Learning: the Erlangen Programme of ML
  • Prof. Albert Chern, UCSD: Gauge Theory for Vector Field Design
  • Prof. Bruce Fischl, Harvard/MGH: Geometry and the Human Brain
  • Dr. Fernando de Goes, Pixar: Geometry Processing at Pixar
  • Prof. Rana Hanocka, University of Chicago: Deep Learning on Meshes
  • Prof. Alexandra Ion, Carnegie-Mellon University: Interactive Structures – Materials that can move, walk, compute
  • Prof. Chenfanfu Jiang, UCLA: Developments in smooth optimization contact
  • Prof. Theodore Kim, Yale University: Anti-Racist Graphics Research
  • Prof. Yaron Lipman, Weizmann Institute: Implicit Neural Representations
  • Prof. Mina Konaković Luković, MIT: Transforming Design and Fabrication with Computational Discovery
  • Prof. Lakshminarayanan Mahadevan, Harvard: Folding and cutting paper: art and science
  • Prof. Caitlin Mueller, MIT: Geometry of High-Performance Architecture
  • Prof. Matthias Niessner, TU Munich: Learning to Reconstruct 3D Scenes
  • Íñigo Quílez: Intro to SDFs and Examples
  • Dr. Elissa Ross, Metafold: Periodic geometry: from art to math and back again
  • Dr. Ryan Schmidt, Epic Games and Gradientspace: Geometry Processing in Practice
  • Prof. Tamar Shinar, UC Riverside: Partitioned solid-fluid coupling
  • Prof. Emily Whiting, Boston University: Mechanics-Based Design for Computational Fabrication

Last but not least, incoming MIT PhD student Leticia Mattos Da Silva organized a talk and panel discussion on the graduate school application process, including a Q&A with Silvia Sellán, Jiayi Eris Zhang, and Oded Stein.

The cast of thousands above is a testament to the dedication of the geometry research community to developing a diverse, energetic community of young researchers.

SGI 2021 comes to a close as quietly as it began, as our Fellows and mentors close one final Zoom call and return to their lives scattered over the globe. In the months and years to come, we look forward to keeping in touch as our Fellows become colleagues, collaborators, and leaders of future generations of researchers.

SGI research projects

Volume-encoded parameterization

by Alice Mehalek, Marcus Vidaurri, and Zhecheng Wang


UV mapping or UV parameterization is the process of mapping a 3D surface to a 2D plane. A UV map assigns every point on the surface to a point on the plane, so that a 2D texture can be applied to a 3D object. A simple example is the representation of the Earth as a 2D map. Because there is no perfect way to project the surface of a globe onto a plane, many different map projections have been made, each with different types and degrees of distortion. Some of these maps use “cuts” or discontinuities to fragment the Earth’s surface.

In computer graphics, the typical method of UV mapping is by representing a 3D object as a triangle mesh and explicitly assigning each vertex on the mesh to a point on the UV plane. This method is slow and must be repeated often, as the same texture map can’t be used by different meshes of the same object.

For any kind of parametrization or UV mapping, a good UV map must be injective and should be conformal (preserving angles) while having few cuts. Ideally it should also be isometric (preserve relative areas). In general, however, more cuts are needed to achieve less distortion.

Volume-encoded parameterization

Our research mentor, Marco Tarini, developed the method of Volume-encoded UV mapping. In this method, a surface is represented by parametric functions and each point on the surface is mapped to a UV position as a function of its 3D position. This is done by enclosing the surface or portion of the surface in a unit cube, or voxel, and assigning UV coordinates to each of the cube’s eight vertices. All other points on the surface can be mapped by trilinear interpolation. Volume-encoded parametrization has the advantage of only needing to store eight sets of UV coordinates per voxel, instead of unique locations of many mesh vertices, and can be applied to different types of surface representations, not just meshes. 

We spent the first week of our research project learning about volume-encoded parametrization by exploring the 2D equivalent: mapping 2D curves to a one-dimensional line, u. Given a curve enclosed within a unit square, our task was to find the u-value for each corner of the square that optimized injectivity, conformality, and orthogonality. We did this using a Least Squares method to solve a linear system consisting of a set of constraints applied to points on the surface. All other points on the curve could then be mapped to u by bilinear interpolation.

A quarter circle (the red curve on the xy plane of the 3D plot) is mapped to a line, u, which is also represented as the height on the z axis in the 3D plot. The surface in the 3D plot is obtained by bilinear interpolation of the height of each corner of the square, and the red curve along the surface represents the path of the quarter circle mapped to 1D.
Each point on the path has a unique height, indicating an injective mapping.

An injective mapping was not possible for portions of a circle greater than a semicircle.
Each point on the path of u does not have a unique height, indicating loss of injectivity. There is also distortion in the middle where the slope is variable.

In Week 2 of the project, we moved on to 3D and created UV maps for portions of a sphere, using a similar method to the 2D case. Some of the questions we wanted to answer were: For what types of surfaces is volume-encoded parametrization possible? At what level of shape complexity is it necessary to introduce cuts, or split up the surface into smaller voxels? In the 2D case, we found that injectivity could be preserved when mapping curves less than half a circle, but there were overlaps for curves greater than a semicircle.

One dimension up, when we go into the 3D space, we found that uniform sampling on the quarter sphere was challenging. Sampling uniformly in the 2D parametric space will result in a distorted sampling that becomes denser when getting closer to the north pole of the sphere. We tried three different approaches: rewriting the parametric equations, sampling in a unit-sphere, and sampling in the 3D target space and then projecting the sample points back to the surface. Unfortunately, all three methods only worked with a certain parametric sphere equation.

When mapping a portion of a sphere to 2D, the grid allows us to see where the mapping is distorted. This mapping is most accurate at the north pole, while areas and angles both become distorted toward the equator.

In conclusion, over the two weeks, we designed and tested a volume-encoded least-squares conformal mapping in both 2D and 3D. In the future, we plan to rewrite the code in C++ and run more experiments.

Geometry processing Math

When Algorithms Go Wrong

By Natasha Diederen, Sahra Yusuf, and Olga Guțan

I. Introduction

An algorithm is a finite sequence of well-defined, computer-implementable instructions, typically designed to solve a class of specific problems or to perform a computation. Similarly to how we use our brains to (sometimes wrongly) categorize things in order to make sense of the world, a computer needs an algorithm to make sense of what the user is asking it to do. Since an algorithm is a way of communication between two vastly different entities — a human and a machine — some information gets lost along the way. The process is intellectually intriguing to witness, however, problems can arise when we use algorithms to make decisions that influence the lives of humans and other self-aware living beings. 

Algorithms are indispensable for efficiently processing data; therefore they will continue to be part of our programming ecosystem. However, we can (and must) keep some space in our minds for the additional nuances that reality contains. While a programmer may not be able to fully convey these nuances to a computer, we must be aware that no algorithm’s output is final, all-encompassing, and universally true. Furthermore, we must be mindful of the conclusions we draw from the output of algorithms. 

II. Algorithms Gone Wrong

Broadly speaking, potential pitfalls for algorithms are manifold. The issues stem from the nature of the algorithms — a communication tool between a human entity and a nonhuman entity (a computer).  These problems morph into different real-life issues based on what types of algorithms we use and what we use them for. 

Even when playing with algorithms intended to solve toy problems that we already have the answers to, we can notice errors. However, in real life, data is (even) more messy and consequences are far larger. In 2018, Elaine Herzberg was hit and killed by a self-driving taxi. She was jaywalking in the dark with a bicycle by her side and the car alternated between classifying her as a person and a bicycle, thus miscalculating her trajectory and not classifying Elaine as an obstruction. In this case, the safety driver was distracted by a television program, and thus not entirely blameless. However, this serves as an example of how our blind trust in the reliability of algorithms can often result in complacency, with far reaching consequences. 

A further example of error in algorithms is adversarial attacks on neural networks. This is when researchers (or hackers) feed a neural network a modified image, where changes made to the image are imperceptible to the human eye. Nevertheless, this can result in incorrect classification by neural networks, with high confidence to boot. 

Figure 1. A demonstration of fast adversarial example generation applied to GoogLeNet (Szegedy et al., 2014a) on ImageNet. Source.

In Figure 1 we see how by adding an imperceptibly small vector whose elements are equal to the sign of the elements of the gradient of the cost function with respect to the input, we can change GoogLeNet’s classification of the image. Here the \(\epsilon\) of 0.007 corresponds to the magnitude of the smallest bit of an 8 bit image encoding after GoogLeNet’s conversion to real numbers.

Although researchers are working on making neural networks more robust to these sorts of attacks, susceptibility to adversarial attacks seems to be an inherent weakness of deep neural networks. This has serious security implications as computer vision increasingly relies on these models for facial recognition, autonomous vehicles, speech recognition and much more. 

III. Geometry Processing with Imperfect Data

In geometry processing, there is often a need for refinement of geometric data, especially when the data is taken from “real” contexts. Examples of imperfections in 3D models include gaps, self-intersections, and non-manifold geometry, or geometry which cannot exist in reality (e.g. edges included in more than two faces and disconnected vertices).

One common source of imperfect, “real-life” data is 3D object scanning. The resulting models typically include gaps and large self-intersections as a result of incomplete scanning or other error arising from the scanning method. Despite these significant problems, scanned data is still invaluable for rapidly generating assets for use in game development and other applications. During our time at the Summer Geometry Institute, Dr. Matthias Nießner spoke about 3D scene reconstruction with scanned data. His work demonstrated a method of improving the overall quality of reconstructed scenes using bundle adjustment and other techniques. He also worked on solving the problem of incomplete scanned objects using machine learning. Previously, we have written about possible mistakes arising from the weaknesses of machine learning, but Dr. Nießner’s work demonstrates that machine learning is a valuable tool for refining data and eliminating mistakes as well.

Although error in geometry processing is not as frequently discussed, the implications are just as important as those of mistakes arising from the use of machine learning. This is primarily due to the fact that machine learning and geometry processing are not isolated fields or tools, but are often used together, especially in the sorts of situations we described earlier. By researching and developing new methods of data refinement, we can improve the usability of natural data and increase, for example, the safety of systems which rely on visual data.

IV. Human Error and Bias

The errors we have discussed so far exclude human error and bias, which can aggravate existing inequalities in society. For example, a Fellow one of us worked with during SGI mentioned how he worked on a project which used face tracking to animate digital characters. However, the state-of-the-art trackers only worked on him (a white male) and could not track eye or mouth movement for those in his team of black descent. Additionally, as we heard from Theodore Kim in another brilliant SGI talk, animation is focused on designing white skin and non-type 4 hair, further highlighting the systemic racism present in society.

Moreover, the fact that 94.8% of Automatic Gender Recognition (AGR) papers recognize gender as binary has huge implications for the misgendering of trans people. This could lead to issues with AGR based access control for things like public bathrooms, where transgender people may be excluded from both male and female facilities. The combination of machine and human error is especially dangerous, and it is important to recognize this, so that we can mitigate against the worst harm.

V. Conclusion

Algorithms have become a fundamental part of human existence, however our blind faith that algorithms are always (1) correct and (2) unbiased is deeply flawed. The world is a complicated place, with far too much data for computers to handle, placing a strong reliance on simplification algorithms. 

As we have seen from the examples above, these algorithms are far from perfect and can sometimes erase or distort important data. This does not mean we should stop using algorithms entirely. What this does, however, mean is that we must employ a hearty dose of critical thinking and skepticism when analyzing results outputted by an algorithm. We must be especially careful if we use these results for making decisions that would influence the lives of other humans. 

SGI research projects

Intrinsic Parameterization: Weeks 1-2

By Joana Portmann, Tal Rastopchin, and Sahra Yusuf.
Mentored by Professor Keenan Crane.

Intrinsic parameterization

During these last two weeks, we explored intrinsic encoding of triangle meshes. As an introduction to this new topic, we wrote a very simple algorithm that lays out a triangle mesh flat. We then improved this algorithm via line search over a week. In connection with this, we looked into terms like ‘angle defects,’ ‘cotangent weights,’ and the ‘cotangent Laplacian’ in preparation for more current research during the week.

From intrinsic to extrinsic parameterization

To get a short introduction into intrinsic parameterization and its applications, I quote some sentences from Keenan’s course. If you’re interested in the subject I can recommend the course notes “Geometry Processing with Intrinsic Triangulations.” “As triangle meshes are a basic representation for 3D geometry, playing the same central role as pixel arrays in image processing, even seemingly small shifts in the way we think about triangle meshes can have major consequences for a wide variety of applications. Instead of thinking about a triangle as a collection of vertex positions in \(\mathbb{R}^n\) from the intrinsic perspective, it’s a collection of edge lengths associated with edges.” 

Many properties of a surface such as the shortest path do only depend on local measurements such as angles and distances along the surfaces and do not depend on how the surface is embedded in space (e.g. vertex positions in \(\mathbb{R}^n\)), so an intrinsic representation of the mesh works fine. Intrinsic triangulations bring several deep ideas from topology and differential geometry into a discrete, computational setting. And the framework of intrinsic triangulations is particularly useful for improving the robustness of existing algorithms.

Laying out edge lengths in the plane

Our first task was to implement a simple algorithm that uses intrinsic edge lengths and a breadth-first search to flatten a triangle mesh onto the plane. The key idea driving this algorithm is that given just a triangle’s edge lengths we can use the law of cosines to compute its internal angles. Given a triangle in the plane we can use the internal angles to flatten out its neighbors into the plane. We will later use these angles to modify the edge lengths so that we “better” flatten the model. The algorithm works by choosing some root triangle and then performing a breadth-first traversal to flatten each of the adjacent triangles into the plane until we have visited every triangle in the mesh.

Breadth-first search pseudocode

Some initial geometry central pseudocode for this breadth first search looks something like

// Pick and flatten some starting triangle
Face rootTriangle = mesh->face(0);

// Calculate the root triangle’s flattened vertex positions

// Initialize a map encoding visited faces
FaceData<bool> isVisited(*mesh, false);

// Initialize a queue for the BFS
std::queue<Face> visited;

// Mark the root triangle as visited and pop it into the queue
isVisited[rootTriangle] = true;

// While our queue is not empty
    // Pop the current Face off the front of the queue
    Face currentFace = visited.front();

    // Visit the adjacent faces
    For (Face adjFace : currentFace.adjacentFaces()) {
        // If we have not already visited the face
        if (!visited[adjFace]
            // Calculate the triangle’s flattened vertex positions
            // And push it onto the queue

In order to make sure we lay down adjacent triangles with respect to the computed flattened plane coordinates of their parent triangle we need to know exactly how a child triangle connects to its parent triangle. Specifically, we need to know which edge is shared by the parent and child triangle and which point belongs to the child triangle but not the parent triangle. One way we could retrieve this information is by computing the set difference between the vertices belonging to the parent triangle and child triangle, all while carefully keeping track of vertex indices and edge orientation. This certainly works, however, it can be cumbersome to write a brute force combinatorial helper method for each unique mesh element traversal routine.

The halfedge mesh data structure

Professor Keenan Crane explained that a popular mesh data structure that allows a scientist to conveniently implement mesh traversal routines is that of the halfedge mesh. At its core a halfedge mesh encodes the connectivity information of a combinatorial surface by keeping track of a set of halfedges and the two connectivity functions known as twin and next. Here the set of halfedges are none other than the directed edges obtained from an oriented triangle mesh. The twin function takes a halfedge and brings it to its corresponding oppositely oriented twin halfedge. In this fashion, if we apply the twin function to some halfedge twice we get the same initial halfedge back. The next function takes a halfedge and brings it to the next halfedge in the current triangle. In this fashion if we take  the next function and apply it to a halfedge belonging to a triangle three times, we get the same initial halfedge back.

A diagram depicting the halfedge data structure connectivity relationships. Source.

Professor Keenan Crane has a well written introduction to the halfedge data structure in section 2.5 of his course notes on Discrete Differential Geometry.

It turns out that geometry central uses the halfedge mesh data structure and so we can rewrite the traversal of the adjacent faces loop to more easily retrieve our desired connectivity information. In the geometry central implementation, every mesh element (vertex, edge, face, etc.) contains a reference to a halfedge, and vice versa.

// Pop the current Face off the front of the queue
Face currentFace = visited.front();

// Get the face’s halfedge
Halfedge currentHalfedge = currentFace.halfedge();

// Visit the adjacent faces
    // Get the current adjacent face
    Face adjFace = currentHalfedge.twin().face();
    if (!isVisited[adjFace])
        // Retrieve our desired vertices
        Vertex a = currentHalfedge.vertex();
        Vertex b =;
        Vertex c = currentHalfedge.twin().next().next().vertex();

        // Calculate the triangle’s flattened vertex positions
        calculateFlatVertices(a, b, c);
        // And push it onto the queue
    // Iterate to the next halfedge
    currentHalfEdge =;
// Exit the loop when we reach our starting halfedge
while (currentHalfEdge != currentFace.halfedge());

Here’s a diagram illustrating the relationship between the currentHalfedge and vertices a, b, and c.

A diagram illustrating the connectivity relationship between the currentHalfedge and vertices a, b, and c. Note that cH abbreviates currentHalfedge.

Segfaults, debugging, and ghost faces

This all looks great right? Now we need to determine the specifics of calculating the flat vertices? Well, not exactly. When we were running a version of this code in which we attempted to visualize the resulting breadth-first traversal we encountered several random segfaults. When Sahra ran a debugger (shout out to GDB and LLDB 🥰) we learned that the culprit was the isVisited[adjFace]  function call on the line

if (!isVisited[adjFace])

We were very confused as to why we would be getting a segfault while trying to retrieve the value corresponding to the key adjFace contained in the map FaceData<bool> isVisited. Sahra inspected the contents of the adjFace object and observed that it had index 248 whereas the mesh we were testing the breadth-first search on only had 247 faces. Because C++ zero indexes its elements, this means we somehow retrieved a face with index out of range by 2! How could this have happened? How did we retrieve that face in the first place? Looking at the lines

// Get the current adjacent face
Face adjFace = currentHalfedge.twin().face();

we realized that we made an unsafe assumption about currentHalfedge. In particular, we assumed that it was not a boundary edge. What does the twin of a boundary halfedge that has no real twin look like? If the issue we were running to was that the currentHalfedge was a boundary halfedge, why didn’t we get a segfault on simply currentHalfedge.twin()? Doing some research, we found that the geometry central internals documentation explains that

“We can implicitly encode the twin() relationship by storing twinned halfedges adjacent to one another– that is, the twin of an even-numbered halfedge numbered he is he+1, and the twin of and odd-numbered halfedge is he-1.”

Geometry central internals documentation

Aha! This explains exactly why currentHalfedge.twin() still works on a boundary halfedge; behind the scenes, it is just adding or subtracting one to the halfedge’s index. Where did the face index come from? We’re still unsure, but we realized that the face currentHalfedge.twin().face() only makes sense (and hence can only be used as a key for the visited map) when currentHalfedge is not a boundary halfedge. Here is a diagram of the “ghost face” that we think the line Face adjFace = currentHalfedge.twin().face() was producing.

A diagram depicting how taking the face of the twin of a boundary halfedge produces a nonexistent face.
A diagram depicting how taking the face of the twin of a boundary halfedge produces a nonexistent face.

Changing the map access line in the if statement to

if (!currentHalfedge.edge().isBoundary() && !isVisited[adjFace])

resolved the segfaults and produced a working breadth-first traversal.

Conformal parameterization

Here is a picture illustrating applying the flattening algorithm to a model of a cat head.

A picture illustrating the application of the flattening algorithm to a model of a cat head.

You can see that there are many cracks and this is because the model of the cat head is not flat—in particular, it has vertices with nonzero angle defect. A vertex angle defect for a given vertex is equal to the difference between \(2 \pi\) and the sum of the corner angles containing that vertex. This is a good measure of how flat a vertex is because for a perfectly flat vertex, all angles surrounding it will sum to \(2 \pi\).

After laying out the edges on the plane \(z=0\), we began the necessary steps to compute a conformal flattening (an angle-preserving parameterization). In order to complete this step, we needed to find a set of new edge lengths which would both be related to the original ones by a scale factor and minimize the angle defects,

\(l_{ij} := \sqrt{ \phi_i \phi_j} l_{ij}^0, \quad \forall ij \in E\),

where where \(l_{ij}\) is the new intrinsic edge length, \(\phi_i, \phi_j\) are the scale factors at vertices \(i, j\), and \(l_{ij}^0\) is the initial edge length.

Discrete Yamabe flow

At this stage, we have a clear goal: to optimize the scale factors in order to scale the edge lengths and minimize the angle defects across the mesh (i.e. fully flatten the mesh). From here, we use the discrete Yamabe flow to meet both of these requirements. Before implementing this algorithm, we began by substituting the scale factors with their logarithms

\(l_{ij} = e^{(u_i + u_j)/2} l_{ij}^0\),

where \(l_{ij}\) is the new intrinsic edge length, \(u_i, u_j\) are the scale factors at vertices \(i, j\), and \(l_{ij}^0\) is the initial edge length. This ensures the new intrinsic edge lengths are always positive and that the optimization is convex.

To implement the algorithm, we followed this procedure:

1. Calculate the initial edge lengths

2. While all angle defects are below a certain threshold epsilon:

  • Compute the scaled edge lengths
  • Compute the current angle defects based on the new interior angles based on the scaled edge lengths
  • Update the scale factors using the step and the angle defects:

\(u_i \leftarrow u_i – h \Omega _i\),

where \(u_i\) is the scale factor of the \(i\)th vertex, \(h\) is the gradient descent step size, and \(\Omega_i\) is the intrinsic angle defect at the \(i\)th vertex.

After running this algorithm and displaying the result, we found that we were able to obtain a perfect conformal flattening of the input mesh (there were no cracks!). There was one issue, however: we needed to manually choose a step size that would work well for the chosen epsilon value. Our next step was to extend our current algorithm by implementing a backtracking line search which would change the step size based on the energy.


Here are two videos demonstrating the Yamabe flow algorithm. The first video illustrates how each iteration of the flow improves the flattened mesh and the second video illustrates how that translates into UV coordinates for texture mapping. We are really happy with how these turned out!

A video visualizing intermediate 2D parameterizations produced by the Yamabe flow.
A video visualizing the intermediate UV coordinates on the cat head mesh produced by the Yamabe flow

Line search

To implement this, we added a sub-loop to our existing Yamabe flow procedure which repeatedly test smaller step sizes until one is found which will decrease the energy, e.g., backtracking line search. A good resource on this topic is Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004. After resolving a bug in which our step size would stall at very small values with no progress, we were successful in implementing this line search. Now, a large step size could be given without missing the minima.

Newton’s method

Next, we worked on using Newton’s method to improve the descent direction by changing the gradient to more easily reach the minimum. To complete this, we needed to calculate the Hessian — in this case, the Hessian of the discrete conformal energy was the cotan-Laplace matrix \(L\). This matrix had square dimensions (the number of rows and the number of columns was equal to the number of interior vertices) and has off-diagonal entries:

\(L_{ij} = -\frac{1}{2} (\cot \theta_k ^{ij} + \cot \theta _k ^{ji})\)

for each edge \(ij\), as well as diagonal entries

\(L_{ii} = – \sum _{ij} L_{ij}\)

for each edge \(ij\) incident to the \(i\)th vertex.

The newton’s descent algorithm is as follows:

  1. First, build the cotan-Laplace matrix L based on the above definitions
  2. Determine the descent direction \(\dot{u} \in \mathbb{R}^{|V_{\text{int}}|}\), by solving the system \(L \dot{u} = – \Omega\) with \(\Omega \in \mathbb{R}^{|V_{\text{int}}|}\).as the vector containing all of the angle defects at interior vertices.
  3. Run line search again, but with \(\dot{u}\) replacing \(– \Omega\) as the search direction.

This method yielded a completely flat mesh with no cracks. Newton’s method was also significantly faster: on one of our machines, Newton’s method took 3 ms to compute a crack-free parameterization for the cat head model while the original Yamabe flow implementation took 58 ms.

SGI research projects

Branch-and-bound method for calculating Hausdorff distance

This week I worked on the “Robust computation of the Hausdorff distance between triangle meshes” project with our mentor Dr. Leonardo Sacht, TA Erik Amezquita and my team Talant Talipov and Bryce Van Ross.

We started our project by doing some initial reading about the topic. Hausdorff distance from triangle meshes A to B is defined as $$h(A,B) = \max_{p \in A}d(p,B)$$ where d is the Euclidean norm. Finding the Hausdorff distance between two triangle meshes is one way of comparing these meshes. We note that the Hausdorff distance from A to B might be different from the Hausdorff distance from B to A, as you can see in the figure below, so it is important to distinguish which one is being computed.

Figure 1: Hausdorff distance from Mesh A to Mesh B
Figure 2: Hausdorff distance from Mesh B to Mesh A

Finally we define $$H(a,b) = \max{h(A,B), h(B,A)}$$ and use this value that is symmetric when comparing triangle meshes A and B.

Our first task was to implement Branch and Bound methods for calculating this distance. Suppose we want to calculate the Hausdorff distance from mesh A to mesh B. There were three main steps in the algorithm to do this: calculation of the upper and lower bound for the Hausdorff distance, and discarding and subdividing triangles according to the values of the upper and lower bound. 

The upper bound for each triangle in A is calculated by taking the maximum of the distances from the given triangle to every vertex in B.  The lower bound is calculated over A by taking the minimum of the distances from each vertex p in A to the triangle in B that is closest to p. If a triangle in A has an upper bound that is less than the lower bound, the triangle is discarded. After the discarding process, the remaining triangles are subdivided into four triangles and the process is repeated with recalculation of the upper bounds and the lower bound. The algorithm is run until the values for the upper and lower bound are within some ε of each other. Ultimately, we get a triangle region that contains the point that realizes the Hausdorff distance from A to B.

To implement this method, my teammates tackled the upper and lower bound codes while I wrote an algorithm for the discarding and subdividing process:

We ran this algorithm with testing values u = [1;2;3;4;5] and l = 3 and got these results:

Figure 3: Initial mesh
Figure 4: After discarding and subdividing

As expected, the two triangles that had the upper bound of 1 and 2 were discarded and the rest were subdivided.

 The lower bound algorithm turned out to be more challenging than we anticipated, and we worked together to come up with a robust method. Currently, we are working on finishing up this part of the code, so that we can run and test the whole algorithm. After this, we are looking forward to extending this project in different directions, whether it is a more theoretical analysis of the topic or working on improving our algorithm!


Week 5 Talks Summary

Authors: Bryce Van Ross, Deniz Ozbay, Talant Talipov

The following is a summary of three Week 5 talks, presented on behalf of students in the MIT Summer Geometry Institute. These summaries were made by members of the team of the Hausdorff Distance group. Topics include: Partitioned Solid-Fluid Coupling, Graduate School Panel and Q&A, and Developments in Smooth Optimization Contact.

Partitioned Solid-Fluid Coupling

This talk was presented by Dr. Tamar Shinar at University of California, Riverside. I liked this talk for several reasons. From an introductory physics background, we are aware that forces act on multiple objects. Depending on the type of object, certain forces behave differently with respect to the object and we use certain equations to model this interaction. However, as the system of objects becomes mixed, our use of forces becomes more complex. So, the modeling efforts must adapt accordingly. This lecture focused on how to model fluid-structure interaction (FSI), commonly referred to as solid-fluid coupling, via numerical algorithms, and discussed related applications. Examples of applications typically include fluid dynamics or animation. As suggested by Dr. Shinar, more intuitive, concrete examples are modeling the flight of a bird in the air (Figure 1), or the movement of a straight hair flip relative to a medium of water (Figure 2).

Figure 1: A bird (solid) flying in the air (fluid) adjacent to water droplets (another fluid). Image from
Figure 2: A woman (solid) flipping her hair in a pool of water (fluid). Image from

Solid-fluid coupling models the dynamics of fluids and solids. Specifically, these dynamics are computed simultaneously, which can be problematic. The challenge in this field is to make it less so, while ensuring a reflective representation of the physical forces at work. In 1- way coupling, the solid (or fluid) affects the fluid (solid), but not vice versa. This is not the case for 2-way coupling, where both objects affect each other, imposing more constraints upon the system. To account for each object’s phenomena, there exists solid and fluid handlers to model these phenomena, respectively. These solvers utilize a variety of boundary conditions, familiar to us from past math/physics courses, such as the Dirichlet, Navier-Stokes, and Neumann equations. Depending on the object and conditions, desirable properties are optimized. For solids, these would be mass and acceleration, whereas for fluids, these are position and velocity. The solution lies in the mixed-solver, and which numerical methods best model these behaviors.

There are multiple options to go about optimizing coupled-interaction. Partitioning the solver is a popular approach. In this approach, fixed point iteration and Newton’s method are strategies that solve the nonlinear system of equations, linearizing the 2-way coupling in terms of the solid and fluid. The partitioned methods tend to be more modular and easier to implement. In contrast, monolithic methods are more stable, have better performance, and handle closed regions better. In comparison to each other, their strengths are the other’s weaknesses. Common risks are that iterations can be computationally expensive and/or that they won’t converge at all. This is remedied by incorporating a reduced model interface (RMI), which serves as an intermediary lightweight solver between the solid and fluid. In effect, the RMI helps estimate tricky Jacobians in a least-weight sense. Dr. Shinar’s research specifically compares these numerical schemes in reference to an under-relaxation variable, which optimizes solution oscillations and aids in stability. The smaller this value, the slower and less efficient the computational time, in tradeoff for increased stability. The takeaways here are that the partition approach in isolation fails in terms of convergence, and the Partition-Coupling and RMI under low under-relaxation will yield optimal convergence in the shortest time. There still remain certain bounded regions, like (incompressible) Dirichlet regions, which interfere with the dynamics of known solid-fluid coupling techniques, but this has yet to be explored. For more information regarding Dr. Shinar’s research, please refer here.

An extended partitioned method for conservative solid-fluid coupling (Video from YouTube)

Developments in Smooth Optimization Contact

Professor Chenfanfu Jiang’s talk on “Developments in Smooth Optimization Contact” was really interesting for many reasons. Essentially the problem they were trying to solve was “How to model contact as a smooth problem”. This helps us accurately simulate solids which allows us to create realistic animations and benefits industry fields that use these simulations. For example in the fashion industry, these simulations provide feedback to the design of the garment.

Firstly, Dr. Jiang presented different ways to approach the problem and the advantages and shortcomings of these approaches. Then, he gave a formulation of how they came up with a solution using these, which was really interesting as it combined various approaches to come up with a better working algorithm. The research problem was also simple to understand but the tools used to solve the problem were very advanced, which was also very intriguing. 

To solve this problem three main tools were used: nonlinear optimization based integrator which simulated large deformation with large steps, variational frictional contact and its extension to higher dimensions. To start with, the solution is an optimization based time integrator for electrodynamics. It is based on incremented time steps, which creates an entirely new problem at each time step. The ODE consists of an inertia term and elasticity term – which makes the problem harder as the elasticity term can be highly nonlinear and non convex. This simulates deformations based on squashing and stretching. For a more physically accurate solution, contact is added as a constraint. This constraint in the form of an inequality, makes the differential equation even harder to solve. To overcome this challenge, different approaches were developed based on previous ways to solve the problem.

The difficulty arises from the fact that the formulated constraints and ODEs that are non-smooth and non-linear, for example, when the distance between objects is zero, we have a non-smooth constraint. To overcome this, the problem is turned into a smooth and unconstrained problem by using the incremental potential contact method which guarantees intersection and parameter tuning free simulations which are stable and accurate, and the barrier method, which turns the problem into an unconstrained problem where Newton iteration can be applied.

Dr. Jiang also discussed some methods that can improve their algorithm such as clamping the barrier function, and the shortcomings such as the difficulty of implementing friction in the optimization problem. Again, the challenges arise from the function becoming non-smooth, which he mentioned can be solved by approximately defining the variational friction as a smooth function.

The most intriguing part of the talk was how they were able to build on existing approaches to turn the problem into a smooth one. Although the techniques they used were very advanced, the idea behind the solution was simple to understand. The simulations (especially the twisting of the cloth) Dr. Jiang showed were also really interesting; it was fascinating to see the difference between two simulations of the same object and how they become visibly more accurate as the constraints change.

Figure 3: Twisting of a cloth (image from
Figure 4: Rod twist (image from

Graduate School Panel

Last week I attended a talk about applications to PhD programs that was organized by Leticia Mattos Da Silva, Silvia Sellán, Jiayi Eris Zhang and Oded Stein.

One of the most exciting deals for every student is application for graduate programs. Most students are always confused with this procedure: they do not know what to write in a motivation letter and how to contact the prospective referee. All these questions were answered in this lecture. This lecture was filled with a lot of useful information about all the steps of applications to universities: from motivation letters to the choice of prospective scientific advisor. It was really interesting on the grounds that the lecturers experienced this procedure firsthand. Moreover, the speakers showed their own essays as successful examples.

As I know from this talk, the motivation letter has a specific structure. The applicant has to mention his previous research experience, why he wants to study exactly on this program, and why the university has to choose this person among the others. Furthermore, there should not be any information that is not relevant to the program: if you are applying to a PhD program in Mathematics, then you do not have to say about your background in knitting. The best moments were when the lecturer highlighted the most frequent mistakes of applicants and I realized that I did it in my motivation letter.

The recommendation letters are one of the most significant parts of the application. That is why the applicant should connect with his referee in advance. The best choice would be the person who did some research with you. In addition, you should have more referees in the case that the main ones reject.

It was a very valuable experience to listen to this lecture!

SGI research projects

Bayesian Rotation Synchronization

By Adrish Dey, Dorothy Najjuma Kamya, and Lily Kimble

Note: Although this is an ongoing work, this report documents our progress between the official 2 weeks of the project. (August 2, 2021 – August 13, 2021)

The past 2 weeks at SGI, we have been working with David Palmer on investigating a novel Bayesian approach towards the angular synchronization problem. This blog post is written to document our work and share a sneak peek into our research.


Consider a set of unknown absolute orientations \(\{q_1, q_2, \ldots, q_n\}\) with respect to some fixed basis. The problem of angular synchronization deals with the accurate estimation of these orientations from noisy observations of their relative offsets \(O_{i, j}\), up to a constant additive phase. We are particularly interested in estimating these “true” orientations in the presence of many outlier measurements. Our interest in this topic stems from the fact that the angular synchronization problem arises in various avenues of science, including reconstruction problems in computer vision, ranking problems, sensor network localization, structural biology, etc. In our work, we study this problem from a Bayesian perspective, by modelling the observed data as a mixture between noisy observations and outliers. We also investigate the problem of continuous label switching, a global ambiguity that arises from the lack of knowledge about the basis of the absolute orientations \(q_i\). Finally, we experiment on a novel Riemannian gradient descent method for alleviating this continuous label switching problem and provide our observations herein.

Brief Primer on Bayesian Inference

Before going deeper, we’ll briefly discuss Bayesian inference. At the heart of Bayesian inference lies the celebrated Bayes’ rule (read \(a|b\) as “a given b”):

\[\underbrace{P(q | O)}_{\textrm{posterior}} = \frac{\overbrace{P(O|q)}^{\textrm{likelihood}} \cdot \overbrace{P(q)}^{\textrm{prior}}}{\underbrace{\int\limits_q P(O|q)\cdot P(q)}}_{\textrm{evidence}}\]

In our problem, \(q\) and \(O\) denote the true orientations that we are estimating and the noisy observations with outliers respectively. We are interested in finding the posterior distribution (or at least samples from it) over the ground truth \(q\) given the noisy observations \(O\).

The denominator (i.e., the evidence or partition function) is an integral over all \(q\)s. Exactly evaluating this integral is often intractable if \(q\) lies on some continuous manifold, as in our problem. This makes drawing samples from the posterior becomes hard.

One way to avoid computing the partition function is a sampling method called Markov Chain Monte Carlo (MCMC). Intuitively, the posterior is approximated by a markov chain whose transitions can be computed using a simpler distribution called the proposal distribution. Successive samples are then accepted or rejected based on an acceptance probability designed to ensure convergence to the posterior distribution in the limit of infinite samples. Simply put, after enough samples are drawn using MCMC, they will look like the samples from the posterior\(P(q|O) \propto P(O|q) \cdot P(q)\) without requiring us to calculate the intractable normalization \(P(O) = \int\limits_q P(O|q)\cdot P(q)\). In our work, we use Hamiltonian Monte Carlo (HMC), an efficient variant of MCMC, which uses Hamiltonian Dynamics to propose the next sample. From an implementation perspective, we use the built-in HMC sampler in Stan for drawing samples.

Mixture Model

As mentioned before, we model the noisy observation as a mixture model of true distribution and outliers. This is denoted by (Equation 1):

O_{i, j} = \begin{cases}
q_i q_j^T + \eta_{i, j} & \textrm{with prob. } p \\
\textrm{Uniform}(\textrm{SO}(D)) & \textrm{with prob. } 1 – p

where \(\eta_{i j}\) is the additive noise to our true observation, \(q_i q_j^T\) is the relative orientation between the \(i^\textrm{th}\) and \(j^\textrm{th}\) objects, and \(\textrm{Uniform}(\textrm{SO}(D))\) is the uniform distribution over the rotation group \(\textrm{SO}(D)\), representing our outlier distribution. \(\textrm{SO}(D)\) is the space where every element represents a D-dimensional rotation. This mixture model serves as the likelihood \(P(O|q)\) for our Bayesian framework.

Sampled Result

The orientations \(\hat{q}_i\) sampled from the posterior look significantly rotated with respect to the original samples. Notice this is a global rotation since all the samples are rotated equally. This problem of global ambiguity of absolute orientations \(q_i\) arises from the fact that the relative orientations \(q_i q_j^T\) and \(\tilde{q}_i \tilde{q_j}^T\) of two different set of vectors can be the same even if the absolute orientations are different. The following section goes over this and provides a sneak peek into our solution for alleviating this problem.

Continuous Label Switching

A careful observation of our problem formulation (Equation 1) would reveal the problem is invariant to transformation of the absolute orientations as long as the relative orientations are preserved. Consider the 2 pairs of observations in the figure below. (Blue and Red; Yellow and Green) Let the absolute orientations be \(q_1\), \(q_2\), \(\tilde{q}_1\) and \(\tilde{q}_2\) and relative orientations between the pairs be \(R_{12}\) and \(\tilde{R}_{12}\). As the absolute orientations \(q_1\) and \(q_2\) are equally rotated by a rotation matrix \(A\), the relative orientation between them \(R_{12}\) is preserved.

More formally, Let \(A \in \textrm{SO}(D)\) be a random orientation matrix in D-dimensions. The following equation demonstrates how rotating two absolute orientations \(q_i\) and \(q_j\) by a rotation matrix \(A\) preserves the relative orientations — which in turn gives rise to a global ambiguity in our framework.

R_{ij} = q_i q_j = q_i A A^T q_j = (q_i A) (A^T q_j) = \tilde{q_i} \tilde{q_j}

Since our inputs to the model are relative orientations, this ambiguity (known as label switching) causes our Bayesian estimates to come randomly rotated by some rotation \(A\).

Proposed Solution

Based on Monteiller et al., in this project we explored a novel solution for alleviating this problem. The intuition is that we believe the unknown ground truth is close to the posterior samples up to a global rotation. Hence we try to approximate the ground truth by starting out with a random guess and optimizing for the alignment map between the estimate and the ground truth. Using this alignment map, and the posterior samples, we iteratively update the guess, using a custom Riemannian Stochastic Gradient Descent over \(\textrm{SO}(D)\).


  1. Start with a random guess \(\mu \sim \mathrm{Uniform}(\mathrm{SO}(D))\).
  2. Sample \(\hat{q} \sim P(q | O)\), where \(P(q | O)\) is the posterior.
  3. Find the global ambiguity \(R\), between \(\hat{q}\) and \(\mu\). This can be obtained by solving for \(\mathrm{argmin}_R \, \| \mu – R \hat{q}\|_\mathrm{F}\).
  4. Move \(\mu\) along the shortest geodesic toward \(\hat{q}\).
  5. Repeat Steps 2 – 4 until convergence. Convergence is detected by a threshold on  geodesic distance.


We use this method to estimate the mean of the posterior over \(\mathrm{SO}(2)\) and plot the results (i.e. 2D orientations) in the complex plane as shown below.


In conclusion, we study the rotation synchronisation problem from a Bayesian perspective. We explore a custom Riemannian Gradient Descent procedure and perform experiments in the \(\mathrm{SO}(2)\) case. The current method is tested on a simple toy dataset. As future work we are interested in improving our Bayesian model and benchmarking it against the current state-of-the-art. There are certain performance bottlenecks in our current architecture, which constrain us to test only on \(\mathrm{SO}(2)\). In the future, we are also interested in carrying out experiments more thoroughly in various dimensions. While the current MCMC procedure we are using does not account for the non-Euclidean geometry of the space of orientations, \(\mathrm{SO}(D)\), we are looking into replacing it with Riemannian versions of MCMC.

Math SGI research projects

Minimal Surfaces, But Periodic

By Zhecheng Wang, Zeltzyn Guadalupe Montes Rosales, and Olga Guțan

Note: This post describes work that has occurred between August 9 and August 20. The project will continue for a third week; more details to come.

I. Introduction

For the past two weeks we had the pleasure of working with Nicholas Sharp, Etienne Vouga, Josh Vekhter, and Erik Amezquita. We learned about a special type of minimal surfaces: triply-periodic minimal surfaces. Their name stems from their repeating pattern. Broadly speaking, a minimal surface minimizes its surface area. This is equivalent to having zero mean curvature. A triply-periodic minimal surface (TPMS) is a surface in \(\mathbb{R}^{3}\) that is invariant under a rank-3 lattice of translations.

Figure 1. (Left) A Minimal Surface [source] and (right) a TPMS [source].

Let’s talk about nonmanifold surfaces. “Manifold” is a geometry term that means: every local region of the surface looks like the plane (more formally — it is homeomorphic to a subset of Euclidean space). Non-manifold then allows for parts of the surface that do not look like the plane, such as T-junctions. Within the context of triangle meshes, a nonmanifold surface is a surface where more than 3 faces share an edge.

II. What We Did

First, we read and studied the 1993 paper by Pinkhall and Polthier that describes the algorithm for generating minimal surfaces. Our next goal was to generate minimal surfaces. Initially, we used pinned (Dirichlet) boundary conditions and regular manifold shapes. After ensuring that our code worked on manifold surfaces, we tested it on non-manifold input. 

Additionally, our team members learned how to use Blender. It has been a very enjoyable process, and the work was deeply satisfying, because of the embedded mathematical ideas intertwined with the artistic components. 

III. Reading the Pinkall Paper

“Computing Discrete Minimal Surfaces and Their Conjugates,” by Ulrich Pinkall and Konrad Polthier, is the classic paper on this subject; it introduces the iterative scheme we used to find minimal surfaces. Reading this paper was our first step in this project. The algorithm that finds a discrete locally area-minimizing surface is as follows:

  • Take the initial surface \(M_0\) with boundary \(∂M_0\) as first approximation of M. Set i to 0.
  • Compute the next surface \(M_i\) by solving the linear Dirichlet problem $$ \min_{M} \frac{1}{2}\int_{M_{i}}|\nabla (f:M_{i} \to M)|^{2}$$
  • Set i to i+1 and continue with step 2. 

The stopping condition is \(|\text{area}(M_i)-\text{area}(M_{i+1})|<\epsilon\). In our case, we used a maximum number of iterations, set by the user, as a stopping condition.

There are additional subtleties that must be considered (such as “what to do with degenerate triangles?”), but since we did not implement them — their discussion is beyond the scope of this post.

IV. Adding Periodic Boundary Conditions

This is, at its core, an optimization problem. To ensure that the optimization works, the boundary conditions have to be periodic instead of fixed in space. This is because we are enforcing a set of boundary conditions on periodic shapes — that is, tiling in a 3D space. 

IV(a). Matching Vertices

First, we need to check every two pairs of vertices in the mesh. We are looking  to see if they have identical coordinates in two dimensions, but are separated by exactly two units in the third dimension. When we find such pairs of vertices, we classify them to \(G_x\), \(G_y\), or \(G_z\). Note that we only store unique pairs.

IV(b). Laplacian Smoothness at the Boundary Vertices

Instead of using the discrete Laplacian, now we introduce a sparse matrix K to adjust our smoothness term based on the new boundary:

$$\min_{x}x^T(L^TK^TKL)x \text{ s.t. } x[b] = x_0[b].$$

Next, we construct the matrix K, which is a sparse square matrix of dimension #vertices by #vertices. To do so, we set \(K(i,i) = 1\), \(K(i,j) = 1\), and set the \(j\)th row entries to 0 for every pair of unique matched boundary vertices. For every interior vertex \(i\), we set \(K(i,i) = 1\).

IV(c). Aligning Boundary Vertices

Now we no longer want to pin boundary vertices to their original location in space. Instead, we want to allow our vertices to move, while the opposite sides of the boundary still match. To do that, we need to adjust the existing constraint term and to include additional linear constraints \(Ax=b\).

Therefore, we add two sets of linear constraints to our linear system:

  • For any pair of boundary vertices distanced by 2 units in one direction, the new coordinates should differ by 2 units.
  • For any pair of boundary vertices matched in the two other directions, the new coordinates should differ by 0.

We construct a selection matrix \(A\), which is a #pairs of boundary vertices in \(G\) by #vertices sparse matrix, to get the distance between any pair of boundary vertices. For every \(r\)th row, \(A(r,i)=1, A(r,j)=-1.\)$

Then we need to construct 3 \(b\) vectors, each of which is a sparse square matrix of the size [# of vector pairs of boundary vertices in G for a 3D coordinate (x,y,z)]. Based on whether at one given moment we are working with \(G_x\), \(G_y\), or \(G_z\), we enter \(2\) for those selected pairs, and \(0\) for the rest.

V. Correct Outcomes

Below, we can see the algorithm being correctly implemented. Each video represents a different mesh.

VI. Aesthetically Pleasing Bugs

Nothing is perfect, and coding in Matlab is no exception. We went through many iterations of our code before we got a functional version. Below are some examples of cool-looking bugs we encountered along the way, while testing the code on (what has become) one of our favorite shapes. Each video represents a different bug applied onto the same mesh.

VII. Conclusion and Future Work

Further work may include studying the physical properties of nonmanifold TPMS. It may also include additional basic structural simulations for physical properties, and establishing a comparison between the results for nonmanifold surfaces and the existing results for manifold surfaces.

Additional goals may be of computational or algebraic nature. For example, one can write scripts to generate many possible initial conditions, then use code to convert the surfaces with each of these conditions into minimal surfaces. An algebraic goal may be to enumerate all possible possibly-nonmanifold structures and perhaps categorize them based on their properties. This is, in fact, an open problem.

The possibilities are truly endless, and potential directions depend on the interests of the group of researchers undertaking this project further.