diff --git a/src/catmull_clark.cpp b/src/catmull_clark.cpp index de8ce26..9a1d9ac 100644 --- a/src/catmull_clark.cpp +++ b/src/catmull_clark.cpp @@ -2,6 +2,23 @@ #include #include #include +#include + +/** + * @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) { - //////////////////////////////////////////////////////////////////////////// - // Replace with your code here: - SV = V; - SF = F; - //////////////////////////////////////////////////////////////////////////// + Eigen::MatrixXd SVX = V; + Eigen::MatrixXi SFX = 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 n = 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); + n++; + } + } + } + + if (n == 0) { + std::cerr << "Error: A vertex should be shared by at least one face, but it is shared by " << n << " faces." << std::endl; + throw std::runtime_error("A vertex should be shared by at least one face."); + } + + /* Determine the edges that share the vertex: 2R */ + /* NOTE that the edge midpoint point is computed twice for each edge */ + int m = 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 + 3) % 4)) + V.row(F(i, j))) / 2.0; + part_R += (V.row(F(i, j)) + V.row(F(i, (j + 1) % 4))) / 2.0; + m++; + } + } + } + + if (n != m) { + std::cerr << "Error: A vertex should be shared by exactly two edges, but it is shared by " << m << " edges." << std::endl; + throw std::runtime_error("A vertex should be shared by exactly two edges."); + } + + /* (F + 2R + (n - 3)P) / n */ + return (part_F + part_R + (n - 3) * V.row(v)) / n; +} + +struct RowVector2iHash { + std::size_t operator()(const Eigen::RowVector2i & v) const { + return std::hash()(v[0]) ^ std::hash()(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 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)); + + /* 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, + SV.block(V.rows(), 0, F.rows(), 3), + 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, SV.block(V.rows(), 0, F.rows(), 3), 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)); + } }