Wrapper to automatically use Eigen_diagonalize if Eigen is available or Internal_diagonalize otherwise

This commit is contained in:
Simon Giraudot 2015-08-31 09:34:17 +02:00
parent bb0d1c7b32
commit 0875fa17ff
10 changed files with 112 additions and 141 deletions

View File

@ -19,11 +19,7 @@
// Author(s) : Kaspar Fischer <fischerk@inf.ethz.ch>
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_diagonalize_traits.h>
#else
#include <CGAL/Internal_diagonalize_traits.h>
#endif
#include <CGAL/Default_diagonalize_traits.h>
#include <CGAL/Approximate_min_ellipsoid_d.h>
@ -161,13 +157,8 @@ namespace CGAL {
CGAL::cpp11::array<double, 4> eigenvectors; // Note: not neces. normalized.
CGAL::cpp11::array<double, 2> eigenvalues; // Note: sorted ascendent.
#ifdef CGAL_EIGEN3_ENABLED
CGAL::Eigen_diagonalize_traits<double, 2>::diagonalize_selfadjoint_covariance_matrix
CGAL::Default_diagonalize_traits<double, 2>::diagonalize_selfadjoint_covariance_matrix
(matrix, eigenvalues, eigenvectors);
#else
CGAL::Internal_diagonalize_traits<double, 2>::diagonalize_selfadjoint_covariance_matrix
(matrix, eigenvalues, eigenvectors);
#endif
// normalize eigenvectors:
double l1=1.0/std::sqrt(eigenvectors[2]*eigenvectors[2]+
@ -208,13 +199,9 @@ namespace CGAL {
CGAL::cpp11::array<double, 9> eigenvectors; // Note: not necessarily normalized.
CGAL::cpp11::array<double, 3> eigenvalues; // Note: sorted ascendent.
#ifdef CGAL_EIGEN3_ENABLED
CGAL::Eigen_diagonalize_traits<double, 3>::diagonalize_selfadjoint_covariance_matrix
CGAL::Default_diagonalize_traits<double, 3>::diagonalize_selfadjoint_covariance_matrix
(matrix, eigenvalues, eigenvectors);
#else
CGAL::Internal_diagonalize_traits<double, 3>::diagonalize_selfadjoint_covariance_matrix
(matrix, eigenvalues, eigenvectors);
#endif
// normalize eigenvectors:
double l1 = 1.0/std::sqrt(eigenvectors[0] * eigenvectors[0]+ // x^2

View File

@ -22,8 +22,8 @@
#include <CGAL/Simple_cartesian.h>
#include <CGAL/circulator.h>
#include <CGAL/Linear_algebraCd.h>
#include <CGAL/eigen.h>
#include <CGAL/Cartesian_converter.h>
#include <CGAL/Default_diagonalize_traits.h>
#include <math.h>
#include <utility>
#ifdef CGAL_EIGEN3_ENABLED
@ -362,27 +362,30 @@ compute_PCA(InputIterator begin, InputIterator end)
// assemble covariance matrix as a
// semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
// 3 4 5
FT covariance[6] = {xx,xy,yy,xz,yz,zz};
FT eigen_values[3];
FT eigen_vectors[9];
// 0 1 2
// 3 4
// 5
CGAL::cpp11::array<FT, 6> covariance = {{ xx,xy,xz,yy,yz,zz }};
CGAL::cpp11::array<FT, 3> eigen_values = {{ 0., 0., 0. }};
CGAL::cpp11::array<FT, 9> eigen_vectors = {{ 0., 0., 0. }};
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in descending order,
// eigen values are sorted in ascending order,
// eigen vectors are sorted in accordance.
CGAL::internal::eigen_symmetric<FT>(covariance,3,eigen_vectors,eigen_values);
CGAL::Default_diagonalize_traits<FT,3>::diagonalize_selfadjoint_covariance_matrix
(covariance, eigen_values, eigen_vectors);
//store in m_pca_basis
for (int i=0; i<3; i++)
{
m_pca_basis[i].first = eigen_values[i];
m_pca_basis[i].first = eigen_values[2-i];
}
Vector_3 v1(eigen_vectors[0],eigen_vectors[1],eigen_vectors[2]);
Vector_3 v1(eigen_vectors[6],eigen_vectors[7],eigen_vectors[8]);
m_pca_basis[0].second = v1;
Vector_3 v2(eigen_vectors[3],eigen_vectors[4],eigen_vectors[5]);
m_pca_basis[1].second = v2;
Vector_3 v3(eigen_vectors[6],eigen_vectors[7],eigen_vectors[8]);
Vector_3 v3(eigen_vectors[0],eigen_vectors[1],eigen_vectors[2]);
m_pca_basis[2].second = v3;
switch_to_direct_orientation(m_pca_basis[0].second,
m_pca_basis[1].second,
@ -523,15 +526,17 @@ compute_Monge_basis(const FT* A, Monge_form& monge_form)
//in the new orthonormal basis (Y,Z) of the tangent plane :
weingarten = inv *(1/det) * weingarten * change_XuXv2YZ;
//switch to eigen_symmetric algo for diagonalization of weingarten
FT W[3] = {weingarten(0,0), weingarten(1,0), weingarten(1,1)};
FT eval[2];
FT evec[4];
//eval in decreasing order
CGAL::internal::eigen_symmetric<FT>(W,2,evec,eval);
// diagonalization of weingarten
CGAL::cpp11::array<FT,3> W = {{ weingarten(0,0), weingarten(1,0), weingarten(1,1) }};
CGAL::cpp11::array<FT,2> eval = {{ 0., 0. }};
CGAL::cpp11::array<FT,4> evec = {{ 0., 0., 0., 0. }};
Vector_3 d_max = evec[0]*Y + evec[1]*Z,
d_min = evec[2]*Y + evec[3]*Z;
//eval in increasing order
CGAL::Default_diagonalize_traits<FT,2>::diagonalize_selfadjoint_covariance_matrix
(W, eval, evec);
Vector_3 d_max = evec[2]*Y + evec[3]*Z,
d_min = evec[0]*Y + evec[1]*Z;
switch_to_direct_orientation(d_max, d_min, normal);
Aff_transformation change_basis (d_max[0], d_max[1], d_max[2],
@ -547,8 +552,8 @@ compute_Monge_basis(const FT* A, Monge_form& monge_form)
monge_form.maximal_principal_direction() = L2D_converter(this->change_world2fitting.inverse()(d_max));
monge_form.minimal_principal_direction() = L2D_converter(this->change_world2fitting.inverse()(d_min));
monge_form.normal_direction() = L2D_converter(this->change_world2fitting.inverse()(normal));
monge_form.coefficients()[0] = L2D_NTconverter()(eval[0]);
monge_form.coefficients()[1] = L2D_NTconverter()(eval[1]);
monge_form.coefficients()[0] = L2D_NTconverter()(eval[1]);
monge_form.coefficients()[1] = L2D_NTconverter()(eval[0]);
}
//end else
}

View File

@ -29,12 +29,7 @@
#include <CGAL/linear_least_squares_fitting_circles_2.h>
#include <CGAL/linear_least_squares_fitting_rectangles_2.h>
#include <CGAL/Dimension.h>
#include <CGAL/Internal_diagonalize_traits.h>
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_diagonalize_traits.h>
#endif
#include <CGAL/Default_diagonalize_traits.h>
#include <iterator>
namespace CGAL {
@ -77,12 +72,7 @@ linear_least_squares_fitting_2(InputIterator first,
typedef typename Kernel_traits<Value_type>::Kernel Kernel;
return CGAL::linear_least_squares_fitting_2
(first,beyond,line,centroid,tag,Kernel(),
#ifdef CGAL_EIGEN3_ENABLED
Eigen_diagonalize_traits<typename Kernel::FT, 2>()
#else
Internal_diagonalize_traits<typename Kernel::FT, 2>()
#endif
);
Default_diagonalize_traits<typename Kernel::FT, 2>());
}
template < typename InputIterator,
@ -99,12 +89,7 @@ linear_least_squares_fitting_2(InputIterator first,
typedef typename Kernel_traits<Value_type>::Kernel Kernel;
typename Kernel::Point_2 centroid; // unused
return CGAL::linear_least_squares_fitting_2(first,beyond,line,centroid,tag,Kernel(),
#ifdef CGAL_EIGEN3_ENABLED
Eigen_diagonalize_traits<typename Kernel::FT, 2>()
#else
Internal_diagonalize_traits<typename Kernel::FT, 2>()
#endif
);
Default_diagonalize_traits<typename Kernel::FT, 2>());
}

View File

@ -29,10 +29,7 @@
#include <CGAL/linear_least_squares_fitting_tetrahedra_3.h>
#include <CGAL/linear_least_squares_fitting_spheres_3.h>
#include <CGAL/Internal_diagonalize_traits.h>
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_diagonalize_traits.h>
#endif
#include <CGAL/Default_diagonalize_traits.h>
#include <CGAL/Dimension.h>
@ -80,12 +77,7 @@ linear_least_squares_fitting_3(InputIterator first,
typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
typedef typename Kernel_traits<Value_type>::Kernel Kernel;
return CGAL::linear_least_squares_fitting_3(first,beyond,object,centroid,tag,Kernel(),
#ifdef CGAL_EIGEN3_ENABLED
Eigen_diagonalize_traits<typename Kernel::FT, 3>()
#else
Internal_diagonalize_traits<typename Kernel::FT, 3>()
#endif
);
Default_diagonalize_traits<typename Kernel::FT, 3>());
}
@ -105,12 +97,7 @@ linear_least_squares_fitting_3(InputIterator first,
typedef typename Kernel_traits<Value_type>::Kernel Kernel;
typename Kernel::Point_3 centroid; // not used by caller
return CGAL::linear_least_squares_fitting_3(first,beyond,object,centroid,tag,Kernel(),
#ifdef CGAL_EIGEN3_ENABLED
Eigen_diagonalize_traits<typename Kernel::FT, 3>()
#else
Internal_diagonalize_traits<typename Kernel::FT, 3>()
#endif
);
Default_diagonalize_traits<typename Kernel::FT, 3>());
}

View File

@ -22,13 +22,6 @@
#include <CGAL/basic.h>
#include <CGAL/centroid.h>
#include <CGAL/Internal_diagonalize_traits.h>
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_diagonalize_traits.h>
#endif
#include <iterator>
#include <cmath>

View File

@ -3,11 +3,7 @@
#include <CGAL/algorithm.h>
#include <CGAL/linear_least_squares_fitting_2.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/Internal_diagonalize_traits.h>
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_diagonalize_traits.h>
#endif
#include <CGAL/Default_diagonalize_traits.h>
#include <vector>
#include <cassert>
@ -37,12 +33,7 @@ void test_2D()
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>(),k,
#ifdef CGAL_EIGEN3_ENABLED
CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 2>()
#else
CGAL::Internal_diagonalize_traits<typename Kernel::FT, 2>()
#endif
);
CGAL::Default_diagonalize_traits<FT,2>());
std::cout << "done (quality: " << quality << ")" << std::endl;
@ -84,12 +75,8 @@ void test_2D_point_set(const unsigned int nb_points)
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>(),k,
#ifdef CGAL_EIGEN3_ENABLED
CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 2>()
#else
CGAL::Internal_diagonalize_traits<typename Kernel::FT, 2>()
#endif
);
CGAL::Default_diagonalize_traits<FT,2>());
std::cout << "done (quality: " << quality << ")" << std::endl;
if(!line.is_horizontal())

