diff --git a/src/catmull_clark.cpp b/src/catmull_clark.cpp index de8ce26..f46520d 100644 --- a/src/catmull_clark.cpp +++ b/src/catmull_clark.cpp @@ -3,6 +3,22 @@ #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_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()(v[0]) ^ std::hash()(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 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]; + } + } } diff --git a/src/cube.cpp b/src/cube.cpp index 5c98dc6..9462ff1 100644 --- a/src/cube.cpp +++ b/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); }