CS 184 Project 2 Writeup

Tianchen Liu 3034681144, Riley Peterlinz 3036754368

This doc is hosted on https://cal-cs184-student.github.io/project-webpages-sp23-tc-liu/proj2/index.html

Overview

For this project we work with modelling smooth 2D and triangular 3D geometry.

Section I: Bezier Curves and Surfaces

We implemented the lerp function using its definition, as well as implenting the evaluateStep function for each recursive step of the Bezier Curve. Then using the evaluateStep function, we implemented evalute1D that returns a single point on the Bezier curve on a vector of control point inputs.

By calling these functions, we can use De Casteljau’s algorithm to evaluate a Bezier Curve and Surface.

Section II: Triangle Meshes and Half-Edge Data Structure

For triangular meshes, we utilized the half-edge data structure to perform edge flipping and edge splitting in a way that retains the continuity of the mesh. We use the typical half edge data structure and utilize pointer assignments for edge flipping and splitting. We then use a combination of edge splitting and flipping to subdivide the mesh.

Section I: Bezier Curves and Surfaces

Task 1

Briefly explain de Casteljau’s algorithm and how you implemented it in order to evaluate Bezier curves.

De Casteljau’s algorithm is a recursive algorithm that evaluates a Bezier curve given a set of control points.

At each recursive step of the algorithm, it takes in $n$ control points and outputs $n-1$ control points. The $i$ th output control point of the recursive step is the result of calling the LERP function on the $i$ th and the $i+1$ th input control point, where the LERP function is defined as

$\text{lerp}(v_0, v_1, t) = (1-t)v_0 + tv_1$

The final evaluated point after the $n-1$ recursive steps lies on the curve.

The $t$ parameter is an adjustable parameter that signfies “how far the point is towards the destination”.

When $t=0$ , the final evaluated point lies on the beginning endpoint of the Bezier curve - which is where the first control point is.
When $t=1$ , the final evaluated point lies on the ending endpoint of the Bezier curve - which is where the final control point is.

I simply translated the algorithm described above into code. I wrote the LERP function as defined and for each recursive step of the evaluation function:

std::vector<Vector2D> BezierCurve::evaluateStep(std::vector<Vector2D> const &points)
  { 
    vector<Vector2D> ret;
    for (size_t i = 0; i < points.size() - 1; i++) {
      ret.push_back(lerp(points[i], points[i+1], t));
    }
    return ret;
  }

Take a look at the provided .bzc files and create your own Bezier curve with 6 control points of your choosing. Use this Bezier curve for your screenshots below.

Show screenshots of each step / level of the evaluation from the original control points down to the final evaluated point. Press E to step through. Toggle C to show the completed Bezier curve as well.

Show a screenshot of a slightly different Bezier curve by moving the original control points around and modifying the parameter tt via mouse scrolling.

Task 2

Briefly explain how de Casteljau algorithm extends to Bezier surfaces and how you implemented it in order to evaluate Bezier surfaces.

De Casteljau’s algorithm essentially evaluates any Bezier curve given a set of control points. Therefore, in order to extend this to 3D, we simply Break down the surface into a set of $n$ Bezier curves that lies on the surface, then perform De Casteljau’s algorithm once more on $n$ points, each point lies on one of the $n$ curves, with the $t$ param on each curve evaluate to the same value. This second-layer Bezier curve is controlled with another param, let’s say $u$ to avoid notation overload.

By adjusting $t$ and $u$ both from 0 to 1, we can evaluate any points on the Bezier surface.

Here’s a picture that demonstrates it perfectly:

To implement it, I first have a function that outputs a point on a Bezier curve given $n$ control points and the param $t$ .

Vector3D BezierPatch::evaluate1D(std::vector<Vector3D> const &points, double t) const
  {
    if (points.size() == 1) {
      return points[0];
    }
    return evaluate1D((evaluateStep(points, t)), t);
  }

where the evaluateStep function is the same as the one in Task 1 (except this one the $t$ is a param).

Then for every vector of points in the 2D vector controlPoints, I perform De Casteljau at param $u$ to gather n = size(controlPoints) interim points. Then I perform De Casteljau again on the $n$ interim points with param $v$ to get a point on the Bezier surface.

Vector3D BezierPatch::evaluate(double u, double v) const
  {  
    vector<Vector3D> calculated_surface_points;
    for (const auto& points : controlPoints) {
      calculated_surface_points.push_back(evaluate1D(points, u));
    }
    return evaluate1D(calculated_surface_points, v);
  }

Then for each point on the control

Show a screenshot of bez/teapot.bez (not .dae) evaluated by your implementation.

Section II: Triangle Meshes and Half-Edge Data Structure

The goal of this section is to create a function that subdivides our jagged meshes into something smooth and pleasant.

Task 3

Briefly explain how you implemented the area-weighted vertex normals.

