diff --git a/Polyhedron/demo/Polyhedron/CMakeLists.txt b/Polyhedron/demo/Polyhedron/CMakeLists.txt index 0ae3f7a2a2a..f30aaa149d6 100644 --- a/Polyhedron/demo/Polyhedron/CMakeLists.txt +++ b/Polyhedron/demo/Polyhedron/CMakeLists.txt @@ -57,6 +57,12 @@ if(Qt5_FOUND) endif(Qt5_FOUND) +find_package(OpenMP) +if (OPENMP_FOUND) + set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") +endif() + find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater) if (EIGEN3_FOUND) include( ${EIGEN3_USE_FILE} ) diff --git a/Polyhedron/demo/Polyhedron/MainWindow.ui b/Polyhedron/demo/Polyhedron/MainWindow.ui index d19a1c77d0c..9ae0a1f1b4e 100644 --- a/Polyhedron/demo/Polyhedron/MainWindow.ui +++ b/Polyhedron/demo/Polyhedron/MainWindow.ui @@ -37,7 +37,7 @@ 0 0 978 - 21 + 26 @@ -95,6 +95,7 @@ + @@ -652,6 +653,11 @@ &Preferences + + + Least square conformal maps + + diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_parameterization_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_parameterization_plugin.cpp index 495fa8d28e9..b1fb77bc71c 100644 --- a/Polyhedron/demo/Polyhedron/Polyhedron_demo_parameterization_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_parameterization_plugin.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include @@ -35,7 +36,8 @@ public: // used by Polyhedron_demo_plugin_helper QStringList actionsNames() const { return QStringList() << "actionMVC" - << "actionDCP"; + << "actionDCP" + << "actionLSC"; } bool applicable(QAction*) const { @@ -45,9 +47,10 @@ public: public Q_SLOTS: void on_actionMVC_triggered(); void on_actionDCP_triggered(); + void on_actionLSC_triggered(); protected: - enum Parameterization_method { PARAM_MVC, PARAM_DCP }; + enum Parameterization_method { PARAM_MVC, PARAM_DCP, PARAM_LSC }; void parameterize(Parameterization_method method); }; // end Polyhedron_demo_parameterization_plugin @@ -91,6 +94,13 @@ void Polyhedron_demo_parameterization_plugin::parameterize(const Parameterizatio typedef CGAL::Discrete_conformal_map_parameterizer_3 Parameterizer; Parameterizer::Error_code err = CGAL::parameterize(adaptor,Parameterizer()); success = err == Parameterizer::OK; + } + case PARAM_LSC: + { + std::cout << "Parameterize (LSC)..."; + typedef CGAL::LSCM_parameterizer_3 Parameterizer; + Parameterizer::Error_code err = CGAL::parameterize(adaptor,Parameterizer()); + success = err == Parameterizer::OK; } } @@ -148,4 +158,10 @@ void Polyhedron_demo_parameterization_plugin::on_actionDCP_triggered() parameterize(PARAM_DCP); } +void Polyhedron_demo_parameterization_plugin::on_actionLSC_triggered() +{ + std::cerr << "LSC..."; + parameterize(PARAM_LSC); +} + #include "Polyhedron_demo_parameterization_plugin.moc" diff --git a/Solver_interface/include/CGAL/Eigen_least_squares_solver.h b/Solver_interface/include/CGAL/Eigen_least_squares_solver.h new file mode 100644 index 00000000000..a446c551699 --- /dev/null +++ b/Solver_interface/include/CGAL/Eigen_least_squares_solver.h @@ -0,0 +1,335 @@ +// Copyright (c) 2005-2008 Inria Loria (France). +/* + * author: Bruno Levy, INRIA, project ALICE + * website: http://www.loria.fr/~levy/software + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation, either version 3 + * of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + * + * Scientific work that use this software can reference the website and + * the following publication: + * + * @INPROCEEDINGS {levy:NMDGP:05, + * AUTHOR = Bruno Levy, + * TITLE = Numerical Methods for Digital Geometry Processing, + * BOOKTITLE =Israel Korea Bi-National Conference, + * YEAR=November 2005, + * URL=http://www.loria.fr/~levy/php/article.php?pub=../publications/papers/2005/Numerics + * } + * + * Laurent Saboret 2005-2006: Changes for CGAL: + * - Added OpenNL namespace + * - DefaultLinearSolverTraits is now a model of the SparseLinearAlgebraTraits_d concept + * - Added SymmetricLinearSolverTraits + * - copied Jacobi preconditioner from Graphite 1.9 code + */ + + +#ifndef __EIGEN_LEAST_SQUARES_SOLVER__ +#define __EIGEN_LEAST_SQUARES_SOLVER__ + +#include + +#include +#include +#include + + +namespace CGAL { + + +/* + * Solves a linear system or minimizes a quadratic form. + */ +class EigenLeastSquaresSolver +{ +protected: + + enum State { + INITIAL, IN_SYSTEM, IN_ROW, CONSTRUCTED, SOLVED + } ; + +public: + // typedef Eigen_solver_traits::EigenType, + // Eigen::Lower|Eigen::Upper> > Solver; + typedef Eigen_solver_traits::EigenType> > Solver; + typedef typename Solver::Matrix Matrix ; + typedef typename Solver::Vector Vector ; + typedef typename Solver::NT CoeffType ; + + class Variable { + public: + Variable() : x_(0), index_(-1), locked_(false) { } + double value() const { return x_; } + void set_value(double x_in) { x_ = x_in ; } + void lock() { locked_ = true ; } + void unlock() { locked_ = false ; } + bool is_locked() const { return locked_ ; } + unsigned int index() const { + CGAL_assertion(index_ != -1) ; + return (unsigned int)(index_) ; + } + void set_index(unsigned int index_in) { + index_ = index_in ; + } + private: + CoeffType x_ ; + int index_ ; + bool locked_ ; + } ; + + + EigenLeastSquaresSolver(unsigned int nb_variables) { + state_ = INITIAL ; + nb_variables_ = nb_variables ; + variable_ = new Variable[nb_variables] ; + A_ = NULL ; + x_ = NULL ; + b_ = NULL ; + } + + ~EigenLeastSquaresSolver() { + delete[] variable_ ; + delete A_ ; + delete x_ ; + delete b_ ; + } + + // __________________ Parameters ________________________ + + void set_least_squares(bool x) { } + + // __________________ Access ____________________________ + + int nb_variables() const { return nb_variables_ ; } + + Variable& variable(unsigned int idx) { + CGAL_assertion(idx < nb_variables_) ; + return variable_[idx] ; + } + + const Variable& variable(unsigned int idx) const { + CGAL_assertion(idx < nb_variables_) ; + return variable_[idx] ; + } + + // _________________ Construction _______________________ + + void begin_system() { + current_row_ = 0 ; + transition(INITIAL, IN_SYSTEM) ; + // Enumerate free variables. + unsigned int index = 0 ; + for(int ii=0; ii < nb_variables() ; ii++) { + Variable& v = variable(ii) ; + if(!v.is_locked()) { + v.set_index(index) ; + index++ ; + } + } + unsigned int n = index ; + A_ = new Matrix(static_cast(n),static_cast(n)) ; + x_ = new Vector(n) ; + b_ = new Vector(n) ; + for(unsigned int i=0; iset(i, 0.); + b_->set(i, 0.); + } + variables_to_vector() ; + } + + void begin_row() { + transition(IN_SYSTEM, IN_ROW) ; + af_.clear() ; + if_.clear() ; + al_.clear() ; + xl_.clear() ; + bk_ = 0 ; + } + + void set_right_hand_side(double b) { + check_state(IN_ROW) ; + bk_ = b ; + } + + void add_coefficient(unsigned int iv, double a) { + check_state(IN_ROW) ; + Variable& v = variable(iv) ; + if(v.is_locked()) { + al_.push_back(a) ; + xl_.push_back(v.value()) ; + } else { + af_.push_back(a) ; + if_.push_back(v.index()) ; + } + } + + void normalize_row(CoeffType weight = 1) { + check_state(IN_ROW) ; + CoeffType norm = 0.0 ; + unsigned int nf = af_.size() ; + for(unsigned int i=0; i1e-40 ); + scale_row(weight / norm) ; + } + + void scale_row(CoeffType s) { + check_state(IN_ROW) ; + unsigned int nf = af_.size() ; + for(unsigned int i=0; iadd_coef (if_[i], if_[j], af_[i] * af_[j]); + } + } + CoeffType S = - bk_ ; + for(unsigned int j=0; jset(if_[i], b_->get(if_[i]) - af_[i] * S); + } + current_row_++ ; + transition(IN_ROW, IN_SYSTEM) ; + } + + void end_system() { + A_->assemble_matrix (); + transition(IN_SYSTEM, CONSTRUCTED) ; + } + + // ----------------------------- Solver ------------------------------- + + // Solves a linear system or minimizes a quadratic form. + // Return true on success. + // (modified for SparseLinearAlgebraTraits_d concept) + bool solve() + { + check_state(CONSTRUCTED) ; + + // Solve the sparse linear system "A*X = B". On success, the solution is (1/D) * X. + Solver solver_; + CoeffType D; + + bool success = solver_.linear_solver(*A_, *b_, *x_, D) ; + CGAL_assertion(D == 1.0); // WARNING: this library does not support homogeneous coordinates! + + vector_to_variables() ; + + transition(CONSTRUCTED, SOLVED) ; + + delete A_ ; A_ = NULL ; + delete b_ ; b_ = NULL ; + delete x_ ; x_ = NULL ; + + return success; + } + +protected: + + // ----------- Converting between user representation and the internal representation ----- + + void vector_to_variables() { + for(int ii=0; ii < nb_variables(); ii++) { + Variable& v = variable(ii) ; + if(!v.is_locked()) { + v.set_value( x_->get(v.index()) ) ; + } + } + } + + void variables_to_vector() { + for(int ii=0; ii < nb_variables(); ii++) { + Variable& v = variable(ii) ; + if(!v.is_locked()) { + x_->set(v.index(), v.value()); + } + } + } + + // ----------- Finite state automaton (checks that calling sequence is respected) --------- + + std::string state_to_string(State s) { + switch(s) { + case INITIAL: + return "initial" ; + case IN_SYSTEM: + return "in system" ; + case IN_ROW: + return "in row" ; + case CONSTRUCTED: + return "constructed" ; + case SOLVED: + return "solved" ; + } + // Should not go there. + CGAL_error(); + return "undefined" ; + } + + void check_state(State s) { + CGAL_USE(s); + CGAL_assertion(state_ == s) ; + } + + void transition(State from, State to) { + check_state(from) ; + state_ = to ; + } + +private: + + // --------------- user representation -------------- + unsigned int nb_variables_ ; + Variable* variable_ ; + + // --------------- construction ----------------------- + State state_ ; + unsigned int current_row_ ; + std::vector af_ ; + std::vector if_ ; + std::vector al_ ; + std::vector xl_ ; + double bk_ ; + + // --------------- internal representation --------- + Matrix* A_ ; + Vector* x_ ; + Vector* b_ ; + +} ; + + +} // namespace OpenNL + +#endif diff --git a/Solver_interface/include/CGAL/Eigen_solver_traits.h b/Solver_interface/include/CGAL/Eigen_solver_traits.h index e641c121b89..e2d11bd09df 100644 --- a/Solver_interface/include/CGAL/Eigen_solver_traits.h +++ b/Solver_interface/include/CGAL/Eigen_solver_traits.h @@ -49,10 +49,20 @@ namespace internal { }; template - struct Get_eigen_matrix< ::Eigen::ConjugateGradient,FT>{ + struct Get_eigen_matrix< ::Eigen::ConjugateGradient,FT>{ typedef Eigen_sparse_symmetric_matrix type; }; + template + struct Get_eigen_matrix< ::Eigen::SimplicialLDLT,FT>{ + 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; @@ -89,6 +99,7 @@ public: Eigen_solver_traits():m_mat(NULL), m_solver_sptr(new EigenSolverT) { + } EigenSolverT& solver() { return *m_solver_sptr; } @@ -110,6 +121,13 @@ public: X = m_solver_sptr->solve(B); + if (m_solver_sptr->info() == Eigen::NumericalIssue) + std::cerr << "Numerical Issue" << std::endl; + if (m_solver_sptr->info() == Eigen::NoConvergence) + std::cerr << "No Convergence" << std::endl; + if (m_solver_sptr->info() == Eigen::InvalidInput) + std::cerr << "Invalid Input" << std::endl; + return m_solver_sptr->info() == Eigen::Success; } diff --git a/Solver_interface/include/CGAL/Eigen_vector.h b/Solver_interface/include/CGAL/Eigen_vector.h index 6280e289be1..c8acdfe74b5 100644 --- a/Solver_interface/include/CGAL/Eigen_vector.h +++ b/Solver_interface/include/CGAL/Eigen_vector.h @@ -86,6 +86,10 @@ public: void set(std::size_t i,NT value) { this->operator[](static_cast(i))=value; } + + NT get(std::size_t i) const { + return this->operator[](static_cast(i)); + } NT* vector() {return this->data();} diff --git a/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Eigen_parameterization.cpp b/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Eigen_parameterization.cpp index 7c97b3f8823..0d0300fbee9 100644 --- a/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Eigen_parameterization.cpp +++ b/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Eigen_parameterization.cpp @@ -3,7 +3,7 @@ #include #include #include - +#include #include #include @@ -17,7 +17,7 @@ typedef CGAL::Simple_cartesian Kernel; typedef CGAL::Polyhedron_3 Polyhedron; - +typedef CGAL::Timer Timer; // ---------------------------------------------------------------------------- // main() @@ -65,7 +65,8 @@ int main(int argc, char * argv[]) typedef CGAL::Parameterization_polyhedron_adaptor_3 Parameterization_polyhedron_adaptor; - + Timer t; + t.start(); Parameterization_polyhedron_adaptor mesh_adaptor(mesh); //*************************************** @@ -77,16 +78,17 @@ int main(int argc, char * argv[]) typedef CGAL::Circular_border_arc_length_parameterizer_3 Border_parameterizer; // Eigen solver - typedef CGAL::Eigen_solver_traits<> Solver; + typedef CGAL::Eigen_solver_traits::EigenType, Eigen::IncompleteLUT< double > > > Solver; // Floater Mean Value Coordinates parameterization // (circular border) with Eigen solver typedef CGAL::Mean_value_coordinates_parameterizer_3 + Solver> Parameterizer; Parameterizer::Error_code err = CGAL::parameterize(mesh_adaptor, Parameterizer()); + t.stop(); switch(err) { case Parameterizer::OK: // Success @@ -119,6 +121,6 @@ int main(int argc, char * argv[]) double v = mesh_adaptor.info(pVertex->halfedge())->uv().y(); std::cout << "(u,v) = (" << u << "," << v << ")" << std::endl; } - + std::cerr << t.time() << "sec." << std::endl; return EXIT_SUCCESS; } diff --git a/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h index ef820d9b813..b05fde07c7f 100644 --- a/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h @@ -32,6 +32,12 @@ #include #include +#ifdef CGAL_EIGEN3_ENABLED +#include +#endif + +#define DEBUG_TRACE + /// \file LSCM_parameterizer_3.h namespace CGAL { @@ -136,8 +142,7 @@ private: typedef typename Sparse_LA::Vector Vector; typedef typename Sparse_LA::Matrix Matrix; - typedef typename OpenNL::LinearSolver - LeastSquaresSolver ; + typedef CGAL::internal::Eigen_least_squares_solver LeastSquaresSolver; // Public operations public: @@ -289,7 +294,6 @@ 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) ; // Initialize the "A*X = B" linear system after // (at least two) border vertices parameterization @@ -412,13 +416,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()) ; + solver.set_value (2*index, uv.x()); + solver.set_value (2*index + 1, 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() ; + solver.lock (2*index); + solver.lock (2*index + 1); } } } @@ -574,8 +578,8 @@ 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() ; + NT u = solver.value (2*index); + NT v = solver.value (2*index + 1); // Fill vertex (u,v) and mark it as "parameterized" mesh.set_vertex_uv(vertexIt, Point_2(u,v)); diff --git a/Surface_mesh_parameterization/include/CGAL/internal/Eigen_least_squares_solver.h b/Surface_mesh_parameterization/include/CGAL/internal/Eigen_least_squares_solver.h new file mode 100644 index 00000000000..452a3ea1d00 --- /dev/null +++ b/Surface_mesh_parameterization/include/CGAL/internal/Eigen_least_squares_solver.h @@ -0,0 +1,181 @@ +#ifndef __EIGEN_LEAST_SQUARES_SOLVER__ +#define __EIGEN_LEAST_SQUARES_SOLVER__ + +#include +#include + +#include +#include +#include + + +namespace CGAL { + + namespace internal + { + class Eigen_least_squares_solver + { + + public: + typedef Eigen_solver_traits::EigenType> > Solver; + typedef typename Solver::Matrix Matrix ; + typedef typename Solver::Vector Vector ; + typedef typename Solver::NT Scalar ; + typedef typename CGAL::cpp11::tuple Indexed_scalar; + + + private: + + std::vector _indexed_scalars; + + std::vector _af; + std::vector _if; + std::vector _al; + std::vector _xl; + + Matrix* _A ; + Vector* _x ; + Vector* _b ; + + public: + + Eigen_least_squares_solver (unsigned int n_bvariables) + { + _indexed_scalars.resize (n_bvariables); + for (std::size_t i = 0; i < _indexed_scalars.size (); ++ i) + { + _indexed_scalars[i].get<0>() = 0.; + _indexed_scalars[i].get<1>() = -1; + _indexed_scalars[i].get<2>() = false; + } + _A = NULL ; + _x = NULL ; + _b = NULL ; + } + + ~Eigen_least_squares_solver() + { + if (_A != NULL) delete _A ; + if (_x != NULL) delete _x ; + if (_b != NULL) delete _b ; + } + + const Scalar& value (std::size_t index) const + { + return _indexed_scalars[index].get<0>(); + } + + void set_value (std::size_t index, Scalar value) + { + _indexed_scalars[index].get<0>() = value; + } + + void lock (std::size_t index) + { + _indexed_scalars[index].get<2>() = true; + } + + void store_scalars_in_x () + { + for (std::size_t i = 0; i < _indexed_scalars.size (); ++ i) + if (!(_indexed_scalars[i].get<2>())) + _x->set (_indexed_scalars[i].get<1>(), _indexed_scalars[i].get<0>()); + } + + void get_result () + { + for (std::size_t i = 0; i < _indexed_scalars.size (); ++ i) + if (!(_indexed_scalars[i].get<2>())) + _indexed_scalars[i].get<0>() = _x->get(_indexed_scalars[i].get<1>()); + } + + void begin_system () + { + std::size_t index = 0; + for (std::size_t i = 0; i < _indexed_scalars.size (); ++ i) + if(!(_indexed_scalars[i].get<2>())) + _indexed_scalars[i].get<1>() = index ++; + + _A = new Matrix (index, index); + _x = new Vector (index); + _b = new Vector (index); + + for (std::size_t i = 0; i < index; ++ i) + { + _x->set(i, 0.); + _b->set(i, 0.); + } + + store_scalars_in_x() ; + } + + void begin_row () + { + _af.clear() ; + _if.clear() ; + _al.clear() ; + _xl.clear() ; + } + + void add_coefficient (unsigned int iv, double a) + { + if(_indexed_scalars[iv].get<2>()) + { + _al.push_back(a) ; + _xl.push_back(_indexed_scalars[iv].get<0>()); + } + else + { + _af.push_back(a) ; + _if.push_back(_indexed_scalars[iv].get<1>()); + } + } + + void end_row () + { + std::size_t nf = _af.size() ; + std::size_t nl = _al.size() ; + + for (std::size_t i = 0; i < nf; ++ i) + for (std::size_t j = 0; j < nf; ++ j) + _A->add_coef (_if[i], _if[j], _af[i] * _af[j]); + + Scalar S = 0.; + + for (std::size_t j = 0; j < nl; ++ j) + S += _al[j] * _xl[j] ; + + for (std::size_t i = 0; i < nf; ++ i) + _b->set(_if[i], _b->get(_if[i]) - _af[i] * S); + } + + void end_system () + { + _A->assemble_matrix (); + } + + bool solve () + { + Scalar D; + Solver solver; + bool success = solver.linear_solver(*_A, *_b, *_x, D) ; + + get_result() ; + + delete _A ; _A = NULL ; + delete _b ; _b = NULL ; + delete _x ; _x = NULL ; + + return success; + } + + + + } ; + + + } // namespace internal + +} // namespace CGAL + +#endif