chore: Cleanup cube

src/catmull_clark.cpp

@@ -3,6 +3,22 @@
 #include <utility>
 #include <functional>
 
+/**
+ * @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_quad_mesh_catmull_clark(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & F,
+        Eigen::MatrixXd & SV,
+        Eigen::MatrixXi & SF);
+
 void catmull_clark(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
@@ -10,9 +26,172 @@ void catmull_clark(
   Eigen::MatrixXd & SV,
   Eigen::MatrixXi & SF)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Replace with your code here:
-  SV = V;
-  SF = F;
-  ////////////////////////////////////////////////////////////////////////////
+  Eigen::MatrixXd SVX = V; /* Copy of const because we need to substitute */
+  Eigen::MatrixXi SFX = F;
+
+  for (int i = 0; i < num_iters; i++) {
+    subdivide_quad_mesh_catmull_clark(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 */
+  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);
+      }
+    }
+  }
+
+  /* (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 FC #F by 3 list of face centroids
+ * @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;
+
+  /* Determine the edges that share the vertex: 2R */
+  int r = 0;
+  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_quad_mesh_catmull_clark(
+  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
+
+  /* Compute face centroids: FC */
+  for (int i = 0; i < F.rows(); i++) {
+    SV.row(V.rows() + i) = compute_face_centroid(V, F.row(i));
+  }
+
+  Eigen::MatrixXd FC = SV.block(V.rows(), 0, F.rows(), 3);
+
+  /* Compute edge average points */
+  for (int i = 0; i < F.rows(); i++) {
+    for (int j = 0; j < F.cols(); j++) {
+      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))))
+    };
+
+    for (int j = 0; j < 4; j++) {
+      SF(4 * i + j, j) = F(i, j);
+      SF(4 * i + j, (j + 1) % 4) = edge_average_points[j];
+      SF(4 * i + j, (j + 2) % 4) = face_centroid;
+      SF(4 * i + j, (j + 3) % 4) = edge_average_points[(j + 3) % 4];
+    }
+  }
 }

src/cube.cpp

@@ -11,20 +11,6 @@ void cube(
   // Construct default Blender cube.
 
   V.resize(8,3);
-  F.resize(6,4);
-//  for (int i = 0; i < 8; i++) {
-//    V(i,0) =  1 - (i&4) / 2;
-//    V(i,1) =  1 - 2 * (i&1);
-//    V(i,2) = -1 + 1 * (i&2);
-
-//    V.row(i) = Eigen::RowVector3d(1 - (i&4) / 2, 1 - 2 * (i&1), -1 + 1 * (i&2));
-
-    /* Add Vertex to the faces it belongs to */
-//    F(0 + (i&4) / 2, i % 4) = i;
-//    F(1 + 2 * (i&1), i / 2) = i;
-//    F(5 - (i&2) / 2, (i % 4) % 2 + (i / 4) * 2) = i;
-//  }
-
   V.row(0) = Eigen::RowVector3d( 1,  1, -1);
   V.row(1) = Eigen::RowVector3d( 1, -1, -1);
   V.row(2) = Eigen::RowVector3d( 1,  1,  1);
@@ -34,20 +20,21 @@ void cube(
   V.row(6) = Eigen::RowVector3d(-1,  1,  1);
   V.row(7) = Eigen::RowVector3d(-1, -1,  1);
 
-  F.row(0) = Eigen::RowVector4i(1, 2, 4, 3);
-  F.row(1) = Eigen::RowVector4i(5, 1, 3, 7);
-  F.row(2) = Eigen::RowVector4i(6, 5, 7, 8);
-  F.row(3) = Eigen::RowVector4i(2, 6, 8, 4);
-  F.row(4) = Eigen::RowVector4i(3, 4, 8, 7);
-  F.row(5) = Eigen::RowVector4i(5, 6, 2, 1);
+  F.resize(6,4);
+  F.row(0) = Eigen::RowVector4i(1, 3, 4, 2);
+  F.row(1) = Eigen::RowVector4i(5, 7, 3, 1);
+  F.row(2) = Eigen::RowVector4i(6, 8, 7, 5);
+  F.row(3) = Eigen::RowVector4i(2, 4, 8, 6);
+  F.row(4) = Eigen::RowVector4i(3, 7, 8, 4);
+  F.row(5) = Eigen::RowVector4i(5, 1, 2, 6);
 
   NV.resize(6,3);
-  NV.row(0) = Eigen::RowVector3d(1, 0, 0);
-  NV.row(1) = Eigen::RowVector3d(0, 1, 0);
-  NV.row(2) = Eigen::RowVector3d(-1, 0, 0);
-  NV.row(3) = Eigen::RowVector3d(0, -1, 0);
-  NV.row(4) = Eigen::RowVector3d(0, 0, 1);
-  NV.row(5) = Eigen::RowVector3d(0, 0, -1);
+  NV.row(0) = Eigen::RowVector3d( 1,  0,  0);
+  NV.row(1) = Eigen::RowVector3d( 0,  1,  0);
+  NV.row(2) = Eigen::RowVector3d(-1,  0,  0);
+  NV.row(3) = Eigen::RowVector3d( 0, -1,  0);
+  NV.row(4) = Eigen::RowVector3d( 0,  0,  1);
+  NV.row(5) = Eigen::RowVector3d( 0,  0, -1);
 
   NF.resize(6,4);
   for (int i = 1; i <= 6; i++) {
@@ -56,19 +43,19 @@ void cube(
 
   UV.resize(14,2);
   UV.row(0) = Eigen::RowVector2d(0.25, 0);
-  UV.row(1) = Eigen::RowVector2d(0.5, 0);
+  UV.row(1) = Eigen::RowVector2d(0.50, 0);
   for (int i = 0; i < 5; i++) {
     UV.row(2 + i) = Eigen::RowVector2d(0.25 * i, 0.25);
-    UV.row(7 + i) = Eigen::RowVector2d(0.25 * i, 0.5);
+    UV.row(7 + i) = Eigen::RowVector2d(0.25 * i, 0.50);
   }
   UV.row(12) = Eigen::RowVector2d(0.25, 0.75);
-  UV.row(13) = Eigen::RowVector2d(0.5, 0.75);
+  UV.row(13) = Eigen::RowVector2d(0.50, 0.75);
 
   UF.resize(6,4);
-  UF.row(0) = Eigen::RowVector4i(4, 5, 10, 9);
-  UF.row(1) = Eigen::RowVector4i(5, 6, 11, 10);
-  UF.row(2) = Eigen::RowVector4i(6, 7, 12, 11);
-  UF.row(3) = Eigen::RowVector4i(3, 4, 9, 8);
-  UF.row(4) = Eigen::RowVector4i(9, 10, 14, 13);
-  UF.row(5) = Eigen::RowVector4i(1, 2, 5, 4);
+  UF.row(0) = Eigen::RowVector4i( 4,  5, 10,  9);
+  UF.row(1) = Eigen::RowVector4i( 5,  6, 11, 10);
+  UF.row(2) = Eigen::RowVector4i( 6,  7, 12, 11);
+  UF.row(3) = Eigen::RowVector4i( 3,  4,  9,  8);
+  UF.row(4) = Eigen::RowVector4i( 9, 10, 14, 13);
+  UF.row(5) = Eigen::RowVector4i( 1,  2,  5,  4);
 }