mirror of https://github.com/CGAL/cgal
First real setup for the documentation.
This commit is contained in:
parent
1729e9b098
commit
941e01ffa3
|
|
@ -0,0 +1,9 @@
|
|||
@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
|
||||
|
||||
PROJECT_NAME = "CGAL ${CGAL_CREATED_VERSION_NUM} - Scale-space surface reconstruction"
|
||||
|
||||
INPUT = ${CMAKE_SOURCE_DIR}/Scale_space_reconstruction_3/
|
||||
|
||||
EXTRACT_ALL = NO
|
||||
HIDE_UNDOC_CLASSES = YES
|
||||
HIDE_UNDOC_MEMBERS = NO
|
||||
|
|
@ -0,0 +1,55 @@
|
|||
/// \defgroup PkgScaleSpaceReconstruction3 Scale-space surface reconstruction Reference
|
||||
/// \defgroup PkgScaleSpaceReconstruction3Concepts Concepts
|
||||
/// \ingroup PkgScaleSpaceReconstruction3
|
||||
|
||||
|
||||
/// \defgroup PkgScaleSpaceReconstruction3Classes Reconstruction Classes
|
||||
/// \ingroup PkgScaleSpaceReconstruction3
|
||||
|
||||
/// \defgroup PkgScaleSpaceReconstruction3Auxiliary Auxiliary Classes
|
||||
/// \ingroup PkgScaleSpaceReconstruction3
|
||||
|
||||
|
||||
/*!
|
||||
\addtogroup PkgScaleSpaceReconstruction3
|
||||
\todo check generated documentation
|
||||
\cgalPkgDescriptionBegin{Scale-space surface reconstruction,PkgScaleSpaceReconstruction3Summary}
|
||||
\cgalPkgPicture{knot_thumb.png}
|
||||
\cgalPkgSummaryBegin
|
||||
\cgalPkgAuthors{Thijs van Lankveld}
|
||||
\cgalPkgDesc{This method allows to reconstruct a surface in 3D from a set of points. This method provides an efficient alternative to the poisson surface reconstruction method. The main difference in output is that this method reconstructs a surface that interpolates the point set (as opposed to approximating the point set). How the surface connects the points depends on a scale variable, which can be estimated semi-automatically.}
|
||||
\cgalPkgManuals{Chapter_Scale_space_reconstruction,PkgScaleSpaceReconstruction3}
|
||||
\cgalPkgSummaryEnd
|
||||
\cgalPkgShortInfoBegin
|
||||
\cgalPkgSince{4.5}
|
||||
\cgalPkgBib{cgal:ssr3}
|
||||
\cgalPkgLicense{\ref licensesGPL "GPL" }
|
||||
\cgalPkgShortInfoEnd
|
||||
\cgalPkgDescriptionEnd
|
||||
|
||||
The scale-space surface reconstruction method is twofold: first the point set is transformed into a coarse scale while maintaining point references to the original point set, then the mesh is reconstructed at this coarse scale (where reconstruction problem is less ill-posed) and transported back to the original reference points.
|
||||
|
||||
The transformation to the coarse scale is done by projecting each point to the `density`-weighted principal component analysis (PCA) of the local (\f$ \delta \f$-distance) neighborhood. If the high-frequency deformation is significantly smaller than the neighborhood size, the scale is coarse enough for mesh reconstruction after a few iterations of the tranformation step. For the PCA procedure, we use the efficient Eigen libraries. The mesh reconstruction at the coarse scale is performed by filtering a 3D alpha-shape. The result is returned as a collection of triples on the indices of the points. This collection may be sorted per `shell`, where a shell is a collection of connected triangles that are locally oriented towards the same side of the surface.
|
||||
|
||||
This method requires a user-specified parameter related to the resolution of the data, but we may be able to estimate this parameter through statistical analysis. Specifically, we use a kD-tree to estimate the mean distance to the n-th nearest neighbor and we use this distance as an approximator for the resolution. Generally a neighborhood containing 30 points is a good estimate. The method generally works well as long as this estimate is not too small and the amount of iterations necessary to reduce the high-frequency signals does not disturb the topology of the data such that back-transportation to the original scale produces too many self-intersections.
|
||||
|
||||
The method provides access to intermediate results and the user can adjust these to better suit his needs. These intermediate results are the estimate of the scale variable, the smoothed point set, and the final collection of surface triangles.
|
||||
|
||||
\cgalClassifedRefPages
|
||||
|
||||
## Concepts ##
|
||||
|
||||
## Classes ##
|
||||
|
||||
### Main Classes ###
|
||||
|
||||
- `CGAL::Scale_space_surface_reconstructer_3<DelaunayTriangulationTraits_3,Boolean_tags,Boolean_tags>`
|
||||
- `CGAL::Shape_type<DelaunayTriangulationTraits_3,Boolean_tags>`
|
||||
|
||||
### Auxiliary Classes ###
|
||||
|
||||
- `CGAL::internal::Auto_count<T>`
|
||||
- `CGAL::internal::_Env<T>`
|
||||
|
||||
*/
|
||||
|
||||
|
|
@ -0,0 +1,35 @@
|
|||
namespace CGAL {
|
||||
/*!
|
||||
|
||||
\mainpage User Manual
|
||||
\anchor Chapter_Scale_space_reconstruction
|
||||
\anchor chapterScaleSpaceReconstruction3
|
||||
\cgalAutoToc
|
||||
|
||||
\authors Thijs van Lankveld
|
||||
|
||||
General introduction.
|
||||
|
||||
\section ScaleSpaceReconstruction3secintro Representation
|
||||
|
||||
Something.
|
||||
|
||||
\section ScaleSpaceReconstruction3secdesign Software Design
|
||||
|
||||
Something else..
|
||||
|
||||
\section Triangulation3seccomplexity Complexity and Performance
|
||||
|
||||
Should I have this section?
|
||||
|
||||
\subsection Triangulation_3RunningTime Running Time
|
||||
|
||||
Probably, right?
|
||||
|
||||
\section Triangulation_3Design Design and Implementation History
|
||||
|
||||
This method has in fact, actually, and really been implemented at some point.
|
||||
|
||||
*/
|
||||
} /* namespace CGAL */
|
||||
|
||||
|
|
@ -0,0 +1,5 @@
|
|||
Manual
|
||||
Kernel_23
|
||||
STL_Extension
|
||||
Triangulation_3
|
||||
Alpha_shapes_3
|
||||
|
|
@ -0,0 +1,3 @@
|
|||
/*!
|
||||
\example Scale_space_reconstruction_3/scale_space.cpp
|
||||
*/
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 18 KiB |
|
|
@ -110,7 +110,8 @@ int main(int argc, char** argv) {
|
|||
sscanf(argv[4], "%u", &samples);
|
||||
|
||||
std::string input = base + ".off";
|
||||
std::string output = base + "_ss.off";
|
||||
std::string output_ss = base + "_ss.off";
|
||||
std::string output_sm = base + "_sm.off";
|
||||
|
||||
// Read the data.
|
||||
std::cout << "Input: " << input << std::flush;
|
||||
|
|
@ -128,9 +129,18 @@ int main(int argc, char** argv) {
|
|||
std::cout << "Reconstruction done." << std::endl;
|
||||
|
||||
// Write the reconstruction.
|
||||
std::cout << "Output: " << output << std::flush;
|
||||
if( !writeOFF( output, points, reconstruct.surface() ) ) {
|
||||
std::cerr << std::endl << "Error writing " << output << std::endl;
|
||||
std::cout << "Output: " << output_ss << std::flush;
|
||||
if( !writeOFF( output_ss, points, reconstruct.surface() ) ) {
|
||||
std::cerr << std::endl << "Error writing " << output_ss << std::endl;
|
||||
exit(-1);
|
||||
}
|
||||
std::cout << " written." << std::endl;
|
||||
|
||||
// Write the reconstruction.
|
||||
std::cout << "Output: " << output_sm << std::flush;
|
||||
Pointset smoothed( reconstruct.scale_space_begin(), reconstruct.scale_space_end() );
|
||||
if( !writeOFF( output_sm, smoothed, reconstruct.surface() ) ) {
|
||||
std::cerr << std::endl << "Error writing " << output_ss << std::endl;
|
||||
exit(-1);
|
||||
}
|
||||
std::cout << " written." << std::endl;
|
||||
|
|
|
|||
|
|
@ -45,17 +45,19 @@ namespace CGAL {
|
|||
*/
|
||||
template < class Kernel, class Fixed_shape >
|
||||
class Shape_type {
|
||||
typedef typename Kernel::FT Scalar;
|
||||
typedef typename Kernel::FT Scalar;
|
||||
typedef typename Kernel::Point_3 Point;
|
||||
typedef internal::Auto_count<Point> PointIndex;
|
||||
|
||||
typedef CGAL::Alpha_shape_vertex_base_3< Kernel,
|
||||
CGAL::Triangulation_vertex_base_with_info_3< unsigned int, Kernel > > Vb;
|
||||
typedef CGAL::Alpha_shape_cell_base_3< Kernel,
|
||||
CGAL::Triangulation_cell_base_with_info_3< unsigned int, Kernel > > Cb;
|
||||
typedef CGAL::Triangulation_data_structure_3<Vb,Cb> Tds;
|
||||
typedef Triangulation_vertex_base_with_info_3< unsigned int, Kernel > Vb;
|
||||
typedef Alpha_shape_vertex_base_3< Kernel, Vb > aVb;
|
||||
typedef Triangulation_cell_base_with_info_3< unsigned int, Kernel > Cb;
|
||||
typedef Alpha_shape_cell_base_3< Kernel, Cb > aCb;
|
||||
typedef Triangulation_data_structure_3<aVb,aCb> Tds;
|
||||
|
||||
public:
|
||||
typedef CGAL::Delaunay_triangulation_3< Kernel, Tds > Structure; ///< The triangulation that spatially orders the point set.
|
||||
typedef CGAL::Alpha_shape_3< Structure > Shape; ///< The structure that identifies the triangles in the surface.
|
||||
typedef Delaunay_triangulation_3< Kernel, Tds > Structure; ///< The triangulation that spatially orders the point set.
|
||||
typedef Alpha_shape_3< Structure > Shape; ///< The structure that identifies the triangles in the surface.
|
||||
|
||||
/// Default constructor.
|
||||
Shape_type() {}
|
||||
|
|
@ -67,7 +69,7 @@ public:
|
|||
* \param squared_radius the squared scale parameter of the shape.
|
||||
* \return a pointer to the new shape.
|
||||
*
|
||||
* Note: this class does not take responsibility for destroying
|
||||
* Note: Shape_type does not take responsibility for destroying
|
||||
* the object after use.
|
||||
*/
|
||||
Shape* construct( Shape* shape, const Scalar& squared_radius ) {
|
||||
|
|
@ -77,41 +79,43 @@ public:
|
|||
|
||||
/// Construct a new shape.
|
||||
/** \tparam InputIterator an interator over the points of the shape.
|
||||
* It should point to a CGAL::Point_3<Kernel>.
|
||||
* The iterator should point to a model of CGAL::Point_3<Kernel>.
|
||||
* \param start an iterator to the first point of the shape.
|
||||
* \param end a past-the-end iterator for the points of the shape.
|
||||
* \param squared_radius the squared scale parameter of the shape.
|
||||
* \return a pointer to the new shape.
|
||||
*
|
||||
* Note: this class does not take responsibility for destroying
|
||||
* Note: Shape_type does not take responsibility for destroying
|
||||
* the object after use.
|
||||
*/
|
||||
template < class InputIterator >
|
||||
Shape* construct( InputIterator start, InputIterator end, const Scalar& squared_radius ) {
|
||||
return new Shape( boost::make_transform_iterator( start, Auto_count<Point>() ),
|
||||
boost::make_transform_iterator( end, Auto_count<Point>() ),
|
||||
return new Shape( boost::make_transform_iterator( start, PointIndex() ),
|
||||
boost::make_transform_iterator( end, PointIndex() ),
|
||||
squared_radius, Shape::GENERAL );
|
||||
}
|
||||
}; // class Shape_type
|
||||
|
||||
/// The yype for the shape of a set of points with fixed scale.
|
||||
/** \tparam Kernel the geometric traits class. It should be a model of
|
||||
// The type for the shape of a set of points with fixed scale.
|
||||
/* \tparam Kernel the geometric traits class. It should be a model of
|
||||
* DelaunayTriangulationTraits_3.
|
||||
*/
|
||||
template < class Kernel >
|
||||
class Shape_type < Kernel, CGAL::Tag_true > {
|
||||
typedef typename Kernel::FT Scalar;
|
||||
class Shape_type < Kernel, Tag_true > {
|
||||
typedef typename Kernel::FT Scalar;
|
||||
typedef typename Kernel::Point_3 Point;
|
||||
typedef internal::Auto_count<Point> PointIndex;
|
||||
|
||||
typedef CGAL::Fixed_alpha_shape_vertex_base_3< Kernel,
|
||||
CGAL::Triangulation_vertex_base_with_info_3< unsigned int, Kernel > > Vb;
|
||||
typedef CGAL::Fixed_alpha_shape_cell_base_3< Kernel,
|
||||
CGAL::Triangulation_cell_base_with_info_3< unsigned int, Kernel > > Cb;
|
||||
typedef Triangulation_vertex_base_with_info_3< unsigned int, Kernel > Vb;
|
||||
typedef Fixed_alpha_shape_vertex_base_3< Kernel, Vb > aVb;
|
||||
typedef Triangulation_cell_base_with_info_3< unsigned int, Kernel > Cb;
|
||||
typedef Fixed_alpha_shape_cell_base_3< Kernel, Cb > aCb;
|
||||
|
||||
typedef CGAL::Triangulation_data_structure_3<Vb,Cb> Tds;
|
||||
typedef Triangulation_data_structure_3<aVb,aCb> Tds;
|
||||
|
||||
public:
|
||||
typedef CGAL::Delaunay_triangulation_3< Kernel, Tds > Structure; ///< The triangulation that spatially orders the point set.
|
||||
typedef CGAL::Fixed_alpha_shape_3< Structure > Shape; ///< The structure that identifies the triangles in the surface.
|
||||
typedef Delaunay_triangulation_3< Kernel, Tds > Structure; ///< The triangulation that spatially orders the point set.
|
||||
typedef Fixed_alpha_shape_3< Structure > Shape; ///< The structure that identifies the triangles in the surface.
|
||||
|
||||
/// Default constructor.
|
||||
Shape_type() {}
|
||||
|
|
@ -123,7 +127,7 @@ public:
|
|||
* \param squared_radius the squared scale parameter of the shape.
|
||||
* \return a pointer to the new shape.
|
||||
*
|
||||
* Note: this class does not take responsibility for destroying
|
||||
* Note: Shape_type does not take responsibility for destroying
|
||||
* the object after use.
|
||||
*/
|
||||
Shape* construct( Shape* shape, const Scalar& squared_radius ) {
|
||||
|
|
@ -133,19 +137,19 @@ public:
|
|||
|
||||
/// Construct a new shape.
|
||||
/** \tparam InputIterator an interator over the points of the shape.
|
||||
* It should point to a CGAL::Point_3<Kernel>.
|
||||
* The iterator should point to a model of CGAL::Point_3<Kernel>.
|
||||
* \param start an iterator to the first point of the shape.
|
||||
* \param end a past-the-end iterator for the points of the shape.
|
||||
* \param squared_radius the squared scale parameter of the shape.
|
||||
* \return a pointer to the new shape.
|
||||
*
|
||||
* Note: this class does not take responsibility for destroying
|
||||
* Note: Shape_type does not take responsibility for destroying
|
||||
* the object after use.
|
||||
*/
|
||||
template < class InputIterator >
|
||||
Shape* construct( InputIterator start, InputIterator end, const Scalar& squared_radius ) {
|
||||
return new Shape( boost::make_transform_iterator( start, Auto_count<Point>() ),
|
||||
boost::make_transform_iterator( end, Auto_count<Point>() ),
|
||||
return new Shape( boost::make_transform_iterator( start, PointIndex() ),
|
||||
boost::make_transform_iterator( end, PointIndex() ),
|
||||
squared_radius );
|
||||
}
|
||||
}; // class Shape_type
|
||||
|
|
|
|||
|
|
@ -32,7 +32,8 @@ namespace internal {
|
|||
/** \tparam T the type of object to count.
|
||||
*/
|
||||
template < class T >
|
||||
class Auto_count: public std::unary_function< const T&, std::pair< T, unsigned int > > {
|
||||
class Auto_count
|
||||
: public std::unary_function< const T&, std::pair< T, unsigned int > > {
|
||||
mutable unsigned int i;
|
||||
public:
|
||||
///< Start a new count.
|
||||
|
|
|
|||
|
|
@ -25,7 +25,7 @@ namespace CGAL {
|
|||
namespace internal {
|
||||
|
||||
/// A general system parameters definition.
|
||||
template<int sz> struct _EnvDef{
|
||||
template< int sz > struct _EnvDef{
|
||||
typedef void s_ptr_type; ///< The pointer size type.
|
||||
typedef void ptr_type; ///< The pointer type.
|
||||
}; // struct _EnvDef
|
||||
|
|
@ -44,7 +44,7 @@ template<> struct _EnvDef<8> {
|
|||
|
||||
/// The current system configuration.
|
||||
struct _Env {
|
||||
typedef _EnvDef< sizeof(void*) > _DEF; ///< The system definitions.
|
||||
typedef _EnvDef< sizeof( void* ) > _DEF; ///< The system definitions.
|
||||
|
||||
typedef _DEF::s_ptr_type s_ptr_type; ///< The pointer size type.
|
||||
typedef _DEF::ptr_type ptr_type; ///< The pointer type.
|
||||
|
|
@ -54,7 +54,7 @@ struct _Env {
|
|||
static const bool is_x64 = ( ptr_size == 64 ); ///< Whether the system is 64-bit.
|
||||
}; // struct _Env
|
||||
|
||||
typedef _Env _ENV;
|
||||
typedef _Env _ENVIRONMENT; ///< The system environment.
|
||||
|
||||
} // namespace internal
|
||||
|
||||
|
|
|
|||
|
|
@ -71,17 +71,17 @@ namespace CGAL {
|
|||
* \sa Scale_space_transform.
|
||||
* \sa Surface_mesher.
|
||||
*/
|
||||
template < class Kernel, class Fixed_shape = CGAL::Tag_true, class Shells = CGAL::Tag_true >
|
||||
template < class Kernel, class Fixed_shape = Tag_true, class Shells = Tag_true >
|
||||
class Scale_space_surface_reconstructer_3 {
|
||||
// Searching for neighbors.
|
||||
typedef typename CGAL::Search_traits_3< Kernel > Search_traits;
|
||||
typedef typename CGAL::Orthogonal_k_neighbor_search< Search_traits >
|
||||
typedef typename Search_traits_3< Kernel > Search_traits;
|
||||
typedef typename Orthogonal_k_neighbor_search< Search_traits >
|
||||
Static_search;
|
||||
typedef typename CGAL::Orthogonal_incremental_neighbor_search< Search_traits >
|
||||
typedef typename Orthogonal_incremental_neighbor_search< Search_traits >
|
||||
Dynamic_search;
|
||||
typedef typename Dynamic_search::Tree Search_tree;
|
||||
|
||||
typedef CGAL::Random Random;
|
||||
typedef Random Random;
|
||||
|
||||
// Constructing the surface.
|
||||
typedef Shape_type< Kernel, Fixed_shape > Shape_generator;
|
||||
|
|
@ -213,7 +213,7 @@ public:
|
|||
if (samples >= (_tree.size() - handled) || _generator.get_double() < double(samples - checked) / (_tree.size() - handled)) {
|
||||
// The neighborhood should contain the point itself as well.
|
||||
Static_search search( _tree, *it, _mean_neighbors+1 );
|
||||
radius += CGAL::sqrt( CGAL::to_double( squared_distance( *it, ( search.end()-1 )->first ) ) );
|
||||
radius += sqrt( to_double( squared_distance( *it, ( search.end()-1 )->first ) ) );
|
||||
++checked;
|
||||
}
|
||||
++handled;
|
||||
|
|
@ -256,32 +256,35 @@ public:
|
|||
* \param iterations the number of iterations to perform.
|
||||
*/
|
||||
void smooth_scale_space(unsigned int iterations = 1) {
|
||||
|
||||
typedef std::vector<unsigned int> CountVec;
|
||||
typedef std::vector< unsigned int > CountVec;
|
||||
typedef typename std::map<Point, size_t> PIMap;
|
||||
|
||||
typedef Eigen::Matrix<double, 3, Eigen::Dynamic> Matrix3D;
|
||||
typedef Eigen::Array<double, 1, Eigen::Dynamic> Array1D;
|
||||
typedef Eigen::Matrix3d Matrix3;
|
||||
typedef Eigen::Vector3d Vector3;
|
||||
typedef Eigen::SelfAdjointEigenSolver<Matrix3> EigenSolver;
|
||||
typedef Eigen::Matrix<double, 3, Eigen::Dynamic> EMatrix3D;
|
||||
typedef Eigen::Array<double, 1, Eigen::Dynamic> EArray1D;
|
||||
typedef Eigen::Matrix3d EMatrix3;
|
||||
typedef Eigen::Vector3d EVector3;
|
||||
typedef Eigen::SelfAdjointEigenSolver<EMatrix3> EigenSolver;
|
||||
|
||||
typedef _ENV::s_ptr_type p_size_t;
|
||||
typedef internal::_ENVIRONMENT::s_ptr_type p_size_t;
|
||||
|
||||
// This method must be called after filling the point collection.
|
||||
CGAL_assertion(!_tree.empty());
|
||||
|
||||
|
||||
if( !has_squared_mean_neighborhood() )
|
||||
estimate_mean_neighborhood( _mean_neighbors, _samples );
|
||||
|
||||
// In order to reduce memory consumption, we maintain two data structures:
|
||||
// a local tree (because openMP cannot work on global members), and
|
||||
// a vector for the points after smoothing.
|
||||
Pointset points;
|
||||
points.assign( _tree.begin(), _tree.end() );
|
||||
_tree.clear();
|
||||
Search_tree tree;
|
||||
|
||||
// Construct a search tree of the points.
|
||||
// Note that the tree has to be local for openMP.
|
||||
for (unsigned int iter = 0; iter < iterations; ++iter) {
|
||||
Search_tree tree( points.begin(), points.end() );
|
||||
tree.insert( points.begin(), points.end() );
|
||||
if(!tree.is_built())
|
||||
tree.build();
|
||||
|
||||
|
|
@ -295,10 +298,10 @@ public:
|
|||
for (p_size_t i = 0; i < count; ++i) {
|
||||
// Iterate over the neighbors until the first one is found that is too far.
|
||||
Dynamic_search search(tree, points[i]);
|
||||
for (Dynamic_search::iterator nit = search.begin(); nit != search.end() && compare(points[i], nit->first, squared_radius) != CGAL::LARGER; ++nit)
|
||||
for (Dynamic_search::iterator nit = search.begin(); nit != search.end() && compare(points[i], nit->first, squared_radius) != LARGER; ++nit)
|
||||
++neighbors[i];
|
||||
}
|
||||
|
||||
|
||||
// Construct a mapping from each point to its index.
|
||||
PIMap indices;
|
||||
p_size_t index = 0;
|
||||
|
|
@ -315,24 +318,25 @@ public:
|
|||
|
||||
// Collect the vertices within the ball and their weights.
|
||||
Dynamic_search search(tree, points[i]);
|
||||
Matrix3D pts(3, neighbors[i]);
|
||||
Array1D wts(1, neighbors[i]);
|
||||
int column = 0;
|
||||
for (Dynamic_search::iterator nit = search.begin(); nit != search.end() && compare(points[i], nit->first, squared_radius) != CGAL::LARGER; ++nit, ++column) {
|
||||
pts(0, column) = CGAL::to_double(nit->first[0]);
|
||||
pts(1, column) = CGAL::to_double(nit->first[1]);
|
||||
pts(2, column) = CGAL::to_double(nit->first[2]);
|
||||
EMatrix3D pts(3, neighbors[i]);
|
||||
EArray1D wts(1, neighbors[i]);
|
||||
unsigned int column = 0;
|
||||
for (Dynamic_search::iterator nit = search.begin(); nit != search.end() && column < neighbors[i]; ++nit, ++column) {
|
||||
pts(0, column) = to_double(nit->first[0]);
|
||||
pts(1, column) = to_double(nit->first[1]);
|
||||
pts(2, column) = to_double(nit->first[2]);
|
||||
wts(column) = 1.0 / neighbors[indices[nit->first]];
|
||||
}
|
||||
CGAL_assertion( column == neighbors[i] );
|
||||
|
||||
// Construct the barycenter.
|
||||
Vector3 bary = (pts.array().rowwise() * wts).rowwise().sum() / wts.sum();
|
||||
EVector3 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 covariance = Matrix3::Zero();
|
||||
EMatrix3 covariance = EMatrix3::Zero();
|
||||
for (column = 0; column < pts.cols(); ++column)
|
||||
covariance += wts.matrix()(column) * pts.col(column) * pts.col(column).transpose();
|
||||
|
||||
|
|
@ -346,20 +350,22 @@ public:
|
|||
// 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()) {
|
||||
if (solver.eigenvectors().col(index).imag() != EVector3::Zero()) {
|
||||
// This should actually never happen!
|
||||
CGAL_assertion(false);
|
||||
continue;
|
||||
}
|
||||
Vector3 n = solver.eigenvectors().col(index).real();
|
||||
EVector3 n = solver.eigenvectors().col(index).real();
|
||||
|
||||
// The vertex is moved by projecting it onto the plane
|
||||
// through the barycenter and orthogonal to the Eigen vector with smallest Eigen value.
|
||||
Vector3 bv = Vector3(CGAL::to_double(points[i][0]), CGAL::to_double(points[i][1]), CGAL::to_double(points[i][2])) - bary;
|
||||
Vector3 per = bary + bv - (n.dot(bv) * n);
|
||||
EVector3 bv = EVector3(to_double(points[i][0]), to_double(points[i][1]), to_double(points[i][2])) - bary;
|
||||
EVector3 per = bary + bv - (n.dot(bv) * n);
|
||||
|
||||
points[i] = Point(per(0), per(1), per(2));
|
||||
}
|
||||
|
||||
tree.clear();
|
||||
}
|
||||
|
||||
_tree.insert( points.begin(), points.end() );
|
||||
|
|
|
|||
Loading…
Reference in New Issue