View File

@ -3,10 +3,7 @@
#include <CGAL/Cartesian.h>
#include <CGAL/linear_least_squares_fitting_3.h>
#include <CGAL/Internal_diagonalize_traits.h>
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_diagonalize_traits.h>
#endif
#include <CGAL/Default_diagonalize_traits.h>
#include <cassert>
@ -49,12 +46,7 @@ void fit_point_set(std::list<Point>& points,
quality = linear_least_squares_fitting_3(points.begin(),points.end(),line,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_3(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_3(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>(),kernel,
#ifdef CGAL_EIGEN3_ENABLED
CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 3>()
#else
CGAL::Internal_diagonalize_traits<typename Kernel::FT, 3>()
#endif
);
CGAL::Default_diagonalize_traits<FT,3>());
std::cout << "done (quality: " << quality << ")" << std::endl;
@ -62,12 +54,8 @@ void fit_point_set(std::list<Point>& points,
quality = linear_least_squares_fitting_3(points.begin(),points.end(),plane,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_3(points.begin(),points.end(),plane,centroid,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_3(points.begin(),points.end(),plane,centroid,CGAL::Dimension_tag<0>(),kernel,
#ifdef CGAL_EIGEN3_ENABLED
CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 3>()
#else
CGAL::Internal_diagonalize_traits<typename Kernel::FT, 3>()
#endif
);
CGAL::Default_diagonalize_traits<FT,3>());
std::cout << "done (quality: " << quality << ")" << std::endl;
}

View File

@ -5,10 +5,7 @@
#include <CGAL/algorithm.h>
#include <CGAL/linear_least_squares_fitting_3.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/Internal_diagonalize_traits.h>
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_diagonalize_traits.h>
#endif
#include <CGAL/Default_diagonalize_traits.h>
#include <vector>
#include <cassert>
@ -39,26 +36,14 @@ FT fit_set(std::list<Segment>& segments,
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),line,CGAL::Dimension_tag<1>());
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),line,centroid,CGAL::Dimension_tag<1>());
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),line,centroid,CGAL::Dimension_tag<1>(),kernel,
#ifdef CGAL_EIGEN3_ENABLED
CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 3>()
#else
CGAL::Internal_diagonalize_traits<typename Kernel::FT, 3>()
#endif
);
CGAL::Default_diagonalize_traits<FT,3>());
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,CGAL::Dimension_tag<1>());
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,centroid,CGAL::Dimension_tag<1>());
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,centroid,CGAL::Dimension_tag<1>(),kernel,
#ifdef CGAL_EIGEN3_ENABLED
CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 3>()
#else
CGAL::Internal_diagonalize_traits<typename Kernel::FT, 3>()
#endif
);
CGAL::Default_diagonalize_traits<FT,3>());
return quality;
return quality;
}

