From a3487cc05facb06dfb453c86b9a6bf2b5a31af3e Mon Sep 17 00:00:00 2001
From: Tibo De Peuter <tibo.depeuter@ugent.be>
Date: Tue, 12 Nov 2024 12:11:43 +0100
Subject: [PATCH] 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 */