diff --git a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Doxyfile.in b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Doxyfile.in new file mode 100644 index 00000000000..38aaf28d221 --- /dev/null +++ b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Doxyfile.in @@ -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 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 new file mode 100644 index 00000000000..a49f134a7f9 --- /dev/null +++ b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/PackageDescription.txt @@ -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` +- `CGAL::Shape_type` + +### Auxiliary Classes ### + +- `CGAL::internal::Auto_count` +- `CGAL::internal::_Env` + +*/ + 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 new file mode 100644 index 00000000000..61800a4639e --- /dev/null +++ b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Scale_space_reconstruction_3.txt @@ -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 */ + diff --git a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/dependencies b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/dependencies new file mode 100644 index 00000000000..6fff58bf110 --- /dev/null +++ b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/dependencies @@ -0,0 +1,5 @@ +Manual +Kernel_23 +STL_Extension +Triangulation_3 +Alpha_shapes_3 diff --git a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/examples.txt b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/examples.txt new file mode 100644 index 00000000000..0c36fa93497 --- /dev/null +++ b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/examples.txt @@ -0,0 +1,3 @@ +/*! +\example Scale_space_reconstruction_3/scale_space.cpp +*/ diff --git a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/fig/knot_thumb.png b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/fig/knot_thumb.png new file mode 100644 index 00000000000..9242306e251 Binary files /dev/null and b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/fig/knot_thumb.png differ diff --git a/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/scale_space.cpp b/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/scale_space.cpp index 53a49352912..43e1390775d 100644 --- a/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/scale_space.cpp +++ b/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/scale_space.cpp @@ -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; diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Shape_type.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Shape_type.h index c6e079c8d94..a5231a53762 100644 --- a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Shape_type.h +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Shape_type.h @@ -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 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 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 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. + * The iterator should point to a model of CGAL::Point_3. * \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() ), - boost::make_transform_iterator( end, Auto_count() ), + 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 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 Tds; + typedef Triangulation_data_structure_3 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. + * The iterator should point to a model of CGAL::Point_3. * \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() ), - boost::make_transform_iterator( end, Auto_count() ), + return new Shape( boost::make_transform_iterator( start, PointIndex() ), + boost::make_transform_iterator( end, PointIndex() ), squared_radius ); } }; // class Shape_type diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/internal/Auto_count.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/internal/Auto_count.h index 3ff59254885..f549d16fbf1 100644 --- a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/internal/Auto_count.h +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/internal/Auto_count.h @@ -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. diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/internal/check3264.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/internal/check3264.h index 97a3739b3d0..fb99ebee4e2 100644 --- a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/internal/check3264.h +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/internal/check3264.h @@ -25,7 +25,7 @@ namespace CGAL { namespace internal { /// A general system parameters definition. -template 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 diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstructer_3.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstructer_3.h index 7c60c9ba00c..879025ca7c2 100644 --- a/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstructer_3.h +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstructer_3.h @@ -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 CountVec; + typedef std::vector< unsigned int > CountVec; typedef typename std::map PIMap; - typedef Eigen::Matrix Matrix3D; - typedef Eigen::Array Array1D; - typedef Eigen::Matrix3d Matrix3; - typedef Eigen::Vector3d Vector3; - typedef Eigen::SelfAdjointEigenSolver EigenSolver; + typedef Eigen::Matrix EMatrix3D; + typedef Eigen::Array EArray1D; + typedef Eigen::Matrix3d EMatrix3; + typedef Eigen::Vector3d EVector3; + typedef Eigen::SelfAdjointEigenSolver 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() );