View File

@ -2,6 +2,7 @@
#include <CGAL/Cartesian.h>
#include <CGAL/linear_least_squares_fitting_3.h>
#include <CGAL/Default_diagonalize_traits.h>
#include <list>
@ -34,28 +35,23 @@ int main(void)
Kernel kernel;
Point centroid;
#ifdef CGAL_EIGEN3_ENABLED
CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 3> diagonalize_traits;
#else
CGAL::Internal_diagonalize_traits<typename Kernel::FT, 3> diagonalize_traits;
#endif
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,CGAL::Dimension_tag<3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,CGAL::Dimension_tag<2>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,centroid,CGAL::Dimension_tag<3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,centroid,CGAL::Dimension_tag<2>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,centroid,CGAL::Dimension_tag<3>(),kernel,
diagonalize_traits);
CGAL::Default_diagonalize_traits<FT,3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,centroid,CGAL::Dimension_tag<2>(),kernel,
diagonalize_traits);
CGAL::Default_diagonalize_traits<FT,3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,CGAL::Dimension_tag<3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,CGAL::Dimension_tag<2>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,centroid,CGAL::Dimension_tag<3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,centroid,CGAL::Dimension_tag<2>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,centroid,CGAL::Dimension_tag<3>(),kernel,
diagonalize_traits);
CGAL::Default_diagonalize_traits<FT,3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,centroid,CGAL::Dimension_tag<2>(),kernel,
diagonalize_traits);
CGAL::Default_diagonalize_traits<FT,3>());
return 0;
}

