Checkpoint

This commit is contained in:
Tibo De Peuter 2024-11-12 11:50:34 +01:00
parent 98c06721d3
commit 8656f1edb3
Signed by: tdpeuter
GPG key ID: 38297DE43F75FFE2

View file

@ -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)
{
////////////////////////////////////////////////////////////////////////////
// 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<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));
/* 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));
}
}