I’m super proud of this task because I spent around 2 hours coming up with a solution and the code turned out to be seven-ish lines of code.

To explain it most clearly, I wrote my explanation of my implementation as comments in my code. Please see below.

  Vector3D Vertex::normal( void ) const
  {
    // Returns an approximate unit normal at this vertex, computed by
    // taking the area-weighted average of the normals of neighboring
    // triangles, then normalizing.

    Vector3D normal_vector;
    auto h = this->halfedge();
    // Loop through every half edge sprouting from this vertex
    do {
      // Obtain two vectors representing the two edges of a triangle neighboring the vertex.
      auto v1 = h->next()->vertex()->position - h->vertex()->position;
      auto v2 = h->next()->next()->vertex()->position - h->next()->vertex()->position;
      // The cross product of the two edges is
      // 1) perpendicular to the two edges, thus normal (perpendicular) to the neighboring triangle
      // 2) of length (Area of triangle * 2)
      // Hence we only need to sum up all the cross product vectors and then normalize it to get the area-weighted average of the normals of neighboring triangles.
      normal_vector += cross(v1, v2);
      // To get the next half edge, get its twin edge and the half edge that sprouts from there
      h = h->twin()->next();
    } while (h != this->halfedge());

    // Normalize by calling the .unit() method
    return normal_vector.unit();
  }

Show screenshots of dae/teapot.dae (not .bez) comparing teapot shading with and without vertex normals. Use Q to toggle default flat shading and Phong shading.

With (Phong shading):

Without (flat shading):
Without

Task 4

Briefly explain how you implemented the edge flip operation and describe any interesting implementation / debugging tricks you have used.

Firstly, let’s check if the edge is a boundary and name all the edges and half edges for convenience and ease to debug:

    if (e0->isBoundary()) return e0;

    auto h = e0->halfedge();
    auto h_twin = h->twin();
    auto BC = h;
    auto CB = h_twin;
    auto BD = CB->next();
    auto DC = BD->next();
    auto CA = BC->next();
    auto AB = CA->next();
    auto B = BC->vertex();
    auto D = DC->vertex();
    auto C = CB->vertex();
    auto A = AB->vertex();

Next, let’s reassign pointers of edges, faces, and vertices. See comments for explanation of implementation

    // Edges no need to change
    // Faces:
    // Make sure the half edges associated with the two faces are the two half edges flipped, respectivelly
    auto face = h->face();
    face->halfedge() = h;
    auto face_twin = h_twin->face();
    face_twin->halfedge() = h_twin;
    // Vertices:
    // For the vertex where the half-edge sprouts from pre-flip, make sure that the new half-edge associated with the vertex is any *other* half edge that sprouts from this vertex.
    // In this way, we make sure that the half-edges assoicated with the vertices are not affected by the edge flip.
    // Simply aquire another half-edge that sprouts from the vertex by calling the .next() method
    auto vertex_root = h->vertex();
    vertex_root->halfedge() = h_twin->next();
    auto vertex_root_twin = h_twin->vertex();
    vertex_root_twin->halfedge() = h->next();

Finally, reassign pointers of the half edges. After flipping, BC (h) is now DA and CB (h_twin) is now AD

    // The next of AB and DC need to change
    AB->setNeighbors(BD, AB->twin(), AB->vertex(), AB->edge(), AB->face());
    DC->setNeighbors(CA, DC->twin(), DC->vertex(), DC->edge(), DC->face());

    // The next and face of BD and CA need to change
    BD->setNeighbors(h, BD->twin(), BD->vertex(), BD->edge(), face);
    CA->setNeighbors(h_twin, CA->twin(), CA->vertex(), CA->edge(), face_twin);

    // h is now DA and h_twin is now AD
    h->setNeighbors(AB, h_twin, D, h->edge(), h->face());
    h_twin->setNeighbors(DC, h, A, h_twin->edge(), h_twin->face());

The debugging trick that I believe helped me the most is honestly simply drawing the picture out, giving a name to all half edges and verticese. Once I drew the diagram I just reassigned all the pointers using the named pointers. This makes it way easier to debug too.

Of course, printing out info using the check_for function is a classic debugging technique that is incredibly useful here.

Show screenshots of the teapot before and after some edge flips.

Before edge flips:

After edge flips:

Write about your eventful debugging journey, if you have experienced one.

I was getting holes after I flip an edge:
https://edstem.org/us/courses/33060/discussion/2594393?comment=5984155

So I debugged by printing out half edge info of the two half edge associated with the edge that I flipped. Turns out the bug is caused by me not accounting for all reassignments of half-edges. I only reassigned pointers of half-edges on the flipped edge but not the other four.

