diff --git a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Concepts/WeightedApproximation_3.h b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Concepts/WeightedApproximation_3.h deleted file mode 100644 index 58c303a5add..00000000000 --- a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Concepts/WeightedApproximation_3.h +++ /dev/null @@ -1,136 +0,0 @@ -//Orthogonal projection of a point onto a weighted PCA plane. -//Copyright (C) 2014 INRIA - Sophia Antipolis -// -//This program is free software: you can redistribute it and/or modify -//it under the terms of the GNU General Public License as published by -//the Free Software Foundation, either version 3 of the License, or -//(at your option) any later version. -// -//This program 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 General Public License for more details. -// -//You should have received a copy of the GNU General Public License -//along with this program. If not, see . -// -// Author(s): Thijs van Lankveld - - -/// A concept for computing an approximation of a weighted point set. -/** \ingroup PkgScaleSpaceReconstruction3Concepts - * \cgalConcept - * This approximation can be used to fit other points to the point set. - * - * \cgalHasModel `CGAL::Weighted_PCA_approximation_3` - */ -class WeightedApproximation_3 { -public: -/// \name Types -/// \{ - typedef unspecified_type FT; ///< defines the field number type. - typedef unspecified_type Point; ///< defines the point type. - -/// \} - -public: -/// \name Constructors -/// \{ - /// constructs an approximation of a undefined point set. - /** The point set holds a fixed number of points with undefined coordinates. - * - * \param size is the size of the points set. - * - * \note this does not compute the approximation. - */ - WeightedApproximation_3( unsigned int size ); - -/// \} - - // computes the approximation of a point set. - /* Similar to constructing a approximation and calling - * [set_points(points_begin, points_end, weights_begin)](\ref WeightedApproximation_3::set_points ) - * - * \tparam PointIterator is an input iterator over the point collection. - * The value type of the iterator must be a `Point`. - * \tparam WeightIterator is an input iterator over the collection of - * weights. The value type of the iterator must be an `FT`. - * - * \param points_begin is an iterator to the first point. - * \param points_end is a past-the-end oterator for the points. - * \param weights_begin is an iterator to the weight of the first point. - * - * \pre The number of weights must be at least equal to the number of - * points. - */ - template < typename PointIterator, typename WeightIterator > - WeightedApproximation_3( PointIterator points_begin, PointIterator points_end, WeightIterator weights_begin ); - -public: -/// \name Point Set. -/// \{ - /// changes a weighted point in the set. - /** This invalidates the approximation. `compute()` should be called - * after all points have been set. - * \pre i must be smaller than the total size of the point set. - */ - void set_point( unsigned int i, const Point& p, const FT& w ); - - /// gives the size of the weighted point set. - std::size_t size() const; - -/// \} - - // changes the weighted point set. - /* After these points are set, the approximation is immediately computed. - * - * \tparam PointIterator is an input iterator over the point collection. - * The value type of the iterator must be a `Point`. - * \tparam WeightIterator is an input iterator over the collection of - * weights. The value type of the iterator must be an `FT`. - * - * \param points_begin is an iterator to the first point. - * \param points_end is a past-the-end iterator for the points. - * \param weights_begin is an iterator to the weight of the first point. - * - * \return whether the approximation converges. If the approximation does - * not converge this may indicate that the point set is too small, or the - * affine hull of the points cannot contain the approximation. - * - * \pre The points must fit in the approximation. - * \pre The number of weights must be at least equal to the number of - * points. - */ - template < typename PointIterator, typename WeightIterator > - bool set_points( PointIterator points_begin, PointIterator points_end, - WeightIterator weights_begin ); - -public: -/// \name Approximation -/// \{ - - /// computes the approximation. - /** \return whether the approximation converges. If the approximation does - * not converge this may indicate that the point set is too small, or the - * affine hull of the points cannot contain the approximation. - */ - bool compute(); - - /// checks whether the approximation has been computed successfully. - bool is_computed() const; -/// \} - -public: -/// \name Fitting -/// \{ - - /// fits a point to the approximation. - /** \param p is the point to fit. - * - * \return the point on the approximation closest to `p`. - * - * \pre The approximation must have been computed. - */ - Point fit( const Point& p ); -/// \} -}; // class WeightedApproximation_3 diff --git a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/PackageDescription.txt b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/PackageDescription.txt index 1ae18cbf86a..8351a47f4f8 100644 --- a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/PackageDescription.txt +++ b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/PackageDescription.txt @@ -26,14 +26,9 @@ \cgalClassifedRefPages -## Concepts ## - -- `WeightedApproximation_3` - ## Classes ## - `CGAL::Scale_space_surface_reconstruction_3` -- `CGAL::Weighted_PCA_approximation_3` */ diff --git a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Scale_space_reconstruction_3.txt b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Scale_space_reconstruction_3.txt index c16422509a8..beb6acb9d4a 100644 --- a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Scale_space_reconstruction_3.txt +++ b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Scale_space_reconstruction_3.txt @@ -90,7 +90,7 @@ The main class `Scale_space_surface_reconstruction_3` contains all the functiona The neighborhood size is estimated using `Orthogonal_k_neighbor_search`. The point set is generally stored in a `Orthogonal_k_neighbor_search::Tree`. When the neighborhood size is estimated, this tree is searched for nearest neighbors. -The scale-space is constructed at the original scale of the points. An iteration of increasing the scale is computed using a weighted PCA procedure. As described by Digne et al. \cgalCite{cgal:dmsl-ssmrp-11}, unlike similar methods this procedure does not lead to an undesirable clustering effect. By default the efficient \ref thirdpartyEigen library is used for this procedure. It is also possible to provide your own model for the `WeightedApproximation_3` concept. The weighted PCA procedure is performed locally per point, so it can be performed with parallel computing (linking with Intel TBB and passing the `Parallel_tag` to the reconstruction class is required). +The scale-space is constructed at the original scale of the points. An iteration of increasing the scale is computed using a weighted PCA procedure. As described by Digne et al. \cgalCite{cgal:dmsl-ssmrp-11}, unlike similar methods this procedure does not lead to an undesirable clustering effect. By default the efficient \ref thirdpartyEigen library is used for this procedure if available. Otherwise, the internal fallback `Internal_diagonalize_traits` is used. It is also possible to provide your own model for the `DiagonalizeTraits` concept. The weighted PCA procedure is performed locally per point, so it can be performed with parallel computing (linking with Intel TBB and passing the `Parallel_tag` to the reconstruction class is required). The mesh reconstruction is performed by filtering a \ref Chapter_3D_Alpha_Shapes "3D alpha shape" of the point set at a fixed scale. This filtering constructs a triangle for each regular facet; each singular facet results in two triangles facing opposite directions. diff --git a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/dependencies b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/dependencies index 73d4d9f43cf..9d74ca1bbde 100644 --- a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/dependencies +++ b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/dependencies @@ -4,3 +4,4 @@ STL_Extension Triangulation_3 Alpha_shapes_3 Spatial_searching +Solver_interface diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Scale_space_surface_reconstruction_3_impl.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Scale_space_surface_reconstruction_3_impl.h index dbeae87ca29..cfb763a5a85 100644 --- a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Scale_space_surface_reconstruction_3_impl.h +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Scale_space_surface_reconstruction_3_impl.h @@ -18,8 +18,6 @@ #ifndef CGAL_SCALE_SPACE_RECONSTRUCTION_3_SCALE_SPACE_SURFACE_RECONSTRUCTION_3_IMPL_H #define CGAL_SCALE_SPACE_RECONSTRUCTION_3_SCALE_SPACE_SURFACE_RECONSTRUCTION_3_IMPL_H -#include - //#include #ifdef CGAL_LINKED_WITH_TBB #include "tbb/blocked_range.h" @@ -204,7 +202,7 @@ public: // Compute the weighted least-squares planar approximation of the point set. CGAL::cpp11::array vnorm = {{ 0., 0., 0. }}; - CGAL::Default_diagonalize_traits::extract_largest_eigenvector_of_covariance_matrix + wA::extract_smallest_eigenvector_of_covariance_matrix (covariance, vnorm); // The vertex is moved by projecting it onto the plane diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Weighted_PCA_approximation_3.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Weighted_PCA_approximation_3.h deleted file mode 100644 index afeef08ed56..00000000000 --- a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Weighted_PCA_approximation_3.h +++ /dev/null @@ -1,177 +0,0 @@ -//Copyright (C) 2014 INRIA - Sophia Antipolis -// -//This program is free software: you can redistribute it and/or modify -//it under the terms of the GNU General Public License as published by -//the Free Software Foundation, either version 3 of the License, or -//(at your option) any later version. -// -//This program 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 General Public License for more details. -// -//You should have received a copy of the GNU General Public License -//along with this program. If not, see . -// -// Author(s): Thijs van Lankveld - - -#ifndef CGAL_SCALE_SPACE_RECONSTRUCTION_3_WEIGHTED_PCA_PROJECTION_3_H -#define CGAL_SCALE_SPACE_RECONSTRUCTION_3_WEIGHTED_PCA_PROJECTION_3_H - -#include - - -namespace CGAL { - -/// approximates a point set using a weighted least-squares plane. -/** \ingroup PkgScaleSpaceReconstruction3Classes - * The weighted least-squares planar approximation contains the barycenter - * of the points and is orthogonal to the eigenvector corresponding to the - * smallest eigenvalue. - * - * Point are fitted to this approximation by projecting them orthogonally onto - * the weighted least-squares plane. - * - * This class requires the eigenvector solvers of \ref thirdpartyEigen version - * 3.1.2 (or greater). - * - * \tparam Kernel the geometric traits class of the input and output. It - * must have a `RealEmbeddable` field number type. - * - * \note Irrespective of the geometric traits class, the approximation is - * always estimated up to double precision. - * - * \cgalModels `WeightedApproximation_3` - */ -template < class Kernel > -class Weighted_PCA_approximation_3 { -public: - typedef typename Kernel::FT FT; - typedef typename Kernel::Point_3 Point; - -private: - typedef Eigen::Matrix Matrix3D; // 3-by-dynamic double-value matrix. - typedef Eigen::Array Array1D; // 1-by-dynamic double-value array. - typedef Eigen::Matrix3d Matrix3; // 3-by-3 double-value matrix. - typedef Eigen::Vector3d Vector3; // 3-by-1 double-value matrix. - typedef Eigen::SelfAdjointEigenSolver< Matrix3 > EigenSolver; // Eigen value decomposition solver. - -private: - bool _comp; // Whether the eigenvalues are computed. - - Matrix3D _pts; // points. - Array1D _wts; // weights. - - Vector3 _bary; // barycenter. - Vector3 _norm; // normal. - -public: - // constructs an default approximation to hold the points. - /* \param size is the number of points that will be added. - */ - Weighted_PCA_approximation_3( unsigned int size ): _comp(false), _pts(3,size), _wts(1,size) {} - - // computes the weighted least-squares planar approximation of a point set. - /* Similar to constructing an empty projection and calling - * [set_points(points_begin, points_end, weights_begin)](\ref WeightedApproximation_3::set_points ) - * - * \tparam PointIterator is an input iterator over the point collection. - * The value type of the iterator must be a `Point`. - * \tparam WeightIterator is an input iterator over the collection of - * weights. The value type of the iterator must be an `FT`. - * - * \param points_begin is an iterator to the first point. - * \param points_end is a past-the-end oterator for the points. - * \param weights_begin is an iterator to the weight of the first point. - */ - template < typename PointIterator, typename WeightIterator > - Weighted_PCA_approximation_3( PointIterator points_begin, PointIterator points_end, WeightIterator weights_begin ) - : _comp(false) { - std::size_t size = std::distance( points_begin, points_end ); - _pts = Matrix3D(3,size); - _wts = Array1D(1,size); - set_points( points_begin, points_end, weights_begin ); - } - -public: - // sets a weighted point in the collection. - void set_point( unsigned int i, const Point& p, const FT& w ) { - _pts( 0, i ) = CGAL::to_double( p[0] ); - _pts( 1, i ) = CGAL::to_double( p[1] ); - _pts( 2, i ) = CGAL::to_double( p[2] ); - _wts( i ) = CGAL::to_double( w ); - _comp = false; - } - - // sets the weighted points and compute PCA. - template < typename PointIterator, typename WeightIterator > - bool set_points( PointIterator points_begin, PointIterator points_end, WeightIterator weights_begin ) { - for( unsigned int column = 0; points_begin != points_end; ++column, ++points_begin, ++weights_begin ) { - _pts( 0, column ) = CGAL::to_double( (*points_begin)[0] ); - _pts( 1, column ) = CGAL::to_double( (*points_begin)[1] ); - _pts( 2, column ) = CGAL::to_double( (*points_begin)[2] ); - _wts( column ) = CGAL::to_double( *weights_begin ); - } - return compute(); - } - - // gives the size of the weighted point set. - std::size_t size() const { return _pts.cols(); } - - // compute weighted PCA. - bool compute() { - // Construct the barycenter. - _bary = ( _pts.array().rowwise() * _wts ).rowwise().sum() / _wts.sum(); - - // Replace the points by their difference with the barycenter. - _pts = ( _pts.colwise() - _bary ).array().rowwise() * _wts; - - // Construct the weighted covariance matrix. - Matrix3 covar = Matrix3::Zero(); - for( int column = 0; column < _pts.cols(); ++column ) - covar += _wts.matrix()(column) * _pts.col(column) * _pts.col(column).transpose(); - - // Construct the Eigen system. - EigenSolver solver( covar ); - - // If the Eigen system does not converge, the vertex is not moved. - if( solver.info() != Eigen::Success ) - return false; - - // Find the Eigen vector with the smallest Eigen value. - std::ptrdiff_t index; - solver.eigenvalues().minCoeff( &index ); - if( solver.eigenvectors().col( index ).imag() != Vector3::Zero() ) { - // This should actually never happen, - // because the covariance matrix is symmetric! - CGAL_assertion( false ); - return false; - } - - // The normal is the Eigen vector with smallest Eigen value. - _norm = solver.eigenvectors().col( index ).real(); - _comp = true; - return true; - } - - // checks whether the weighted least-squares approximating plane has been computed. - bool is_computed() const { return _comp; } - -public: - // projects a point onto the weighted PCA plane. - Point fit( const Point& p ) { - CGAL_assertion( _comp ); - // The point is moved by projecting it onto the plane through the - // barycenter and orthogonal to the normal. - Vector3 to_p = Vector3( to_double( p[0] ), - to_double( p[1] ), - to_double( p[2] ) ) - _bary; - Vector3 proj = _bary + to_p - ( _norm.dot(to_p) * _norm ); - return Point( proj(0), proj(1), proj(2) ); - } -}; // class Weighted_PCA_approximation_3 - -} // namespace CGAL - -#endif // CGAL_SCALE_SPACE_RECONSTRUCTION_3_WEIGHTED_PCA_PROJECTION_3_H diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h index dae57b68e83..3462a7c92fb 100644 --- a/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h @@ -36,13 +36,10 @@ #include #include #include +#include #include -#ifdef CGAL_EIGEN3_ENABLED -#include -#endif // CGAL_EIGEN3_ENABLED - #include namespace CGAL { @@ -99,22 +96,18 @@ namespace CGAL { * only has an impact on the run-time. * \tparam Sh determines whether to collect the surface per shell. It * must be a `Boolean_tag` type. The default value is `Tag_true`. - * \tparam wA must be a model of `WeightedPCAProjection_3` and determines how - * to approximate a weighted point set. If \ref thirdpartyEigen 3.1.2 (or - * greater) is available and CGAL_EIGEN3_ENABLED is defined, then - * `Weighted_PCA_approximation_3` is used by default. + * \tparam wA must be a model of `DiagonalizeTraits` and determines + * how to diagonalize a weighted covariance matrix to approximate a + * weighted point set. It can be omitted: if Eigen 3 (or greater) is + * available and `CGAL_EIGEN3_ENABLED` is defined then an overload + * using `Eigen_diagonalize_traits` is provided. Otherwise, the + * internal implementation `Internal_diagonalize_traits` is used. * \tparam Ct indicates whether to use concurrent processing. It must be * either `Sequential_tag` or `Parallel_tag` (the default value). */ -template < class Gt, class FS = Tag_true, class Sh = Tag_true, class wA = Default, class Ct = Parallel_tag > +template < class Gt, class FS = Tag_true, class Sh = Tag_true, + class wA = Default_diagonalize_traits, class Ct = Parallel_tag > class Scale_space_surface_reconstruction_3 { - typedef typename Default::Get< wA, -#ifdef CGAL_EIGEN3_ENABLED - Weighted_PCA_approximation_3 -#else // CGAL_EIGEN3_ENABLED - void -#endif // CGAL_EIGEN3_ENABLED - >::type Approximation; public: typedef typename Gt::Point_3 Point; ///< defines the point type.