diff --git a/Solver_interface/include/CGAL/Eigen_solver_traits.h b/Solver_interface/include/CGAL/Eigen_solver_traits.h index e641c121b89..9e065c65558 100644 --- a/Solver_interface/include/CGAL/Eigen_solver_traits.h +++ b/Solver_interface/include/CGAL/Eigen_solver_traits.h @@ -53,6 +53,11 @@ namespace internal { typedef Eigen_sparse_symmetric_matrix type; }; + template + struct Get_eigen_matrix< ::Eigen::LeastSquaresConjugateGradient,FT>{ + typedef Eigen_sparse_matrix type; + }; + template struct Get_eigen_matrix< ::Eigen::SimplicialCholesky,FT>{ typedef Eigen_sparse_symmetric_matrix type; diff --git a/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h index ef820d9b813..47d42c98f34 100644 --- a/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h @@ -24,14 +24,22 @@ #include #include + +#ifdef CGAL_EIGEN3_ENABLED +#include +#endif + #include + #include #include #include #include #include +#define DEBUG_TRACE + /// \file LSCM_parameterizer_3.h namespace CGAL { @@ -136,8 +144,12 @@ private: typedef typename Sparse_LA::Vector Vector; typedef typename Sparse_LA::Matrix Matrix; - typedef typename OpenNL::LinearSolver - LeastSquaresSolver ; + typedef Eigen_solver_traits::EigenType> > + LeastSquaresSolver ; + + typedef typename LeastSquaresSolver::Matrix LSS_Matrix; + typedef typename LeastSquaresSolver::Vector LSS_Vector; + typedef typename LeastSquaresSolver::NT LSS_Scalar; // Public operations public: @@ -173,7 +185,8 @@ private: /// \pre Vertices must be indexed. /// \pre X and B must be allocated and empty. /// \pre At least 2 border vertices must be parameterized. - void initialize_system_from_mesh_border(LeastSquaresSolver& solver, + void initialize_system_from_mesh_border(LSS_Matrix& A, + LSS_Vector& B, const Adaptor& mesh) ; /// Utility for setup_triangle_relations(): @@ -185,7 +198,7 @@ private: /// Create two lines in the linear system per triangle (one for u, one for v). /// /// \pre vertices must be indexed. - Error_code setup_triangle_relations(LeastSquaresSolver& solver, + Error_code setup_triangle_relations(LSS_Matrix& A, LSS_Vector& B, const Adaptor& mesh, Facet_const_handle facet) ; @@ -288,34 +301,39 @@ parameterize(Adaptor& mesh) // Create sparse linear system "A*X = B" of size 2*nbVertices x 2*nbVertices // (in fact, we need only 2 lines per triangle x 1 column per vertex) - LeastSquaresSolver solver(2*nbVertices); - solver.set_least_squares(true) ; - + LSS_Matrix A (2 * nbVertices, 2 * nbVertices); + LSS_Vector B (2 * nbVertices); + // Initialize the "A*X = B" linear system after // (at least two) border vertices parameterization - initialize_system_from_mesh_border(solver, mesh); + initialize_system_from_mesh_border(A, B, mesh); // Fill the matrix for the other vertices - solver.begin_system() ; for (Facet_iterator facetIt = mesh.mesh_facets_begin(); facetIt != mesh.mesh_facets_end(); facetIt++) { // Create two lines in the linear system per triangle (one for u, one for v) - status = setup_triangle_relations(solver, mesh, facetIt); + status = setup_triangle_relations(A, B, mesh, facetIt); if (status != Base::OK) return status; } - solver.end_system() ; + + A.assemble_matrix(); + #ifdef DEBUG_TRACE std::cerr << " matrix filling (" << 2*mesh.count_mesh_facets() << " x " << nbVertices << "): " << timer.time() << " seconds." << std::endl; timer.reset(); #endif - + // Solve the "A*X = B" linear system in the least squares sense - if ( ! solver.solve() ) - status = Base::ERROR_CANNOT_SOLVE_LINEAR_SYSTEM; + LSS_Vector X (2 * nbVertices); + LSS_Scalar D; + LeastSquaresSolver solver; + if (!(solver.linear_solver (A, B, X, D))) + status = Base::ERROR_CANNOT_SOLVE_LINEAR_SYSTEM; + #ifdef DEBUG_TRACE std::cerr << " solving linear system: " << timer.time() << " seconds." << std::endl; @@ -323,7 +341,7 @@ parameterize(Adaptor& mesh) #endif if (status != Base::OK) return status; - + return status; // Copy X coordinates into the (u,v) pair of each vertex set_mesh_uv_from_system(mesh, solver); #ifdef DEBUG_TRACE @@ -396,7 +414,7 @@ check_parameterize_preconditions(Adaptor& mesh) template inline void LSCM_parameterizer_3:: -initialize_system_from_mesh_border(LeastSquaresSolver& solver, +initialize_system_from_mesh_border(LSS_Matrix& A, LSS_Vector& B, const Adaptor& mesh) { for (Vertex_const_iterator it = mesh.mesh_vertices_begin(); @@ -412,14 +430,13 @@ initialize_system_from_mesh_border(LeastSquaresSolver& solver, // Write (u,v) in X (meaningless if vertex is not parameterized) // Note : 2*index --> u // 2*index + 1 --> v - solver.variable(2*index ).set_value(uv.x()) ; - solver.variable(2*index + 1).set_value(uv.y()) ; // Copy (u,v) in B if vertex is parameterized - if (mesh.is_vertex_parameterized(it)) { - solver.variable(2*index ).lock() ; - solver.variable(2*index + 1).lock() ; - } + if (mesh.is_vertex_parameterized(it)) + { + B.set (2 * index, uv.x ()); + B.set (2 * index + 1, uv.y ()); + } } } @@ -474,7 +491,7 @@ template inline typename LSCM_parameterizer_3::Error_code LSCM_parameterizer_3:: -setup_triangle_relations(LeastSquaresSolver& solver, +setup_triangle_relations(LSS_Matrix& A, LSS_Vector& B, const Adaptor& mesh, Facet_const_handle facet) { @@ -537,24 +554,15 @@ setup_triangle_relations(LeastSquaresSolver& solver, // // Real part // Note: b = 0 - solver.begin_row() ; - solver.add_coefficient(u0_id, -a+c) ; - solver.add_coefficient(v0_id, b-d) ; - solver.add_coefficient(u1_id, -c) ; - solver.add_coefficient(v1_id, d) ; - solver.add_coefficient(u2_id, a) ; - solver.end_row() ; + + // TODO + // // Imaginary part // Note: b = 0 - solver.begin_row() ; - solver.add_coefficient(u0_id, -b+d) ; - solver.add_coefficient(v0_id, -a+c) ; - solver.add_coefficient(u1_id, -d) ; - solver.add_coefficient(v1_id, -c) ; - solver.add_coefficient(v2_id, a) ; - solver.end_row() ; + // TODO + return Base::OK; } @@ -574,8 +582,10 @@ set_mesh_uv_from_system(Adaptor& mesh, // Note : 2*index --> u // 2*index + 1 --> v - NT u = solver.variable(2*index ).value() ; - NT v = solver.variable(2*index + 1).value() ; + + double u, v; + // TODO + // Fill vertex (u,v) and mark it as "parameterized" mesh.set_vertex_uv(vertexIt, Point_2(u,v));