From 8656f1edb3d15dd01498d0e15b872ba391e001a4 Mon Sep 17 00:00:00 2001 From: Tibo De Peuter Date: Tue, 12 Nov 2024 11:50:34 +0100 Subject: [PATCH 1/2] Checkpoint --- src/catmull_clark.cpp | 219 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 214 insertions(+), 5 deletions(-) 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)); + } } From a3487cc05facb06dfb453c86b9a6bf2b5a31af3e Mon Sep 17 00:00:00 2001 From: Tibo De Peuter Date: Tue, 12 Nov 2024 12:11:43 +0100 Subject: [PATCH 2/2] Fixed it --- src/catmull_clark.cpp | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/catmull_clark.cpp b/src/catmull_clark.cpp index 9a1d9ac..4942c74 100644 --- a/src/catmull_clark.cpp +++ b/src/catmull_clark.cpp @@ -112,41 +112,37 @@ Eigen::RowVector3d compute_vertex_moved_position( Eigen::RowVector3d part_R = Eigen::RowVector3d::Zero(); /* Determine the faces that share the vertex: F */ - int n = 0; + 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); - n++; + f++; } } } + part_F /= f; - 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; + 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 */ - 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++; + r++; } } } - - 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."); - } + part_R /= r; /* (F + 2R + (n - 3)P) / n */ - return (part_F + part_R + (n - 3) * V.row(v)) / n; + return (part_F + 2 * part_R + (f - 3) * V.row(v)) / f; } struct RowVector2iHash { @@ -169,7 +165,11 @@ void subdivide_quadmesh_catmullclark( 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)); @@ -180,7 +180,7 @@ void subdivide_quadmesh_catmullclark( SV.row(V.rows() + F.rows() + EAP.size()) = compute_edge_average_point( V, F, - SV.block(V.rows(), 0, F.rows(), 3), + FC, Eigen::RowVector2i(v1, v2)); EAP[Eigen::RowVector2i(v1, v2)] = V.rows() + F.rows() + EAP.size(); } @@ -190,7 +190,7 @@ void subdivide_quadmesh_catmullclark( /* 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); + SV.row(i) = compute_vertex_moved_position(V, F, FC, i); } /* Compute the faces */