mirror of https://github.com/CGAL/cgal
Merge branch 'Surface_mesh_parametrisation-OpenNL_Eigen-GF-old' into Surface_mesh_parametrisation-OpenNL_Eigen-GF
This commit is contained in:
Simon Giraudot
commit
b9bbadf631
|
|
@ -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} )
|
||||
|
|
|
|||
|
|
@ -37,7 +37,7 @@
|
|||
<x>0</x>
|
||||
<y>0</y>
|
||||
<width>978</width>
|
||||
<height>21</height>
|
||||
<height>26</height>
|
||||
</rect>
|
||||
</property>
|
||||
<widget class="QMenu" name="menuFile">
|
||||
|
|
@ -95,6 +95,7 @@
|
|||
</property>
|
||||
<addaction name="actionMVC"/>
|
||||
<addaction name="actionDCP"/>
|
||||
<addaction name="actionLSC"/>
|
||||
</widget>
|
||||
<widget class="QMenu" name="menuPCA">
|
||||
<property name="title">
|
||||
|
|
@ -652,6 +653,11 @@
|
|||
<string>&Preferences</string>
|
||||
</property>
|
||||
</action>
|
||||
<action name="actionLSC">
|
||||
<property name="text">
|
||||
<string>Least square conformal maps</string>
|
||||
</property>
|
||||
</action>
|
||||
</widget>
|
||||
<customwidgets>
|
||||
<customwidget>
|
||||
|
|
|
|||
|
|
@ -12,6 +12,7 @@
|
|||
#include <CGAL/Parameterization_polyhedron_adaptor_3.h>
|
||||
#include <CGAL/parameterize.h>
|
||||
#include <CGAL/Discrete_conformal_map_parameterizer_3.h>
|
||||
#include <CGAL/LSCM_parameterizer_3.h>
|
||||
#include <CGAL/Two_vertices_parameterizer_3.h>
|
||||
|
||||
#include <CGAL/Textured_polyhedron_builder.h>
|
||||
|
|
@ -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<Adaptor> 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<Adaptor> 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"
|
||||
|
|
|
|||
|
|
@ -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 <CGAL/Eigen_solver_traits.h>
|
||||
|
||||
#include <vector>
|
||||
#include <iostream>
|
||||
#include <cstdlib>
|
||||
|
||||
|
||||
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<Eigen::ConjugateGradient<Eigen_sparse_symmetric_matrix<double>::EigenType,
|
||||
// Eigen::Lower|Eigen::Upper> > Solver;
|
||||
typedef Eigen_solver_traits<Eigen::SimplicialLDLT<Eigen_sparse_symmetric_matrix<double>::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<int>(n),static_cast<int>(n)) ;
|
||||
x_ = new Vector(n) ;
|
||||
b_ = new Vector(n) ;
|
||||
for(unsigned int i=0; i<n; i++) {
|
||||
x_->set(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; i<nf; i++) {
|
||||
norm += af_[i] * af_[i] ;
|
||||
}
|
||||
unsigned int nl = al_.size() ;
|
||||
for(unsigned int i=0; i<nl; i++) {
|
||||
norm += al_[i] * al_[i] ;
|
||||
}
|
||||
norm = sqrt(norm) ;
|
||||
CGAL_assertion( fabs(norm)>1e-40 );
|
||||
scale_row(weight / norm) ;
|
||||
}
|
||||
|
||||
void scale_row(CoeffType s) {
|
||||
check_state(IN_ROW) ;
|
||||
unsigned int nf = af_.size() ;
|
||||
for(unsigned int i=0; i<nf; i++) {
|
||||
af_[i] *= s ;
|
||||
}
|
||||
unsigned int nl = al_.size() ;
|
||||
for(unsigned int i=0; i<nl; i++) {
|
||||
al_[i] *= s ;
|
||||
}
|
||||
bk_ *= s ;
|
||||
}
|
||||
|
||||
void end_row() {
|
||||
unsigned int nf = af_.size() ;
|
||||
unsigned int nl = al_.size() ;
|
||||
for(unsigned int i=0; i<nf; i++) {
|
||||
for(unsigned int j=0; j<nf; j++) {
|
||||
A_->add_coef (if_[i], if_[j], af_[i] * af_[j]);
|
||||
}
|
||||
}
|
||||
CoeffType S = - bk_ ;
|
||||
for(unsigned int j=0; j<nl; j++) {
|
||||
S += al_[j] * xl_[j] ;
|
||||
}
|
||||
for(unsigned int i=0; i<nf; i++) {
|
||||
b_->set(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<CoeffType> af_ ;
|
||||
std::vector<unsigned int> if_ ;
|
||||
std::vector<CoeffType> al_ ;
|
||||
std::vector<CoeffType> xl_ ;
|
||||
double bk_ ;
|
||||
|
||||
// --------------- internal representation ---------
|
||||
Matrix* A_ ;
|
||||
Vector* x_ ;
|
||||
Vector* b_ ;
|
||||
|
||||
} ;
|
||||
|
||||
|
||||
} // namespace OpenNL
|
||||
|
||||
#endif
|
||||
|
|
@ -49,10 +49,20 @@ namespace internal {
|
|||
};
|
||||
|
||||
template <class FT,class EigenMatrix>
|
||||
struct Get_eigen_matrix< ::Eigen::ConjugateGradient<EigenMatrix>,FT>{
|
||||
struct Get_eigen_matrix< ::Eigen::ConjugateGradient<EigenMatrix, ::Eigen::Lower|::Eigen::Upper>,FT>{
|
||||
typedef Eigen_sparse_symmetric_matrix<FT> type;
|
||||
};
|
||||
|
||||
template <class FT,class EigenMatrix>
|
||||
struct Get_eigen_matrix< ::Eigen::SimplicialLDLT<EigenMatrix>,FT>{
|
||||
typedef Eigen_sparse_symmetric_matrix<FT> type;
|
||||
};
|
||||
|
||||
template <class FT,class EigenMatrix>
|
||||
struct Get_eigen_matrix< ::Eigen::LeastSquaresConjugateGradient<EigenMatrix>,FT>{
|
||||
typedef Eigen_sparse_matrix<FT> type;
|
||||
};
|
||||
|
||||
template <class FT,class EigenMatrix>
|
||||
struct Get_eigen_matrix< ::Eigen::SimplicialCholesky<EigenMatrix>,FT>{
|
||||
typedef Eigen_sparse_symmetric_matrix<FT> 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;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -86,6 +86,10 @@ public:
|
|||
void set(std::size_t i,NT value) {
|
||||
this->operator[](static_cast<int>(i))=value;
|
||||
}
|
||||
|
||||
NT get(std::size_t i) const {
|
||||
return this->operator[](static_cast<int>(i));
|
||||
}
|
||||
|
||||
NT* vector() {return this->data();}
|
||||
|
||||
|
|
|
|||
|
|
@ -3,7 +3,7 @@
|
|||
#include <CGAL/IO/Polyhedron_iostream.h>
|
||||
#include <CGAL/Parameterization_polyhedron_adaptor_3.h>
|
||||
#include <CGAL/parameterize.h>
|
||||
|
||||
#include <CGAL/Timer.h>
|
||||
#include <CGAL/Eigen_solver_traits.h>
|
||||
|
||||
#include <iostream>
|
||||
|
|
@ -17,7 +17,7 @@
|
|||
|
||||
typedef CGAL::Simple_cartesian<double> Kernel;
|
||||
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
|
||||
|
||||
typedef CGAL::Timer Timer;
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
// main()
|
||||
|
|
@ -65,7 +65,8 @@ int main(int argc, char * argv[])
|
|||
|
||||
typedef CGAL::Parameterization_polyhedron_adaptor_3<Polyhedron>
|
||||
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<Parameterization_polyhedron_adaptor>
|
||||
Border_parameterizer;
|
||||
// Eigen solver
|
||||
typedef CGAL::Eigen_solver_traits<> Solver;
|
||||
typedef CGAL::Eigen_solver_traits<Eigen::BiCGSTAB<CGAL::Eigen_sparse_matrix<double>::EigenType, Eigen::IncompleteLUT< double > > > Solver;
|
||||
|
||||
// Floater Mean Value Coordinates parameterization
|
||||
// (circular border) with Eigen solver
|
||||
typedef CGAL::Mean_value_coordinates_parameterizer_3<Parameterization_polyhedron_adaptor,
|
||||
Border_parameterizer,
|
||||
Solver>
|
||||
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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -32,6 +32,12 @@
|
|||
#include <CGAL/Parameterization_mesh_feature_extractor.h>
|
||||
#include <iostream>
|
||||
|
||||
#ifdef CGAL_EIGEN3_ENABLED
|
||||
#include <CGAL/internal/Eigen_least_squares_solver.h>
|
||||
#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<Sparse_LA>
|
||||
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));
|
||||
|
|
|
|||
|
|
@ -0,0 +1,181 @@
|
|||
#ifndef __EIGEN_LEAST_SQUARES_SOLVER__
|
||||
#define __EIGEN_LEAST_SQUARES_SOLVER__
|
||||
|
||||
#include <CGAL/Eigen_solver_traits.h>
|
||||
#include <CGAL/tuple.h>
|
||||
|
||||
#include <vector>
|
||||
#include <iostream>
|
||||
#include <cstdlib>
|
||||
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
namespace internal
|
||||
{
|
||||
class Eigen_least_squares_solver
|
||||
{
|
||||
|
||||
public:
|
||||
typedef Eigen_solver_traits<Eigen::SimplicialLDLT<Eigen_sparse_symmetric_matrix<double>::EigenType> > Solver;
|
||||
typedef typename Solver::Matrix Matrix ;
|
||||
typedef typename Solver::Vector Vector ;
|
||||
typedef typename Solver::NT Scalar ;
|
||||
typedef typename CGAL::cpp11::tuple<Scalar, std::size_t, bool> Indexed_scalar;
|
||||
|
||||
|
||||
private:
|
||||
|
||||
std::vector<Indexed_scalar> _indexed_scalars;
|
||||
|
||||
std::vector<Scalar> _af;
|
||||
std::vector<unsigned int> _if;
|
||||
std::vector<Scalar> _al;
|
||||
std::vector<Scalar> _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
|
||||
Loading…
Reference in New Issue