Use eigen wrapper

This commit is contained in:
Julian Stahl 2024-01-02 19:46:12 +01:00
parent 86ac3311e8
commit ea95b6c3b6
1 changed files with 24 additions and 13 deletions

View File

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