diff --git a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Mean_neighborhood_ball.h b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Mean_neighborhood_ball.h deleted file mode 100644 index b55d585f069..00000000000 --- a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Mean_neighborhood_ball.h +++ /dev/null @@ -1,113 +0,0 @@ -//A method that computes the size of a mean neighborhood. -//Copyright (C) 2013 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 MEAN_NEIGHBORHOOD_BALL -#define MEAN_NEIGHBORHOOD_BALL - -#include -#include -#include - -/// Compute the mean neighborhood size. -/** Each neighborhood size is expressed as the radius of the smallest - * ball centered on a sample point such that the ball contains at least - * a specified number of points. - * - * The mean neighborhood size is then the mean of these radii, taken over - * a number of point samples. - * \tparam Kernel the geometric traits class. This class specifies, amongst - * other things, the number types and predicates used. - */ -template < typename Kernel > -class Mean_neighborhood_ball { - typedef typename Kernel::FT Scalar; ///< The number type. - - typedef typename Kernel::Point_3 Point; ///< The point type. - - typedef Scalar result_type; ///< The type for expressing the mean neigborhood size. - -public: - /// The constructor. - /** \param nn the number of nearest neighbors that the neigborhood should contain. - * \param smp the number of smaple points used to estimate the neighborhood size. - */ - Mean_neighborhood_ball(unsigned int nn = 30, unsigned int smp = 200); - - /// Accessor for the number of neighbors required in the neighborhood. - /** \return the number of neighbors the neighborhood must contain. - */ - unsigned int get_neighbors() const; - - /// Accessor for the number of sample points used to estimate the neighborhood size. - /** \return the number of sample points used to estimate the neighborhood size. - */ - unsigned int get_samples() const; - - /// Mutator for the number of neighbors required in the neighborhood. - /** \param nn the number of neighbors required in the neighborhood. - */ - void set_neighbors( unsigned int nn ); - - /// Mutator for the number of sample points used to estimate the neighborhood size. - /** \param smp the number of sample points used to estimate the neighborhood size. - */ - void set_samples( unsigned int smp ); - - /// Construct a search tree from a collection of points. - /** This tree must be constructed before handling sample points. - * \tparam InputIterator an iterator over a collection of points. - * The iterator must point to a Point type. - * \param start an iterator to the first point of the collection. - * \param end a past-the-end iterator for the point collection. - * \sa compute(). - * \sa operator()(InputIterator start, InputIterator end). - */ - template < class InputIterator > - void constructTree(InputIterator start, InputIterator end); - - /// Incorporate the neighborhood of a sample point into the estimate. - /** Note that a total of get_samples() points should be handled. - * \param p the point to incorporate in the estimate. - * \sa compute(). - */ - void handlePoint(const Point& p); - - /// Estimate the mean neighborhood size based on a number of sample points. - /** It is assumed that the seach tree has already been constructed. - * \return the estimated mean neighborhood size. - * \sa handlePoint(const Point& p). - * \sa operator()(InputIterator start, InputIterator end). - */ - result_type compute(); - - /// Estimate the mean neighborhood size of a collection of points. - /** This method is equivalent to running [constructTree(start, end)](\ref constructTree) - * followed by compute(). - * \tparam InputIterator an iterator over a collection of points. - * The iterator must point to a Point type. - * \param start an iterator to the first point of the collection. - * \param end a past-the-end iterator for the point collection. - * \sa constructTree(InputIterator start, InputIterator end). - * \sa compute(). - */ - template < class InputIterator > - result_type operator()(InputIterator start, InputIterator end); -}; // class Mean_neighborhood_ball - -#endif MEAN_NEIGHBORHOOD_BALL diff --git a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Scale_space_surface_constructer.h b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Scale_space_surface_constructer.h deleted file mode 100644 index 0b3e4feb532..00000000000 --- a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Scale_space_surface_constructer.h +++ /dev/null @@ -1,110 +0,0 @@ -//A method to construct a surface. -//Copyright (C) 2013 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 SCALE_SPACE_SURFACE_CONSTRUCTER -#define SCALE_SPACE_SURFACE_CONSTRUCTER - -#include -#include -#include - -#include "Mean_neighborhood_ball.h" -#include "Scale_space_transform.h" -#include "Surface_mesher.h" - -/// Compute a smoothed surface mesh from a collection of points. -/** This is a convenience class for using Mean_neighborhood_ball, - * Scale_space_transform, and Surface_mesher. - * - * An appropriate neighborhood size is estimated, followed by a - * number of smoothing iterations. Finally, the surface of the - * smoothed point set is computed. - * - * The order of the point set remains the same, meaning that - * the smoothed surface can be transposed back to its unsmoothed - * version by overwriting the smoothed point collection with its - * unsmoothed version. - * \tparam Kernel the geometric traits class. This class - * specifies, amongst other things, the number types and - * predicates used. - * \tparam Vb the vertex base used in the internal data - * structure that spatially orders the point set. - * \tparam Cb the cell base used in the internal data - * structure that spatially orders the point set. - * \sa Mean_neighborhood_ball. - * \sa Scale_space_transform. - * \sa Surface_mesher. - */ -template < typename Kernel, - typename Vb = CGAL::Triangulation_vertex_base_3, - typename Cb = CGAL::Triangulation_cell_base_3 > -class Scale_space_surface_constructer { -public: - typedef typename Kernel::FT Scalar; ///< The number type. - - typedef typename Kernel::Point_3 Point; ///< The point type. - typedef typename Kernel::Triangle_3 Triangle; ///< The triangle type. - - typedef typename Surface_mesher::Vertex_handle Vertex_handle; ///< A handle to access the vertices. - typedef typename Surface_mesher::Facet Facet; ///< A facet of the triangulation. - - typedef typename std::vector PointCollection; ///< A collection of points. - -public: - /// The constructor. - /** \param iterations the number of smoothing iterations to perform. - */ - Scale_space_surface_constructer(unsigned int iterations = 1); - - /// Mutator for the number of smoothing iterations. - /** \param iterations the number of smoothing iterations to perform. - */ - void set_iterations(unsigned int iterations); - - /// Accessor for the smoothed points. - /** \return the current collection of smoothed points. - * Note that this collection may change as smoothing - * iterations are performed. - */ - const PointCollection& moved() const; - - /// Compute a smoothed surface mesh from a collection of points. - /** \tparam InputIterator an iterator over the point collection. - * The iterator must point to a Point type. - * \tparam OutputIterator an output iterator for the surface - * triangles. The iterator must point to a Triangle type. - * \param start an iterator to the first point of the collection. - * \param end a past-the-end iterator for the point collection. - * \param out an iterator to place to output the triangles. - */ - template < class InputIterator, class OutputIterator > - OutputIterator operator()(InputIterator start, InputIterator end, OutputIterator out); - - /// Get the vertices of a facet ordered to point towards the outside of the surface. - /** This orientation is expressed using the 'right-hand rule' - * on the ordered vertices of the facet. - * \param f a facet of the data structure as seen from the outside. - * \param v0 the first vertex of the facet. - * \param v1 the second vertex of the facet. - * \param v2 the third vertex of the facet. - */ - void ordered_vertices(const Facet& f, Vertex_handle& v0, Vertex_handle& v1, Vertex_handle& v2); -}; // class Scale_space_surface_constructer - -#endif // SCALE_SPACE_SURFACE_CONSTRUCTER \ No newline at end of file diff --git a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Scale_space_transform.h b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Scale_space_transform.h deleted file mode 100644 index f046f21564d..00000000000 --- a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Scale_space_transform.h +++ /dev/null @@ -1,139 +0,0 @@ -//A method to transform a point set. -//Copyright (C) 2013 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 SCALE_SPACE_PERTURB -#define SCALE_SPACE_PERTURB - -#include - -#include -#include - -#include -#include - -#include - -#include "check3264.h" - - -/// Transform a point set to a rougher scale. -/** This transfor is achieved by projecting each point onto - * the best planar fit of the weighted local neighborhood. - * \tparam Kernel the geometric traits class. This class specifies, amongst - * other things, the number types and predicates used. - */ -template < typename Kernel > -class Scale_space_transform { - typedef typename Kernel::FT Scalar; ///< The number type. - - typedef typename Kernel::Point_3 Point; ///< The point type. - typedef typename std::vector Collection; ///< A collection of points. - -public: - /// The constructor. - /** \param r the radius of the local smoothing neighborhood. - */ - Scale_space_transform(const Scalar& r); - - /// Accessor for the radius of the local smoothing neighborhood. - /** \return the radius of the local smoothing neighborhood. - */ - Scalar get_radius2() const; - - /// Mutator for the radius of the local smoothing neighborhood. - /** \param r the radius of the local smoothing neighborhood. - */ - void set_radius2(const Scalar& r); - - /// Assign a collection of points to be smoothed. - /** \tparam InputIterator an iterator over a collection of points. - * The iterator must point to a Point type. - * \param start an iterator to the first point of the collection. - * \param end a past-the-end iterator for the point collection. - * \sa operator()(InputIterator start, InputIterator end). - * \sa operator()(InputIterator start, InputIterator end, unsigned int iterations). - */ - template < class InputIterator > - void assign(InputIterator start, InputIterator end); - - /// Accessor for the collection of points to be smoothed. - /** \return the collection of points to be smoothed. - */ - Collection& get_collection(); - - /// Accessor for the collection of points to be smoothed. - /** \return the collection of points to be smoothed. - */ - const Collection& get_collection() const; - - /// Compute one iteration of scale-space transforming. - /** The result of this iteration is stored in the point collection. - * If earlier iterations have been computed, calling iterate() - * will add the next iteration. - * \sa iterate(unsigned int iterations). - * \sa get_collection(). - * \sa get_collection() const. - */ - void iterate(); - - /// Compute a number of iterations of scale-space transforming. - /** The result of this iteration is stored in the point collection. - * \param iterations the number of iterations to perform. - * \sa iterate(). - * \sa get_collection(). - * \sa get_collection() const. - */ - void iterate(unsigned int iterations); - - /// Compute one transform iteration on a collection of points. - /** This method is equivalent to running [assign(start, end)](\ref assign) - * followed by iterate(). - * \tparam InputIterator an iterator over a collection of points. - * The iterator must point to a Point type. - * \param start an iterator to the first point of the collection. - * \param end a past-the-end iterator for the point collection. - * \sa operator()(InputIterator start, InputIterator end, unsigned int iterations). - * \sa assign(InputIterator start, InputIterator end). - * \sa iterate(). - */ - template < class InputIterator > - void operator()(InputIterator start, InputIterator end); - - /// Compute a number of transform iterations on a collection of points. - /** This method is equivalent to running [assign(start, end)](\ref assign) - * followed by [iterate(iterations)](\ref iterate). - * \tparam InputIterator an iterator over a collection of points. - * The iterator must point to a Point type. - * \param start an iterator to the first point of the collection. - * \param end a past-the-end iterator for the point collection. - * \param iterations the number of iterations to perform. - * \sa operator()(InputIterator start, InputIterator end). - * \sa assign(InputIterator start, InputIterator end). - * \sa iterate(unsigned int iterations). - */ - template < class InputIterator > - void operator()(InputIterator start, InputIterator end, unsigned int iterations) { - // Perturb the points for a number of iterations. - _points.assign(start, end); - iterate(iterations); - } -}; // class Scale_space_transform - -#endif // SCALE_SPACE_PERTURB \ No newline at end of file diff --git a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Surface_mesher.h b/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Surface_mesher.h deleted file mode 100644 index 20de1b37c36..00000000000 --- a/Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Surface_mesher.h +++ /dev/null @@ -1,276 +0,0 @@ -//The surface meshing algorithm implemented using the alpha-shape. -//Copyright (C) 2013 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 BALL_PIVOTING -#define BALL_PIVOTING - -#include - -#define BALL_PIVOTING_FIXED -#define BALL_PIVOTING_CONNECTED - -#ifdef BALL_PIVOTING_FIXED -#include -#include -#include -#else -#include -#include -#include -#endif - -#include "Flagged_triangulation_cell_base_3.h" - -/// Compute a surface mesh from a collection of points. -/** The surface mesh indicates the boundary of a sampled - * object. The point sample must be dense enough, such - * that any region that fits an empty ball with a given - * radius can be assumed to be ouside the object. - * - * The resulting surface need not be manifold and in - * many cases where the points sample the surface of - * an object, the surface computed will contain both - * an 'outward-facing' and an 'inward-facing' surface. - * - * However, this surface cannot have holes and its - * triangles are all oriented towards the same side of - * the surface. This orientation is expressed using the - * 'right-hand rule' on the ordered vertices of each - * triangle. - * - * The computed triangles will be returned ordered such - * that a connected 'shell' contains a set of consecutive - * triangles in the output. For this purpose, a 'shell' - * is a maximal collection of triangles that are connected - * and locally facing towards the same side of the surface. - * \tparam Kernel the geometric traits class. This class - * specifies, amongst other things, the number types and - * predicates used. - * \tparam Vb the vertex base used in the internal data - * structure that spatially orders the point set. - * \tparam Cb the cell base used in the internal data - * structure that spatially orders the point set. - */ -template < typename Kernel, - typename Vb = CGAL::Triangulation_vertex_base_3, - typename Cb = CGAL::Triangulation_cell_base_3 > -class Surface_mesher { -public: - typedef typename Kernel::FT Scalar; ///< The number type. - - // Convert the DT types to be in an alpha-shape. - /// The vertex type. -#ifdef BALL_PIVOTING_FIXED - typedef CGAL::Fixed_alpha_shape_vertex_base_3 AVb; -#else - typedef CGAL::Alpha_shape_vertex_base_3 AVb; -#endif - - /// The cell type. - typedef CGAL::Flagged_triangulation_cell_base_3 FCb; -#ifdef BALL_PIVOTING_FIXED - typedef CGAL::Fixed_alpha_shape_cell_base_3 ACb; -#else - typedef CGAL::Alpha_shape_cell_base_3 ACb; -#endif - - // Triangulation. - typedef CGAL::Triangulation_data_structure_3 Tds; ///< The data structure that stores the point set. - typedef CGAL::Delaunay_triangulation_3 ADT; ///< The base triangulation that spatially orders the point set. - /// The structure that identifies the triangles in the surface. -#ifdef BALL_PIVOTING_FIXED - typedef CGAL::Fixed_alpha_shape_3 AS; -#else - typedef CGAL::Alpha_shape_3 AS; -#endif - -public: - typedef ADT Triangulation; ///< The triangulation that spatially orders the point set. - typedef AS Alpha_shape; ///< The structure that identifies the triangles in the surface. - - typedef typename AS::Vertex_handle Vertex_handle; ///< A handle to access the vertices. - typedef typename AS::Cell_handle Cell_handle; ///< A handle to access the cells. - typedef typename AS::Facet Facet; ///< A facet of the triangulation. - - typedef typename Kernel::Point_3 Point; ///< The point type. - typedef typename Kernel::Triangle_3 Triangle; ///< The triangle type. - - /// The constructor that starts the structure from scratch. - /** \param r2 the squared radius of the ball used to indicate - * regions outside the shape. - * \sa Surface_mesher(const Triangulation& tr, const Scalar& r2). - */ - Surface_mesher(const Scalar& r2); - - /// The constructor that uses an existing spatial structure. - /** \param tr the spatial structure on the point sample. - * \param r2 the squared radius of the ball used to indicate - * regions outside the shape. - * \sa Surface_mesher(const Scalar& r2). - */ - Surface_mesher(const Triangulation& tr, const Scalar& r2); - - /// Construct a new spatial structure on a set of sample points. - /** \tparam InputIterator an iterator over the point sample. - * The iterator must point to a Point type. - * \param start an iterator to the first point of the collection. - * \param end a past-the-end iterator for the point collection. - * \sa is_constructed(). - */ - template < class InputIterator > - void construct_shape(InputIterator start, InputIterator end); - - /// Check whether the spatial structure has been constructed. - /** \return true if the structure exists and false otherwise. - * \sa construct_shape(InputIterator start, InputIterator end). - */ - bool is_constructed() const; - - /// Clear the spatial data structure. - void clear(); - - /// Accessor for the structure indicating the triangles of the mesh. - /** \return the structure indicating the triangles of the mesh. - */ - Alpha_shape* shape(); - - /// Accessor for the ball indicating regions outside the shape. - /** \return the squared radius of the ball used to indicate - * regions outside the shape. - */ - Scalar get_radius2() const; - - /// Mutator for the ball indicating regions outside the shape. - /** \param r2 the squared radius of the ball used to indicate - * regions outside the shape. - */ - void set_radius2(const Scalar& r2); - - /// Collect the triangles of one shell of the surface. - /** A shell is a maximal collection of triangles that are - * connected and locally facing towards the same side of the - * surface. - * \tparam OutputIterator an output iterator for a collection - * of triangles. The iterator must point to a Triangle type. - * \param c a cell touching a triangle of the shell from the - * outside of the object. - * \param li the index of the vertex of the cell opposite to - * the triangle touched. - * \param out an iterator to place to output the triangles. - * \sa collect_shell(const Facet& f, OutputIterator out). - * \sa collect_surface(OutputIterator out). - */ - template < class OutputIterator > - OutputIterator collect_shell(Cell_handle c, unsigned int li, OutputIterator out); - - /// Collect the triangles of one shell of the surface. - /** A shell is a maximal collection of triangles that are - * connected and locally facing towards the same side of the - * surface. - * \tparam OutputIterator an output iterator for a collection - * of triangles. The iterator must point to a Triangle type. - * \param f a facet of the shell as seen from the outside. - * \param out an iterator to place to output the triangles. - * \sa collect_shell(Cell_handle c, unsigned int li, OutputIterator out). - * \sa collect_surface(OutputIterator out). - */ - template < class OutputIterator > - OutputIterator collect_shell(const Facet& f, OutputIterator out); - - /// Collect the triangles of the complete surface. - /** These triangles are given ordered per shell. - * A shell is a maximal collection of triangles that are - * connected and locally facing towards the same side of the - * surface. - * \tparam OutputIterator an output iterator for a collection - * of triangles. The iterator must point to a Triangle type. - * \param out an iterator to place to output the triangles. - * \sa collect_shell(const Facet& f, OutputIterator out). - * \sa collect_shell(Cell_handle c, unsigned int li, OutputIterator out). - */ - template < class OutputIterator > - OutputIterator collect_surface(OutputIterator out); - - /// Construct the surface trinagles from a set of sample points. - /** These triangles are given ordered per shell. - * A shell is a maximal collection of triangles that are - * connected and locally facing towards the same side of the - * surface. - * - * This method is equivalent to running [construct_shape(start, end)](\ref construct_shape) - * followed by [collect_surface(out)](\ref collect_surface). - * \tparam InputIterator an iterator over the point sample. - * The iterator must point to a Point type. - * \tparam OutputIterator an output iterator for a collection - * of triangles. The iterator must point to a Triangle type. - * \param start an iterator to the first point of the collection. - * \param end a past-the-end iterator for the point collection. - * \param out an iterator to place to output the triangles. - * \sa construct_shape(InputIterator start, InputIterator end). - * \sa collect_surface(OutputIterator out). - */ - template < class InputIterator, class OutputIterator > - OutputIterator operator()(InputIterator start, InputIterator end, OutputIterator out); - - /// Locate a vertex in the spatial data structure. - /** \param p the point of the vertex. - * \param v the vertex in the data structure. - * \param hint where to start looking for the vertex. - * Provinding a hint near the vertex will greatly speed - * up location. - * \return whether the point actually has a vertex in the structure. - */ - bool locate_vertex(const Point& p, Vertex_handle& v, Cell_handle hint = Cell_handle()) const; - - /// Get the vertices of a facet ordered to point towards the outside of the surface. - /** This orientation is expressed using the 'right-hand rule' - * on the ordered vertices of the facet. - * \param f a facet of the data structure as seen from the outside. - * \param v0 the first vertex of the facet. - * \param v1 the second vertex of the facet. - * \param v2 the third vertex of the facet. - */ - void ordered_vertices(const Facet& f, Vertex_handle& v0, Vertex_handle& v1, Vertex_handle& v2); - - /// Get a triangle oriented towards the outside of the surface. - /** This orientation is expressed using the 'right-hand rule' - * on the ordered vertices of the triangle. - * \param c a cell touching the triangle from the outside of - * the object. - * \param li the index of the vertex of the cell opposite to - * the triangle touched. - * \return the triangle oriented towards the outside of the - * surface. - * \sa oriented_triangle(const Facet& f) const. - */ - Triangle oriented_triangle(Cell_handle c, unsigned int li) const; - - /// Get a triangle oriented towards the outside of the surface. - /** This orientation is expressed using the 'right-hand rule' - * on the ordered vertices of the triangle. - * \param f a facet of the data structure indicating the triangle - * as seen from the outside. - * \return the triangle oriented towards the outside of the - * surface. - * \sa oriented_triangle(Cell_handle c, unsigned int li) const. - */ - Triangle oriented_triangle(const Facet& f) const; -}; // class Surface_mesher - -#endif // BALL_PIVOTING diff --git a/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/CMakeLists.txt b/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/CMakeLists.txt new file mode 100644 index 00000000000..5a60bf81c35 --- /dev/null +++ b/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/CMakeLists.txt @@ -0,0 +1,45 @@ +project( Scale_space ) + +# Check CMake version +cmake_minimum_required( VERSION 2.6.2 ) +if( "${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}" VERSION_GREATER 2.6 ) + if( "${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}.${CMAKE_PATCH_VERSION}" VERSION_GREATER 2.8.3 ) + cmake_policy( VERSION 2.8.4 ) + else() + cmake_policy( VERSION 2.6 ) + endif() +endif() + + +# Preprocessor directives +if( WIN32 ) + add_definitions( "/W3 /D_CRT_SECURE_NO_WARNINGS /wd4005 /wd4307 /wd4308 /nologo" ) +endif() +if( UNIX ) + add_definitions( "-fpermissive" ) +endif() + +# Enable openMP for parallel processing +if( WIN32 ) + set( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /openmp" ) + set( CMAKE_C_FLAGS "${CMAKE_C_FLAGS} /openmp" ) +endif() + + +find_package(CGAL QUIET COMPONENTS Core ) + +if ( CGAL_FOUND ) + + include( ${CGAL_USE_FILE} ) + + include( CGAL_CreateSingleSourceCGALProgram ) + + include_directories (BEFORE "../../include" $ENV{EIGEN3_DIR}/ ) + + create_single_source_cgal_program( "scale_space.cpp" ) + +else() + + message(STATUS "This program requires the CGAL library, and will not be compiled.") + +endif() diff --git a/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/quad.h b/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/quad.h deleted file mode 100644 index d44585c6bb2..00000000000 --- a/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/quad.h +++ /dev/null @@ -1,60 +0,0 @@ -//A struct with four values. -//Copyright (C) 2013 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 SCALE_SPACE_QUAD -#define SCALE_SPACE_QUAD - -// A class that holds four values, possibly a color. -template < typename Scalar > -struct Quad { - Scalar first, second, third, fourth; - Quad(): first(0), second(0), third(0), fourth(0) {} - Quad(const Quad& q): first(q.first), second(q.second), third(q.third), fourth(q.fourth) {} - Quad(const Scalar& s1, const Scalar& s2, const Scalar& s3, const Scalar& s4): first(s1), second(s2), third(s3), fourth(s4) {} - inline Scalar x() const {return first;} - inline Scalar y() const {return second;} - inline Scalar z() const {return third;} - inline Scalar w() const {return fourth;} - inline Scalar r() const {return first;} - inline Scalar g() const {return second;} - inline Scalar b() const {return third;} - inline Scalar a() const {return fourth;} - inline bool operator<(const Quad& q) const { - if (first < q.first) return true; - if (q.first < first) return false; - if (second < q.second) return true; - if (q.second < second) return false; - if (third < q.third) return true; - if (q.third < third) return false; - return fourth < q.fourth; - } -}; // struct Clock - -template < typename Scalar > -std::ostream& operator<<(std::ostream& os, const Quad& q) { - os << q.first << " " << q.second << " " << q.third << " " << q.fourth; - return os; -} - -typedef Quad int4; -typedef Quad long4; -typedef Quad double4; -typedef Quad float4; - -#endif // SCALE_SPACE_QUAD \ No newline at end of file diff --git a/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/readFile.h b/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/readFile.h deleted file mode 100644 index 700fbb4d606..00000000000 --- a/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/readFile.h +++ /dev/null @@ -1,178 +0,0 @@ -//A general purpose file reader. -//Note that if we define more file formats, -//it may be beneficial to split this into multiple plugin classes/files. -//Copyright (C) 2013 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 READ_FILE -#define READ_FILE - -#include -#include - -/*#include -#include */ - -#include "quad.h" - -const unsigned int LINE_SIZE = 1024; - -void readLine(std::ifstream& fin, char* line) { - fin.getline(line, LINE_SIZE); - while (fin && line[0] == '#') { - // Comment line. - std::cout << "Comment: " << line << std::endl; - fin.getline(line, LINE_SIZE); - } -} - - -template < class Kernel, class PointOutput, class NormalOutput > -bool loadXYZN(char* file, PointOutput po, NormalOutput no) { - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::Vector_3 Normal; - - char line[LINE_SIZE]; - - std::ifstream fin(file); - if (!fin) return false; - - unsigned int num = 0; - - // Collect the points. - while (fin) { - readLine(fin, line); - - if (strlen(line) > 0) { - ++num; - - // Collect the point and normal. - double x, y, z; - double nx, ny, nz; - int a = sscanf(line,"%lf %lf %lf %lf %lf %lf", &x, &y, &z, &nx, &ny, &nz); - - if (a > 2) - *po++ = Point(x, y, z); - if (a > 5) - *no++ = Normal(nx, ny, nz); - } - } - - fin.close(); - - return true; -} - - -/*class VectorCollecter: public osg::ConstValueVisitor { - inline bool comp(double d1, double d2) const {return abs(d1 - d2) < 0.0001;} - -public: - VectorCollecter(): osg::ConstValueVisitor() {} - - inline void apply(const osg::Vec3d& vec) {_vec3 = osg::Vec3d(vec[0],vec[1],vec[2]);} - inline void apply(const osg::Vec3& vec) {_vec3 = osg::Vec3d(vec[0],vec[1],vec[2]);} - inline void apply(const osg::Vec4d& vec) {_vec4 = osg::Vec4d(vec[0],vec[1],vec[2],vec[3]);} - inline void apply(const osg::Vec4& vec) {_vec4 = osg::Vec4d(vec[0],vec[1],vec[2],vec[3]);} - inline void apply(const osg::Vec4ub& vec) {_vec4 = osg::Vec4d(vec[0]/256.f, vec[1]/256.f, vec[2]/256.f, vec[3]/256.f);} - - inline bool compareColor(const osg::Vec4d& c) {return comp(_vec4[0], c[0]) && comp(_vec4[1], c[1]) && comp(_vec4[2], c[2]);} - - osg::Vec3d _vec3; - osg::Vec4d _vec4; - - template < class ToType > - ToType convert3() const {return ToType(_vec3[0], _vec3[1], _vec3[2]);} - template < class ToType > - ToType convert4() const {return ToType(_vec4[0], _vec4[1], _vec4[2], _vec4[3]);} -}; - -template < class Kernel, typename Color, class PointOutput, class NormalOutput, class ColorOutput > -class PointCollecter: public osg::NodeVisitor { - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::Vector_3 Normal; - - PointOutput _po; - NormalOutput _no; - ColorOutput _co; -public: - PointCollecter(PointOutput po, NormalOutput no, ColorOutput co): osg::NodeVisitor(osg::NodeVisitor::TRAVERSE_ALL_CHILDREN), _po(po), _no(no), _co(co) {} - - void apply(osg::Geode& geode) { - VectorCollecter vc; - for (unsigned int i = 0; i < geode.getNumDrawables(); ++i) { - osg::Geometry* geom = geode.getDrawable(i)->asGeometry(); - - // Only collect points. - if (geom && geom->getNumPrimitiveSets() > 0 && geom->getPrimitiveSet(0)->getMode() == osg::PrimitiveSet::POINTS) { - // Collect locations. - osg::Array* vertices = geom->getVertexArray(); - for (unsigned int j = 0; j < vertices->getNumElements(); ++j) { - vertices->accept(j, vc); - *_po++ = vc.convert3(); - } - - // Collect normals. - osg::Array* normals = geom->getNormalArray(); - if (geom->getNormalBinding() == osg::Geometry::BIND_OVERALL) { - normals->accept(0, vc); - Normal n = vc.convert3(); - for (unsigned int j = 0; j < vertices->getNumElements(); ++j) - *_no++ = n; - } - else if (normals) { - for (unsigned int j = 0; j < normals->getNumElements() && j < vertices->getNumElements(); ++j) { - normals->accept(j, vc); - *_no++ = vc.convert3(); - } - for (unsigned int j = 0; j < (vertices->getNumElements() - normals->getNumElements()); ++j) { - *_no++ = Normal(0, 0, 0); - } - } - - // Collect colors. - osg::Array* colors = geom->getColorArray(); - if (geom->getColorBinding() == osg::Geometry::BIND_OVERALL) { - colors->accept(0, vc); - Color c = vc.convert4(); - for (unsigned int j = 0; j < vertices->getNumElements(); ++j) - *_co++ = c; - } - else if (colors) { - for (unsigned int j = 0; j < colors->getNumElements() && j < vertices->getNumElements(); ++j) { - colors->accept(j, vc); - *_co++ = vc.convert4(); - } - for (unsigned int j = 0; j < (vertices->getNumElements() - colors->getNumElements()); ++j) { - *_co++ = Color(0, 0, 0, 0); - } - } - } - } - } -}; - -template < class Kernel, typename Color, class PointOutput, class NormalOutput, class ColorOutput > -bool loadOSG(char* file, PointOutput po, NormalOutput no, ColorOutput co) { - osg::ref_ptr scene = osgDB::readNodeFile(file); - PointCollecter pc(po, no, co); - scene->accept(pc); - return scene.valid(); -}*/ - -#endif // READ_FILE \ No newline at end of file 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 53b9b606b1d..53a49352912 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 @@ -18,136 +18,120 @@ #include -#include +#include +#include +#include #include +#include + +#include //#include #include -#include "readFile.h" -#include "writeFile.h" -#include "quad.h" +//typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel; +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; -#include "Scale_space_surface_constructer.h" +typedef CGAL::Scale_space_surface_reconstructer_3< Kernel > Reconstructer; +typedef Reconstructer::Point Point; +typedef std::vector< Point > Pointset; + +typedef Reconstructer::Tripleset Tripleset; + +const unsigned int LINE_SIZE = 1024; + + +void readLine( std::ifstream& fin, char* line ) { + fin.getline( line, LINE_SIZE ); + while( fin && line[0] == '#' ) { + // Comment line. + std::cout << "Comment: " << line << std::endl; + fin.getline( line, LINE_SIZE ); + } +} + +bool readOFF( std::string input, Pointset& points ) { + // Read the input file. + std::ifstream fin( input.c_str() ); + + // Any .off file should start with the line "OFF" + char line[LINE_SIZE]; + readLine( fin, line ); + if( std::strncmp(line, "OFF", 3) != 0 ) return false; + + // The next line contains the object counts. + unsigned int n_points, n_polygons, n_segments; + fin >> n_points >> n_polygons >> n_segments; + + // Collect the points. + points.reserve( n_points ); + for( unsigned int n = 0; fin && n < n_points; ++n ) { + Point p; + fin >> p; + points.push_back( p ); + } + if( !fin ) return false; + + // Ignore the rest + fin.close(); + return true; +} + +bool writeOFF( std::string output, const Pointset& points, const Tripleset& triples ) { + // Write the output file. + std::ofstream fout( output ); + fout << "OFF" << std::endl; + fout << points.size() << " " << triples.size() << " 0" << std::endl; + + // Write the points. + for( Pointset::const_iterator pit = points.begin(); pit != points.end(); ++pit) + fout << *pit << std::endl; + + // Write the triples. + for( Tripleset::const_iterator tit = triples.begin(); tit != triples.end(); ++tit ) + fout << "3 " << *tit << std::endl; + + return true; +} int main(int argc, char** argv) { - if (argc > 1 && (!std::strcmp(argv[1], "-help") || !std::strcmp(argv[1], "--help"))) { - std::cout << "usage:\nDigne.exe [input].xyzn [outfile] [#neighbors] [#iterations] [#samples for r-estimate]" << std::endl; - exit(0); - } - - char* input = "kitten.xyzn"; - char* outFile = "kitten_o"; - unsigned int nn = 30; - unsigned int iter = 4; - unsigned int samples = 400; + std::string base = "kitten"; + unsigned int neighbors = 30; + unsigned int iterations = 4; + unsigned int samples = 200; if (argc > 1) - input = argv[1]; + base = argv[1]; if (argc > 2) - outFile = argv[2]; + sscanf(argv[2], "%u", &neighbors); if (argc > 3) - sscanf(argv[3], "%u", &nn); + sscanf(argv[3], "%u", &iterations); if (argc > 4) - sscanf(argv[4], "%u", &iter); - if (argc > 5) - sscanf(argv[5], "%u", &samples); - - //typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel; - typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; - - typedef Kernel::Point_3 Point_3; - typedef Kernel::Vector_3 Vector_3; - typedef Kernel::Triangle_3 Triangle_3; - - typedef std::vector PointCollection; - typedef std::vector VectorCollection; - typedef float4 Color; - typedef std::vector ColorCollection; - - // Read the input file. - PointCollection points; - VectorCollection normals; - ColorCollection colors; - std::string is(input); - std::transform(is.begin(), is.end(), is.begin(), std::tolower); - bool loaded = false; - if (is.rfind(".xyzn") == is.length()-5) - loaded = loadXYZN(input, std::back_inserter(points), std::back_inserter(normals)); -/* else if (is.rfind(".ive") == is.length()-4 || is.rfind(".osg") == is.length()-4) - loaded = loadOSG(input, std::back_inserter(points), std::back_inserter(normals), std::back_inserter(colors));*/ - if (!loaded) { - std::cout << "Error loading input." << std::endl; - exit(1); - } - std::cout << "Input loaded:" << std::endl; - std::cout << points.size() << " points, " << normals.size() << " normals, " << colors.size() << " colors" << std::endl; - if (points.size() != normals.size()) - normals.resize(points.size()); + sscanf(argv[4], "%u", &samples); + + std::string input = base + ".off"; + std::string output = base + "_ss.off"; + + // Read the data. + std::cout << "Input: " << input << std::flush; + Pointset points; + if( !readOFF( input, points ) ) { + std::cerr << std::endl << "Error reading " << input << std::endl; + exit(-1); + } + std::cout << " loaded: " << points.size() << " points." << std::endl; // Construct the mesh in a scale space. - typedef Scale_space_surface_constructer Scale_space_surface_constructer; - typedef Scale_space_surface_constructer::Triangle Triangle; - std::list triangles; - Scale_space_surface_constructer sssc(iter); - sssc(points.begin(), points.end(), std::back_inserter(triangles)); + Reconstructer reconstruct; - std::cout << "Save to: " << outFile << std::endl; - std::string fs(outFile); + reconstruct.compute_surface( points.begin(), points.end(), neighbors, iterations, samples ); + std::cout << "Reconstruction done." << std::endl; - // Save the points. - std::cout << "Save as xyzn" << std::endl; - std::ofstream fout; - fout.open((fs+".xyzn").c_str()); - if (fout) { - VectorCollection::const_iterator nit = normals.begin(); - for (Scale_space_surface_constructer::PointCollection::const_iterator pit = sssc.moved().begin(); pit != sssc.moved().end(); ++pit, ++nit) { - fout << *pit << " " << *nit << std::endl; - } - fout.close(); - } - else - std::cout << "Error writing to .xyzn output." << std::endl; - - - // Construct the normals of the triangles. - Kernel::Construct_cross_product_vector_3 cross; - VectorCollection tri_norm; - tri_norm.reserve(triangles.size()); - for (std::list::const_iterator tit = triangles.begin(); tit != triangles.end(); ++tit) { - const Point_3& p1 = tit->vertex(0); - const Point_3& p2 = tit->vertex(1); - const Point_3& p3 = tit->vertex(1); - - Vector_3 n = cross((p2 - p1), (p3 - p1)); - tri_norm.push_back(n / CGAL::sqrt(CGAL::to_double(n.squared_length()))); - } - - // Link the colors of the points. - std::map toColor; - PointCollection::const_iterator pit = points.begin(); - for (ColorCollection::const_iterator cit = colors.begin(); cit != colors.end(); ++cit) - toColor[*pit++] = *cit; - ColorCollection tri_colors; - tri_colors.reserve(3*triangles.size()); - if (!toColor.empty()) - for (std::list::const_iterator tit = triangles.begin(); tit != triangles.end(); ++tit) - for (unsigned int i = 0; i < 3; ++i) - tri_colors.push_back(toColor[tit->vertex(i)]); - - // Save the shapes (OSG style). -/* bool saveIve = false; - std::cout << "Save as IVE" << std::endl; - osg::ref_ptr shape = new osg::Switch; - if (savePointsIVE(*shape, sssc.moved().begin(), sssc.moved().end(), normals.begin(), normals.end(), colors.begin(), colors.end())) - if (saveTrianglesIVE(*shape, facets.begin(), facets.end(), tri_norm.begin(), tri_norm.end(), tri_colors.begin(), tri_colors.end())) - if (saveIVE(shape.get(), (fs+".ive").c_str())) - saveIve = true; - if (!saveIve) - std::cout << "Error writing to .ive output." << std::endl;*/ - - // Save the shapes (OFF style). - std::cout << "Save as OFF" << std::endl; - if (!saveOFF((fs+".off").c_str(), triangles.begin(), triangles.end(), tri_norm.begin(), tri_norm.end(), colors.begin(), colors.end())) - std::cout << "Error writing to .off output." << 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; + exit(-1); + } + std::cout << " written." << std::endl; } \ No newline at end of file diff --git a/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/writeFile.h b/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/writeFile.h deleted file mode 100644 index e56eb86733c..00000000000 --- a/Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/writeFile.h +++ /dev/null @@ -1,185 +0,0 @@ -//A general purpose file writer. -//Note that if we define more file formats, -//it may be beneficial to split this into multiple plugin classes/files. -//Copyright (C) 2013 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 WRITE_FILE -#define WRITE_FILE - -#include - -/*#include -#include -#include - -#include -#include - -template < class Point > -osg::Vec3 toOSG(const Point& p) { - return osg::Vec3(CGAL::to_double(p.x()), CGAL::to_double(p.y()), CGAL::to_double(p.z())); -} - -template < class PointIterator, class NormalIterator, class ColorIterator > -bool savePointsIVE(osg::Switch& shape, - PointIterator points_start, PointIterator points_end, - NormalIterator normals_start, NormalIterator normals_end, - ColorIterator colors_start, ColorIterator colors_end) { - osg::ref_ptr geode = new osg::Geode; - osg::ref_ptr geom = new osg::Geometry; - - osg::ref_ptr vert = new osg::Vec3Array; - for (PointIterator pit = points_start; pit != points_end; ++pit) - vert->push_back(toOSG(*pit)); - geom->setVertexArray(vert.get()); - - osg::ref_ptr norm = new osg::Vec3Array; - norm->reserve(vert->getNumElements()); - for (NormalIterator nit = normals_start; nit != normals_end; ++nit) - norm->push_back(toOSG(*nit)); - if (norm->getNumElements() != vert->getNumElements()) - return false; - geom->setNormalBinding(osg::Geometry::BIND_PER_VERTEX); - geom->setNormalArray(norm.get()); - - osg::ref_ptr colors = new osg::Vec4Array; - if (colors_start == colors_end) { - colors->reserve(1); - colors->push_back(osg::Vec4(0.5,0.5,0.5,1.0)); - geom->setColorBinding(osg::Geometry::BIND_OVERALL); - } - else { - colors->reserve(vert->getNumElements()); - for (ColorIterator cit = colors_start; cit != colors_end; ++cit) - colors->push_back(osg::Vec4(cit->r(), cit->g(), cit->b(), cit->a())); - geom->setColorBinding(osg::Geometry::BIND_PER_VERTEX); - } - geom->setColorArray(colors.get()); - - geom->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::POINTS, 0, vert->size())); - - geode->addDrawable(geom.get()); - shape.addChild(geode.get()); - - return true; -} - -template < class FacetIterator, class NormalIterator, class ColorIterator > -bool saveTrianglesIVE(osg::Switch& shape, - FacetIterator facets_start, FacetIterator facets_end, - NormalIterator normals_start, NormalIterator normals_end, - ColorIterator colors_start, ColorIterator colors_end) { - size_t size = std::distance(facets_start, facets_end); - if (size != std::distance(normals_start, normals_end)) - return false; - - osg::ref_ptr geode = new osg::Geode; - osg::ref_ptr geom = new osg::Geometry; - - osg::ref_ptr vert = new osg::Vec3Array; - vert->reserve(3 * size); - osg::ref_ptr norm = new osg::Vec3Array; - norm->reserve(size); - NormalIterator nit = normals_start; - for (FacetIterator tit = facets_start; tit != facets_end; ++tit, ++nit) { - vert->push_back(toOSG(tit->first->vertex((tit->second+1)&3)->point())); - vert->push_back(toOSG(tit->first->vertex((tit->second+2)&3)->point())); - vert->push_back(toOSG(tit->first->vertex((tit->second+3)&3)->point())); - norm->push_back(toOSG(*nit)); - } - - osg::ref_ptr colors = new osg::Vec4Array; - if (colors_start == colors_end) { - colors->reserve(1); - colors->push_back(osg::Vec4(0.5,0.5,0.5,0.5)); - geom->setColorBinding(osg::Geometry::BIND_OVERALL); - } - else { - colors->reserve(vert->getNumElements()); - for (ColorIterator cit = colors_start; cit != colors_end; ++cit) - colors->push_back(osg::Vec4(cit->r(), cit->g(), cit->b(), cit->a())); - geom->setColorBinding(osg::Geometry::BIND_PER_VERTEX); - } - - geom->setVertexArray(vert.get()); - geom->setNormalArray(norm.get()); - geom->setNormalBinding(osg::Geometry::BIND_PER_PRIMITIVE); - geom->setColorArray(colors.get()); - geom->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::TRIANGLES, 0, vert->size())); - geode->addDrawable(geom.get()); - - osg::StateSet* state = geode->getOrCreateStateSet(); - state->setAttributeAndModes(new osg::CullFace(osg::CullFace::BACK)); -// osg::Material* material = new osg::Material; -// material->setDiffuse(osg::Material::FRONT_AND_BACK, osg::Vec4(0.2, 0.2, 0.2, 0.25)); -// state->setAttributeAndModes(material); - shape.addChild(geode.get()); - - return true; -} - -bool saveIVE(osg::Node* node, const char* file) { - osg::ref_ptr options = new osgDB::ReaderWriter::Options("precision 10"); - return osgDB::writeNodeFile(*node, file, options.get()); -}*/ - -template < class TriangleIterator, class NormalIterator, class ColorIterator > -bool saveOFF(const char* file, - TriangleIterator triangles_start, TriangleIterator triangles_end, - NormalIterator normals_start, NormalIterator normals_end, - ColorIterator colors_start, ColorIterator colors_end) { - typedef CGAL::Point_3 Point; - - // Construct a mapping from the points to their indices. - std::map map; - unsigned int ind = 0, tri = 0; - for (TriangleIterator tit = triangles_start; tit != triangles_end; ++tit) { - for (unsigned int i = 0; i < 4; ++i) - if (map.find(tit->vertex(i)) == map.end()) - map[tit->vertex(i)] = ind++; - ++tri; - } - - if (tri != std::distance(normals_start, normals_end)) - return false; - - std::ofstream fout(file); - fout << "OFF" << std::endl; - fout << map.size() << " " << tri << " 0" << std::endl; - - // Write the points. - ColorIterator cit = colors_start; - for (std::map::const_iterator pit = map.begin(); pit != map.end(); ++pit) { - fout << pit->first; - if (cit != colors_end) - fout << " " << *cit++; - fout << std::endl; - } - - // Write the facets. - for (TriangleIterator tit = triangles_start; tit != triangles_end; ++tit) { - fout << "3 " << map[tit->vertex(0)] - << " " << map[tit->vertex(1)] - << " " << map[tit->vertex(2)] << std::endl; - } - - return true; -} - -#endif // WRITE_FILE \ No newline at end of file diff --git a/Scale_space_reconstruction_3/include/CGAL/Mean_neighborhood_ball.h b/Scale_space_reconstruction_3/include/CGAL/Mean_neighborhood_ball.h deleted file mode 100644 index 1b1157f6ad0..00000000000 --- a/Scale_space_reconstruction_3/include/CGAL/Mean_neighborhood_ball.h +++ /dev/null @@ -1,104 +0,0 @@ -//A method that computes the size of a mean neighborhood. -//Copyright (C) 2013 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 MEAN_NEIGHBORHOOD_BALL -#define MEAN_NEIGHBORHOOD_BALL - -#include -#include -#include - -// Compute the mean radius of the spheres containing the specified number of nearest neighbors from a point sample. -template < typename Kernel > -class Mean_neighborhood_ball { - typedef typename Kernel::FT Scalar; - - typedef typename CGAL::Search_traits_3 Tree_traits; - - typedef typename CGAL::Orthogonal_k_neighbor_search Neighbor_search; - typedef typename Neighbor_search::iterator Neighbor_iterator; - typedef typename Neighbor_search::Tree Neighbor_tree; - - typedef CGAL::Random Random; - - unsigned int handled, checked; - Scalar dist_far; - - Neighbor_tree tree; - - Random generator; - -public: - typedef typename Kernel::Point_3 Point; - - typedef Scalar result_type; - -private: - unsigned int neighbors, samples; - -public: - Mean_neighborhood_ball(unsigned int nn = 30, unsigned int smp = 200): neighbors(nn), samples(smp) {} - - unsigned int get_neighbors() const { return neighbors; } - void set_neighbors( unsigned int nn ) { neighbors = nn; } - unsigned int get_samples() const { return samples; } - void set_samples( unsigned int smp ) { samples = smp; } - - template < class InputIterator > - void constructTree(InputIterator start, InputIterator end) { - std::cout << "MNB: tree" << std::endl; - tree.clear(); - tree.insert(start, end); - if(!tree.is_built()) - tree.build(); - } - - void handlePoint(const Point& p) { - Neighbor_search search(tree, p, neighbors+1); - - dist_far += CGAL::sqrt(CGAL::to_double(CGAL::squared_distance(p, (search.end()-1)->first))); - } - - result_type compute() { - handled = 0; - checked = 0; - dist_far = 0; - - std::cout << "MNB: sample" << std::endl; - for (typename Neighbor_tree::const_iterator it = tree.begin(); it != tree.end(); ++it) { - if (generator.get_double() < double(samples - checked) / (tree.size() - handled)) { - handlePoint(*it); - ++checked; - } - ++handled; - } - - dist_far /= double(checked); - - return dist_far; - } - - template < class InputIterator > - result_type operator()(InputIterator start, InputIterator end) { - constructTree(start, end); - return compute(); - } -}; // Mean_neighborhood_ball - -#endif MEAN_NEIGHBORHOOD_BALL diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_constructer.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_constructer.h index b07e8f69499..1387902ac61 100644 --- a/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_constructer.h +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_constructer.h @@ -20,103 +20,687 @@ #ifndef SCALE_SPACE_SURFACE_CONSTRUCTER #define SCALE_SPACE_SURFACE_CONSTRUCTER +#include + #include #include #include +#include -#include "Mean_neighborhood_ball.h" -#include "Scale_space_transform.h" -#include "Surface_mesher.h" +#include -// Smooth a set of points and then construct a surface mesh on the smoothed set; -// Finally transport the surface back to the original point locations. -template < typename Kernel, - typename Vb = CGAL::Triangulation_vertex_base_3, - typename Cb = CGAL::Triangulation_cell_base_3 > -class Scale_space_surface_constructer { - typedef typename Mean_neighborhood_ball ComputeRadius; - typedef typename Scale_space_transform Scale_space_transform; - typedef typename Surface_mesher Surface_mesher; +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include + +#include + + +namespace CGAL { + +//a functor that returns a std::pair. +//the unsigned integer is incremented at each call to +//operator() +template < class T > +class Auto_count: public std::unary_function< const T&, std::pair< T, unsigned int > > { + mutable unsigned int i; +public: + Auto_count(): i(0) {} + std::pair< T,unsigned int > operator()( const T& t ) const { return std::make_pair( t, i++ ); } +}; + +// Struct to choose between the alpha-shape and fixed alpha-shape +template < class Kernel, class Fixed_shape > +class _Shape { + typedef typename Kernel::FT Scalar; + + 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; public: - typedef typename Kernel::FT Scalar; + typedef CGAL::Delaunay_triangulation_3< Kernel, Tds > Structure; + typedef CGAL::Alpha_shape_3< Structure > Shape; + + _Shape() {} + Shape* construct( Shape* shape, const Scalar& squared_radius ) { + if( shape ) return new Shape( *shape, squared_radius, Shape::GENERAL ); + else return new Shape( squared_radius, Shape::GENERAL ); + } + 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() ), + squared_radius, Shape::GENERAL ); + } +}; // class _Shape - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::Triangle_3 Triangle; - - typedef typename Surface_mesher::Vertex_handle Vertex_handle; - typedef typename Surface_mesher::Facet Facet; +// Struct to choose between the alpha-shape and fixed alpha-shape +// specialization: fixed.. +template < class Kernel > +class _Shape < Kernel, CGAL::Tag_true > { + typedef typename Kernel::FT Scalar; - typedef typename std::vector PointCollection; + 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 CGAL::Triangulation_data_structure_3 Tds; + +public: + typedef CGAL::Delaunay_triangulation_3< Kernel, Tds > Structure; + typedef CGAL::Fixed_alpha_shape_3< Structure > Shape; + + _Shape() {} + Shape* construct( Shape* shape, const Scalar& squared_radius ) { + if( shape ) return new Shape( *shape, squared_radius ); + else return new Shape( squared_radius ); + } + 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() ), + squared_radius ); + } +}; // class _Shape + + +/// Compute a smoothed surface mesh from a collection of points. +/** An appropriate neighborhood size is estimated, followed by a + * number of smoothing iterations. Finally, the surface of the + * smoothed point set is computed. + * + * The order of the point set remains the same, meaning that + * the smoothed surface can be transposed back to its unsmoothed + * version by overwriting the smoothed point collection with its + * unsmoothed version. + * \tparam Kernel the geometric traits class. This class + * specifies, amongst other things, the number types and + * predicates used. + * \tparam Fixed_shape indicates whether shape of the object + * should be constructed for a fixed neighborhood size. + * + * Generally, constructing for a fixed neighborhood size is more + * efficient. This is not the case if the surface should be + * constructed for different neighborhood sizes without changing + * the point set or recomputing the scale-space. + * \tparam Shells indicates whether to collect the surface per shell. + * + * A shell is a connected component of the surface where connected + * facets are locally oriented towards the same side of the surface. + * \sa Mean_neighborhood_ball. + * \sa Scale_space_transform. + * \sa Surface_mesher. + */ +template < class Kernel, class Fixed_shape = CGAL::Tag_true, class Shells = CGAL::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 > + Static_search; + typedef typename CGAL::Orthogonal_incremental_neighbor_search< Search_traits > + Dynamic_search; + typedef typename Dynamic_search::Tree Search_tree; + + typedef CGAL::Random Random; + + // Constructing the surface. + typedef _Shape< Kernel, Fixed_shape > Shape_generator; + + typedef typename Shape_generator::Shape Shape; + typedef typename Shape::Vertex_handle Vertex_handle; + typedef typename Shape::Cell_handle Cell_handle; + typedef typename Shape::Facet Facet; + + typedef typename Shape::All_cells_iterator All_cells_iterator; + typedef typename Shape::Finite_facets_iterator Finite_facets_iterator; + typedef typename Shape::Classification_type Classification_type; + +public: + typedef Shells Collect_per_shell; ///< Whether to collect the surface per shell. + + typedef typename Kernel::FT Scalar; ///< The number type. + + typedef typename Kernel::Point_3 Point; ///< The point type. + typedef typename Kernel::Triangle_3 Triangle; ///< The triangle type. + + typedef Triple< unsigned int, unsigned int, unsigned int > + Triple; ///< A triangle of the surface. + typedef std::list< Triple > Tripleset; ///< A collection of triples. + + typedef typename Search_tree::iterator Point_iterator; ///< An iterator over the points. + typedef typename Search_tree::const_iterator Point_const_iterator; ///< A constant iterator over the points. private: - PointCollection _moved; + typedef typename std::vector Pointset; ///< A collection of points. - ComputeRadius mnb; - unsigned int iterations; - Surface_mesher bp; +private: + Search_tree _tree; // To quickly search for nearest neighbors. + + Random _generator; // For sampling random points. + + unsigned int _mean_neighbors; // The number of nearest neighbors in the mean neighborhood. + unsigned int _samples; // The number of sample points for estimating the mean neighborhood. + + Scalar _squared_radius; // The squared mean neighborhood radius. + + // The shape must be a pointer, because the alpha of + // a Fixed_alpha_shape_3 can only be set at + // construction and its assignment operator is private. + Shape* _shape; + + // The surface. If the surface is collected per shell, + // consecutive triples belong to the same shell and + // different shells are separated by a (0,0,0) triple. + Tripleset _surface; + +private: + void clear_tree() { _tree.clear(); } public: - Scale_space_surface_constructer(unsigned int iter = 1): mnb(), iterations(iter), bp(0) {} + /// Default constructor. + /** \param sq_radius the squared radius of the + * neighborhood size. If this value is negative when + * the point set is smoothed or when the surface is computed, + * the neighborhood size will be computed automatically. + */ + Scale_space_surface_reconstructer_3(unsigned int neighbors = 30, unsigned int samples = 200, Scalar sq_radius = -1 ): _mean_neighbors(neighbors), _samples(samples), _squared_radius( sq_radius ), _shape(0) {} + ~Scale_space_surface_reconstructer_3() { deinit_shape(); } - void set_iterations(unsigned int iter) {iterations = iter;} +private: + void deinit_shape() { if( _shape != 0 ) { delete _shape; _shape = 0; } } - const PointCollection& moved() const {return _moved;} +public: + Point_iterator scale_space_begin() { return _tree.begin(); } + Point_iterator scale_space_end() { return _tree.end(); } + + Point_const_iterator scale_space_begin() const { return _tree.begin(); } + Point_const_iterator scale_space_end() const { return _tree.begin(); } - // Input: iterators over Point, output: collection of Triangles. - template < class InputIterator, class OutputIterator > - OutputIterator operator()(InputIterator start, InputIterator end, OutputIterator out) { - typedef std::map Map; - typedef std::list List; - // Compute the radius for which the mean ball would contain 30 points. - Scalar radius2 = mnb(start, end); - radius2 *= radius2; + bool has_surface() const { return _shape != 0; } - // Perturb the points. - Scale_space_transform ssp(radius2); - ssp(start, end, iterations); - // Note that the transform has constructed a novel vector to hold the tranformed points. - _moved = PointCollection(ssp.points().begin(), ssp.points().end()); + void clear_surface() { + if( has_surface() ) { + _shape->clear(); + } + } + + /// Insert a collection of points. + /** \tparam InputIterator an iterator over a collection of points. + * The iterator must point to a Point type. + * \param start an iterator to the first point of the collection. + * \param end a past-the-end iterator for the point collection. + * \sa compute_surface(InputIterator start, InputIterator end). + */ + template < class InputIterator > + void insert_points( InputIterator start, InputIterator end ) { + _tree.insert( start, end ); + } - // Perform ball-pivoting on the perturbed points. - std::list facets; - bp = Surface_mesher(radius2); - bp(_moved.begin(), _moved.end(), std::back_inserter(facets)); + void clear() { + clear_tree(); + clear_surface(); + } + + /// Estimate the mean neighborhood size based on a number of sample points. + /** A neighborhood size is expressed as the radius of the smallest + * ball centered on a point such that the ball contains at least + * a specified number of points. + * + * The mean neighborhood size is then the mean of these radii, + * taken over a number of point samples. + * \return the estimated mean neighborhood size. + * \sa handlePoint(const Point& p). + * \sa operator()(InputIterator start, InputIterator end). + */ + Scalar estimate_mean_neighborhood( unsigned int neighbors = 30, unsigned int samples = 200 ) { + Kernel::Compute_squared_distance_3 squared_distance = Kernel().compute_squared_distance_3_object(); - // Transport the points back to their original position. - // Note that if you don't need the vertex base to have info, - // you can replace the vertices by CGAL::Triangulation_vertex_base_with_info_3 - // to reduce the computation time for retrieving the original point locations - // (just set the info of each vertex to its index in the original collection). - Map map; - InputIterator oit = start; - for (PointCollection::iterator pit = _moved.begin(); pit != _moved.end(); ++pit, ++oit) - map[*pit] = *oit; + _mean_neighbors = neighbors; + _samples = samples; - Triangle t; - for (List::const_iterator fit = facets.begin(); fit != facets.end(); ++fit) { - t = bp.oriented_triangle(*fit); - *out++ = Triangle(map[t.vertex(0)], - map[t.vertex(1)], - map[t.vertex(2)]); + unsigned int handled = 0; + unsigned int checked = 0; + Scalar radius = 0; + + if( !_tree.is_built() ) + _tree.build(); + + for (typename Search_tree::const_iterator it = _tree.begin(); it != _tree.end(); ++it) { + 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 ) ) ); + ++checked; + } + ++handled; + } + radius /= double(checked); + + set_mean_neighborhood( radius ); + return radius; + } + + + template < class InputIterator > + Scalar estimate_mean_neighborhood(InputIterator start, InputIterator end, unsigned int neighbors = 30, unsigned int samples = 200) { + insert_points(start, end); + return estimate_mean_neighborhood(neighbors, samples); + } + + // -1 if not yet set. + Scalar get_squared_mean_neighborhood() const { return _squared_radius; } + void set_mean_neighborhood( const Scalar& radius ) { + _squared_radius = radius * radius; + if( has_surface() ) + _shape = Shape_generator().construct( _shape, _squared_radius ); + } + bool has_squared_mean_neighborhood() const { + return sign( _squared_radius ) == POSITIVE; + } + + unsigned int get_mean_neighbors() const { return _mean_neighbors; } + unsigned int get_number_samples() const { return _samples; } + void set_mean_neighbors(unsigned int neighbors) { _mean_neighbors = neighbors;} + void set_number_samples(unsigned int samples) { _samples = samples; } + + + /// Compute a number of iterations of scale-space transforming. + /** If earlier iterations have been computed, calling smooth_scale_space() + * will add more iterations. + * + * If the mean neighborhood is negative, it will be computed first. + * \param iterations the number of iterations to perform. + */ + void smooth_scale_space(unsigned int iterations = 1) { + + typedef std::vector 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 _ENV::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 ); + + Pointset smoothed_points; + smoothed_points.assign( _tree.begin(), _tree.end() ); + _tree.clear(); + + // 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( smoothed_points.begin(), smoothed_points.end() ); + if(!tree.is_built()) + tree.build(); + + // Collect the number of neighbors of each point. + // This can be done parallel. + CountVec neighbors(_tree.size(), 0); + Kernel::Compare_squared_distance_3 compare; + p_size_t count = _tree.size(); // openMP can only use signed variables. + const Scalar squared_radius = _squared_radius; // openMP can only use local variables. +#pragma omp parallel for shared(count,tree,smoothed_points,squared_radius,neighbors) firstprivate(compare) + 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, tree[i]); + for (Dynamic_search::iterator nit = search.begin(); nit != search.end() && compare(tree[i], nit->first, squared_radius) != CGAL::LARGER; ++nit) + ++neighbors[i]; + } + + // Construct a mapping from each point to its index. + PIMap indices; + p_size_t index = 0; + for( Search_tree::const_iterator tit = tree.begin(); tit != tree.end(); ++tit, ++index) + indices[ *tit ] = index; + + // Compute the tranformed point locations. + // This can be done parallel. +#pragma omp parallel for shared(count,neighbors,tree,squared_radius) firstprivate(compare) + for (p_size_t i = 0; i < count; ++i) { + // If the neighborhood is too small, the vertex is not moved. + if (neighbors[i] < 4) + continue; + + // Collect the vertices within the ball and their weights. + Dynamic_search search(tree, tree[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(tree[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]); + wts(column) = 1.0 / neighbors[indices[nit->first]]; + } + + // Construct the barycenter. + Vector3 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(); + for (column = 0; column < pts.cols(); ++column) + covariance += wts.matrix()(column) * pts.col(column) * pts.col(column).transpose(); + + // Construct the Eigen system. + EigenSolver solver(covariance); + + // If the Eigen system does not converge, the vertex is not moved. + if (solver.info() != Eigen::Success) + continue; + + // 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! + CGAL_assertion(false); + continue; + } + Vector3 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(tree[i][0]), CGAL::to_double(tree[i][1]), CGAL::to_double(tree[i][2])) - bary; + Vector3 per = bary + bv - (n.dot(bv) * n); + + smoothed_points[i] = Point(per(0), per(1), per(2)); + } + } + _tree.insert( smoothed_points.begin(), smoothed_points.end() ); + } + + + /// Compute a number of transform iterations on a collection of points. + /** This method is equivalent to running [insert_points(start, end)](\ref insert_points) + * followed by [smooth_scale_space(iterations)](\ref smooth_scale_space). + * \tparam InputIterator an iterator over a collection of points. + * The iterator must point to a Point type. + * \param start an iterator to the first point of the collection. + * \param end a past-the-end iterator for the point collection. + * \param iterations the number of iterations to perform. + * \sa operator()(InputIterator start, InputIterator end). + * \sa assign(InputIterator start, InputIterator end). + * \sa iterate(unsigned int iterations). + */ + template < class InputIterator > + void smooth_scale_space(InputIterator start, InputIterator end, unsigned int iterations = 1) { + insert_points(start, end); + smooth_scale_space(iterations); + } + +private: + // Once a facet is added to the surface, it is marked as handled. + inline bool is_handled( Cell_handle c, unsigned int li ) const { + switch( li ) { + case 0: return ( c->info()&1 ) != 0; + case 1: return ( c->info()&2 ) != 0; + case 2: return ( c->info()&4 ) != 0; + case 3: return ( c->info()&8 ) != 0; + } + return false; + } + inline bool is_handled( const Facet& f ) const { return is_handled( f.first, f.second ); } + inline void set_handled( Cell_handle c, unsigned int li ) { + switch( li ) { + case 0: c->info() |= 1; return; + case 1: c->info() |= 2; return; + case 2: c->info() |= 4; return; + case 3: c->info() |= 8; return; + } + } + inline void set_handled( Facet f ) { set_handled( f.first, f.second ); } + +public: + + // make new construct_shape() method for when pts already known... + void construct_shape() { + construct_shape( scale_space_begin(), scale_space_end() ); + } + + /*// For if you already have a Delaunay triangulation of the points. + void construct_shape(Shape_generator::Structure& tr ) { + deinit_shape(); + if( !has_squared_mean_neighborhood() ) + estimate_mean_neighborhood( _mean_neighbors, _samples ); + _shape = Shape_generator().construct( *tr, r2 ); + }*/ + + // If you don't want to smooth the point set. + template < class InputIterator > + void construct_shape(InputIterator start, InputIterator end) { + deinit_shape(); + if( !has_squared_mean_neighborhood() ) + estimate_mean_neighborhood( _mean_neighbors, _samples ); + _shape = Shape_generator().construct( start, end, _squared_radius ); + } + +private: + Triple ordered_facet_indices( const Facet& f ) const { + if( (f.second&1) == 0 ) + return Triple( f.first->vertex( (f.second+2)&3 )->info(), + f.first->vertex( (f.second+1)&3 )->info(), + f.first->vertex( (f.second+3)&3 )->info() ); + else + return Triple( f.first->vertex( (f.second+1)&3 )->info(), + f.first->vertex( (f.second+2)&3 )->info(), + f.first->vertex( (f.second+3)&3 )->info() ); + } + + void collect_shell( Cell_handle c, unsigned int li ) { + // Collect one surface mesh from the alpha-shape in a fashion similar to ball-pivoting. + // Invariant: the facet is regular or singular. + + // To stop stack overflows: use own stack. + std::stack stack; + stack.push( Facet(c, li) ); + + Facet f; + Cell_handle n, p; + int ni, pi; + Vertex_handle a; + Classification_type cl; + bool processed = false; + while( !stack.empty() ) { + f = stack.top(); + stack.pop(); + + // Check if the cell was already handled. + // Note that this is an extra check that in many cases is not necessary. + if( is_handled(f) ) + continue; + + // The facet is part of the surface. + CGAL_triangulation_assertion( !_shape->is_infinite(f) ); + set_handled(f); + + // Output the facet as a triple of indices. + _surface.push_back( ordered_facet_indices(f) ); + + if( Shells::value ) { + // Pivot over each of the facet's edges and continue the surface at the next regular or singular facet. + for( unsigned int i = 0; i < 4; ++i ) { + // Skip the current facet. + if( i == f.second || is_handled(f.first, i) ) + continue; + + // Rotate around the edge (starting from the shared facet in the current cell) until a regular or singular facet is found. + n = f.first; + ni = i; + a = f.first->vertex(f.second); + cl = _shape->classify( Facet(n, ni) ); + while( cl != Shape::REGULAR && cl != Shape::SINGULAR ) { + p = n; + n = n->neighbor(ni); + ni = n->index(a); + pi = n->index(p); + a = n->vertex(pi); + cl = _shape->classify( Facet(n, ni) ); + } + + // Continue the surface at the next regular or singular facet. + stack.push( Facet(n, ni) ); + } + processed = true; + } + } + + // We indicate the end of a shell by the (0,0,0) triple. + if( Shells::value && processed ) + _surface.push_back( Triple(0,0,0) ); + } + + void collect_shell( const Facet& f ) { + collect_shell( f.first, f.second ); + } + +public: + template < class OutputIterator > + OutputIterator collect_surface( OutputIterator out ) { + if( !has_squared_mean_neighborhood() ) + estimate_mean_neighborhood( _mean_neighbors, _samples ); + + // Collect all surface meshes from the alpha-shape in a fashion similar to ball-pivoting. + // Reset the facet handled markers. + for( All_cells_iterator cit = _shape->all_cells_begin(); cit != _shape->all_cells_end(); ++cit ) + cit->info() = 0; + + // We check each of the facets: if it is not handled and either regular or singular, + // we start collecting the next surface from here. + Facet m; + int ns = 0; + for( Finite_facets_iterator fit = _shape->finite_facets_begin(); fit != _shape->finite_facets_end(); ++fit ) { + m = _shape->mirror_facet( *fit ); + switch( _shape->classify( *fit ) ) { + case Shape::REGULAR: + if( !is_handled(*fit) && !is_handled(m) ) + ++ns; + // Build a surface from the outer cell. + if( _shape->classify(fit->first) == Shape::EXTERIOR ) + collect_shell( *fit ); + else + collect_shell( m ); + break; + case Shape::SINGULAR: + if( !is_handled( *fit ) ) + ++ns; + // Build a surface from both incident cells. + collect_shell( *fit ); + if( !is_handled(m) ) + ++ns; + collect_shell( m ); + break; + } } return out; } - void ordered_vertices(const Facet& f, Vertex_handle& v0, Vertex_handle& v1, Vertex_handle& v2) { - if ((f.second&1) == 0) { - v0 = f.first->vertex((f.second+2)&3); - v1 = f.first->vertex((f.second+1)&3); - v2 = f.first->vertex((f.second+3)&3); - } - else { - v0 = f.first->vertex((f.second+1)&3); - v1 = f.first->vertex((f.second+2)&3); - v2 = f.first->vertex((f.second+3)&3); - } + template < class InputIterator, class OutputIterator > + OutputIterator collect_surface( InputIterator start, InputIterator end, OutputIterator out ) { + construct_shape( start, end ); + return collect_surface( out ); } -}; // class Scale_space_surface_constructer + + template < class InputIterator > + void collect_surface( InputIterator start, InputIterator end ) { + construct_shape( start, end ); + collect_surface( std::back_inserter(_surface) ); + } + + void collect_surface() { + collect_surface( scale_space_begin(), scale_space_end() ); + } + + const Tripleset& surface() const { return _surface; } + Tripleset& surface() { return _surface; } + + + + /// Compute a smoothed surface mesh from a collection of points. + /** This is equivalent to calling insert_points(), + * iterate(), and compute_surface(). + * + * After this computation, the surface can be accessed using + * surface(). + * \tparam InputIterator an iterator over the point collection. + * The iterator must point to a Point type. + * \param start an iterator to the first point of the collection. + * \param end a past-the-end iterator for the point collection. + * + * \sa surface(). + * \sa surface() const. + */ + template < class InputIterator > + void compute_surface(InputIterator start, InputIterator end, unsigned int neighbors = 30, unsigned int iterations = 4, unsigned int samples = 200) { + typedef std::map Map; + + + // Compute the radius for which the mean ball would contain the required number of neighbors. + clear(); + insert_points( start, end ); + + _mean_neighbors = neighbors; + _samples = samples; + + + if( !has_squared_mean_neighborhood() ) + estimate_mean_neighborhood( _mean_neighbors, samples ); // TMP: move this into transform() and compute_surface() + + + // Smooth the scale space. + smooth_scale_space(iterations); + + // Mesh the perturbed points. + collect_surface(); + } +}; // class Scale_space_surface_reconstructer_3 + + +template< typename T1, typename T2, typename T3 > +std::ostream& +operator<<( std::ostream& os, const Triple< T1, T2, T3 >& t ) { + return os << t.first << " " << t.second << " " << t.third; +} + +template< typename T1, typename T2, typename T3 > +std::istream& +operator>>( std::istream& is, Triple< T1, T2, T3 >& t ) { + return is >> t.first >> t.second >> t.third; +} + +} // namespace CGAL #endif // SCALE_SPACE_SURFACE_CONSTRUCTER \ No newline at end of file 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 new file mode 100644 index 00000000000..39938a608f5 --- /dev/null +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstructer_3.h @@ -0,0 +1,794 @@ +//A method to construct a surface. +//Copyright (C) 2013 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 SCALE_SPACE_SURFACE_CONSTRUCTER +#define SCALE_SPACE_SURFACE_CONSTRUCTER + +#include + +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include + +#include + + +namespace CGAL { + +//a functor that returns a std::pair. +//the unsigned integer is incremented at each call to +//operator() +template < class T > +class Auto_count: public std::unary_function< const T&, std::pair< T, unsigned int > > { + mutable unsigned int i; +public: + Auto_count(): i(0) {} + std::pair< T,unsigned int > operator()( const T& t ) const { return std::make_pair( t, i++ ); } +}; + +// Struct to choose between the alpha-shape and fixed alpha-shape +template < class Kernel, class Fixed_shape > +class _Shape { + typedef typename Kernel::FT Scalar; + + 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; + +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. + + _Shape() {} + Shape* construct( Shape* shape, const Scalar& squared_radius ) { + if( shape ) return new Shape( *shape, squared_radius, Shape::GENERAL ); + else return new Shape( squared_radius, Shape::GENERAL ); + } + 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() ), + squared_radius, Shape::GENERAL ); + } +}; // class _Shape + +// Struct to choose between the alpha-shape and fixed alpha-shape +// specialization: fixed.. +template < class Kernel > +class _Shape < Kernel, CGAL::Tag_true > { + typedef typename Kernel::FT Scalar; + + 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 CGAL::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. + + _Shape() {} + Shape* construct( Shape* shape, const Scalar& squared_radius ) { + if( shape ) return new Shape( *shape, squared_radius ); + else return new Shape( squared_radius ); + } + 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() ), + squared_radius ); + } +}; // class _Shape + + +/// Compute a smoothed surface mesh from a collection of points. +/** An appropriate neighborhood size is estimated, followed by a + * number of smoothing iterations. Finally, the surface of the + * smoothed point set is computed. + * + * The order of the point set remains the same, meaning that + * the smoothed surface can be transposed back to its unsmoothed + * version by overwriting the smoothed point collection with its + * unsmoothed version. + * \tparam Kernel the geometric traits class. This class + * specifies, amongst other things, the number types and + * predicates used. + * \tparam Fixed_shape indicates whether shape of the object + * should be constructed for a fixed neighborhood size. + * + * Generally, constructing for a fixed neighborhood size is more + * efficient. This is not the case if the surface should be + * constructed for different neighborhood sizes without changing + * the point set or recomputing the scale-space. + * \tparam Shells indicates whether to collect the surface per shell. + * + * A shell is a connected component of the surface where connected + * facets are locally oriented towards the same side of the surface. + * \sa Mean_neighborhood_ball. + * \sa Scale_space_transform. + * \sa Surface_mesher. + */ +template < class Kernel, class Fixed_shape = CGAL::Tag_true, class Shells = CGAL::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 > + Static_search; + typedef typename CGAL::Orthogonal_incremental_neighbor_search< Search_traits > + Dynamic_search; + typedef typename Dynamic_search::Tree Search_tree; + + typedef CGAL::Random Random; + + // Constructing the surface. + typedef _Shape< Kernel, Fixed_shape > Shape_generator; + + typedef typename Shape_generator::Shape Shape; + typedef typename Shape::Vertex_handle Vertex_handle; + typedef typename Shape::Cell_handle Cell_handle; + typedef typename Shape::Facet Facet; + + typedef typename Shape::All_cells_iterator All_cells_iterator; + typedef typename Shape::Finite_facets_iterator Finite_facets_iterator; + typedef typename Shape::Classification_type Classification_type; + +public: + typedef Shells Collect_per_shell; ///< Whether to collect the surface per shell. + + typedef typename Kernel::FT Scalar; ///< The number type. + + typedef typename Kernel::Point_3 Point; ///< The point type. + typedef typename Kernel::Triangle_3 Triangle; ///< The triangle type. + + typedef Triple< unsigned int, unsigned int, unsigned int > + Triple; ///< A triangle of the surface. + typedef std::list< Triple > Tripleset; ///< A collection of triples. + + typedef typename Search_tree::iterator Point_iterator; ///< An iterator over the points. + typedef typename Search_tree::const_iterator Point_const_iterator; ///< A constant iterator over the points. + +private: + typedef typename std::vector Pointset; ///< A collection of points. + +private: + Search_tree _tree; // To quickly search for nearest neighbors. + + Random _generator; // For sampling random points. + + unsigned int _mean_neighbors; // The number of nearest neighbors in the mean neighborhood. + unsigned int _samples; // The number of sample points for estimating the mean neighborhood. + + Scalar _squared_radius; // The squared mean neighborhood radius. + + // The shape must be a pointer, because the alpha of + // a Fixed_alpha_shape_3 can only be set at + // construction and its assignment operator is private. + Shape* _shape; + + // The surface. If the surface is collected per shell, + // consecutive triples belong to the same shell and + // different shells are separated by a (0,0,0) triple. + Tripleset _surface; + +private: + void clear_tree() { _tree.clear(); } + +public: + /// Default constructor. + /** \param sq_radius the squared radius of the + * neighborhood size. If this value is negative when + * the point set is smoothed or when the surface is computed, + * the neighborhood size will be computed automatically. + */ + Scale_space_surface_reconstructer_3(unsigned int neighbors = 30, unsigned int samples = 200, Scalar sq_radius = -1 ): _mean_neighbors(neighbors), _samples(samples), _squared_radius( sq_radius ), _shape(0) {} + ~Scale_space_surface_reconstructer_3() { deinit_shape(); } + +private: + void deinit_shape() { if( _shape != 0 ) { delete _shape; _shape = 0; } } + +public: + Point_iterator scale_space_begin() { return _tree.begin(); } + Point_iterator scale_space_end() { return _tree.end(); } + + Point_const_iterator scale_space_begin() const { return _tree.begin(); } + Point_const_iterator scale_space_end() const { return _tree.begin(); } + + /// Check whether the spatial structure has been constructed. + /** \return true if the structure exists and false otherwise. + * \sa construct_shape(InputIterator start, InputIterator end). + */ + bool has_surface() const { return _shape != 0; } + + void clear_surface() { + if( has_surface() ) { + _shape->clear(); + } + } + + /// Insert a collection of points. + /** \tparam InputIterator an iterator over a collection of points. + * The iterator must point to a Point type. + * \param start an iterator to the first point of the collection. + * \param end a past-the-end iterator for the point collection. + * \sa compute_surface(InputIterator start, InputIterator end). + */ + template < class InputIterator > + void insert_points( InputIterator start, InputIterator end ) { + _tree.insert( start, end ); + } + + void clear() { + clear_tree(); + clear_surface(); + } + + /// Estimate the mean neighborhood size based on a number of sample points. + /** A neighborhood size is expressed as the radius of the smallest + * ball centered on a point such that the ball contains at least + * a specified number of points. + * + * The mean neighborhood size is then the mean of these radii, + * taken over a number of point samples. + * \return the estimated mean neighborhood size. + * \sa handlePoint(const Point& p). + * \sa operator()(InputIterator start, InputIterator end). + */ + Scalar estimate_mean_neighborhood( unsigned int neighbors = 30, unsigned int samples = 200 ) { + Kernel::Compute_squared_distance_3 squared_distance = Kernel().compute_squared_distance_3_object(); + + _mean_neighbors = neighbors; + _samples = samples; + + unsigned int handled = 0; + unsigned int checked = 0; + Scalar radius = 0; + + if( !_tree.is_built() ) + _tree.build(); + + for (typename Search_tree::const_iterator it = _tree.begin(); it != _tree.end(); ++it) { + 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 ) ) ); + ++checked; + } + ++handled; + } + radius /= double(checked); + + set_mean_neighborhood( radius ); + return radius; + } + + + template < class InputIterator > + Scalar estimate_mean_neighborhood(InputIterator start, InputIterator end, unsigned int neighbors = 30, unsigned int samples = 200) { + insert_points(start, end); + return estimate_mean_neighborhood(neighbors, samples); + } + + // -1 if not yet set. + Scalar get_squared_mean_neighborhood() const { return _squared_radius; } + void set_mean_neighborhood( const Scalar& radius ) { + _squared_radius = radius * radius; + if( has_surface() ) + _shape = Shape_generator().construct( _shape, _squared_radius ); + } + bool has_squared_mean_neighborhood() const { + return sign( _squared_radius ) == POSITIVE; + } + + unsigned int get_mean_neighbors() const { return _mean_neighbors; } + unsigned int get_number_samples() const { return _samples; } + void set_mean_neighbors(unsigned int neighbors) { _mean_neighbors = neighbors;} + void set_number_samples(unsigned int samples) { _samples = samples; } + + + /// Compute a number of iterations of scale-space transforming. + /** If earlier iterations have been computed, calling smooth_scale_space() + * will add more iterations. + * + * If the mean neighborhood is negative, it will be computed first. + * \param iterations the number of iterations to perform. + */ + void smooth_scale_space(unsigned int iterations = 1) { + + typedef std::vector 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 _ENV::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 ); + + Pointset points; + points.assign( _tree.begin(), _tree.end() ); + _tree.clear(); + + // 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() ); + if(!tree.is_built()) + tree.build(); + + // Collect the number of neighbors of each point. + // This can be done parallel. + CountVec neighbors(tree.size(), 0); + Kernel::Compare_squared_distance_3 compare; + p_size_t count = tree.size(); // openMP can only use signed variables. + const Scalar squared_radius = _squared_radius; // openMP can only use local variables. +#pragma omp parallel for shared(count,tree,points,squared_radius,neighbors) firstprivate(compare) + 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) + ++neighbors[i]; + } + + // Construct a mapping from each point to its index. + PIMap indices; + p_size_t index = 0; + for( Pointset::const_iterator pit = points.begin(); pit != points.end(); ++pit, ++index) + indices[ *pit ] = index; + + // Compute the tranformed point locations. + // This can be done parallel. +#pragma omp parallel for shared(count,neighbors,tree,squared_radius) firstprivate(compare) + for (p_size_t i = 0; i < count; ++i) { + // If the neighborhood is too small, the vertex is not moved. + if (neighbors[i] < 4) + continue; + + // 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]); + wts(column) = 1.0 / neighbors[indices[nit->first]]; + } + + // Construct the barycenter. + Vector3 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(); + for (column = 0; column < pts.cols(); ++column) + covariance += wts.matrix()(column) * pts.col(column) * pts.col(column).transpose(); + + // Construct the Eigen system. + EigenSolver solver(covariance); + + // If the Eigen system does not converge, the vertex is not moved. + if (solver.info() != Eigen::Success) + continue; + + // 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! + CGAL_assertion(false); + continue; + } + Vector3 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); + + points[i] = Point(per(0), per(1), per(2)); + } + } + + _tree.insert( points.begin(), points.end() ); + } + + + /// Compute a number of transform iterations on a collection of points. + /** This method is equivalent to running [insert_points(start, end)](\ref insert_points) + * followed by [smooth_scale_space(iterations)](\ref smooth_scale_space). + * \tparam InputIterator an iterator over a collection of points. + * The iterator must point to a Point type. + * \param start an iterator to the first point of the collection. + * \param end a past-the-end iterator for the point collection. + * \param iterations the number of iterations to perform. + * \sa operator()(InputIterator start, InputIterator end). + * \sa assign(InputIterator start, InputIterator end). + * \sa iterate(unsigned int iterations). + */ + template < class InputIterator > + void smooth_scale_space(InputIterator start, InputIterator end, unsigned int iterations = 1) { + insert_points(start, end); + smooth_scale_space(iterations); + } + +private: + // Once a facet is added to the surface, it is marked as handled. + inline bool is_handled( Cell_handle c, unsigned int li ) const { + switch( li ) { + case 0: return ( c->info()&1 ) != 0; + case 1: return ( c->info()&2 ) != 0; + case 2: return ( c->info()&4 ) != 0; + case 3: return ( c->info()&8 ) != 0; + } + return false; + } + inline bool is_handled( const Facet& f ) const { return is_handled( f.first, f.second ); } + inline void set_handled( Cell_handle c, unsigned int li ) { + switch( li ) { + case 0: c->info() |= 1; return; + case 1: c->info() |= 2; return; + case 2: c->info() |= 4; return; + case 3: c->info() |= 8; return; + } + } + inline void set_handled( Facet f ) { set_handled( f.first, f.second ); } + +public: + + // make new construct_shape() method for when pts already known... + void construct_shape() { + construct_shape( scale_space_begin(), scale_space_end() ); + } + + /*// For if you already have a Delaunay triangulation of the points. + void construct_shape(Shape_generator::Structure& tr ) { + deinit_shape(); + if( !has_squared_mean_neighborhood() ) + estimate_mean_neighborhood( _mean_neighbors, _samples ); + _shape = Shape_generator().construct( *tr, r2 ); + }*/ + + // If you don't want to smooth the point set. + /// Construct a new spatial structure on a set of sample points. + /** \tparam InputIterator an iterator over the point sample. + * The iterator must point to a Point type. + * \param start an iterator to the first point of the collection. + * \param end a past-the-end iterator for the point collection. + * \sa is_constructed(). + */ + template < class InputIterator > + void construct_shape(InputIterator start, InputIterator end) { + deinit_shape(); + if( !has_squared_mean_neighborhood() ) + estimate_mean_neighborhood( _mean_neighbors, _samples ); + _shape = Shape_generator().construct( start, end, _squared_radius ); + } + +private: + Triple ordered_facet_indices( const Facet& f ) const { + if( (f.second&1) == 0 ) + return Triple( f.first->vertex( (f.second+2)&3 )->info(), + f.first->vertex( (f.second+1)&3 )->info(), + f.first->vertex( (f.second+3)&3 )->info() ); + else + return Triple( f.first->vertex( (f.second+1)&3 )->info(), + f.first->vertex( (f.second+2)&3 )->info(), + f.first->vertex( (f.second+3)&3 )->info() ); + } + + /// Collect the triangles of one shell of the surface. + /** A shell is a maximal collection of triangles that are + * connected and locally facing towards the same side of the + * surface. + * \tparam OutputIterator an output iterator for a collection + * of triangles. The iterator must point to a Triangle type. + * \param c a cell touching a triangle of the shell from the + * outside of the object. + * \param li the index of the vertex of the cell opposite to + * the triangle touched. + * \param out an iterator to place to output the triangles. + * \sa collect_shell(const Facet& f, OutputIterator out). + * \sa collect_surface(OutputIterator out). + */ + void collect_shell( Cell_handle c, unsigned int li ) { + // Collect one surface mesh from the alpha-shape in a fashion similar to ball-pivoting. + // Invariant: the facet is regular or singular. + + // To stop stack overflows: use own stack. + std::stack stack; + stack.push( Facet(c, li) ); + + Facet f; + Cell_handle n, p; + int ni, pi; + Vertex_handle a; + Classification_type cl; + bool processed = false; + while( !stack.empty() ) { + f = stack.top(); + stack.pop(); + + // Check if the cell was already handled. + // Note that this is an extra check that in many cases is not necessary. + if( is_handled(f) ) + continue; + + // The facet is part of the surface. + CGAL_triangulation_assertion( !_shape->is_infinite(f) ); + set_handled(f); + + // Output the facet as a triple of indices. + _surface.push_back( ordered_facet_indices(f) ); + + if( Shells::value ) { + // Pivot over each of the facet's edges and continue the surface at the next regular or singular facet. + for( unsigned int i = 0; i < 4; ++i ) { + // Skip the current facet. + if( i == f.second || is_handled(f.first, i) ) + continue; + + // Rotate around the edge (starting from the shared facet in the current cell) until a regular or singular facet is found. + n = f.first; + ni = i; + a = f.first->vertex(f.second); + cl = _shape->classify( Facet(n, ni) ); + while( cl != Shape::REGULAR && cl != Shape::SINGULAR ) { + p = n; + n = n->neighbor(ni); + ni = n->index(a); + pi = n->index(p); + a = n->vertex(pi); + cl = _shape->classify( Facet(n, ni) ); + } + + // Continue the surface at the next regular or singular facet. + stack.push( Facet(n, ni) ); + } + processed = true; + } + } + + // We indicate the end of a shell by the (0,0,0) triple. + if( Shells::value && processed ) + _surface.push_back( Triple(0,0,0) ); + } + + /// Collect the triangles of one shell of the surface. + /** A shell is a maximal collection of triangles that are + * connected and locally facing towards the same side of the + * surface. + * \tparam OutputIterator an output iterator for a collection + * of triangles. The iterator must point to a Triangle type. + * \param f a facet of the shell as seen from the outside. + * \param out an iterator to place to output the triangles. + * \sa collect_shell(Cell_handle c, unsigned int li, OutputIterator out). + * \sa collect_surface(OutputIterator out). + */ + void collect_shell( const Facet& f ) { + collect_shell( f.first, f.second ); + } + +public: + /// Compute a surface mesh from a collection of points. + /** The surface mesh indicates the boundary of a sampled + * object. The point sample must be dense enough, such + * that any region that fits an empty ball with a given + * radius can be assumed to be ouside the object. + * + * The resulting surface need not be manifold and in + * many cases where the points sample the surface of + * an object, the surface computed will contain both + * an 'outward-facing' and an 'inward-facing' surface. + * + * However, this surface cannot have holes and its + * triangles are all oriented towards the same side of + * the surface. This orientation is expressed using the + * 'right-hand rule' on the ordered vertices of each + * triangle. + * + * The computed triangles will be returned ordered such + * that a connected 'shell' contains a set of consecutive + * triangles in the output. For this purpose, a 'shell' + * is a maximal collection of triangles that are connected + * and locally facing towards the same side of the surface. + * \tparam Kernel the geometric traits class. This class + * specifies, amongst other things, the number types and + * predicates used. + * \tparam Vb the vertex base used in the internal data + * structure that spatially orders the point set. + * \tparam Cb the cell base used in the internal data + * structure that spatially orders the point set. + */ + /// Collect the triangles of the complete surface. + /** These triangles are given ordered per shell. + * A shell is a maximal collection of triangles that are + * connected and locally facing towards the same side of the + * surface. + * \tparam OutputIterator an output iterator for a collection + * of triangles. The iterator must point to a Triangle type. + * \param out an iterator to place to output the triangles. + * \sa collect_shell(const Facet& f, OutputIterator out). + * \sa collect_shell(Cell_handle c, unsigned int li, OutputIterator out). + */ + template < class OutputIterator > + OutputIterator collect_surface( OutputIterator out ) { + if( !has_squared_mean_neighborhood() ) + estimate_mean_neighborhood( _mean_neighbors, _samples ); + + // Collect all surface meshes from the alpha-shape in a fashion similar to ball-pivoting. + // Reset the facet handled markers. + for( All_cells_iterator cit = _shape->all_cells_begin(); cit != _shape->all_cells_end(); ++cit ) + cit->info() = 0; + + // We check each of the facets: if it is not handled and either regular or singular, + // we start collecting the next surface from here. + Facet m; + int ns = 0; + for( Finite_facets_iterator fit = _shape->finite_facets_begin(); fit != _shape->finite_facets_end(); ++fit ) { + m = _shape->mirror_facet( *fit ); + switch( _shape->classify( *fit ) ) { + case Shape::REGULAR: + if( !is_handled(*fit) && !is_handled(m) ) + ++ns; + // Build a surface from the outer cell. + if( _shape->classify(fit->first) == Shape::EXTERIOR ) + collect_shell( *fit ); + else + collect_shell( m ); + break; + case Shape::SINGULAR: + if( !is_handled( *fit ) ) + ++ns; + // Build a surface from both incident cells. + collect_shell( *fit ); + if( !is_handled(m) ) + ++ns; + collect_shell( m ); + break; + } + } + + return out; + } + + /// Construct the surface triangles from a set of sample points. + /** These triangles are given ordered per shell. + * A shell is a maximal collection of triangles that are + * connected and locally facing towards the same side of the + * surface. + * + * This method is equivalent to running [construct_shape(start, end)](\ref construct_shape) + * followed by [collect_surface(out)](\ref collect_surface). + * \tparam InputIterator an iterator over the point sample. + * The iterator must point to a Point type. + * \tparam OutputIterator an output iterator for a collection + * of triangles. The iterator must point to a Triangle type. + * \param start an iterator to the first point of the collection. + * \param end a past-the-end iterator for the point collection. + * \param out an iterator to place to output the triangles. + * \sa construct_shape(InputIterator start, InputIterator end). + * \sa collect_surface(OutputIterator out). + */ + template < class InputIterator, class OutputIterator > + OutputIterator collect_surface( InputIterator start, InputIterator end, OutputIterator out ) { + construct_shape( start, end ); + return collect_surface( out ); + } + + template < class InputIterator > + void collect_surface( InputIterator start, InputIterator end ) { + construct_shape( start, end ); + collect_surface( std::back_inserter(_surface) ); + } + + void collect_surface() { + construct_shape(); + collect_surface( std::back_inserter(_surface) ); + } + + const Tripleset& surface() const { return _surface; } + Tripleset& surface() { return _surface; } + + + + /// Compute a smoothed surface mesh from a collection of points. + /** This is equivalent to calling insert_points(), + * iterate(), and compute_surface(). + * + * After this computation, the surface can be accessed using + * surface(). + * \tparam InputIterator an iterator over the point collection. + * The iterator must point to a Point type. + * \param start an iterator to the first point of the collection. + * \param end a past-the-end iterator for the point collection. + * + * \sa surface(). + * \sa surface() const. + */ + template < class InputIterator > + void compute_surface(InputIterator start, InputIterator end, unsigned int neighbors = 30, unsigned int iterations = 4, unsigned int samples = 200) { + // Compute the radius for which the mean ball would contain the required number of neighbors. + clear(); + insert_points( start, end ); + + _mean_neighbors = neighbors; + _samples = samples; + + // Smooth the scale space. + smooth_scale_space( iterations ); + + // Mesh the perturbed points. + collect_surface(); + } +}; // class Scale_space_surface_reconstructer_3 + + +template< typename T1, typename T2, typename T3 > +std::ostream& +operator<<( std::ostream& os, const Triple< T1, T2, T3 >& t ) { + return os << t.first << " " << t.second << " " << t.third; +} + +template< typename T1, typename T2, typename T3 > +std::istream& +operator>>( std::istream& is, Triple< T1, T2, T3 >& t ) { + return is >> t.first >> t.second >> t.third; +} + +} // namespace CGAL + +#endif // SCALE_SPACE_SURFACE_CONSTRUCTER \ No newline at end of file diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_transform.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_transform.h deleted file mode 100644 index 8f6e011adac..00000000000 --- a/Scale_space_reconstruction_3/include/CGAL/Scale_space_transform.h +++ /dev/null @@ -1,192 +0,0 @@ -//A method to transform a point set. -//Copyright (C) 2013 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 SCALE_SPACE_PERTURB -#define SCALE_SPACE_PERTURB - -#include - -#include -#include - -#include -#include - -#include - -#include "check3264.h" - - -// Transform a point set by projecting each point onto the best planar fit of the weighted local neighborhood. -template < typename Kernel > -class Scale_space_transform { - typedef typename Kernel::FT Scalar; - -public: - typedef typename Kernel::Point_3 Point; - typedef typename std::vector Collection; - -private: - typedef CGAL::Search_traits_3 Tree_traits; - - typedef typename CGAL::Orthogonal_incremental_neighbor_search Neighbor_search; - typedef typename Neighbor_search::iterator Neighbor_iterator; - typedef typename Neighbor_search::Tree Neighbor_Tree; - - typedef typename Kernel::Compare_squared_distance_3 Compare_squared_distance; - - Scalar _radius2; - Collection _points; - -public: - Scale_space_transform(const Scalar& r): _radius2(r) {} - - void set_radius2(const Scalar& r) {_radius2 = r;} - Scalar get_radius2() const {return _radius2;} - - template < class InputIterator > - void assign(InputIterator start, InputIterator end) { - std::cout << "SST: Assign" << std::endl; - _points.assign(start, end); - } - - Collection& get_collection() {return _points;} - const Collection& get_collection() const {return _points;} - - void iterate() { - std::cout << "SST: Iterate" << std::endl; - typedef std::vector 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 _ENV::s_ptr_type p_size_t; - - // This method must be called after filling the point collection. - CGAL_assertion(!_points.empty()); - Collection points; - points.swap(_points); - Scalar radius = _radius2; - - // Construct a search tree of the points. - Neighbor_Tree tree(points.begin(), points.end()); - if(!tree.is_built()) - tree.build(); - - // Collect the number of neighbors of each point. - // This can be done parallel. - CountVec neighbors(points.size(), 0); - Compare_squared_distance compare; - p_size_t count = points.size(); // openMP can only use signed variable -#pragma omp parallel for shared(count,tree,points,radius,neighbors) firstprivate(compare) - for (p_size_t i = 0; i < count; ++i) { - // Iterate over the neighbors until the first one is found that is too far. - Neighbor_search search(tree, points[i]); - for (Neighbor_iterator nit = search.begin(); nit != search.end() && compare(points[i], nit->first, radius) != CGAL::LARGER; ++nit) - ++neighbors[i]; - } - - // Construct a mapping from each point to its index. - PIMap indices; - for (p_size_t i = 0; i < count; ++i) - indices[points[i]] = i; - - // Compute the tranformed point locations. - // This can be done parallel. -#pragma omp parallel for shared(count,neighbors,tree,points,radius) firstprivate(compare) - for (p_size_t i = 0; i < count; ++i) { - // If the neighborhood is too small, the vertex is not moved. - if (neighbors[i] < 4) - continue; - - // Collect the vertices within the ball and their weights. - Neighbor_search search(tree, points[i]); - Matrix3D pts(3, neighbors[i]); - Array1D wts(1, neighbors[i]); - int column = 0; - for (Neighbor_iterator nit = search.begin(); nit != search.end() && compare(points[i], nit->first, 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]); - wts(column) = 1.0 / neighbors[indices[nit->first]]; - } - - // Construct the barycenter. - Vector3 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(); - for (column = 0; column < pts.cols(); ++column) - covariance += wts.matrix()(column) * pts.col(column) * pts.col(column).transpose(); - - // Construct the Eigen system. - EigenSolver solver(covariance); - - // If the Eigen system does not converge, the vertex is not moved. - if (solver.info() != Eigen::Success) - continue; - - // 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! - CGAL_assertion(false); - continue; - } - Vector3 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); - - points[i] = Point(per(0), per(1), per(2)); - } - - _points.swap(points); - } - - void iterate(unsigned int iterations) { - for (unsigned int iter = 0; iter < iterations; ++iter) - iterate(); - } - - template < class InputIterator > - void operator()(InputIterator start, InputIterator end) { - _points.assign(start, end); - iterate(); - } - - template < class ForwardIterator > - void operator()(ForwardIterator start, ForwardIterator end, unsigned int iterations) { - // Perturb the points for a number of iterations. - _points.assign(start, end); - iterate(iterations); - } -}; // class Scale_space_transform - -#endif // SCALE_SPACE_PERTURB \ No newline at end of file diff --git a/Scale_space_reconstruction_3/include/CGAL/Surface_mesher.h b/Scale_space_reconstruction_3/include/CGAL/Surface_mesher.h deleted file mode 100644 index b71ce51bff3..00000000000 --- a/Scale_space_reconstruction_3/include/CGAL/Surface_mesher.h +++ /dev/null @@ -1,295 +0,0 @@ -//The ball pivoting algorithm implemented using the alpha-shape. -//Copyright (C) 2013 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 BALL_PIVOTING -#define BALL_PIVOTING - -#include - -#define BALL_PIVOTING_FIXED -#define BALL_PIVOTING_CONNECTED - -#ifdef BALL_PIVOTING_FIXED -#include -#include -#include -#else -#include -#include -#include -#endif - -#include "Flagged_triangulation_cell_base_3.h" - -// Collect the vertices within a ball of fixed radius. -template < typename Kernel, - typename Vb = CGAL::Triangulation_vertex_base_3, - typename Cb = CGAL::Triangulation_cell_base_3 > -class Surface_mesher { - typedef typename Kernel::FT Scalar; - - // Convert the DT types to be in an alpha-shape. - // Vertex. -#ifdef BALL_PIVOTING_FIXED - typedef CGAL::Fixed_alpha_shape_vertex_base_3 AVb; -#else - typedef CGAL::Alpha_shape_vertex_base_3 AVb; -#endif - - // Cell. - typedef CGAL::Flagged_triangulation_cell_base_3 FCb; -#ifdef BALL_PIVOTING_FIXED - typedef CGAL::Fixed_alpha_shape_cell_base_3 ACb; -#else - typedef CGAL::Alpha_shape_cell_base_3 ACb; -#endif - - // Triangulation. - typedef CGAL::Triangulation_data_structure_3 Tds; - typedef CGAL::Delaunay_triangulation_3 ADT; -#ifdef BALL_PIVOTING_FIXED - typedef CGAL::Fixed_alpha_shape_3 AS; -#else - typedef CGAL::Alpha_shape_3 AS; -#endif - -public: - typedef ADT Triangulation; - typedef AS Alpha_shape; - - typedef typename AS::Vertex_handle Vertex_handle; - typedef typename AS::Cell_handle Cell_handle; - typedef typename AS::Facet Facet; - -private: - typedef typename AS::All_cells_iterator All_cells_iterator; - typedef typename AS::Finite_facets_iterator Finite_facets_iterator; - typedef typename AS::Classification_type Classification_type; - - AS* _shape; - - Scalar radius2; - - // Once a facet is outputted, it is marked as handled. - inline bool isHandled(Cell_handle c, unsigned int li) const { - switch (li) { - case 0: return (c->flag()&1) != 0; - case 1: return (c->flag()&2) != 0; - case 2: return (c->flag()&4) != 0; - case 3: return (c->flag()&8) != 0; - } - return false; - } - inline bool isHandled(const Facet& f) const {return isHandled(f.first, f.second);} - - inline void setHandled(Cell_handle c, unsigned int li) { - switch (li) { - case 0: c->flag() |= 1; return; - case 1: c->flag() |= 2; return; - case 2: c->flag() |= 4; return; - case 3: c->flag() |= 8; return; - } - } - inline void setHandled(Facet f) {setHandled(f.first, f.second);} - -public: - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::Triangle_3 Triangle; - - Surface_mesher(const Scalar& r2): _shape(0), radius2(r2) {} - Surface_mesher(const Triangulation& tr, const Scalar& r2): radius2(r2) { -#ifdef BALL_PIVOTING_FIXED - _shape = new AS(tr, radius2); -#else - _shape = new AS(tr, radius2, AS::GENERAL); -#endif - } - ~Surface_mesher() {if (_shape != 0) {delete _shape; _shape = 0;}} - - template < class InputIterator > - void construct_shape(InputIterator start, InputIterator end) { - std::cout << "BP: Construct shape"; - ADT tr(start, end); - - if (_shape != 0) - delete _shape; -#ifdef BALL_PIVOTING_FIXED - _shape = new AS(tr, radius2); -#else - _shape = new AS(tr, radius2, AS::GENERAL); -#endif - std::cout << ": " << _shape->number_of_vertices() << " points" << std::endl; - } - bool is_constructed() const {return _shape != 0;} - void clear() {if (_shape != 0) {delete _shape; _shape = 0;}} - - AS* shape() {return _shape;} - - Scalar get_radius2() const {return radius2;} - void set_radius2(const Scalar& r2) { - radius2 = r2; - if (_shape != 0) { - ADT tr; - tr.swap(*_shape); - - delete _shape; -#ifdef BALL_PIVOTING_FIXED - _shape = new AS(tr, radius2); -#else - _shape = new AS(tr, radius2, AS::GENERAL); -#endif - } - } - - template < class OutputIterator > - OutputIterator collect_shell(Cell_handle c, unsigned int li, OutputIterator out) { - // Collect one surface mesh from the alpha-shape in a fashion similar to ball-pivoting. - // Invariant: the facet is regular or singular. - - // To stop stack overflows: use own stack. - std::stack stack; - stack.push(Facet(c, li)); - - Facet f; - Cell_handle n, p; int ni, pi; - Vertex_handle a; - Classification_type cl; - while (!stack.empty()) { - f = stack.top(); - stack.pop(); - - // Check if the cell was already handled. - // Note that this is an extra check that in many cases is not necessary. - if (isHandled(f)) - continue; - - // The facet is part of the surface. - CGAL_triangulation_assertion(!_shape->is_infinite(f)); - *out++ = f; - setHandled(f); - -#ifdef BALL_PIVOTING_CONNECTED - // Pivot over each of the facet's edges and continue the surface at the next regular or singular facet. - for (unsigned int i = 0; i < 4; ++i) { - // Skip the current facet. - if (i == f.second || isHandled(f.first, i)) - continue; - - // Rotate around the edge (starting from the shared facet in the current cell) until a regular or singular facet is found. - n = f.first; - ni = i; - a = f.first->vertex(f.second); - cl = _shape->classify(Facet(n, ni)); - while (cl != AS::REGULAR && cl != AS::SINGULAR) { - p = n; - n = n->neighbor(ni); - ni = n->index(a); - pi = n->index(p); - a = n->vertex(pi); - cl = _shape->classify(Facet(n, ni)); - } - - // Continue the surface at the next regular or singular facet. - stack.push(Facet(n, ni)); - } -#endif - } - - return out; - } - - template < class OutputIterator > - OutputIterator collect_shell(const Facet& f, OutputIterator out) { - return collect_shell(f.first, f.second, out); - } - - template < class OutputIterator > - OutputIterator collect_surface(OutputIterator out) { - std::cout << "BP: Collect surfaces" << std::endl; - // Collect all surface meshes from the alpha-shape in a fashion similar to ball-pivoting. - // Reset the facet handled markers. - for (All_cells_iterator cit = _shape->all_cells_begin(); cit != _shape->all_cells_end(); ++cit) - cit->flag() = 0; - - // We check each of the facets: if it is not handled and either regular or singular, - // we start collecting the next surface from here. - Facet m; - int ns = 0; - for (Finite_facets_iterator fit = _shape->finite_facets_begin(); fit != _shape->finite_facets_end(); ++fit) { - m = _shape->mirror_facet(*fit); - switch (_shape->classify(*fit)) { - case AS::REGULAR: - if (!isHandled(*fit) && !isHandled(m)) - ++ns; - // Build a surface from the outer cell. - if (_shape->classify(fit->first) == AS::EXTERIOR) - collect_shell(*fit, out); - else - collect_shell(m, out); - break; - case AS::SINGULAR: - if (!isHandled(*fit)) - ++ns; - // Build a surface from both incident cells. - collect_shell(*fit, out); - if (!isHandled(m)) - ++ns; - collect_shell(m, out); - break; - } - } - - std::cout << ": " << ns << std::flush; - - return out; - } - - template < class InputIterator, class OutputIterator > - OutputIterator operator()(InputIterator start, InputIterator end, OutputIterator out) { - construct_shape(start, end); - return collect_surface(out); - } - - bool locate_vertex(const Point& p, Vertex_handle& v, Cell_handle hint = Cell_handle()) const { - typename AS::Locate_type lt; int li, lj; - Cell_handle c = _shape->locate(p, lt, li, lj, hint); - if (lt != AS::VERTEX) - return false; - v = c->vertex(li); - return true; - } - - void ordered_vertices(const Facet& f, Vertex_handle& v0, Vertex_handle& v1, Vertex_handle& v2) { - if ((f.second&1) == 0) { - v0 = f.first->vertex((f.second+2)&3); - v1 = f.first->vertex((f.second+1)&3); - v2 = f.first->vertex((f.second+3)&3); - } - else { - v0 = f.first->vertex((f.second+1)&3); - v1 = f.first->vertex((f.second+2)&3); - v2 = f.first->vertex((f.second+3)&3); - } - } - - Triangle oriented_triangle(Cell_handle c, unsigned int li) const {return _shape->triangle(c, li);} - Triangle oriented_triangle(const Facet& f) const {return _shape->triangle(f);} -}; // class Surface_mesher - -#endif // BALL_PIVOTING diff --git a/Scale_space_reconstruction_3/include/CGAL/check3264.h b/Scale_space_reconstruction_3/include/CGAL/check3264.h index c025ae17679..1d837a4d893 100644 --- a/Scale_space_reconstruction_3/include/CGAL/check3264.h +++ b/Scale_space_reconstruction_3/include/CGAL/check3264.h @@ -16,9 +16,9 @@ // // Author(s): Thijs van Lankveld - -template struct _Env{ typedef void s_ptr_type; typedef void ptr_type; static const int ptr_size = sizeof(ptr_type)*8; }; + +template struct _Env{ typedef void s_ptr_type; typedef void ptr_type; static const int ptr_size = sizeof(ptr_type)*8; static const bool x64 = false; }; template<> struct _Env<4> { typedef int s_ptr_type; typedef unsigned int ptr_type; static const int ptr_size = sizeof(ptr_type)*8; }; -template<> struct _Env<8> { typedef long long s_ptr_type; typedef unsigned long long ptr_type; static const int ptr_size = sizeof(ptr_type)*8; }; +template<> struct _Env<8> { typedef long long s_ptr_type; typedef unsigned long long ptr_type; static const int ptr_size = sizeof(ptr_type)*8; static const bool x64 = true;}; typedef _Env _ENV; \ No newline at end of file