diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/LSCM_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/LSCM_parameterizer_3.h index 2dcaab67b80..dae71439143 100644 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/LSCM_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/LSCM_parameterizer_3.h @@ -29,7 +29,6 @@ #ifdef CGAL_EIGEN3_ENABLED #include #endif -#include #include @@ -100,9 +99,6 @@ public: typedef typename Default::Get< SolverTraits_, #if defined(CGAL_EIGEN3_ENABLED) - // WARNING: at the moment, the choice of SolverTraits_ is completely - // ignored (see LeastSquaresSolver typedef) and `OpenNL::LinearSolver` - // is always used... CGAL::Eigen_solver_traits< Eigen::SimplicialLDLT::EigenType> > #else @@ -143,7 +139,13 @@ private: typedef typename Solver_traits::Vector Vector; typedef typename Solver_traits::Matrix Matrix; - typedef typename OpenNL::LinearSolver LeastSquaresSolver; +// Fields +private: + // %Object that maps (at least two) border vertices onto a 2D space + Border_parameterizer m_borderParameterizer; + + // Traits object to solve a sparse linear system + Solver_traits m_linearAlgebra; // Public operations public: @@ -157,6 +159,15 @@ public: // Default copy constructor and operator =() are fine +// Private accessors +private: + // Get the object that maps the surface's border onto a 2D space. + Border_parameterizer& get_border_parameterizer() { return m_borderParameterizer; } + + // Get the sparse linear algebra (traits object to access the linear system). + Solver_traits& get_linear_algebra_traits() { return m_linearAlgebra; } + +public: /// returns whether the 3D -> 2D mapping is one-to-one. template bool is_one_to_one_mapping(const Triangle_mesh& mesh, @@ -220,88 +231,59 @@ public: // Compute (u,v) for (at least two) border vertices // and mark them as "parameterized" - Error_code status = - get_border_parameterizer().parameterize(mesh, bhd, uvmap, vimap, vpmap); + Error_code status = get_border_parameterizer().parameterize(mesh, bhd, uvmap, vimap, vpmap); if(status != OK) return status; - // 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) ; + // Vertices that are already parameterized are not part of the system + std::vector fvi(nbVertices); // vimap to free vertex index + int index = 0; + for(vertex_descriptor vd : ccvertices) { + if(get(vpmap,vd)) { + fvi[get(vimap, vd)] = -1; + } else { + fvi[get(vimap, vd)] = index++; + } + } - // Initialize the "A*X = B" linear system after - // (at least two) border vertices parameterization - initialize_system_from_mesh_border(solver, ccvertices, uvmap, vimap, vpmap); - - // Fill the matrix for the other vertices - solver.begin_system(); + // Create sparse linear system "A*X = B" + int nbFreeVertices = index; + int nbVariables = 2 * nbFreeVertices; + Matrix A(nbVariables, nbVariables); // the constant matrix used in the linear system A*X = B + Vector X(nbVariables), B(nbVariables); for(face_descriptor fd : ccfaces) { // Create two lines in the linear system per triangle (one for u, one for v) - status = setup_triangle_relations(solver, mesh, fd, vimap); + status = setup_triangle_relations(mesh, fd, uvmap, vimap, vpmap, fvi, A, B); if (status != OK) return status; } - solver.end_system(); - // Solve the "A*X = B" linear system in the least squares sense - if(!solver.solve()) + double D; + if(!get_linear_algebra_traits().linear_solver(A, B, X, D)) { + std::cerr << "Could not solve linear system" << std::endl; status = ERROR_CANNOT_SOLVE_LINEAR_SYSTEM; - - if(status != OK) return status; + } // Copy X coordinates into the (u,v) pair of each vertex - //set_mesh_uv_from_system(mesh, solver, uvmap); - for(vertex_descriptor vd : ccvertices) { + if(get(vpmap,vd)) + continue; + int index = get(vimap,vd); - NT u = solver.variable(2 * index).value(); - NT v = solver.variable(2 * index + 1).value(); + NT u = X[2 * fvi[index]]; + NT v = X[2 * fvi[index] + 1]; put(uvmap, vd, Point_2(u,v)); } + return status; } // Private operations private: - // Initialize "A*X = B" linear system after - // (at least two) border vertices are parameterized. - // - // \pre Vertices must be indexed. - // \pre X and B must be allocated and empty. - // \pre At least 2 border vertices must be parameterized. - template - void initialize_system_from_mesh_border(LeastSquaresSolver& solver, - const std::unordered_set& ccvertices, - UVmap uvmap, - VertexIndexMap vimap, - VertexParameterizedMap vpmap) const - { - for(vertex_descriptor v : ccvertices) { - // Get vertex index in sparse linear system - int index = get(vimap, v); - - // Get vertex (u,v) (meaningless if vertex is not parameterized) - const Point_2& uv = get(uvmap, v); - // TODO: it is meaningless but must it be called for non-border vertices?? - // 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 (get(vpmap,v)) { - solver.variable(2 * index ).lock(); - solver.variable(2 * index + 1).lock(); - } - } - } - // Utility for setup_triangle_relations(): // Computes the coordinates of the vertices of a triangle // in a local 2D orthonormal basis of the triangle's plane. @@ -344,11 +326,18 @@ private: // Zk = xk + i.yk is the complex number corresponding to local (x,y) coords // cool: no divide with this expression; makes it more numerically stable // in presence of degenerate triangles - template - Error_code setup_triangle_relations(LeastSquaresSolver& solver, - const Triangle_mesh& mesh, + template + Error_code setup_triangle_relations(const Triangle_mesh& mesh, face_descriptor facet, - VertexIndexMap vimap) const + UVmap uvmap, + VertexIndexMap vimap, + VertexParameterizedMap vpmap, + const FreeVertexIndices& fvi, + Matrix& A, + Vector& B) const { const PPM ppmap = get(vertex_point, mesh); @@ -361,11 +350,6 @@ private: halfedge_descriptor h2 = next(h1, mesh); v2 = target(h2, mesh); - // Get the vertices index - int id0 = get(vimap, v0); - int id1 = get(vimap, v1); - int id2 = get(vimap, v2); - // Get the vertices position const Point_3& p0 = get(ppmap, v0); const Point_3& p1 = get(ppmap, v1); @@ -392,51 +376,75 @@ private: // // Note : 2*index --> u // 2*index + 1 --> v - int u0_id = 2*id0 ; - int v0_id = 2*id0 + 1; - int u1_id = 2*id1 ; - int v1_id = 2*id1 + 1; - int u2_id = 2*id2 ; - int v2_id = 2*id2 + 1; + + std::vector vals; + std::vector ids; + std::vector fvals; + std::vector pvals; + + auto add_coefficient = [&](const vertex_descriptor v, const double a, const bool is_u) + { + if(get(vpmap, v)) { + const Point_2& uv = get(uvmap, v); + fvals.push_back(a); + pvals.push_back(is_u ? uv.x() : uv.y()); + } else { + const int index = 2 * fvi[get(vimap, v)] + !is_u; + CGAL_assertion(index >= 0); + vals.push_back(a); + ids.push_back(index); + } + }; + + auto add_u_coefficient = [&](const vertex_descriptor v, const double a) { + return add_coefficient(v, a, true /*is u*/); + }; + auto add_v_coefficient = [&](const vertex_descriptor v, const double a) { + return add_coefficient(v, a, false /*is not u*/); + }; + + auto accumulate = [&]() + { + const std::size_t nf = vals.size(); + const std::size_t nl = fvals.size(); + + for(std::size_t i=0; i