Overview

Learned about the neat (although sometimes painful) halfedge data structure, which gives us a way to navigate any triangular mesh quite nicely (except maybe when you're trying to implement split or flip code lol). We learned smoothing algorithms that allowed us to interpolate in 3D space, as well as subdivide to smooth out more irregular meshes, and in general became more familiar with 3D mesh manipulation techniques that are cool to see in action (subdivision for example is big in 3D modeling).

Section I: Bezier Curves and Surfaces

Part 1: Bezier curves with 1D de Casteljau subdivision

De Casteljau's algorithm is an interpolation algorithm that in a sense parameterizes a smooth (except for maybe endpoints) curve. It's pretty ingenious, taking the line segments between points, using the interpolating parameter t to find a point on each line segment to use for the next set of line segments, and recursing until we have a single point left, which we use as a point on the curve. Moving t from 0 to 1 traces the entire smooth curve.

Below is a neat gif demo, where changing t moves the red point along the curve, and moving control points changes that point in our list of points/line segments (which changes calculations of the blue lines, and with it, the red point).

Part 2: Bezier surfaces with separable 1D de Casteljau subdivision

We extended de Casteljau to Bezier surfaces by basically interpolating between first the u "rows" to get a 1D "column" of positions, then doing a final de Casteljau to interpolate across the 1D u "column" using v to get our final point. This was done for every u,v position on the teapot.

Essentially, this works because we interpolate across 1 dimension at a time, which allows us to use our regular de Casteljau on one dimension at a time, reducing dimensions until we finally reach our single point, kinda like map reduce being able to reduce multidimensional arrays into a single value. The image below is a nice visualization:

Section II: Sampling

Part 3: Average normals for half-edge meshes

We implemented Area-weighted vertex normals by iterating through all triangles incident (touching) the vertex, where at each triangle we calculating the cross product of the two vectors pointing away from the vertex, which gave us the normal vector, then using applied math skillz we knew the length of the cross product divided by 2 gave us the triangle area.

We then weighted the normal by the area and added it to a running sum.

To traverse to subsequent triangles, we went from our current halfedge, which pointed from the vertex to a point within our old triangle, to its twin's next to get an identical halfedge, except translated one triangle over.

Afterwards, we took the total sum of all the weighted normals and normalized it, which gave us the vertex normal.

Part 4: Half-edge flip

We implemented edge flip by treating the input edge as CB in the picture above and creating variables for each of halfedges, edges, vertices, and faces. From there, we simply assigned everything to be the same, except we switched cb to be ad and bc to be da, making sure next was set to be the next counterclockwise edge (eg: bc's next became ab), the vertex was correct (eg: bc's vertex became d since we wanted bc to be da, which started from d).

In terms of debugging, there was a lot of stress testing, basically clicking random edges and flipping until we got black faces or other visual artifacts. Also, pressing Q to enable Phong shading revealed flipping created black splotches, which was solved by setting all our references correctly. Making sure that 2 faces, 5 edges, 4 vertices, and 6 halfedges were set was also key, since these were all the elements involved (originally we only setup halfedges). The real debugging for this happened in part 6 :')

To get into technical details, the vertices accessed from the half-edges after calling next() in both faces two times was set as the source vertices for the twin and half-edge of the flipped edge. Every element pointer for half-edges in the triangle was set to the correct pointer according to the flipped edge. This meant that the vertices on the flipped edge became the source of just one half-edge instead of two before the flip.

Additionally, the half-edge pointers for each vertex and face was set to the correct half-edge counter-clockwise. In the beginning of the function, we added a check to prevent a boundary edge from getting flipped.

When debugging, we also used the check_for() debugging tool to check if the element pointers for the half-edges had changed.

Part 5: Half-edge split

Similar to before, we stored all the relevant edges, faces, etc. as variables, then created a new vertex, 3 new edges (not 4 since we can reuse an old edge) and 2 new faces (not 4 since we can reuse 2). We set all the relevant pointers according to the diagram above.

For debugging, splitting then flipping sometimes made halfedges point from the origin to its place on the mesh. This was solved by going through our implementation and checking to see if we set the correct pointers for everything, which we eventually did. Again, stress testing was key, clicking random spots and splitting over and over to see if something broke. The real debugging for this happened in part 6 :')

At a more detailed, technical level, to implement split, we first accessed and stored the half edges, vertices, edges, and faces that make up the two triangles in which the edge lies between. Then, we created a new vertex on the split edge with its position set as the average of the positions of the vertices it's positioned between. The isNew class variable of the new vertex is set to true. Then, we create three more edges (while the original edge is used as the fourth) to set up the new edges and old edges. Six new half edges are then created for the diagonal edges in reference to the diagram shown in class, with the two edges having their isNew flag set to true. Then we assign the appropriate half edges for the faces, edges, and vertices. At the end of our implementation, we set up all of the half-edges appropriately using the function setNeighbours(). We also implemented boundary checks so that we return the split edge's vertex if the edge is located on a boundary.

Part 6: Loop subdivision for mesh upsampling

At a high level, we first calculate the "new positions" of old and new vertices based off an interpolation formula, store them in the old_vertex->newPosition and old_edge->newPosition respectively, then split all old edges, then flip edges that connect an old and new vertex, then using our precalculated "new positions" move the old and new vertices to their interpolated positions.

In terms of debugging, this was rich. We had to store newPosition in old vertices, which was fine since they didn't really get modified, but for storing it in edges, we had to go back and fix up part 5 to use that old edge before it split, because otherwise it would've gotten rid of the old edge's newPosition value. Other debugging included:

Printing out positions (They were 0 in step 5 but nonzero in step 1/2, indicating the newPositions were calculated but thrown away).

Checking newPosition and isNew were set correctly

Making everything a float/double (except n/degree when necessary)

"Unit testing:" First testing just splitting everything, then splitting and flipping, then setting positions

Manually splitting and flipping to test if they work together (if not usually a bug with part 4 or 5)

Manually calculating positions to see if they matched up (see image below)

We see that the icosahedron when put through iterations of loop subdivision results in a sphere. This is because loop subdivision smooths out the sharp corners and edges due to our interpolation of vertices. This effect could be reduced by pre-splitting some edges near the sharp corners/edges, because as a vertex has more neighbors near that corner/edge, the original vertices get moved less and less by subdivision, because u grows smaller and smaller.

The cube when put through several loop subdivisions becomes asymmetric, which could be lessened by preprocessing that makes edges symmetric. The asymmetry is because the mesh itself started asymmetric; the diagonal line is only symmetric across one axis, not both axes in the viewing plane (below the asymmetric cube is a symmetric version where splitting the diagonal made it symmetric) Making it symmetric allows the mesh to subdivide itself symmetrically; visually the two versions better explain this.

Here is a more symmetric one:

Going more into the details, loop subdivision required a total of five parts to implement. The first thing we implemented was computing the updated positions of the old vertices by iterating through all vertices, getting all neighboring vertices for each vertex, and storing the updated position of the vertex in its newPosition variable. The updated position for an old vertex was calculated using the equation (1 - n * u) * original_position + u * original_neighbor_position_sum where n refers to the degree of the vertex while u refers to a constant set by the value of n. u was 3/16 if n was equal to 3, otherwise, it was calculated as 3/(8n). Then, the positions of new vertices were calculated by iterating through all edges in the mesh, getting the vertex associated with that edge, and using 3/8 * (A + B) + 1/8 * (C + D) as shown in lecture. In addition, when iterating through all the edges, all old edges were stored in a vector to be iterated through in the next part of the implementation where we split all old edges. The next part of our implementation iterated through all edges again to flip all new edges. The two vertices that make up any new edges we find were accessed and if one was a new vertex while the other an old vertex, we flipped that edge. We accessed the isNew variable for each vertex (set before in splitEdge and in this part) to check if they were new or old. Then finally, we iterated through all vertices in the mesh to update their positions. This was done by setting a vertex's position to its newPosition.