CS184/284A Summer 2025 Homework 2 Write-Up

De Casteljau's algorithm is a recursive algorithm used generally for Bezier curves to find points on the curve or split it into smaller curves. It repeatedly performs linear interpolations between consecutive points.

I implemented the evaluateStep function with a for loop running through points using the de Cateljau formula.

• Immediate output: given control points and curve
• Output 1: linear interpolations between the consecutive pairs
• Output 2: linear interpolations between the points from Output 1
• Output 3+: linear interpolations between the points from previous Output
• Final Output: one evaluated point on the curve

std::vector<&Vector2D> BezierCurve::evaluateStep(std::vector<&Vector2D> const &points)
  {
    std::vector<&Vector2D> nextPoints;

    // don't run if less than two points
    if (points.size() < 2) {
      return nextPoints;
    }

    for (size_t i = 0; i < points.size() - 1; ++i) {
      Vector2D interpolated = (1 - t) * points[i] + t * points[i + 1];
      nextPoints.push_back(interpolated);
    }
    return nextPoints;
  }
				

Bezier surfaces are made up of 1D Bezier curves in both u and v directions. So, the de Casteljau algorithm can be extended to surfaces by applying it twice: first across the rows (1st dimension) and then across the columns (2nd dimension).

The evaluateStep function evaluates just one step of de Casteljau on 3D points

  std::vector BezierPatch::evaluateStep(std::vector const &points, double t) const
  {
    std::vector nextPoints;

    // we need at least 2 points to interpolate
    if (points.size() < 2) {
      return nextPoints;
    }

    // interpolate between every pair of adjacent points
    for (size_t i = 0; i < points.size() - 1; ++i) {
      Vector3D interpolated = (1 - t) * points[i] + t * points[i + 1];
      nextPoints.push_back(interpolated);
    }

    return nextPoints;
  }
				

The evaluate1D function recursively evaluates a 1D Bezier curve to a single point.

  Vector3D BezierPatch::evaluate1D(std::vector const &points, double t) const
  {
    std::vector currPoints = points;

    // call recursively until only 1 point left
    while (currPoints.size() > 1) {
      currPoints = evaluateStep(currPoints, t);
    }

    // if empty, return default 3D vector
    if (currPoints.empty()) {
      return Vector3D(0, 0, 0);
    }

    return currPoints[0];
  }
				

The evaluate function evaluates the full Bezier surface at (u, v)

  Vector3D BezierPatch::evaluate(double u, double v) const
  {
    std::vector evaluatedPoints;

    // evaluate each row at u
    for (const std::vector& row : controlPoints) {
      // result is one point on that row's curve
      Vector3D pointOnRowCurve = evaluate1D(row, u);
      evaluatedPoints.push_back(pointOnRowCurve);
    }

    // evaluate vertical curve at v
    Vector3D surfacePoint = evaluate1D(evaluatedPoints, v);
    return surfacePoint;
  }
  					

The function Vertex::normal() approximates the unit normal at a vertex with area-weighted vertex normals. We weight each face normal by its area to make the rendered result more accurate.

Starting from a half-edge at the vertex, we compute the face normal of the current triangle using the cross product and add this to normalSum. We move to the next triangle around the vertex using twin() to get to the twin edge and then the next edge. We loop until we end up at the starting half-edge again. To make it a unit normal, we normalize the final vector and return.

Vector3D Vertex::normal( void ) const
  {
	// start with half-edge at vertex
    Vector3D normalSum(0, 0, 0);
    HalfedgeCIter start = halfedge();
    HalfedgeCIter curr = start;

    do {
      Vector3D p0 = curr->vertex()->position;
      Vector3D p1 = curr->next()->vertex()->position;
      Vector3D p2 = curr->next()->next()->vertex()->position;

		// take area-weighted avg of the normals of neighboring triangles
      Vector3D faceNormal = cross(p1 - p0, p2 - p0);
      normalSum += faceNormal;

      curr = curr->twin()->next();
    } while (curr != start);

	// normalize entire sum of normals
    normalSum.normalize();

    return normalSum;
  }
    

We implemented the edge flip operation by flipping the shared edge on two adjacent triangles, which remeshes the figure. The function checks that the adjacent faces we are looking at are not boundaries, or else exits. It then re-assigns half-edge, vertex, edge, and face pointers to reflect the new connectivity.

