Sync

src/catmull_clark.cpp

@@ -2,6 +2,22 @@
 #include <unordered_map>
 #include <utility>
 #include <functional>
+#include <vector>
+
+/**
+ * Compute the centroids of the faces of a mesh.
+ *
+ * Inputs:
+ * @param V  #V by 3 list of vertex positions
+ * @param F  #F by 3 list of triangle mesh indices into V
+ * Outputs:
+ * @param FC  #F by 3 list of face centroids
+ */
+void compute_face_centroids(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  Eigen::MatrixXd & FC
+);
 
 void catmull_clark(
   const Eigen::MatrixXd & V,
@@ -10,9 +26,59 @@ void catmull_clark(
   Eigen::MatrixXd & SV,
   Eigen::MatrixXi & SF)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Replace with your code here:
+  Eigen::MatrixXd VF;
+
+  /* Calculate new vertices. */
+
+  // Calculate face centroids into VF.
+  compute_face_centroids(V, F, VF);
+
+  // Calculate edge points into VE.
+  // TODO
+
+  // Move original vertices to new positions into VV.
+  // TODO
+
+  /* Form edges and faces in the new mesh. */
+
+  // Connect each new face point (VF) to the new edge points (VE) of all original edges defining the original face.
+  // TODO
+
+  // Connect each new vertex point to the new edge points of all original edges incident on the original vertex.
+  // TODO
+
+  // Define the new faces as enclosed by edges.
+  for (int i = 0; )
+
+  // TODO Replace me
   SV = V;
   SF = F;
-  ////////////////////////////////////////////////////////////////////////////
+}
+
+void extract_edges(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  std::unordered_map<std::pair<int, int>, std::pair<int, int>, std::function<size_t(const std::pair<int, int> &)>> & edge_map
+) {
+  for (int i = 0; i < F.rows(); i++) {
+    for (int j = 0; j < F.cols(); j++) {
+      std::pair<int, int> edge(F(i, j), F(i, (j + 1) % F.cols()));
+      if (edge_map.find(edge) == edge_map.end()) {
+        edge_map[edge] = std::make_pair(i, -1);
+      } else {
+        edge_map[edge].second = i;
+      }
+    }
+  }
+}
+
+void compute_face_centroids(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  Eigen::MatrixXd & FC
+) {
+  FC.resize(F.rows(), 3);
+  for (int i = 0; i < F.rows(); i++) {
+    FC.row(i) = (V.row(F(i, 0)) + V.row(F(i, 1)) + V.row(F(i, 2))) / 3;
+  }
 }

src/per_corner_normals.cpp

@@ -2,7 +2,6 @@
 #include "triangle_area_normal.h"
 // Hint:
 #include "vertex_triangle_adjacency.h"
