songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Combined into one class. 

This did mess up the documentation a bit..
This commit is contained in:
Thijs van Lankveld 2014-03-27 17:13:34 +01:00
parent e016b1360b
commit 273d4264d8
15 changed files with 1595 additions and 1840 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

113
Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Mean_neighborhood_ball.h
View File

@ -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 <http://www.gnu.org/licenses/>.
//
// Author(s): Thijs van Lankveld
#ifndef MEAN_NEIGHBORHOOD_BALL
#define MEAN_NEIGHBORHOOD_BALL
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Random.h>
/// 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

110
Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Scale_space_surface_constructer.h
View File

@ -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 <http://www.gnu.org/licenses/>.
//
// Author(s): Thijs van Lankveld
#ifndef SCALE_SPACE_SURFACE_CONSTRUCTER
#define SCALE_SPACE_SURFACE_CONSTRUCTER
#include <iostream>
#include <list>
#include <map>
#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<Kernel>,
typename Cb = CGAL::Triangulation_cell_base_3<Kernel> >
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<Point> 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

139
Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Scale_space_transform.h
View File

@ -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 <http://www.gnu.org/licenses/>.
//
// Author(s): Thijs van Lankveld
#ifndef SCALE_SPACE_PERTURB
#define SCALE_SPACE_PERTURB
#include <omp.h>
#include <map>
#include <vector>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_incremental_neighbor_search.h>
#include <Eigen/Dense>
#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<Point> 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

276
Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/CGAL/Surface_mesher.h
View File

@ -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 <http://www.gnu.org/licenses/>.
//
// Author(s): Thijs van Lankveld
#ifndef BALL_PIVOTING
#define BALL_PIVOTING
#include <CGAL/Delaunay_triangulation_3.h>
#define BALL_PIVOTING_FIXED
#define BALL_PIVOTING_CONNECTED
#ifdef BALL_PIVOTING_FIXED
#include <CGAL/Fixed_alpha_shape_3.h>
#include <CGAL/Fixed_alpha_shape_vertex_base_3.h>
#include <CGAL/Fixed_alpha_shape_cell_base_3.h>
#else
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>
#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<Kernel>,
typename Cb = CGAL::Triangulation_cell_base_3<Kernel> >
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<Kernel, Vb> AVb;
#else
typedef CGAL::Alpha_shape_vertex_base_3<Kernel, Vb> AVb;
#endif
/// The cell type.
typedef CGAL::Flagged_triangulation_cell_base_3<Kernel, Cb> FCb;
#ifdef BALL_PIVOTING_FIXED
typedef CGAL::Fixed_alpha_shape_cell_base_3<Kernel, FCb> ACb;
#else
typedef CGAL::Alpha_shape_cell_base_3<Kernel, FCb> ACb;
#endif
// Triangulation.
typedef CGAL::Triangulation_data_structure_3<AVb, ACb> Tds; ///< The data structure that stores the point set.
typedef CGAL::Delaunay_triangulation_3<Kernel, Tds> 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<ADT> AS;
#else
typedef CGAL::Alpha_shape_3<ADT> 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

45
Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/CMakeLists.txt Normal file
View File

@ -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()

60
Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/quad.h
View File

@ -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 <http://www.gnu.org/licenses/>.
//
// 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<Scalar>& q) {
os << q.first << " " << q.second << " " << q.third << " " << q.fourth;
return os;
}
typedef Quad<int> int4;
typedef Quad<long> long4;
typedef Quad<double> double4;
typedef Quad<float> float4;
#endif // SCALE_SPACE_QUAD

178
Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/readFile.h
View File

@ -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 <http://www.gnu.org/licenses/>.
//
// Author(s): Thijs van Lankveld
#ifndef READ_FILE
#define READ_FILE
#include <fstream>
#include <iostream>
/*#include <osgDB/ReaderWriter>
#include <osgDB/ReadFile>*/
#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<Point>();
}
// Collect normals.
osg::Array* normals = geom->getNormalArray();
if (geom->getNormalBinding() == osg::Geometry::BIND_OVERALL) {
normals->accept(0, vc);
Normal n = vc.convert3<Normal>();
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<Normal>();
}
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<Color>();
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<Color>();
}
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<osg::Node> scene = osgDB::readNodeFile(file);
PointCollecter<Kernel, Color, PointOutput, NormalOutput, ColorOutput> pc(po, no, co);
scene->accept(pc);
return scene.valid();
}*/
#endif // READ_FILE