We made sure to never flip a boundary edge, or add/delete any mesh elements. Pointers were re-assigned instead of creating new mesh elements to maintain the same number of elements before and after the edge flip.

      if (e0->halfedge()->face()->isBoundary() ||
          e0->halfedge()->twin()->isBoundary()) {
          return e0;
      }

      HalfedgeIter BC = e0->halfedge();
      HalfedgeIter CA = BC->next();
      HalfedgeIter AB = CA->next();

      HalfedgeIter CB = BC->twin();
      HalfedgeIter BD = CB->next();
      HalfedgeIter DC = BD->next();

      VertexIter A = AB->vertex();
      VertexIter B = BC->vertex();
      VertexIter C = CA->vertex();
      VertexIter D = DC->vertex();

      FaceIter ABC = BC->face();
      FaceIter CBD = CB->face();

      BC->vertex() = A;
      CB->vertex() = D;
      
      HalfedgeIter AD = BC;
      HalfedgeIter DA = CB;

      FaceIter ADC = ABC;
      FaceIter DAB = CBD;

      AD->next() = DC;
      DC->next() = CA;
      CA->next() = AD;

      DA->next() = AB;
      AB->next() = BD;
      BD->next() = DA;

      DC->face() = ADC;
      AB->face() = DAB;

      ADC->halfedge() = AD;
      DAB->halfedge() = DA;
      
      A->halfedge() = AD;
      B->halfedge() = BD;
      C->halfedge() = CA;
      D->halfedge() = DA;
    			

We implemented the edge split operation by adding a vertex at the midpoint of an edge, using that point to create 4 smaller triangles that replace the original 2 triangles there.

For implementation, I didn't run into any particular debugging problems. I kept track of all edges and vertices, making sure to never create new ones. Instead, I only re-assigned these variables to avoid any errors typically caused by mistakenly creating new edges and vertices. The edges and vertices in the mesh were defined the same way as in part 4, so I will skip that part below.

Using the formula given in the spec, we calculate updated old vertex positions, and save them temporarily in v->newPosition.

for (VertexIter v = mesh.verticesBegin(); v != mesh.verticesEnd(); ++v) {
    v->isNew = false;

    int deg = v->degree();

    float u;
    if (deg == 3) {
        u = 3.0f / 16.0f;
    }
    else {
        u = 3.0f / (8.0f * deg);
    }

    Vector3D sum;

    // iterate thru all the neighborign vertices (from discussion!)
    HalfedgeIter curr = v->halfedge();
    HalfedgeIter end = curr;
    do {
        sum += curr->next()->vertex()->position;
        curr = curr->next()->next()->twin();
    } while (curr != end);

    v->newPosition = (1 - deg * u) * v->position + u * sum;
}
    			

Using the other formula given in the spec, we compute the new vertex positions, and save them temporarily in e->newPosition, since the vertex will soon be created on that edge

for (EdgeIter e : oldEdges) {
    e->isNew = false;

    HalfedgeIter AB = e->halfedge();
    HalfedgeIter BA = AB->twin();
    Vector3D newPos = (3.0f / 8.0f) * (AB->vertex()->position + BA->vertex()->position);
    newPos += (1.0f / 8.0f) * (AB->next()->next()->vertex()->position + BA->next()->next()->vertex()->position);

    e->newPosition = newPos;
}
    			

We'll go through and split all old edges, then go through again and check if any edge connects an old vertex to a new vertex, in which case we will flip it. note that we also assign the new vertex positions calculated in step 2.

// 3.
for (EdgeIter e : oldEdges) {
    VertexIter v = mesh.splitEdge(e);
    v->isNew = true;
    v->newPosition = e->newPosition;
}

// 4.
for (EdgeIter e = mesh.edgesBegin(); e != mesh.edgesEnd(); ++e) {
    if (e->isNew) {
        if (e->halfedge()->vertex()->isNew != e->halfedge()->twin()->vertex()->isNew) {
            mesh.flipEdge(e);
        }
    }
}
    			

To top it off, simply take the stored new vertex positions and assign them.

for (VertexIter v = mesh.verticesBegin(); v != mesh.verticesEnd(); ++v) {
    v->position = v->newPosition;
}