-#include <iostream>
 
 void per_corner_normals(
   const Eigen::MatrixXd & V,
@@ -11,7 +10,30 @@ void per_corner_normals(
   Eigen::MatrixXd & N)
 {
   N = Eigen::MatrixXd::Zero(F.rows()*3,3);
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here:
-  ////////////////////////////////////////////////////////////////////////////
+  Eigen::VectorXi count = Eigen::VectorXi::Zero(F.rows()*3);
+  std::vector<std::vector<int>> VF;
+
+  vertex_triangle_adjacency(F, V.rows(), VF);
+
+  double cos_threshold = std::cos(corner_threshold * M_PI / 180.0);
+
+  /* Compute the area-weighted average of normals of faces with normals that deviate less than threshold. */
+  for (int i = 0; i < F.rows(); i++) { /* For each face */
+    Eigen::RowVector3d face_area_normal = triangle_area_normal(V.row(F(i,0)), V.row(F(i,1)), V.row(F(i,2)));
+    for (int j = 0; j < 3; j++) { /* For each corner of the face */
+      for (int g = 0; g < VF[F(i,j)].size(); g++) { /* For each face neighboring the corner */
+        int face = VF[F(i, j)][g];
+        Eigen::RowVector3d neighbor_area_normal = triangle_area_normal(V.row(F(face,0)), V.row(F(face,1)), V.row(F(face,2)));
+        double dot = neighbor_area_normal.normalized().dot(face_area_normal.normalized());
+        if (std::max(-1.0, std::min(1.0, dot)) > cos_threshold) {
+          N.row(i * 3 + j) += neighbor_area_normal;
+          count(i * 3 + j)++;
+        }
+      }
+    }
+  }
+  /* Normalize the sum of the normals by the number of faces neighboring each corner. */
+  for (int i = 0; i < N.rows(); i++) {
+    N.row(i) = (N.row(i) / count(i)).normalized();
+  }
 }

src/per_face_normals.cpp

@@ -6,8 +6,9 @@ void per_face_normals(
   const Eigen::MatrixXi & F,
   Eigen::MatrixXd & N)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Replace with your code:
-  N = Eigen::MatrixXd::Zero(F.rows(),3);
-  ////////////////////////////////////////////////////////////////////////////
+  N = Eigen::MatrixXd::Zero(F.rows(), 3);
+
+  for (int i = 0; i < F.rows(); i++) {
+    N.row(i) = triangle_area_normal(V.row(F(i, 0)), V.row(F(i, 1)), V.row(F(i, 2))).normalized();
+  }
 }

src/per_vertex_normals.cpp

@@ -1,5 +1,6 @@
 #include "per_vertex_normals.h"
 #include "triangle_area_normal.h"
+#include "vertex_triangle_adjacency.h" /* Why not use the method that does exactly what we want to do (but faster)? */
 
 void per_vertex_normals(
   const Eigen::MatrixXd & V,
@@ -7,7 +8,20 @@ void per_vertex_normals(
   Eigen::MatrixXd & N)
 {
   N = Eigen::MatrixXd::Zero(V.rows(),3);
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here:
-  ////////////////////////////////////////////////////////////////////////////
+  Eigen::VectorXi count = Eigen::VectorXi::Zero(V.rows());
+  std::vector<std::vector<int>> VF;
+
+  vertex_triangle_adjacency(F, V.rows(), VF);
+
+  /* Calculate the product of the area and the normal of each face neighboring each vertex. */
+  for (int i = 0; i < V.rows(); i++) {
+    for (int j = 0; j < VF[i].size(); j++) {
+      N.row(i) += triangle_area_normal(V.row(F(VF[i][j],0)), V.row(F(VF[i][j],1)), V.row(F(VF[i][j],2)));
+      count(i)++;
+    }
+  }
+  /* Normalize the sum of the normals by the number of faces neighboring each vertex. */
+  for (int i = 0; i < V.rows(); i++) {
+    N.row(i) = (N.row(i) / count (i)).normalized();
+  }
 }

src/triangle_area_normal.cpp

@@ -6,8 +6,7 @@ Eigen::RowVector3d triangle_area_normal(
   const Eigen::RowVector3d & b, 
   const Eigen::RowVector3d & c)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Replace with your code:
-  ////////////////////////////////////////////////////////////////////////////
-  return Eigen::RowVector3d(0,0,0);
+  Eigen::RowVector3d normal = (b - a).cross(c - a);
+  double area = 0.5 * normal.norm(); /* https://en.wikipedia.org/wiki/Cross_product#Geometric_meaning */
+  return normal.normalized() * area;
 }

src/vertex_triangle_adjacency.cpp

@@ -6,8 +6,10 @@ void vertex_triangle_adjacency(
   std::vector<std::vector<int> > & VF)
 {
   VF.resize(num_vertices);
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here:
-  ////////////////////////////////////////////////////////////////////////////
-}
 
+  for (int i = 0; i < F.rows(); i++) {
+    for (int j = 0; j < F.cols(); j++) {
+      VF[F(i, j)].emplace_back(i);
+    }
+  }
+}