216
Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/scale_space.cpp
View File

@ -18,136 +18,120 @@
#include <algorithm>
#include <vector>
#include <fstream>
#include <iostream>
#include <map>
#include <string>
#include <vector>
#include <CGAL/Scale_space_surface_reconstructer_3.h>
//#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include "readFile.h"
#include "writeFile.h"
#include "quad.h"
#include "Scale_space_surface_constructer.h"
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;
if (argc > 1)
input = argv[1];
if (argc > 2)
outFile = argv[2];
if (argc > 3)
sscanf(argv[3], "%u", &nn);
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 CGAL::Scale_space_surface_reconstructer_3< Kernel > Reconstructer;
typedef std::vector<Point_3> PointCollection;
typedef std::vector<Vector_3> VectorCollection;
typedef float4 Color;
typedef std::vector<Color> ColorCollection;
typedef Reconstructer::Point Point;
typedef std::vector< Point > Pointset;
// 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<Kernel>(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<Kernel, Color>(input, std::back_inserter(points), std::back_inserter(normals), std::back_inserter(colors));*/
if (!loaded) {
std::cout << "Error loading input." << std::endl;
exit(1);
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 );
}
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());
}
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) {
std::string base = "kitten";
unsigned int neighbors = 30;
unsigned int iterations = 4;
unsigned int samples = 200;
if (argc > 1)
base = argv[1];
if (argc > 2)
sscanf(argv[2], "%u", &neighbors);
if (argc > 3)
sscanf(argv[3], "%u", &iterations);
if (argc > 4)
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<Kernel> Scale_space_surface_constructer;
typedef Scale_space_surface_constructer::Triangle Triangle;
std::list<Triangle> 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;
// 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);
}
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<Triangle>::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<Point_3, Color> 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<Triangle>::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<osg::Switch> 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;
std::cout << " written." << std::endl;
}

185
Scale_space_reconstruction_3/examples/Scale_space_reconstruction_3/writeFile.h
View File

