mirror of https://github.com/CGAL/cgal
Update scale space (code + doc) using the DiagonalizeTraits concept
This commit is contained in:
Simon Giraudot
parent
22eea088db
commit
a107ae4fff
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
//
|
||||
// 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
|
||||
* <code>[set_points(points_begin, points_end, weights_begin)](\ref WeightedApproximation_3::set_points )</code>
|
||||
*
|
||||
* \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
|
||||
|
|
@ -26,14 +26,9 @@
|
|||
|
||||
\cgalClassifedRefPages
|
||||
|
||||
## Concepts ##
|
||||
|
||||
- `WeightedApproximation_3`
|
||||
|
||||
## Classes ##
|
||||
|
||||
- `CGAL::Scale_space_surface_reconstruction_3<Gt,FS,Sh,wA,Ct>`
|
||||
- `CGAL::Weighted_PCA_approximation_3<Kernel>`
|
||||
|
||||
*/
|
||||
|
||||
|
|
|
|||
|
|
@ -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 <i>et al.</i> \cgalCite{cgal:dmsl-ssmrp-11}, unlike similar methods this procedure does not lead to an <em>undesirable clustering effect</em>. 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 <i>et al.</i> \cgalCite{cgal:dmsl-ssmrp-11}, unlike similar methods this procedure does not lead to an <em>undesirable clustering effect</em>. 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.
|
||||
|
||||
|
|
|
|||
|
|
@ -4,3 +4,4 @@ STL_Extension
|
|||
Triangulation_3
|
||||
Alpha_shapes_3
|
||||
Spatial_searching
|
||||
Solver_interface
|
||||
|
|
|
|||
|
|
@ -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 <CGAL/Default_diagonalize_traits.h>
|
||||
|
||||
//#include <omp.h>
|
||||
#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<FT, 3> vnorm = {{ 0., 0., 0. }};
|
||||
CGAL::Default_diagonalize_traits<double,3>::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
|
||||
|
|
|
|||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
//
|
||||
// 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 <Eigen/Dense>
|
||||
|
||||
|
||||
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<double, 3, Eigen::Dynamic> Matrix3D; // 3-by-dynamic double-value matrix.
|
||||
typedef Eigen::Array<double, 1, Eigen::Dynamic> 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
|
||||
* <code>[set_points(points_begin, points_end, weights_begin)](\ref WeightedApproximation_3::set_points )</code>
|
||||
*
|
||||
* \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
|
||||
|
|
@ -36,13 +36,10 @@
|
|||
#include <CGAL/Orthogonal_k_neighbor_search.h>
|
||||
#include <CGAL/Fuzzy_sphere.h>
|
||||
#include <CGAL/Random.h>
|
||||
#include <CGAL/Default_diagonalize_traits.h>
|
||||
|
||||
#include <CGAL/Scale_space_reconstruction_3/Shape_construction_3.h>
|
||||
|
||||
#ifdef CGAL_EIGEN3_ENABLED
|
||||
#include <CGAL/Scale_space_reconstruction_3/Weighted_PCA_approximation_3.h>
|
||||
#endif // CGAL_EIGEN3_ENABLED
|
||||
|
||||
#include <boost/mpl/and.hpp>
|
||||
|
||||
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<DelaunayTriangulationTraits_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<typename Gt::FT, 3>, class Ct = Parallel_tag >
|
||||
class Scale_space_surface_reconstruction_3 {
|
||||
typedef typename Default::Get< wA,
|
||||
#ifdef CGAL_EIGEN3_ENABLED
|
||||
Weighted_PCA_approximation_3<Gt>
|
||||
#else // CGAL_EIGEN3_ENABLED
|
||||
void
|
||||
#endif // CGAL_EIGEN3_ENABLED
|
||||
>::type Approximation;
|
||||
|
||||
public:
|
||||
typedef typename Gt::Point_3 Point; ///< defines the point type.
|
||||
|
|
|
|||
Loading…
Reference in New Issue