From ea95b6c3b64a25bfdfdcb592e34d7e30b645e96e Mon Sep 17 00:00:00 2001
From: Julian Stahl <julianstahl.dev@gmail.com>
Date: Tue, 2 Jan 2024 19:46:12 +0100
Subject: [PATCH] Use eigen wrapper

---
 .../internal/dual_contouring_functors.h       | 37 ++++++++++++-------
 1 file changed, 24 insertions(+), 13 deletions(-)

diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/dual_contouring_functors.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/dual_contouring_functors.h
index 90c36f17278..fbfbd9e3ac1 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/dual_contouring_functors.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/dual_contouring_functors.h
@@ -19,7 +19,12 @@
 #include <CGAL/centroid.h>
 #include <CGAL/Origin.h>
 
+#ifdef CGAL_EIGEN3_ENABLED
+#include <CGAL/Eigen_matrix.h>
+#include <CGAL/Eigen_vector.h>
+
 #include <Eigen/SVD>
+#endif
 
 #include <array>
 #include <map>
@@ -66,10 +71,10 @@ public:
     using Point_3 = typename Geom_traits::Point_3;
     using Vector_3 = typename Geom_traits::Vector_3;
 
-    using Eigen_vector_3 = Eigen::Vector<FT, 3>;
-    using Eigen_vector_x = Eigen::Vector<FT, Eigen::Dynamic>;
-    using Eigen_matrix_3 = Eigen::Matrix<FT, 3, 3>;
-    using Eigen_matrix_x = Eigen::Matrix<FT, Eigen::Dynamic, Eigen::Dynamic>;
+    using Eigen_vector_3 = Eigen_vector<FT, 3>;
+    using Eigen_vector_x = Eigen_vector<FT>;
+    using Eigen_matrix_3 = Eigen_matrix<FT, 3, 3>;
+    using Eigen_matrix_x = Eigen_matrix<FT>;
 
     using Vertex_descriptor = typename Domain::Vertex_descriptor;
 
@@ -127,21 +132,26 @@ public:
     rhs.setZero();
     for(std::size_t i=0; i<edge_intersections.size(); ++i)
     {
-      Eigen_vector_3 n_k{ x_coord(edge_intersection_normals[i]),
-                            y_coord(edge_intersection_normals[i]),
-                            z_coord(edge_intersection_normals[i]) };
-      Eigen_vector_3 p_k{ x_coord(edge_intersections[i]),
-                            y_coord(edge_intersections[i]),
-                            z_coord(edge_intersections[i]) };
+      Eigen_vector_3 n_k;
+      n_k.set(0, x_coord(edge_intersection_normals[i]));
+      n_k.set(1, y_coord(edge_intersection_normals[i]));
+      n_k.set(2, z_coord(edge_intersection_normals[i]));
+
+      Eigen_vector_3 p_k;
+      p_k.set(0, x_coord(edge_intersections[i]));
+      p_k.set(1, y_coord(edge_intersections[i]));
+      p_k.set(2, z_coord(edge_intersections[i]));
+
       const FT d_k = n_k.transpose() * p_k;
 
       Eigen_matrix_3 A_k = n_k * n_k.transpose();
-      Eigen_vector_3 b_k = d_k * n_k;
+      Eigen_vector_3 b_k;
+      b_k = d_k * n_k;
       A += A_k;
       rhs += b_k;
     }
 
-    Eigen::JacobiSVD<Eigen_matrix_x> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
+    Eigen::JacobiSVD<Eigen_matrix_x::EigenType> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
 
     // set threshold as in Peter Lindstrom's paper, "Out-of-Core Simplification of Large Polygonal Models"
     svd.setThreshold(1e-3);
@@ -151,7 +161,8 @@ public:
     x_hat << x_coord(p), y_coord(p), z_coord(p);
 
     // Lindstrom formula for QEM new position for singular matrices
-    Eigen_vector_x v_svd = x_hat + svd.solve(rhs - A * x_hat);
+    Eigen_vector_x v_svd;
+    v_svd = x_hat + svd.solve(rhs - A * x_hat);
     p = point(v_svd[0], v_svd[1], v_svd[2]);
 
     // bounding box