@ -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 <http://www.gnu.org/licenses/>.
//
// Author(s): Thijs van Lankveld
#ifndef WRITE_FILE
#define WRITE_FILE
#include <map>
/*#include <osg/CullFace>
#include <osg/Material>
#include <osg/Switch>
#include <osgDB/ReaderWriter>
#include <osgDB/WriteFile>
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<osg::Geode> geode = new osg::Geode;
osg::ref_ptr<osg::Geometry> geom = new osg::Geometry;
osg::ref_ptr<osg::Vec3Array> 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<osg::Vec3Array> 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<osg::Vec4Array> 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<osg::Geode> geode = new osg::Geode;
osg::ref_ptr<osg::Geometry> geom = new osg::Geometry;
osg::ref_ptr<osg::Vec3Array> vert = new osg::Vec3Array;
vert->reserve(3 * size);
osg::ref_ptr<osg::Vec3Array> 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<osg::Vec4Array> 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<osgDB::ReaderWriter::Options> 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<typename NormalIterator::value_type::R> Point;
// Construct a mapping from the points to their indices.
std::map<Point, unsigned int> 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<Point, unsigned int>::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

104
Scale_space_reconstruction_3/include/CGAL/Mean_neighborhood_ball.h
View File

@ -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 <http://www.gnu.org/licenses/>.
//
// Author(s): Thijs van Lankveld
#ifndef MEAN_NEIGHBORHOOD_BALL
#define MEAN_NEIGHBORHOOD_BALL
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Random.h>
// 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<Kernel> Tree_traits;
typedef typename CGAL::Orthogonal_k_neighbor_search<Tree_traits> 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

726
Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_constructer.h
View File

@ -20,103 +20,687 @@
#ifndef SCALE_SPACE_SURFACE_CONSTRUCTER
#define SCALE_SPACE_SURFACE_CONSTRUCTER
#include <omp.h>
#include <iostream>
#include <list>
#include <map>
#include <vector>
#include "Mean_neighborhood_ball.h"
#include "Scale_space_transform.h"
#include "Surface_mesher.h"
#include <boost/iterator/transform_iterator.hpp>
// 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<Kernel>,
typename Cb = CGAL::Triangulation_cell_base_3<Kernel> >
class Scale_space_surface_constructer {
typedef typename Mean_neighborhood_ball<Kernel> ComputeRadius;
typedef typename Scale_space_transform<Kernel> Scale_space_transform;
typedef typename Surface_mesher<Kernel, Vb, Cb> Surface_mesher;
#include <CGAL/utility.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_incremental_neighbor_search.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Random.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Triangulation_cell_base_with_info_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>
#include <CGAL/Fixed_alpha_shape_3.h>
#include <CGAL/Fixed_alpha_shape_vertex_base_3.h>
#include <CGAL/Fixed_alpha_shape_cell_base_3.h>
#include <Eigen/Dense>
#include <CGAL/check3264.h>
namespace CGAL {
//a functor that returns a std::pair<Point,unsigned>.
//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 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;
typedef typename std::vector<Point> PointCollection;
private:
PointCollection _moved;
ComputeRadius mnb;
unsigned int iterations;
Surface_mesher bp;
typedef CGAL::Alpha_shape_vertex_base_3< Kernel,
CGAL::Triangulation_vertex_base_with_info_3< unsigned int, Kernel > > Vb;
typedef CGAL::Alpha_shape_cell_base_3< Kernel,
CGAL::Triangulation_cell_base_with_info_3< unsigned int, Kernel > > Cb;
typedef CGAL::Triangulation_data_structure_3<Vb,Cb> Tds;
public:
Scale_space_surface_constructer(unsigned int iter = 1): mnb(), iterations(iter), bp(0) {}
typedef CGAL::Delaunay_triangulation_3< Kernel, Tds > Structure;
typedef CGAL::Alpha_shape_3< Structure > Shape;
void set_iterations(unsigned int iter) {iterations = iter;}
_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<Point>() ),
boost::make_transform_iterator( end, Auto_count<Point>() ),
squared_radius, Shape::GENERAL );
}
}; // class _Shape
const PointCollection& moved() const {return _moved;}
// 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;
// Input: iterators over Point, output: collection of Triangles.
template < class InputIterator, class OutputIterator >
OutputIterator operator()(InputIterator start, InputIterator end, OutputIterator out) {
typedef std::map<Point, Point> Map;
typedef std::list<Facet> List;
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;
// Compute the radius for which the mean ball would contain 30 points.
Scalar radius2 = mnb(start, end);
radius2 *= radius2;
typedef CGAL::Triangulation_data_structure_3<Vb,Cb> Tds;
// 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());
public:
typedef CGAL::Delaunay_triangulation_3< Kernel, Tds > Structure;
typedef CGAL::Fixed_alpha_shape_3< Structure > Shape;
// Perform ball-pivoting on the perturbed points.
std::list<Facet> facets;
bp = Surface_mesher(radius2);
bp(_moved.begin(), _moved.end(), std::back_inserter(facets));
_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<Point>() ),
boost::make_transform_iterator( end, Auto_count<Point>() ),
squared_radius );
}
}; // class _Shape
// 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<size_t, Vb>
// 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;
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)]);
/// 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<Point> 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(); }
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<unsigned int> CountVec;
typedef typename std::map<Point, size_t> PIMap;
typedef Eigen::Matrix<double, 3, Eigen::Dynamic> Matrix3D;
typedef Eigen::Array<double, 1, Eigen::Dynamic> Array1D;
typedef Eigen::Matrix3d Matrix3;
typedef Eigen::Vector3d Vector3;
typedef Eigen::SelfAdjointEigenSolver<Matrix3> EigenSolver;
typedef _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<Facet> 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);
template < class InputIterator, class OutputIterator >
OutputIterator collect_surface( InputIterator start, InputIterator end, OutputIterator out ) {
construct_shape( start, end );
return collect_surface( out );
}
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 >
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() );
}
}; // class Scale_space_surface_constructer
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<Point, unsigned int> 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

794
Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstructer_3.h Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
//
// Author(s): Thijs van Lankveld
#ifndef SCALE_SPACE_SURFACE_CONSTRUCTER
#define SCALE_SPACE_SURFACE_CONSTRUCTER
#include <omp.h>
#include <iostream>
#include <list>
#include <map>
#include <vector>
#include <boost/iterator/transform_iterator.hpp>
#include <CGAL/utility.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_incremental_neighbor_search.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Random.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Triangulation_cell_base_with_info_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>
#include <CGAL/Fixed_alpha_shape_3.h>
#include <CGAL/Fixed_alpha_shape_vertex_base_3.h>
#include <CGAL/Fixed_alpha_shape_cell_base_3.h>
#include <Eigen/Dense>
#include <CGAL/check3264.h>
namespace CGAL {
//a functor that returns a std::pair<Point,unsigned>.
//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<Vb,Cb> 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<Point>() ),
boost::make_transform_iterator( end, Auto_count<Point>() ),
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<Vb,Cb> 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<Point>() ),
boost::make_transform_iterator( end, Auto_count<Point>() ),
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<Point> 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<unsigned int> CountVec;
typedef typename std::map<Point, size_t> PIMap;
typedef Eigen::Matrix<double, 3, Eigen::Dynamic> Matrix3D;
typedef Eigen::Array<double, 1, Eigen::Dynamic> Array1D;
typedef Eigen::Matrix3d Matrix3;
typedef Eigen::Vector3d Vector3;
typedef Eigen::SelfAdjointEigenSolver<Matrix3> EigenSolver;
typedef _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<Facet> 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

192
Scale_space_reconstruction_3/include/CGAL/Scale_space_transform.h
View File

@ -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 <http://www.gnu.org/licenses/>.
//
// Author(s): Thijs van Lankveld
#ifndef SCALE_SPACE_PERTURB
#define SCALE_SPACE_PERTURB
#include <omp.h>
#include <map>
#include <vector>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_incremental_neighbor_search.h>
#include <Eigen/Dense>
#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<Point> Collection;
private:
typedef CGAL::Search_traits_3<Kernel> Tree_traits;
typedef typename CGAL::Orthogonal_incremental_neighbor_search<Tree_traits> 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<unsigned int> CountVec;
typedef typename std::map<Point, size_t> PIMap;
typedef Eigen::Matrix<double, 3, Eigen::Dynamic> Matrix3D;
typedef Eigen::Array<double, 1, Eigen::Dynamic> Array1D;
typedef Eigen::Matrix3d Matrix3;
typedef Eigen::Vector3d Vector3;
typedef Eigen::SelfAdjointEigenSolver<Matrix3> EigenSolver;
typedef _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

295
Scale_space_reconstruction_3/include/CGAL/Surface_mesher.h
View File

@ -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 <http://www.gnu.org/licenses/>.
//
// Author(s): Thijs van Lankveld
#ifndef BALL_PIVOTING
#define BALL_PIVOTING
#include <CGAL/Delaunay_triangulation_3.h>
#define BALL_PIVOTING_FIXED
#define BALL_PIVOTING_CONNECTED
#ifdef BALL_PIVOTING_FIXED
#include <CGAL/Fixed_alpha_shape_3.h>
#include <CGAL/Fixed_alpha_shape_vertex_base_3.h>
#include <CGAL/Fixed_alpha_shape_cell_base_3.h>
#else
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>
#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<Kernel>,
typename Cb = CGAL::Triangulation_cell_base_3<Kernel> >
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<Kernel, Vb> AVb;
#else
typedef CGAL::Alpha_shape_vertex_base_3<Kernel, Vb> AVb;
#endif
// Cell.
typedef CGAL::Flagged_triangulation_cell_base_3<Kernel, Cb> FCb;
#ifdef BALL_PIVOTING_FIXED
typedef CGAL::Fixed_alpha_shape_cell_base_3<Kernel, FCb> ACb;
#else
typedef CGAL::Alpha_shape_cell_base_3<Kernel, FCb> ACb;
#endif
// Triangulation.
typedef CGAL::Triangulation_data_structure_3<AVb, ACb> Tds;
typedef CGAL::Delaunay_triangulation_3<Kernel, Tds> ADT;
#ifdef BALL_PIVOTING_FIXED
typedef CGAL::Fixed_alpha_shape_3<ADT> AS;
#else
typedef CGAL::Alpha_shape_3<ADT> 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<Facet> 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

4
Scale_space_reconstruction_3/include/CGAL/check3264.h
View File

@ -17,8 +17,8 @@
// Author(s): Thijs van Lankveld
template<int sz> struct _Env{ typedef void s_ptr_type; typedef void ptr_type; static const int ptr_size = sizeof(ptr_type)*8; };
template<int sz> 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<sizeof(void*)> _ENV;