Fixed it

src/catmull_clark.cpp

@@ -2,6 +2,23 @@
 #include <unordered_map>
 #include <utility>
 #include <functional>
+#include <iostream>
+
+/**
+ * @brief Compute the vertices and faces after Catmull-Clark subdivision. Perform one pass of the refinement scheme.
+ *
+ * Inputs:
+ * @param V  #V by 3 list of vertex positions
+ * @param F  #F by 4 list of quad mesh indices into V
+ * Outputs:
+ * @param SV  #SV by 3 list of vertex positions after subdivision
+ * @param SF  #SF by 4 list of quad mesh indices into SV after subdivision
+ */
+void subdivide_quadmesh_catmullclark(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & F,
+        Eigen::MatrixXd & SV,
+        Eigen::MatrixXi & SF);
 
 void catmull_clark(
   const Eigen::MatrixXd & V,
@@ -10,9 +27,201 @@ void catmull_clark(
   Eigen::MatrixXd & SV,
   Eigen::MatrixXi & SF)
 {
-  ////////////////////////////////////////////////////////////////////////////
+  Eigen::MatrixXd SVX = V;
-  // Replace with your code here:
+  Eigen::MatrixXi SFX = F;
-  SV = V;
+
-  SF = F;
+  for (int i = 0; i < num_iters; i++) {
-  ////////////////////////////////////////////////////////////////////////////
+    subdivide_quadmesh_catmullclark(SVX, SFX, SV, SF);
+    SVX = SV;
+    SFX = SF;
+  }
+}
+
+/**
+ * @brief Compute the centroid of a face.
+ *
+ * Inputs:
+ * @param V  #V by 3 list of vertex positions
+ * @param f  4 by 1 list of face indices into V
+ * Outputs:
+ * @return  1 by 3 list of the face centroid
+ */
+Eigen::RowVector3d compute_face_centroid(
+  const Eigen::MatrixXd & V,
+  const Eigen::RowVector4i & f)
+{
+  return (V.row(f[0]) + V.row(f[1]) + V.row(f[2]) + V.row(f[3])) / 4.0;
+}
+
+/**
+ * @brief Compute the average point of an edge.
+ *
+ * Inputs:
+ * @param V  #V by 3 list of vertex positions
+ * @param F  #F by 4 list of quad mesh indices into V
+ * @param FC #F by 3 list of face centroids
+ * @param e  2 by 1 list of edge indices into V, which make up the edge
+ * Outputs:
+ * @return  1 by 3 list of the edge average point
+ */
+Eigen::RowVector3d compute_edge_average_point(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & FC,
+  const Eigen::RowVector2i & e)
+{
+  /* M + E */
+  Eigen::RowVector3d edge_average_point = V.row(e[0]) + V.row(e[1]);
+
+  /* Determine the two faces that share the edge: A + F */
+  int counter = 0;
+  for (int i = 0; i < F.rows(); i++) {
+    for (int j = 0; j < F.cols(); j++) {
+      if (F(i, j) == e[0] && (F(i, (j + 1) % 4) == e[1] || F(i, (j + 3) % 4) == e[1])) {
+        edge_average_point += FC.row(i);
+        counter++;
+      }
+    }
+  }
+
+  if (counter != 2) {
+    throw std::runtime_error("An edge should be shared by exactly two faces.");
+  }
+
+  /* (A + F + M + E) / 4 */
+  return edge_average_point / 4.0;
+}
+
+/**
+ * @brief Compute the vertex position after Catmull-Clark subdivision.
+ *
+ * Inputs:
+ * @param V  #V by 3 list of vertex positions
+ * @param F  #F by 4 list of quad mesh indices into V
+ * @param v  index of the vertex to be moved
+ * Outputs:
+ * @return  1 by 3 list of the new vertex position
+ */
+Eigen::RowVector3d compute_vertex_moved_position(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & FC,
+  const int v)
+{
+  Eigen::RowVector3d part_F = Eigen::RowVector3d::Zero();
+  Eigen::RowVector3d part_R = Eigen::RowVector3d::Zero();
+
+  /* Determine the faces that share the vertex: F */
+  int f = 0;
+  for (int i = 0; i < F.rows(); i++) {
+    for (int j = 0; j < F.cols(); j++) {
+      if (F(i, j) == v) {
+        part_F += FC.row(i);
+        f++;
+      }
+    }
+  }
+  part_F /= f;
+
+  if (f == 0) {
+    std::cerr << "Error: A vertex should be shared by at least one face, but it is shared by " << f << " faces." << std::endl;
+    throw std::runtime_error("A vertex should be shared by at least one face.");
+  }
+
+  int r = 0;
+  /* Determine the edges that share the vertex: 2R */
+  /* NOTE that the edge midpoint point is computed twice for each edge */
+  for (int i = 0; i < F.rows(); i++) {
+    for (int j = 0; j < F.cols(); j++) {
+      if (F(i, j) == v) {
+        part_R += (V.row(F(i, j)) + V.row(F(i, (j + 1) % 4))) / 2.0;
+        r++;
+      }
+    }
+  }
+  part_R /= r;
+
+  /* (F + 2R + (n - 3)P) / n */
+  return (part_F + 2 * part_R + (f - 3) * V.row(v)) / f;
+}
+
+struct RowVector2iHash {
+  std::size_t operator()(const Eigen::RowVector2i & v) const {
+    return std::hash<int>()(v[0]) ^ std::hash<int>()(v[1]);
+  }
+};
+
+void subdivide_quadmesh_catmullclark(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  Eigen::MatrixXd & SV,
+  Eigen::MatrixXi & SF)
+{
+  SV.resize(V.rows() + 3 * F.rows(), 3); // Moved vertices + face centroids + edge average points
+  SF.resize(F.rows() * 4, 4);
+
+  std::unordered_map<Eigen::RowVector2i, int, RowVector2iHash> EAP; // Edge average points
+
+  for (int i = 0; i < F.rows(); i++) { /* For each face */
+    /* Compute face centroids: FC */
+    SV.row(V.rows() + i) = compute_face_centroid(V, F.row(i));
+  }
+
+  Eigen::MatrixXd FC = SV.block(V.rows(), 0, F.rows(), 3);
+
+  for (int i = 0; i < F.rows(); i++) {
+    /* Compute edge average points */
+    for (int j = 0; j < F.cols(); j++) { /* For each edge of that face */
+      int v1 = std::min(F(i, j), F(i, (j + 1) % 4));
+      int v2 = std::max(F(i, j), F(i, (j + 1) % 4));
+
+      /* Check if the edge has been computed before */
+      if (EAP.find(Eigen::RowVector2i(v1, v2)) == EAP.end()) {
+        SV.row(V.rows() + F.rows() + EAP.size()) = compute_edge_average_point(
+                V,
+                F,
+                FC,
+                Eigen::RowVector2i(v1, v2));
+        EAP[Eigen::RowVector2i(v1, v2)] = V.rows() + F.rows() + EAP.size();
+      }
+    }
+  }
+
+  /* Compute moved vertices: MV */
+  /* Vertices get moved but retain their original indices. */
+  for (int i = 0; i < V.rows(); i++) {
+    SV.row(i) = compute_vertex_moved_position(V, F, FC, i);
+  }
+
+  /* Compute the faces */
+  for (int i = 0; i < F.rows(); i++) { /* For each face */
+    int face_centroid = V.rows() + i;
+    int edge_average_points[4] = {
+      EAP.at(Eigen::RowVector2i(std::min(F(i, 0), F(i, 1)), std::max(F(i, 0), F(i, 1)))),
+      EAP.at(Eigen::RowVector2i(std::min(F(i, 1), F(i, 2)), std::max(F(i, 1), F(i, 2)))),
+      EAP.at(Eigen::RowVector2i(std::min(F(i, 2), F(i, 3)), std::max(F(i, 2), F(i, 3)))),
+      EAP.at(Eigen::RowVector2i(std::min(F(i, 3), F(i, 0)), std::max(F(i, 3), F(i, 0))))
+    };
+
+    SF.row(4 * i + 0) = Eigen::RowVector4i(
+      F(i, 0),
+      edge_average_points[0],
+      face_centroid,
+      edge_average_points[3]);
+    SF.row(4 * i + 1) = Eigen::RowVector4i(
+      edge_average_points[0],
+      F(i, 1),
+      edge_average_points[1],
+      face_centroid);
+    SF.row(4 * i + 2) = Eigen::RowVector4i(
+      face_centroid,
+      edge_average_points[1],
+      F(i, 2),
+      edge_average_points[2]);
+    SF.row(4 * i + 3) = Eigen::RowVector4i(
+      edge_average_points[3],
+      face_centroid,
+      edge_average_points[2],
+      F(i, 3));
+  }
 }