So I added these lines and it fixed the bug:

    // The next of AB and DC need to change
    AB->setNeighbors(BD, AB->twin(), AB->vertex(), AB->edge(), AB->face());
    DC->setNeighbors(CA, DC->twin(), DC->vertex(), DC->edge(), DC->face());

    // The next and face of BD and CA need to change
    BD->setNeighbors(h, BD->twin(), BD->vertex(), BD->edge(), face);
    CA->setNeighbors(h_twin, CA->twin(), CA->vertex(), CA->edge(), face_twin);

Task 5

Implementation

There are three steps to splitting an edge in our program:

Naming the geometry
Creating New Geometric Elements
Assigning the Relationship between Elements

1. Naming the Geometry

We start with labelling the edge we want to split BC and its components like so:

Splitting the edge means placing a vertex M at the midpoint of BC and then also connecting M to A and D

2. Creating New Geometric Elements

Excluding gory details for the code, we did not delete any geometric elements. Rather repurposing existing elements. This was to help retain relationships that already exist in the original geometry and reduce the probability of bugs.

Half Edges:
BC (left) -> BM (lower left), CB(right) -> CM (upper right)
Edges:
BCCB (center) -> BMMB (lower center)
Faces:
BCA(left) -> BMA (lower left), CBD(right) -> CMD(upper right)

All other elements were newly created.

3. Assigning the Relationship between Elements

Once all the elements were created we used the setNeighbors method for HalfedgeIter to assign the relationships from halfedges to faces, vertices, and edges.

Debugging

Surprisingly, we did not run into many issues with debugging. I think just keeping the picture helped. However, like any well-thought-out piece of software, the first run failed. After splitting, one of the triangle faces went missing. Since it didn’t crash, I figured that we had just not assigned the face correctly rather than missed creating one entirely. Sure enough there was a typo in the face assignment that, once fixed, gave us beautifully split geometry.

Results

Let’s start with the cow’s face

First, I split all the triangles of the square in the middle

Further splitting!

Now let’s flip the edges to get a sort of diamond pattern <>

Even more Splitting

Extra Credit

We can also split the boundary edges!

This was achieved by doing a similar split to normal edges. The only real difference is that for boundary edges, their twin points to the next boundary edge and its face is considered as the hole left by the boundaries on the mesh.

Task 6

A culmination of our work in flipping and splitting, now we can subdivide a whole mesh!

Implementation

We followed the steps recommended in the spec

Compute the positions of both new and old vertices using the original mesh
Subdivide the original mesh via edge splits and flips as described
Update all vertex positions in the subdivided mesh using the values already computed

1. Compute the positions of both new and old vertices using the original mesh

First, we iterated through all the edges of the mesh to calculate where the new point will be placed.

This new position is calculated from the vertices that attach the edge A and B, and the vertices it bisects: C and D. The position of this new vertex is:
3.0 / 8.0 * (B->position + C->position) + 1.0 / 8.0 * (A->position + D->position)
we store the new position in e->newPosition until we are ready to place the vertex there.

The old vertices have a slightly more complicated new location. The location depends on where all the vertices it is connected to are before subdivision:

The position for old vertices in a newly subdivided mesh is:
(1 - n * u) * original_position + u * original_neighbor_position_sum

Finding the variables to the above equation requires looping through the existing vertices v of our mesh. Luckily, n is already an attribute of a vertex so u is trivially calculated. The original_position is also easy.

Getting original_neighbor_position_sum is tricky, as we need to find all the connected vertices. We wrote a loop that traverses the vertices by considering the half edge of our current vertex v. If we take the twin of the half edge connected to v, the vertex of the twin will be u. To cycle around the u’s we simply rotate to the next() half edge and repeat the proceess until we are back at v.

We also store the new position of existing vertices in v->newPosition

2. Subdivide the original mesh via edge splits and flips as described

Once we have our positions, we neeed to commence splitting the mesh then flipping select edges.

First, we split every edge down the middle. We iterate over every edge of a mesh and only split an edge if the vertices that connect the edge are labelled as isNew = 1. This is to avoid over-subdividing and infinite loops. We also assign the newPosition of the edge we split e to the newPosition of the vertex created v.

v->newPosition = e->newPosition;

Next, we flip only new edges e that are created (blue in the figure) that arre connected to exactly one old vertex and one new vertex.

3. Update all vertex positions in the subdivided mesh using the values already computed

Finally, for each vertex in the mesh, we place it into position using

v->position = v->newPosition;

Debugging

No debugging was needed for this part :D (what a relief)

Results/Analysis

For objects with no protruding corners, like a bean, subdivision does very to the overall geometry:

After one subdivision:

It went from bean to more beany

However, an object with protruding corners like a cube becomes smaller since the averaging of subdivision smooths out the corners and edges like digital sandpaper:

However, notice that the cube is oblong after subdivision! This is because the cube is made of triangle that are anisotropic on a face. We can make these triangles isotropic on a face by splitting them before subdividing.

Now the cube is symmetric with respect to its original 6 planar faces. Notice that it is not a sphere because it still retains spatial preference towards its extraneous edges.