View File

@ -0,0 +1,58 @@
#ifndef CGAL_DEFAULT_DIAGONALIZE_TRAITS_H
#define CGAL_DEFAULT_DIAGONALIZE_TRAITS_H
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_diagonalize_traits.h>
#else
#include <CGAL/Internal_diagonalize_traits.h>
#endif
namespace CGAL {
// Wrapper class designed to automatically use
// Eigen_diagonalize_traits if Eigen is available and otherwise use
// the fallback Internal_diagonalize_traits class.
template <typename FT, unsigned int dim = 3>
class Default_diagonalize_traits{
#ifdef CGAL_EIGEN3_ENABLED
typedef Eigen_diagonalize_traits<FT, dim> Base;
#else
typedef Internal_diagonalize_traits<FT, dim> Base;
#endif
public:
static bool
diagonalize_selfadjoint_covariance_matrix(
const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT, dim>& eigenvalues)
{
return Base::diagonalize_selfadjoint_covariance_matrix (cov, eigenvalues);
}
static bool
diagonalize_selfadjoint_covariance_matrix(
const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT, dim>& eigenvalues,
cpp11::array<FT, dim * dim>& eigenvectors)
{
return Base::diagonalize_selfadjoint_covariance_matrix (cov, eigenvalues, eigenvectors);
}
// Extract the eigenvector associated to the largest eigenvalue
static bool
extract_largest_eigenvector_of_covariance_matrix (
const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT,dim> &normal)
{
return Base::extract_largest_eigenvector_of_covariance_matrix (cov, normal);
}
};
} // namespace CGAL
#endif // CGAL_DEFAULT_DIAGONALIZE_TRAITS_H