{"id":998,"date":"2021-08-25T19:44:48","date_gmt":"2021-08-25T19:44:48","guid":{"rendered":"http:\/\/summergeometry.org\/sgi2021\/?p=998"},"modified":"2021-08-25T19:44:49","modified_gmt":"2021-08-25T19:44:49","slug":"intrinsic-parameterization-weeks-1-2","status":"publish","type":"post","link":"https:\/\/summergeometry.org\/sgi2021\/intrinsic-parameterization-weeks-1-2\/","title":{"rendered":"Intrinsic Parameterization: Weeks 1-2"},"content":{"rendered":"\n

By Joana Portmann, Tal Rastopchin, and Sahra Yusuf.<\/em>
Mentored by Professor Keenan Crane.<\/em><\/p>\n\n\n\n

Intrinsic parameterization<\/h2>\n\n\n\n

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.<\/p>\n\n\n\n

From intrinsic to extrinsic parameterization<\/h4>\n\n\n\n

To get a short introduction into intrinsic parameterization and its applications, I quote some sentences from Keenan\u2019s course. If you\u2019re interested in the subject I can recommend the course notes “Geometry Processing with Intrinsic Triangulations<\/a>.” “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.” <\/p>\n\n\n\n

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.<\/p>\n\n\n\n

Laying out edge lengths in the plane<\/h2>\n\n\n\n

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\u2019s 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.<\/p>\n\n\n\n

Breadth-first search pseudocode<\/h4>\n\n\n\n

Some initial geometry central pseudocode for this breadth first search looks something like<\/p>\n\n\n

\n\/\/ Pick and flatten some starting triangle\nFace rootTriangle = mesh->face(0);\n\n\/\/ Calculate the root triangle\u2019s flattened vertex positions\ncalculateFlatVertices(rootTriangle);\n\n\/\/ Initialize a map encoding visited faces\nFaceData<bool> isVisited(*mesh, false);\n\n\/\/ Initialize a queue for the BFS\nstd::queue<Face> visited;\n\n\/\/ Mark the root triangle as visited and pop it into the queue\nisVisited[rootTriangle] = true;\nvisited.push(rootTriangle);\n\n\/\/ While our queue is not empty\nwhile(!visited.empty())\n{\n    \/\/ Pop the current Face off the front of the queue\n    Face currentFace = visited.front();\n    visited.pop();\n\n    \/\/ Visit the adjacent faces\n    For (Face adjFace : currentFace.adjacentFaces()) {\n        \/\/ If we have not already visited the face\n        if (!visited[adjFace]\n            \/\/ Calculate the triangle\u2019s flattened vertex positions\n            calculateFlatVertices(adjFace);\n            \/\/ And push it onto the queue\n            visited.push(adjFace);\n        }\n    }\n}\n\n<\/pre><\/div>\n\n\n
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.<\/p>\n\n\n\n
The halfedge mesh data structure<\/h4>\n\n\n\n
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<\/em> and next<\/em>. 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.<\/p>\n\n\n\n
\"\"
A diagram depicting the halfedge data structure connectivity relationships. Source<\/a>.<\/figcaption><\/figure><\/div>\n\n\n\n
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<\/a>.<\/p>\n\n\n\n
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.<\/p>\n\n\n
\n\/\/ Pop the current Face off the front of the queue\nFace currentFace = visited.front();\nvisited.pop();\n\n\/\/ Get the face\u2019s halfedge\nHalfedge currentHalfedge = currentFace.halfedge();\n\n\/\/ Visit the adjacent faces\ndo\n{\n    \/\/ Get the current adjacent face\n    Face adjFace = currentHalfedge.twin().face();\n    if (!isVisited[adjFace])\n    {\n        \/\/ Retrieve our desired vertices\n        Vertex a = currentHalfedge.vertex();\n        Vertex b = currentHalfedge.next().vertex();\n        Vertex c = currentHalfedge.twin().next().next().vertex();\n\n        \/\/ Calculate the triangle\u2019s flattened vertex positions\n        calculateFlatVertices(a, b, c);\n        \/\/ And push it onto the queue\n        visited.push(adjFace);\n    }\n    \/\/ Iterate to the next halfedge\n    currentHalfEdge = currentHalfEdge.next();\n}\n\/\/ Exit the loop when we reach our starting halfedge\nwhile (currentHalfEdge != currentFace.halfedge());\n\n<\/pre><\/div>\n\n\n
Here\u2019s a diagram illustrating the relationship between the currentHalfedge<\/code> and vertices a<\/code>, b<\/code>, and c<\/code>.<\/p>\n\n\n\n
\"\"
A diagram illustrating the connectivity relationship between the currentHalfedge<\/code> and vertices a<\/code>, b<\/code>, and c<\/code>. Note that  cH<\/code> abbreviates currentHalfedge.<\/code><\/figcaption><\/figure><\/div>\n\n\n\n
Segfaults, debugging, and ghost faces<\/h4>\n\n\n\n
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 \ud83e\udd70) we learned that the culprit was the isVisited[adjFace]<\/code>  function call on the line<\/p>\n\n\n
\nif (!isVisited[adjFace])\n<\/pre><\/div>\n\n\n
We were very confused as to why we would be getting a segfault while trying to retrieve the value corresponding to the key adjFace<\/code> contained in the map FaceData<bool> isVisited<\/code>. Sahra inspected the contents of the adjFace<\/code> 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<\/p>\n\n\n
\n\/\/ Get the current adjacent face\nFace adjFace = currentHalfedge.twin().face();\n<\/pre><\/div>\n\n\n
we realized that we made an unsafe assumption about currentHalfedge<\/code>. 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<\/code> was a boundary halfedge, why didn\u2019t we get a segfault on simply currentHalfedge.twin()<\/code>? Doing some research, we found that the geometry central internals documentation explains that<\/p>\n\n\n\n
“We can implicitly encode the twin()<\/code> relationship by storing twinned halfedges adjacent to one another\u2013 that is, the twin of an even-numbered halfedge numbered he<\/code> is he+1<\/code>, and the twin of and odd-numbered halfedge is he-1<\/code>.”<\/p>Geometry central internals documentation<\/a><\/cite><\/blockquote>\n\n\n\n
Aha! This explains exactly why currentHalfedge.twin()<\/code> still works on a boundary halfedge; behind the scenes, it is just adding or subtracting one to the halfedge\u2019s index. Where did the face index come from? We\u2019re still unsure, but we realized that the face currentHalfedge.twin().face()<\/code> only makes sense (and hence can only be used as a key for the visited map) when currentHalfedge<\/code> is not a boundary halfedge. Here is a diagram of the “ghost face” that we think the line Face adjFace = currentHalfedge.twin().face()<\/code> was producing.<\/p>\n\n\n\n
\"A
A diagram depicting how taking the face of the twin of a boundary halfedge produces a nonexistent face.<\/figcaption><\/figure><\/div>\n\n\n\n
Changing the map access line in the if statement to<\/p>\n\n\n
\nif (!currentHalfedge.edge().isBoundary() && !isVisited[adjFace])\n<\/pre><\/div>\n\n\n
resolved the segfaults and produced a working breadth-first traversal. <\/p>\n\n\n\n
Conformal parameterization<\/h2>\n\n\n\n
Here is a picture illustrating applying the flattening algorithm to a model of a cat head.<\/p>\n\n\n\n
\"\"
A picture illustrating the application of the flattening algorithm to a model of a cat head.<\/figcaption><\/figure><\/div>\n\n\n\n
You can see that there are many cracks and this is because the model of the cat head is not flat\u2014in 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\\).<\/p>\n\n\n\n
After laying out the edges on the plane \\(z=0\\), we began the necessary steps to compute a conformal flattening<\/em> (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,<\/p>\n\n\n\n
\\(l_{ij} := \\sqrt{ \\phi_i \\phi_j} l_{ij}^0, \\quad \\forall ij \\in E\\),<\/p>\n\n\n\n
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.<\/p>\n\n\n\n
Discrete Yamabe flow<\/h3>\n\n\n\n
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<\/em> Yamabe flow <\/em>to meet both of these requirements. Before implementing this algorithm, we began by substituting the scale factors with their logarithms<\/p>\n\n\n\n
\\(l_{ij} = e^{(u_i + u_j)\/2} l_{ij}^0\\),<\/p>\n\n\n\n
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.<\/p>\n\n\n\n
To implement the algorithm, we followed this procedure:<\/p>\n\n\n\n
1. Calculate the initial edge lengths<\/p>\n\n\n\n
2. While all angle defects are below a certain threshold epsilon:<\/p>\n\n\n\n
  • Compute the scaled edge lengths<\/li>
  • Compute the current angle defects based on the new interior angles based on the scaled edge lengths<\/li>
  • Update the scale factors using the step and the angle defects:<\/li><\/ul>\n\n\n\n
    \\(u_i \\leftarrow u_i – h \\Omega _i\\),<\/p>\n\n\n\n
    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.<\/p>\n\n\n\n
    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<\/em> which would change the step size based on the energy.<\/p>\n\n\n\n
    Results<\/h2>\n\n\n\n
    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!<\/p>\n\n\n\n