diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Alpha_shape_mesher.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Alpha_shape_mesher.h new file mode 100644 index 00000000000..9ee89b4129e --- /dev/null +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Alpha_shape_mesher.h @@ -0,0 +1,914 @@ +#ifndef CGAL_SCALE_SPACE_RECONSTRUCTION_3_ALPHA_SHAPE_MESHER_H +#define CGAL_SCALE_SPACE_RECONSTRUCTION_3_ALPHA_SHAPE_MESHER_H + +#include + +#include + +namespace CGAL +{ + +namespace Scale_space_reconstruction_3 +{ + +template +class Alpha_shape_mesher +{ +public: + typedef typename Geom_traits::FT FT; + typedef typename Geom_traits::Point_3 Point; ///< defines the point type. + + typedef CGAL::cpp11::array< unsigned int, 3 > Triple; ///< defines a triple of point indices indicating a triangle of the surface. +private: + typedef std::list< Triple > Tripleset; ///< defines a collection of triples. + // Note that this is a list for two reasons: iterator validity for the shell iterators, and memory requirements for the expected huge collections. + +public: +#ifdef DOXYGEN_RUNNING + typedef unspecified_type Triple_iterator; ///< defines an iterator over the triples. + typedef const unspecified_type Triple_const_iterator; ///< defines a constant iterator over the triples. +#else // DOXYGEN_RUNNING + typedef Tripleset::iterator Triple_iterator; + typedef Tripleset::const_iterator Triple_const_iterator; +#endif // DOXYGEN_RUNNING + +private: + typedef std::vector< Triple_iterator > TripleIterSet; + +private: + + // Constructing the surface. + typedef CGAL::Shape_construction_3 Shape_construction_3; + + typedef typename Shape_construction_3::Shape Shape; + typedef typename Shape_construction_3::Triangulation Triangulation; + + typedef typename Shape::Vertex_handle Vertex_handle; + typedef typename Shape::Cell_handle Cell_handle; + typedef typename Shape::Facet Facet; + typedef typename Shape::Edge Edge; + typedef std::pair VEdge; + + typedef typename Shape::Vertex_iterator Vertex_iterator; + typedef typename Shape::Cell_iterator Cell_iterator; + typedef typename Shape::Facet_iterator Facet_iterator; + typedef typename Shape::Edge_iterator Edge_iterator; + + typedef typename Shape::Finite_cells_iterator Finite_cells_iterator; + typedef typename Shape::Finite_facets_iterator Finite_facets_iterator; + typedef typename Shape::Finite_edges_iterator Finite_edges_iterator; + typedef typename Shape::Finite_vertices_iterator Finite_vertices_iterator; + + typedef typename Shape::Facet_circulator Facet_circulator; + + typedef typename Shape::All_cells_iterator All_cells_iterator; + + typedef typename Shape::Classification_type Classification_type; + + typedef std::map Map_facet_to_shell; + typedef typename CGAL::cpp11::array, 2 > Bubble; + + bool _separate_shells; + bool _force_manifold; + FT _border_angle; + + // 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. + // We want to be able to set the alpha after constructing the scale-space + // reconstructer object. + Shape* _shape; + + // The surface. If the surface is collected per shell, the triples of the + // same shell are stored consecutively. + Tripleset _surface; + + // The shells can be accessed through iterators to the surface. + TripleIterSet _shells; + + // If the surface is forced to be manifold, removed facets are stored + Tripleset _garbage; + + // Map TDS facets to shells + Map_facet_to_shell _map_f2s; + // std::map _map_f2s; + unsigned int _index; + + std::vector _bubbles; + std::map _map_f2b; + + FT _squared_radius; + +public: + + Alpha_shape_mesher (bool separate_shells = false, bool force_manifold = false, + FT border_angle = 45.) + : _separate_shells (separate_shells), + _force_manifold (force_manifold), + _border_angle (border_angle), + _shape (NULL), + _squared_radius (0.) + { + + } + + template + void operator() (InputIterator begin, InputIterator end, OutputIterator output) + { + clear_surface(); + + _shape = Shape_construction_3().construct (begin, end, _squared_radius); + + // If shells are not separated and no manifold constraint is given, + // then the quick collection of facets can be applied + if (!_separate_shells && !_force_manifold) + { + collect_facets_quick (); + _shells.push_back (_surface.begin ()); + } + else + { + // Init shell index + _index = 0; + + // 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; + + if (_force_manifold) + { + // If manifold surface is wanted, small volumes (bubbles) are first detected in + // order to know which layer to ignore when collecting facets + detect_bubbles(); + } + + collect_facets (); + + if (_force_manifold) + { + // Even when taking into account facets, some nonmanifold features might remain + fix_nonmanifold_edges(); + fix_nonmanifold_vertices(); + } + } + + for (Triple_iterator it = _surface.begin(); it != _surface.end(); ++ it) + { + Triple t = *it; + *(output ++) = {{ (*it)[0], (*it)[1], (*it)[2] }}; + } +// std::copy (_surface.begin(), _surface.end(), output); + } + + +private: + + void deinit_shape() + { + if (_shape != NULL) + { + delete _shape; + _shape = NULL; + } + } + + void clear_surface() + { + _shells.clear(); + _surface.clear(); + _garbage.clear(); + deinit_shape(); + } + + void collect_facets () + { + // We check each of the facets: if it is not handled and either regular or singular, + // we start collecting the next surface from here. + for (Finite_facets_iterator fit = _shape->finite_facets_begin(); fit != _shape->finite_facets_end(); ++fit) + { + switch(_shape->classify (*fit)) + { + case Shape::REGULAR: + // Build a surface from the outer cell. + if (_shape->classify(fit->first) == Shape::EXTERIOR) + collect_shell (*fit); + else + collect_shell (_shape->mirror_facet (*fit)); + break; + case Shape::SINGULAR: + // Build a surface from both incident cells. + collect_shell (*fit); + collect_shell (_shape->mirror_facet (*fit)); + break; + default: + break; + } + } + } + + void collect_facets_quick () + { + // Collect all facets from the alpha-shape in an unordered fashion. + for (Finite_facets_iterator fit = _shape->finite_facets_begin(); fit != _shape->finite_facets_end(); ++fit) + { + switch (_shape->classify(*fit)) + { + case Shape::REGULAR: + // Collect the outer cell. + if (_shape->classify(fit->first) == Shape::EXTERIOR) + _surface.push_back (ordered_facet_indices (*fit)); + else + _surface.push_back (ordered_facet_indices (_shape->mirror_facet(*fit))); + break; + case Shape::SINGULAR: + // Collect both incident cells. + _surface.push_back (ordered_facet_indices (*fit)); + _surface.push_back (ordered_facet_indices (_shape->mirror_facet (*fit))); + break; + default: + break; + } + } + } + + 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 mark_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 mark_handled( Facet f ) { mark_handled( f.first, f.second ); } + + inline void mark_opposite_handled( Facet f ) + { + + Classification_type cl = _shape->classify (f); + + // If cell is singular, simply mark mirror facet as handled + if (cl == Shape::SINGULAR) + { + Facet mirror = _shape->mirror_facet (f); + mark_handled (mirror); + } + // If cell is regular, get corresponding bubble and mark + // facets of the other layer of the bubble as handled + else if (cl == Shape::REGULAR) + { + Facet fac = (_shape->classify (f.first) == Shape::EXTERIOR) + ? f + : _shape->mirror_facet (f); + + typename std::map::iterator + search = _map_f2b.find (fac); + + if (search == _map_f2b.end ()) + return; + + unsigned int layer = (_bubbles[search->second][0].find (fac) == _bubbles[search->second][0].end ()) + ? 0 : 1; + + typename std::set::iterator it = _bubbles[search->second][layer].begin (); + + // If bubble has already been handled, no need to do it again + if (is_handled (*it)) + return; + + for (;it != _bubbles[search->second][layer].end (); ++ it) + { + _garbage.push_back (ordered_facet_indices (*it)); + mark_handled (*it); + } + + } + + + } + + + inline Triple ordered_facet_indices( const Facet& f ) const + { + if( (f.second&1) == 0 ) + return make_array( 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 make_array( 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( const Facet& f) + { + collect_shell (f.first, f.second); + } + void collect_shell( Cell_handle c, unsigned int li) + { + // Collect one surface mesh from the alpha-shape in a fashion similar to ball-pivoting. + // Invariant: the facet is regular or singular. + + + // To stop stack overflows: use own stack. + std::stack stack; + stack.push( Facet(c, li) ); + + Facet f; + Cell_handle n, p; + int ni, pi; + Vertex_handle a; + Classification_type cl; + bool processed = false; + + while( !stack.empty() ) { + f = stack.top(); + stack.pop(); + + // Check if the cell was already handled. + // Note that this is an extra check that in many cases is not necessary. + if( is_handled(f) ) + continue; + + // The facet is part of the surface. + CGAL_assertion( !_shape->is_infinite(f) ); + mark_handled(f); + // Output the facet as a triple of indices. + _surface.push_back( ordered_facet_indices(f) ); + if( !processed ) { + if (_separate_shells || _shells.size () == 0) + { + _shells.push_back( --_surface.end() ); + _index ++; + } + processed = true; + } + + // Save in which shell the facet is stored + _map_f2s[f] = _index-1; + + // If the surface is required to be manifold, + // the opposite layer should be ignored + if (_force_manifold) + mark_opposite_handled (f); + + // Pivot over each of the facet's edges and continue the surface at the next regular or singular facet. + for( 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) ); + } + + } + + } + + void detect_bubbles () + { + std::set done; + + unsigned int nb_facets_removed = 0; + + unsigned int nb_skipped = 0; + for (Cell_iterator cit = _shape->cells_begin (); cit != _shape->cells_end (); ++ cit) + { + if (_shape->is_infinite (cit)) + continue; + if (_shape->classify (cit) != Shape::INTERIOR) + continue; + if (done.find (cit) != done.end ()) + continue; + + std::set borders; + std::vector cells; + std::stack todo; + todo.push (cit); + + // Get all cells of volume and all borders + while (!(todo.empty ())) + { + Cell_handle c = todo.top (); + todo.pop (); + + if (!(done.insert (c).second)) + continue; + + cells.push_back (c); + + for (unsigned int i = 0; i < 4; ++ i) + { + if (_shape->classify (c->neighbor (i)) == Shape::INTERIOR) + todo.push (c->neighbor (i)); + else + { + // Test if edge is proper border + for (unsigned int j = 0; j < 3; ++ j) + { + unsigned int i0 = (i + j + 1)%4; + unsigned int i1 = (i + (j+1)%3 + 1)%4; + CGAL_assertion (i0 != i && i1 != i); + Edge edge (c, i0, i1); + + if (_shape->classify (edge) != Shape::REGULAR) + continue; + + VEdge vedge = (c->vertex (i0) < c->vertex (i1)) + ? std::make_pair (c->vertex (i0), c->vertex (i1)) + : std::make_pair (c->vertex (i1), c->vertex (i0)); + + if (borders.find (vedge) != borders.end ()) + continue; + + Facet_circulator start = _shape->incident_facets (edge); + Facet_circulator circ = start; + unsigned int cnt = 0; + do + { + if (_shape->classify (*circ) == Shape::SINGULAR + || _shape->classify (*circ) == Shape::REGULAR) + ++ cnt; + ++ circ; + } + while (circ != start); + + // If edge is non-manifold, use as border + if (cnt > 2) + { + borders.insert (vedge); + continue; + } + + // Else, if facets in cell are regular and angle is + // under _border_angle limit, use as border + Facet f0 (c, i); + Facet f1 (c, (i + (j+2)%3 + 1)%4); + + if (_shape->classify (f0) != Shape::REGULAR + || _shape->classify (f1) != Shape::REGULAR) + continue; + + double angle = Geom_traits().compute_approximate_dihedral_angle_3_object()(vedge.first->point (), + vedge.second->point (), + c->vertex (i)->point (), + c->vertex ((i + (j+2)%3 + 1)%4)->point ()); + + if (-_border_angle < angle && angle < _border_angle) + { + borders.insert (vedge); + } + } + + } + } + } + + int layer = -1; + + // Try to generate bubble from the volume found + _bubbles.push_back (Bubble()); + std::set done; + for (unsigned int c = 0; c < cells.size (); ++ c) + { + for (unsigned int ii = 0; ii < 4; ++ ii) + { + Facet start = _shape->mirror_facet (Facet (cells[c], ii)); + + if (_shape->classify (start) != Shape::REGULAR) + continue; + + if (done.find (start) != done.end ()) + continue; + + ++ layer; + + std::stack stack; + stack.push (start); + + Facet f; + Cell_handle n, p; + int ni, pi; + Vertex_handle a; + Classification_type cl; + + // A bubble is well formed is the border contains one loop that + // separates two layers. + // If the number of layers is different than 2, the volume is completely ignored. + while( !stack.empty() ) + { + f = stack.top(); + stack.pop(); + + if (!(done.insert (f).second)) + continue; + + if (_shape->classify (f.first) == Shape::EXTERIOR) + { + if (layer < 2) + { + _bubbles.back ()[layer].insert (f); + _map_f2b[f] = _bubbles.size () - 1; + } + else + { + nb_facets_removed ++; + mark_handled (f); + _garbage.push_back (ordered_facet_indices (f)); + } + } + else + { + if (layer < 2) + { + _bubbles.back ()[layer].insert (_shape->mirror_facet (f)); + _map_f2b[_shape->mirror_facet (f)] = _bubbles.size () - 1; + } + else + { + nb_facets_removed ++; + mark_handled (_shape->mirror_facet (f)); + _garbage.push_back (ordered_facet_indices (_shape->mirror_facet (f))); + } + } + + + for( int i = 0; i < 4; ++i ) + { + // Skip the current facet. + if( i == f.second) + continue; + + n = f.first; + ni = i; + a = f.first->vertex( f.second ); + cl = _shape->classify( Facet(n, ni) ); + + int n0 = -1, n1 = -1; + bool n0found = false; + for (int j = 0; j < 4; ++ j) + { + if (j != ni && j != f.second) + { + if (n0found) + { + n1 = j; + break; + } + else + { + n0 = j; + n0found = true; + } + } + } + + VEdge vedge = (n->vertex (n0) < n->vertex (n1)) + ? std::make_pair (n->vertex (n0), n->vertex (n1)) + : std::make_pair (n->vertex (n1), n->vertex (n0)); + + // If the edge is a border, propagation stops in this direction. + if (borders.find (vedge) != borders.end ()) + continue; + + 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) ); + } + + stack.push (Facet (n, ni)); + + } + + } + } + + } + + // If number of layers is != 2, ignore volume and discard bubble + if (layer != 1) + { + nb_skipped ++; + for (unsigned int i = 0; i < 2; ++ i) + for (typename std::set::iterator fit = _bubbles.back()[i].begin (); + fit != _bubbles.back()[i].end (); ++ fit) + { + mark_handled (*fit); + _map_f2b.erase (*fit); + _garbage.push_back (ordered_facet_indices (*fit)); + nb_facets_removed ++; + } + _bubbles.pop_back (); + } + + } + } + + + void fix_nonmanifold_edges() + { + + typedef std::map, std::set > Edge_shell_map_triples; + typedef typename Edge_shell_map_triples::iterator Edge_shell_map_triples_iterator; + + unsigned int nb_facets_removed = 0; + + unsigned int nb_nm_edges = 0; + + // Store for each pair edge/shell the incident facets + Edge_shell_map_triples eshell_triples; + std::map map_t2f; + + for (typename Map_facet_to_shell::iterator fit = _map_f2s.begin (); + fit != _map_f2s.end (); ++ fit) + { + Facet f = fit->first; + Triple t = ordered_facet_indices (f); + map_t2f[t] = f; + + for (unsigned int k = 0; k < 3; ++ k) + { + Vertex_handle v0 = f.first->vertex ((f.second + k + 1)%4); + Vertex_handle v1 = f.first->vertex ((f.second + (k+1)%3 + 1)%4); + VEdge vedge = (v0 < v1) ? std::make_pair (v0, v1) : std::make_pair (v1, v0); + + std::pair + search = eshell_triples.insert (std::make_pair (std::make_pair (vedge, fit->second), + std::set())); + + search.first->second.insert (t); + } + } + + for (Edge_shell_map_triples_iterator eit = eshell_triples.begin (); + eit != eshell_triples.end (); ++ eit) + { + // If an edge has more than 2 incident facets for one shell, it is non-manifold + if (eit->second.size () < 3) + continue; + + ++ nb_nm_edges; + + Triple_iterator tit = _shells[eit->first.second]; + Triple_iterator end = (eit->first.second == _shells.size () - 1) + ? _surface.end () : _shells[eit->first.second + 1]; + + // Remove facets until the edge is manifold in this shell + while (tit != end && eit->second.size () > 2) + { + Triple_iterator current = tit ++; + + typename std::set::iterator search = eit->second.find (*current); + + if (search != eit->second.end ()) + { + if (current == _shells[eit->first.second]) + _shells[eit->first.second] = tit; + + _garbage.push_back (*current); + _map_f2s.erase (map_t2f[*current]); + _surface.erase (current); + + ++ nb_facets_removed; + eit->second.erase (search); + } + + } + + } + } + + template + struct operator_less + { + bool operator() (const T& a, const T& b) const + { + return &*a < &*b; + } + }; + + + void find_two_other_vertices(const Facet& f, Vertex_handle v, + Vertex_handle& v1, Vertex_handle& v2) + { + Vertex_handle vother = f.first->vertex (f.second); + bool v1found = false; + + for (unsigned int i = 0; i < 4; ++ i) + { + Vertex_handle vi = f.first->vertex (i); + if (vi != v && vi != vother) + { + if (v1found) + { + v2 = vi; + return; + } + else + { + v1 = vi; + v1found = true; + } + } + } + } + + void fix_nonmanifold_vertices() + { + + typedef ::CGAL::Union_find UF; + typedef typename UF::handle UF_handle; + + + typedef std::map, std::vector > Vertex_shell_map_facets; + typedef typename Vertex_shell_map_facets::iterator Vertex_shell_map_facet_iterator; + + // For faster facet removal, we sort the triples of each shell as a preprocessing + for (unsigned int i = 0; i < _shells.size (); ++ i) + { + Triple_iterator begin = _shells[i]; + Triple_iterator end = (i+1 == _shells.size ()) ? _surface.end () : _shells[i+1]; + + Tripleset tmp; + tmp.splice (tmp.end(), _surface, begin, end); + + tmp.sort(); + _shells[i] = tmp.begin (); + _surface.splice(end, tmp, tmp.begin(), tmp.end()); + } + + unsigned int nb_facets_removed = 0; + unsigned int nb_nm_vertices = 0; + // Removing facets to fix non-manifold vertices might make some other vertices + // become non-manifold, therefore we iterate until no facet needs to be removed. + do + { + nb_nm_vertices = 0; + nb_facets_removed = 0; + + // Store for each pair vertex/shell the incident facets + Vertex_shell_map_facets vshell_facets; + + for (typename Map_facet_to_shell::iterator fit = _map_f2s.begin (); + fit != _map_f2s.end (); ++ fit) + { + Facet f = fit->first; + + for (unsigned int k = 0; k < 3; ++ k) + { + Vertex_handle v = f.first->vertex ((f.second+k+1)%4); + + std::pair + search = vshell_facets.insert (std::make_pair (std::make_pair (v, fit->second), + std::vector())); + search.first->second.push_back (f); + + } + + } + + for (Vertex_shell_map_facet_iterator fit = vshell_facets.begin (); + fit != vshell_facets.end (); ++ fit) + { + if (fit->second.size () < 2) + continue; + + Vertex_handle vit = fit->first.first; + unsigned int shell = fit->first.second; + + UF uf; + std::map map_f2h; + + for (unsigned int i = 0; i < fit->second.size (); ++ i) + map_f2h.insert (std::make_pair (fit->second[i], uf.make_set (fit->second[i]))); + + std::map map_v2f; + + for (unsigned int i = 0; i < fit->second.size (); ++ i) + { + Vertex_handle v1, v2; + find_two_other_vertices (fit->second[i], vit, v1, v2); + std::pair::iterator, bool> + insertion1 = map_v2f.insert (std::make_pair (v1, fit->second[i])); + if (!(insertion1.second)) + uf.unify_sets (map_f2h[fit->second[i]], map_f2h[insertion1.first->second]); + std::pair::iterator, bool> + insertion2 = map_v2f.insert (std::make_pair (v2, fit->second[i])); + if (!(insertion2.second)) + uf.unify_sets (map_f2h[fit->second[i]], map_f2h[insertion2.first->second]); + } + + if (uf.number_of_sets () > 1) + { + ++ nb_nm_vertices; + + typedef std::map, operator_less > Map_uf_sets; + Map_uf_sets map_h2f; + for (unsigned int i = 0; i < fit->second.size (); ++ i) + { + UF_handle handle = uf.find (map_f2h[fit->second[i]]); + + std::pair + insertion = map_h2f.insert (std::make_pair (handle, std::vector())); + + insertion.first->second.push_back (fit->second[i]); + } + + typename Map_uf_sets::iterator largest = map_h2f.end (); + std::size_t nb_largest = 0; + for (typename Map_uf_sets::iterator ufit = map_h2f.begin (); ufit != map_h2f.end (); ++ ufit) + { + std::size_t size = ufit->second.size (); + if (size > nb_largest) + { + nb_largest = size; + largest = ufit; + } + } + + std::vector triples; + + for (typename Map_uf_sets::iterator ufit = map_h2f.begin (); ufit != map_h2f.end (); ++ ufit) + { + if (ufit == largest) + continue; + for (unsigned int i = 0; i < ufit->second.size (); ++ i) + { + _map_f2s.erase (ufit->second[i]); + triples.push_back (ordered_facet_indices (ufit->second[i])); + } + } + std::sort (triples.begin (), triples.end ()); + + Triple_iterator tit = _shells[shell]; + Triple_iterator end = (shell == _shells.size () - 1) + ? _surface.end () : _shells[shell + 1]; + + unsigned int tindex = 0; + + while (tit != end && tindex < triples.size ()) + { + Triple_iterator current = tit ++; + + if (*current == triples[tindex]) + { + if (current == _shells[shell]) + _shells[shell] = tit; + + _garbage.push_back (*current); + _surface.erase (current); + + ++ nb_facets_removed; + ++ tindex; + } + + } + } + + } + + } + while (nb_nm_vertices != 0); + + } + +}; + + +} // namespace Scale_space_reconstruction_3 + +} // namespace CGAL + +#endif // CGAL_SCALE_SPACE_RECONSTRUCTION_3_ALPHA_SHAPE_MESHER_H diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Weighted_PCA_smoother.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Weighted_PCA_smoother.h new file mode 100644 index 00000000000..eeab2b87832 --- /dev/null +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Weighted_PCA_smoother.h @@ -0,0 +1,280 @@ +#ifndef CGAL_SCALE_SPACE_RECONSTRUCTION_3_WEIGHTED_PCA_SMOOTHER_H +#define CGAL_SCALE_SPACE_RECONSTRUCTION_3_WEIGHTED_PCA_SMOOTHER_H + +#include +#include +#include +#include +#include +#include +#include + +#include + +#ifdef CGAL_LINKED_WITH_TBB +#include "tbb/blocked_range.h" +#include "tbb/parallel_for.h" +#endif // CGAL_LINKED_WITH_TBB + +namespace CGAL +{ + +namespace Scale_space_reconstruction_3 +{ + +template , +#ifdef CGAL_LINKED_WITH_TBB + typename ConcurrencyTag = CGAL::Parallel_tag> +#else + typename ConcurrencyTag = CGAL::Sequential_tag> +#endif +class Weighted_PCA_smoother +{ +public: + typedef typename Geom_traits::FT FT; ///< defines the point type. + typedef typename Geom_traits::Point_3 Point; ///< defines the point typ.e + typedef typename Geom_traits::Vector_3 Vector; ///< defines the vector type. +private: + + + typedef boost::tuple Point_and_size_t; + typedef std::vector CountVec; + typedef std::vector Pointset; + + typedef Search_traits_3 Traits_base; + typedef CGAL::Search_traits_adapter, + Traits_base> Search_traits; + typedef Orthogonal_k_neighbor_search Static_search; + typedef Orthogonal_incremental_neighbor_search Dynamic_search; + typedef typename Dynamic_search::Tree Search_tree; + typedef Fuzzy_sphere< Search_traits > Sphere; + + Random _generator; + unsigned int _mean_neighbors; + unsigned int _samples; + FT _squared_radius; + Search_tree _tree; // To quickly search for nearest neighbors. + Pointset _points; + + +public: + + Weighted_PCA_smoother (unsigned int neighbors = 12, unsigned int samples = 300) + : _generator(0), _mean_neighbors (neighbors), _samples (samples),_squared_radius (-1.) { } + + Weighted_PCA_smoother (FT squared_radius) + : _generator(0), _mean_neighbors (0), _samples (0),_squared_radius (squared_radius) { } + + template + void operator() (InputIterator begin, InputIterator end) + { + _tree.clear(); + _points.clear(); + + std::size_t i = 0; + std::size_t size = std::size_t(end - begin); + _tree.reserve(size); + _points.reserve(size); + while(begin!=end) + { + _points.push_back (*begin); + _tree.insert( boost::make_tuple(*(begin ++),i ++) ); + } + + _tree.build(); + + if (_squared_radius == -1) + estimate_neighborhood_squared_radius(); + + // Collect the number of neighbors of each point. + // This can be done concurrently. + CountVec neighbors (_tree.size(), 0); + try_parallel (ComputeNN (_points, _tree, _squared_radius, neighbors), 0, _tree.size()); + + // Compute the transformed point locations. + // This can be done concurrently. + try_parallel (AdvanceSS (_tree, neighbors, _points), 0, _tree.size()); + + } + +private: + + void estimate_neighborhood_squared_radius () + { + typename Geom_traits::Compute_squared_distance_3 squared_distance + = Geom_traits().compute_squared_distance_3_object(); + + unsigned int handled = 0; + unsigned int checked = 0; + FT radius = 0.; + + for (typename Search_tree::const_iterator it = _tree.begin(); it != _tree.end(); ++it ) + { + unsigned int left = (unsigned int)(_tree.size() - handled); + if (_samples >= left || _generator.get_double() < double(_samples - checked) / left) + { + // The neighborhood should contain the point itself as well. + Static_search search (_tree, boost::get<0>(*it), _mean_neighbors + 1); + radius += std::sqrt (to_double (squared_distance( boost::get<0>(*it), boost::get<0>((search.end()-1)->first)))); + ++checked; + } + ++handled; + } + radius /= double(checked); + + _squared_radius = radius * radius; + } + + + template + void try_parallel (const F& func, std::size_t begin, std::size_t end) + { +#ifndef CGAL_LINKED_WITH_TBB + CGAL_static_assertion_msg (!(boost::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else + if (boost::is_convertible::value) + tbb::parallel_for(tbb::blocked_range(begin, end), func); + else +#endif + for (std::size_t i = begin; i < end; ++i) + func(i); + } + + struct Inc + { + unsigned int * i; + + Inc(unsigned int& i) + : i(&i) + {} + + template + void operator()(const T&) const + { + ++(*i); + } + + }; + + // Compute the number of neighbors of a point that lie within a fixed radius. + class ComputeNN + { + private: + typename Geom_traits::Compare_squared_distance_3 compare; + + const Pointset& _pts; + const Search_tree& _tree; + const FT _sq_rd; + CountVec& _nn; + + public: + ComputeNN(const Pointset& points, const Search_tree& tree, + const FT& sq_radius, CountVec& nn) + : _pts(points), _tree(tree), _sq_rd(sqrt(sq_radius)), _nn(nn) {} + +#ifdef CGAL_LINKED_WITH_TBB + void operator()( const tbb::blocked_range< std::size_t >& range ) const { + for( std::size_t i = range.begin(); i != range.end(); ++i ) + (*this)( i ); + } +#endif // CGAL_LINKED_WITH_TBB + void operator()( const std::size_t& i ) const { + + Sphere sp(_pts[i], _sq_rd); + + Inc inc(_nn[i]); + _tree.search(boost::make_function_output_iterator(inc),sp); + } + }; // class ComputeNN + + +// Advance a point to a coarser scale. + class AdvanceSS + { + private: + const Search_tree& _tree; + const CountVec& _nn; + Pointset& _pts; + + public: + AdvanceSS(const Search_tree& tree, const CountVec& nn, Pointset& points) + : _tree(tree), _nn(nn),_pts(points) {} + +#ifdef CGAL_LINKED_WITH_TBB + void operator()( const tbb::blocked_range< std::size_t >& range ) const { + for( std::size_t i = range.begin(); i != range.end(); ++i ) + (*this)( i ); + } +#endif // CGAL_LINKED_WITH_TBB + void operator()( const std::size_t& i ) const { + // If the neighborhood is too small, the vertex is not moved. + if( _nn[i] < 4 ) + return; + + Static_search search(_tree, _pts[i], _nn[i]); + + Point barycenter (0., 0., 0.); + FT weight_sum = 0.; + unsigned int column = 0; + // Compute total weight + for( typename Static_search::iterator nit = search.begin(); + nit != search.end() && column < _nn[i]; + ++nit, ++column ) + weight_sum += (1.0 / _nn[boost::get<1>(nit->first)]); + + column = 0; + // Compute barycenter + for( typename Static_search::iterator nit = search.begin(); + nit != search.end() && column < _nn[i]; + ++nit, ++column ) + { + Vector v (CGAL::ORIGIN, boost::get<0>(nit->first)); + barycenter = barycenter + ((1.0 / _nn[boost::get<1>(nit->first)]) / weight_sum) * v; + } + + CGAL::cpp11::array covariance = {{ 0., 0., 0., 0., 0., 0. }}; + column = 0; + // Compute covariance matrix of Weighted PCA + for( typename Static_search::iterator nit = search.begin(); + nit != search.end() && column < _nn[i]; + ++nit, ++column ) + { + Vector v (barycenter, boost::get<0>(nit->first)); + FT w = (1.0 / _nn[boost::get<1>(nit->first)]); + v = w*v; + covariance[0] += w * v.x () * v.x (); + covariance[1] += w * v.x () * v.y (); + covariance[2] += w * v.x () * v.z (); + covariance[3] += w * v.y () * v.y (); + covariance[4] += w * v.y () * v.z (); + covariance[5] += w * v.z () * v.z (); + } + + // Compute the weighted least-squares planar approximation of the point set. + CGAL::cpp11::array eigenvectors = {{ 0., 0., 0., + 0., 0., 0., + 0., 0., 0. }}; + CGAL::cpp11::array eigenvalues = {{ 0., 0., 0. }}; + DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix + (covariance, eigenvalues, eigenvectors); + + // The vertex is moved by projecting it onto the plane + // through the barycenter and orthogonal to the Eigen vector with smallest Eigen value. + Vector norm (eigenvectors[0], eigenvectors[1], eigenvectors[2]); + Vector b2p (barycenter, _pts[i]); + + _pts[i] = barycenter + b2p - ((norm * b2p) * norm); + } + }; // class AdvanceSS + +}; + + +} // namespace Scale_space_reconstruction_3 + +} // namespace CGAL + +#endif // CGAL_SCALE_SPACE_RECONSTRUCTION_3_WEIGHTED_PCA_SMOOTHER_H diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h index 32d63d7ab84..d15a3ed7a99 100644 --- a/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h @@ -1,968 +1,210 @@ -//Copyright (C) 2014 INRIA - Sophia Antipolis -// -//This program is free software: you can redistribute it and/or modify -//it under the terms of the GNU General Public License as published by -//the Free Software Foundation, either version 3 of the License, or -//(at your option) any later version. -// -//This program is distributed in the hope that it will be useful, -//but WITHOUT ANY WARRANTY; without even the implied warranty of -//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -//GNU General Public License for more details. -// -//You should have received a copy of the GNU General Public License -//along with this program. If not, see . -// -// Author(s): Thijs van Lankveld - #ifndef CGAL_SCALE_SPACE_SURFACE_RECONSTRUCTION_3_H #define CGAL_SCALE_SPACE_SURFACE_RECONSTRUCTION_3_H #include +#include +#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +namespace CGAL +{ -#include -#include +template +class Scale_space_surface_reconstruction_3 +{ +public: -namespace CGAL { + typedef typename Geom_traits::FT FT; ///< defines the field number type. + typedef typename Geom_traits::Point_3 Point; ///< defines the point type. + typedef cpp11::array Facet; ///< defines a facet of the surface (triple of point indices). + +#ifdef DOXYGEN_RUNNING + typedef unspecified_type Point_iterator; ///< defines an iterator over the points. + typedef const unspecified_type Point_const_iterator; ///< defines a constant iterator over the points. +#else + typedef typename std::vector Point_vector; + typedef typename Point_vector::iterator Point_iterator; + typedef typename Point_vector::const_iterator Point_const_iterator; +#endif + +#ifdef DOXYGEN_RUNNING + typedef unspecified_type Facet_iterator; ///< defines an iterator over the points. + typedef const unspecified_type Facet_const_iterator; ///< defines a constant iterator over the points. +#else + typedef typename std::vector Facet_vector; + typedef typename Facet_vector::iterator Facet_iterator; + typedef typename Facet_vector::const_iterator Facet_const_iterator; +#endif + + // Default algorithms used (same as in old API) + typedef Scale_space_reconstruction_3::Weighted_PCA_smoother Weighted_PCA_smoother; + typedef Scale_space_reconstruction_3::Alpha_shape_mesher Alpha_shape_mesher; +private: -/// computes a triangulated surface mesh interpolating a point set. -/** \ingroup PkgScaleSpaceReconstruction3Classes - * This class stores several of the (intermediate) results. This makes it - * easier and more efficient to adjust the parameter settings based on - * preliminary results, or to further increase the scale to improve the - * results. The class stores the point set at the current scale and the - * reconstructed surface, possibly with iterators over the shells. - * - * The class also stores the parameters for estimating the optimal - * neighborhood radius and either the lastest estimate or the manually set - * radius. This way, the radius can be estimated (again) whenever necessary. - * Also note that both increasing the scale and reconstructing the - * surface use this radius. By changing or re-estimating the radius between - * these operations, they can use separate parameter settings. - * - * The surface can be constructed either for a fixed neighborhood radius, or - * for a dynamic radius. When constructing the surface for exactly one - * neighborhood radius, it is faster to set `FS` to `Tag_true`. If - * the correct neighborhood radius should be changed or estimated multiple - * times, it is faster to set `FS` to `Tag_false`. - * - * It is undefined whether a surface with fixed radius may have its radius - * changed, but if so, this will likely require more computation time than - * changing the radius of a dynamic surface. In either case, it is possible to - * change the point set while maintaining the same radius. - * - * The surface can be stored either as an unordered collection of triangles, - * or as a collection ordered by shells. A shell is a maximally connected - * component of the surface where connected facets are locally oriented - * towards the same side of the surface. - * - * \tparam Gt is the geometric traits class. It must be a model of - * `DelaunayTriangulationTraits_3`. It must have a `RealEmbeddable` field - * number type. Generally, `Exact_predicates_inexact_constructions_kernel` is - * preferred. - * \tparam FS determines whether the surface is expected to be constructed - * for a fixed neighborhood radius. It must be a `Boolean_tag` type. The default value is - * `Tag_true`. Note that the value of this parameter does not change the result but - * only has an impact on the run-time. - * \tparam wA must be a model of `DiagonalizeTraits` and determines - * how to diagonalize a weighted covariance matrix to approximate a - * weighted point set. It can be omitted: if Eigen 3 (or greater) is - * available and `CGAL_EIGEN3_ENABLED` is defined then an overload - * using `Eigen_diagonalize_traits` is provided. Otherwise, the - * internal implementation `Diagonalize_traits` is used. - * \tparam Ct indicates whether to use concurrent processing. It must be - * either `Sequential_tag` or `Parallel_tag` (the default value). - */ -template < class Gt, class FS = Tag_true, - class wA = CGAL::Default_diagonalize_traits, class Ct = Parallel_tag > -class Scale_space_surface_reconstruction_3 { + Point_vector m_points; + Facet_vector m_facets; public: - typedef typename Gt::Point_3 Point; ///< defines the point type. - typedef boost::tuple Point_and_size_t; -private: - // Searching for neighbors. - typedef Search_traits_3< Gt > Traits_base; -typedef CGAL::Search_traits_adapter, - Traits_base> Search_traits; - typedef Orthogonal_k_neighbor_search< Search_traits > - Static_search; - typedef Orthogonal_incremental_neighbor_search< Search_traits > - Dynamic_search; - typedef typename Dynamic_search::Tree Search_tree; - typedef Fuzzy_sphere< Search_traits > Sphere; - typedef CGAL::Random Random; + Scale_space_surface_reconstruction_3 () { } - // Constructing the surface. - typedef CGAL::Shape_construction_3< Gt, FS > Shape_construction_3; - - typedef typename Shape_construction_3::Shape Shape; - typedef typename Shape_construction_3::Triangulation Triangulation; - - typedef typename Shape::Vertex_handle Vertex_handle; - typedef typename Shape::Cell_handle Cell_handle; - typedef typename Shape::Facet Facet; - typedef typename Shape::Edge Edge; - typedef std::pair VEdge; + template + Scale_space_surface_reconstruction_3 (InputIterator begin, InputIterator end) + { + insert (begin, end); + } + /// inserts a point into the scale-space at the current scale. + /** \param p is the point to insert. + * + * \note Inserting the point does not automatically construct or + * update the surface. + * + * \note In order to construct the surface, call + * `#reconstruct_surface()`. + * + * \sa `insert(InputIterator begin, InputIterator end)`. + */ + void insert (const Point& p) + { + m_points.push_back (p); + } - - typedef typename Shape::Vertex_iterator Vertex_iterator; - typedef typename Shape::Cell_iterator Cell_iterator; - typedef typename Shape::Facet_iterator Facet_iterator; - typedef typename Shape::Edge_iterator Edge_iterator; + /// inserts a collection of points into the scale-space at the current scale. + /** \tparam InputIterator is an iterator over the point collection. + * The value type of the iterator must be a `Point`. + * + * \param begin is an iterator to the first point of the collection. + * \param end is a past-the-end iterator for the point collection. + * + * \note Inserting the points does not automatically construct or + * update the surface. + * + * \note In order to construct the surface, call + * `reconstruct_surface()` after inserting the points. + * + * \sa `insert(const Point& p)`. + */ + template + void insert (InputIterator begin, InputIterator end) + { + std::copy (begin, end, std::back_inserter (m_points)); + } - typedef typename Shape::Finite_cells_iterator Finite_cells_iterator; - typedef typename Shape::Finite_facets_iterator Finite_facets_iterator; - typedef typename Shape::Finite_edges_iterator Finite_edges_iterator; - typedef typename Shape::Finite_vertices_iterator Finite_vertices_iterator; + /// increases the scale by a number of iterations. + /** Each iteration the scale is increased, the points set at a higher scale + * is computed. At a higher scale, the points set is smoother. + * + * \param iterations is the number of iterations to perform. If + * `iterations` is 0, nothing happens. + * + * \note This method processes the point set at the current scale. The + * points can be set with [insert(begin, end)](\ref insert). + * + * \note If the surface was already constructed, increasing the scale + * will not automatically adjust the surface. + * + * \sa `reconstruct_surface()`. + */ + template + void increase_scale (std::size_t iterations = 1, const Smoother& smoother = Smoother()) + { + for (std::size_t i = 0; i < iterations; ++ i) + const_cast(smoother) (m_points.begin(), m_points.end()); + } - typedef typename Shape::Facet_circulator Facet_circulator; + void increase_scale (std::size_t iterations = 1) + { + increase_scale (iterations, Weighted_PCA_smoother()); + } + + /// constructs a triangle mesh from the point set at a fixed scale. + /** The order of the points at the current scale is the same as the order + * at the original scale, meaning that the surface can interpolate the + * point set at the original scale by applying the indices of the surface + * triangles to the original point set. + * + * After construction, the facets of the surface can be iterated using + * `facets_begin()` and `facets_end()`. + * + * \note This method processes the point set at the current scale. The + * points can be set with [insert(begin, end)](\ref insert). + * + * \sa `increase_scale()`. + */ + template + void reconstruct_surface (const Mesher& mesher = Mesher()) + { + m_facets.clear(); + const_cast(mesher) (m_points.begin(), m_points.end(), std::back_inserter (m_facets)); + } + + void reconstruct_surface () + { + reconstruct_surface (Alpha_shape_mesher()); + } + + /// gives the number of points of the surface. + std::size_t number_of_points() const { return m_points.size(); } - typedef typename Shape::All_cells_iterator All_cells_iterator; + /// gives an iterator to the first point at the current scale. + /** \warning Changes to the scale-space do not cause an automatic update to + * the surface. + */ + Point_iterator points_begin() { return m_points.begin(); } - typedef typename Shape::Classification_type Classification_type; - - typedef std::map Map_facet_to_shell; - typedef typename CGAL::cpp11::array, 2 > Bubble; + /// gives a past-the-end iterator of the points at the current scale. + /** \warning Changes to the scale-space do not cause an automatic update to + * the surface. + */ + Point_iterator points_end() { return m_points.end(); } -public: -/// \name Types -/// \{ - typedef typename Gt::FT FT; ///< defines the field number type. - typedef typename Gt::Vector_3 Vector; ///< defines the vector type. - typedef typename Gt::Plane_3 Plane; ///< defines the plane type. - typedef typename Gt::Triangle_3 Triangle; ///< defines the triangle type. - -#ifdef DOXYGEN_RUNNING - typedef unspecified_type Point_iterator; ///< defines an iterator over the points. - typedef const unspecified_type Point_const_iterator; ///< defines a constant iterator over the points. -#else // DOXYGEN_RUNNING - - // typedef typename Search_tree::iterator Point_iterator; - // typedef typename Search_tree::const_iterator Point_const_iterator; - typedef typename std::vector::iterator Point_iterator; - typedef typename std::vector::const_iterator Point_const_iterator; -#endif // DOXYGEN_RUNNING - - typedef CGAL::cpp11::array< unsigned int, 3 > Triple; ///< defines a triple of point indices indicating a triangle of the surface. -private: - typedef std::list< Triple > Tripleset; ///< defines a collection of triples. - // Note that this is a list for two reasons: iterator validity for the shell iterators, and memory requirements for the expected huge collections. - -public: -#ifdef DOXYGEN_RUNNING - typedef unspecified_type Triple_iterator; ///< defines an iterator over the triples. - typedef const unspecified_type Triple_const_iterator; ///< defines a constant iterator over the triples. -#else // DOXYGEN_RUNNING - typedef Tripleset::iterator Triple_iterator; - typedef Tripleset::const_iterator Triple_const_iterator; -#endif // DOXYGEN_RUNNING - -/// \} - -private: - typedef std::vector< Triple_iterator > TripleIterSet; - -private: - class Finite_point_iterator; - class In_surface_tester; - typedef Filter_iterator< Facet_iterator, In_surface_tester > - Surface_facets_iterator; - - // Parallel processing functors. - class ComputeNN; - class AdvanceSS; - -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 neighborhood radius. - - FT _squared_radius; // The squared 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. - // We want to be able to set the alpha after constructing the scale-space - // reconstructer object. - Shape* _shape; - - // The surface. If the surface is collected per shell, the triples of the - // same shell are stored consecutively. - Tripleset _surface; - - // The shells can be accessed through iterators to the surface. - TripleIterSet _shells; - - // If the surface is forced to be manifold, removed facets are stored - Tripleset _garbage; - - // Map TDS facets to shells - Map_facet_to_shell _map_f2s; - // std::map _map_f2s; - unsigned int _index; - - std::vector _bubbles; - std::map _map_f2b; + /// gives an iterator to the first point at the current scale. + Point_const_iterator points_begin() const { return m_points.begin(); } - typedef std::vector Pointset; - Pointset _points; + /// gives a past-the-end iterator of the points at the current scale. + Point_const_iterator points_end() const { return m_points.end(); } -public: -/// \name Constructors -/// \{ - /// constructs a surface reconstructor that will automatically estimate the neighborhood radius. - /** \param neighbors is the number of neighbors a point's neighborhood should - * contain on average. - * \param samples is the number of points sampled to estimate the - * neighborhood radius. - */ - Scale_space_surface_reconstruction_3( unsigned int neighbors, unsigned int samples); + /// gives the number of facets of the surface. + std::size_t number_of_facets() const { return m_facets.size(); } - /// constructs a surface reconstructor with a given neighborhood radius. - /** \param sq_radius is the squared radius of the neighborhood. - * - * \pre `sq_radius` is not negative. - */ - Scale_space_surface_reconstruction_3( FT sq_radius ); + /// gives an iterator to the first triple in the surface. + /** \warning Changes to the surface may change its topology. + */ + Facet_iterator facets_begin() { return m_facets.begin(); } -/// \} - ~Scale_space_surface_reconstruction_3() { deinit_shape(); } + /// gives a past-the-end iterator of the triples in the surface. + /** \warning Changes to the surface may change its topology. + */ + Facet_iterator facets_end() { return m_facets.end(); } -private: - void deinit_shape() { if( _shape != 0 ) { delete _shape; _shape = 0; } } + /// gives an iterator to the first triple in the surface. + Facet_const_iterator facets_begin() const { return m_facets.begin(); } - void clear_tree() { _tree.clear(); } - void clear_surface() { _shells.clear(); _surface.clear(); _garbage.clear(); deinit_shape(); } - - // SURFACE COLLECTION - // Once a facet is added to the surface, it is marked as handled. - bool is_handled( Cell_handle c, unsigned int li ) const; - inline bool is_handled( const Facet& f ) const { return is_handled( f.first, f.second ); } - void mark_handled( Cell_handle c, unsigned int li ); - void mark_opposite_handled( Facet f ); - inline void mark_handled( Facet f ) { mark_handled( f.first, f.second ); } - - // Get the indices of the points of the facet ordered to point - // towardds the outside of the shape. - Triple ordered_facet_indices( const Facet& f ) const; - - // Collect the triangles of one shell of the surface. - void collect_shell( Cell_handle c, unsigned int li, bool separate_shells, bool force_manifold); - - // Collect the triangles of one shell of the surface. - void collect_shell( const Facet& f, bool separate_shells, bool force_manifold) { - collect_shell( f.first, f.second, separate_shells, force_manifold ); - } - - // Collect the triangles of the complete surface. - void collect_facets( bool separate_shells, bool force_manifold ); - void collect_facets_quick( ); + /// gives a past-the-end iterator of the triples in the surface. + Facet_const_iterator facets_end() const { return m_facets.end(); } +}; -private: - // Get the shape of the scale space. - /* If the shape does not exist, it is constructed first. - * \return the shape of the scale space. - */ - const Shape& shape() const; - -public: -/// \name Point Set Manipulation -/// \{ - /// inserts a collection of points into the scale-space at the current scale. - /** \tparam InputIterator is an iterator over the point collection. - * The value type of the iterator must be a `Point`. - * - * \param begin is an iterator to the first point of the collection. - * \param end is a past-the-end iterator for the point collection. - * - * \note Inserting the points does not automatically construct or - * update the surface. - * - * \note In order to construct the surface, call - * `reconstruct_surface()` after inserting the points. - * - * \warning Inserting new points may invalidate the neighborhood radius if - * it was previously estimated. - * - * \sa `insert(const Point& p)`. - */ - template < class InputIterator > -#ifdef DOXYGEN_RUNNING - void insert( InputIterator begin, InputIterator end ) { -#else // DOXYGEN_RUNNING - void insert( InputIterator begin, InputIterator end, - typename boost::enable_if< CGAL::is_iterator >::type* = NULL ) { -#endif // DOXYGEN_RUNNING - _tree.reserve (std::distance(begin,end) + _points.size()); - - _points.insert(_points.end(), begin, end); - - while(begin!=end) - _tree.insert( boost::make_tuple(*(begin ++),0) ); - } - - /// inserts a point into the scale-space at the current scale. - /** \param p is the point to insert. - * - * \note Inserting the point does not automatically construct or - * update the surface. - * - * \note In order to construct the surface, call - * `#reconstruct_surface()`. - * - * \warning Inserting a new point may invalidate the neighborhood radius - * if it was previously estimated. - * - * \sa `insert(InputIterator begin, InputIterator end)`. - */ - void insert( const Point& p ) { - _tree.insert( boost::make_tuple(p,0) ); - _points.push_back(p); - } - - /// clears the stored scale-space surface reconstruction data. - /** This includes discarding the surface, the scale-space and all its - * points, and any estimation of the neighborhood radius. - * - * Methods called after this point may have to re-estimate the - * neighborhood radius. This method does not discard the parameters for - * estimating this radius (the mean number of neighbors and the sample - * size). - */ - void clear() { - clear_tree(); - clear_surface(); - _map_f2s.clear (); - _bubbles.clear (); - _map_f2b.clear (); - - _squared_radius = -1; - } - -/// \} - -private: - // checks whether the shape has been constructed. - /* The shape contains the structure of the point cloud. - * - * Until the shape is constructed, the surface is undefined. - * \return true if the shape exists and false otherwise. - * \sa construct_shape(InputIterator begin, InputIterator end). - */ - bool has_shape() const { return _shape != 0; } - -public: -/// \name Neighborhood Size Estimation -/// \{ - /// estimates the neighborhood radius. - /** This method is equivalent to running - * [estimate_neighborhood_squared_radius( mean_number_of_neighbors(), neighborhood_sample_size() )](\ref estimate_neighborhood_squared_radius). - * - * This method can be called by the scale-space and surface construction - * methods if the neighborhood radius is not set when they are called. - * - * \return the estimated neighborhood radius. - * - * \note This method processes the point set at the current scale. The - * points can be set with [insert(begin, end)](\ref insert). - * - * \warning If the surface was already constructed, estimating the - * neighborhood radius will automatically adjust the surface. - * - * \sa `set_mean_number_of_neighbors(unsigned int neighbors)`. - * \sa `set_neighborhood_sample_size(unsigned int samples)`. - * \sa `estimate_neighborhood_squared_radius(unsigned int neighbors, unsigned int samples)`. - * \sa `increase_scale(unsigned int iterations)`. - * \sa `reconstruct_surface()`. - */ - inline FT estimate_neighborhood_squared_radius() { - return estimate_neighborhood_squared_radius( mean_number_of_neighbors(), neighborhood_sample_size() ); - } - - /// estimates the neighborhood radius based on a number of sample points. - /** The neighborhood radius 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, not counting the point itself. - * - * The neighborhood radius is set to the mean of these radii, taken over a - * number of point samples. - * - * \param neighbors is the number of neighbors a point's neighborhood - * should contain on average, not counting the point itself. - * \param samples is the number of points sampled to estimate the - * neighborhood radius. - * \return the estimated neighborhood radius. - * - * \note This method processes the point set at the current scale. The - * points can be set with [insert(begin, end)](\ref insert). - * - * \warning If the surface was already constructed, estimating the - * neighborhood radius will automatically adjust the surface. - * - * \sa `estimate_neighborhood_squared_radius()`. - */ - FT estimate_neighborhood_squared_radius( unsigned int neighbors, unsigned int samples ); - - - /// sets the squared radius of the neighborhood. - /** The neighborhood radius is used by - * `#[increase_scale()` to - * compute the point set at the desired scale and by - * `reconstruct_surface()` to - * construct a surface from the point set at the current scale. - * - * \param sq_radius is the squared radius of the neighborhood. - * - * \note If the neighborhood squared radius is negative when the point set - * is smoothed or when the surface is computed, the neighborhood radius - * will be computed automatically. - * - * \warning If the surface was already constructed, changing the - * neighborhood radius will automatically adjust the surface. - * - * \sa `neighborhood_squared_radius()`. - * \sa `has_neighborhood_squared_radius()`. - * \sa `increase_scale(unsigned int iterations)`. - * \sa `reconstruct_surface()`. - */ - void set_neighborhood_squared_radius( const FT& sq_radius ) { - _squared_radius = sq_radius; - if( has_neighborhood_squared_radius() && has_shape() ) - Shape_construction_3().change_scale( _shape, _squared_radius ); - } - - /// gives the squared radius of the neighborhood. - /** The neighborhood radius is used by - * `#increase_scale()` to - * compute the point set at the desired scale and by - * `#reconstruct_surface()` to - * construct a surface from the point set at the current scale. - * - * \return the squared radius of the neighborhood, or -1 if the - * neighborhood radius has not yet been set. - * - * \sa `increase_scale(unsigned int iterations)`. - * \sa `reconstruct_surface()`. - */ - FT neighborhood_squared_radius() const { return _squared_radius; } - - /// checks whether the neighborhood radius has been set. - /** The radius can be set manually, or estimated automatically. - * - * \return `true` iff the radius has been either set manually or estimated. - * - * \sa `set_neighborhood_squared_radius()`. - * \sa `estimate_neighborhood_squared_radius()`. - */ - bool has_neighborhood_squared_radius() const { - return sign( _squared_radius ) == POSITIVE; - } - - /// \cond internal_doc - /// estimates the neighborhood radius of a collection of points. - /** This method is equivalent to running - * `clear()` followed by - * [insert(begin, end)](\ref insert) and - * finally [estimate_neighborhood_squared_radius( mean_number_of_neighbors(), neighborhood_sample_size() )](\ref estimate_neighborhood_squared_radius). - * - * This method can be called by the scale-space and surface construction - * methods if the neighborhood radius is not set when they are called. - * - * \tparam InputIterator is an iterator over the point collection. - * The value type of the iterator must be a `Point`. - * \param begin is an iterator to the first point of the collection. - * \param end is a past-the-end iterator for the point collection. - * \return the estimated neighborhood radius. - * - * \sa `set_mean_number_of_neighbors(unsigned int neighbors)`. - * \sa `set_neighborhood_sample_size(unsigned int samples)`. - * \sa `estimate_neighborhood_squared_radius()`. - * \sa `insert(InputIterator begin, InputIterator end)`. - */ - template < class InputIterator > -#ifdef DOXYGEN_RUNNING - FT estimate_neighborhood_squared_radius( InputIterator begin, InputIterator end ) { -#else // DOXYGEN_RUNNING - FT estimate_neighborhood_squared_radius( InputIterator begin, InputIterator end, - typename boost::enable_if< CGAL::is_iterator >::type* = NULL ) { -#endif // DOXYGEN_RUNNING - return estimate_neighborhood_squared_radius( begin, end, mean_number_of_neighbors(), neighborhood_sample_size() ); - } - /// \endcond - - /// \cond internal_doc - /// estimates the neighborhood radius of a collection of points based on a number of sample points. - /** The neighborhood radius 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, not counting the point itself. - * - * The neighborhood radius is set to the mean of these radii, taken over a - * number of point samples. - * - * This method is equivalent to running - * `clear()` followed by - * [insert(begin, end)](\ref insert) and finally - * [estimate_neighborhood_squared_radius(neighbors, samples)](\ref estimate_neighborhood_squared_radius). - * - * \tparam InputIterator is an iterator over the point collection. - * The value type of the iterator must be a `Point`. - * - * \param begin is an iterator to the first point of the collection. - * \param end is a past-the-end iterator for the point collection. - * \param neighbors is the number of neighbors a point's neighborhood - * should contain on average, not counting the point itself. - * \param samples is the number of points sampled to estimate the - * neighborhood radius. - * \return the estimated neighborhood radius. - * - * \sa `estimate_neighborhood_squared_radius(unsigned int neighbors, unsigned int samples)`. - */ - template < class InputIterator > -#ifdef DOXYGEN_RUNNING - FT estimate_neighborhood_squared_radius( InputIterator begin, InputIterator end, unsigned int neighbors, unsigned int samples ); -#else // DOXYGEN_RUNNING - FT estimate_neighborhood_squared_radius( InputIterator begin, InputIterator end, unsigned int neighbors, unsigned int samples, - typename boost::enable_if< CGAL::is_iterator >::type* = NULL); -#endif // DOXYGEN_RUNNING - /// \endcond -/// \} - -/// \name Neighborhood Size Estimation Parameters -/// \{ - /// gives the mean number of neighbors an estimated neighborhood should contain. - /** This number is only used if the neighborhood radius has not been set - * manually. - * - * When the neighborhood radius is estimated, it should on average contain - * this many neighbors, not counting the neighborhood center. - * - * \return the number of neighbors a neighborhood ball centered on a point - * should contain on average when the radius is estimated, not counting - * the point itself. - * - * \sa `set_mean_number_of_neighbors(unsigned int neighbors)`. - * \sa `has_neighborhood_squared_radius()`. - * \sa `neighborhood_sample_size()`. - * \sa `estimate_neighborhood_squared_radius()`. - */ - unsigned int mean_number_of_neighbors() const { return _mean_neighbors; } - - /// gives the number of sample points the neighborhood estimation uses. - /** This number is only used if the neighborhood radius has not been set - * manually. - * - * If the number of samples is larger than the point cloud, every point is - * used and the optimal neighborhood radius is computed exactly instead of - * estimated. - * - * \return the number of points sampled for neighborhood estimation. - * - * \sa `set_neighborhood_sample_size(unsigned int samples)`. - * \sa `has_neighborhood_squared_radius()`. - * \sa `mean_number_of_neighbors()`. - * \sa `estimate_neighborhood_squared_radius()`. - */ - unsigned int neighborhood_sample_size() const { return _samples; } - - /// sets the mean number of neighbors an estimated neighborhood should contain. - /** This number is only used if the neighborhood radius has not been set - * manually. - * - * When the neighborhood radius is estimated, it should on average contain - * this many neighbors, not counting the neighborhood center. - * - * \param neighbors is the number of neighbors a neighborhood ball centered on a point - * should contain on average when the radius is estimated, not counting - * the point itself. - * - * \note This does not start the estimation process. - * - * \sa `mean_number_of_neighbors()`. - * \sa `has_neighborhood_squared_radius()`. - * \sa `set_neighborhood_sample_size(unsigned int samples)`. - */ - void set_mean_number_of_neighbors( unsigned int neighbors ) { _mean_neighbors = neighbors; } - - /// sets the number of sample points the neighborhood estimation uses. - /** This number is only used if the neighborhood radius has not been set - * manually. - * - * If the number of samples is larger than the point cloud, every point is - * used and the optimal neighborhood radius is computed exactly instead of - * estimated. - * - * \param samples is the number of points to sample for neighborhood - * estimation. - * - * \note This does not start the estimation process. - * - * \sa `neighborhood_sample_size()`. - * \sa `has_neighborhood_squared_radius()`. - * \sa `set_mean_number_of_neighbors(unsigned int neighbors)`. - * \sa `estimate_neighborhood_squared_radius()`. - */ - void set_neighborhood_sample_size( unsigned int samples ) { _samples = samples; } -/// \} - - -/// \name Scale-Space Manipulation -/// \{ - /// increases the scale by a number of iterations. - /** Each iteration the scale is increased, the points set at a higher scale - * is computed. At a higher scale, the points set is smoother. - * - * If the neighborhood radius has not been set before, it is automatically - * estimated using `estimate_neighborhood_squared_radius()`. - * - * \param iterations is the number of iterations to perform. If - * `iterations` is 0, nothing happens. - * - * \note This method processes the point set at the current scale. The - * points can be set with [insert(begin, end)](\ref insert). - * - * \note If the surface was already constructed, increasing the scale - * will not automatically adjust the surface. - * - * \sa `estimate_neighborhood_squared_radius()`. - * \sa `reconstruct_surface()`. - */ - void increase_scale( unsigned int iterations = 1 ); - - /// \cond internal_doc - /// constructs a scale-space of a collection of points. - /** If the neighborhood radius has not been set before, it is automatically - * estimated using `estimate_neighborhood_squared_radius()`. - * - * This method is equivalent to running - * `clear()` followed by - * [insert(begin, end)](\ref insert) and finally - * [increase_scale(iterations)](\ref increase_scale). - * - * \tparam InputIterator is an iterator over the point collection. - * The value type of the iterator must be a `Point`. - * - * \param begin is an iterator to the first point of the collection. - * \param end is a past-the-end iterator for the point collection. - * \param iterations is the number of iterations to perform. If - * `iterations` is 0, nothing happens. - * - * \sa `insert(InputIterator begin, InputIterator end)`. - * \sa `estimate_neighborhood_squared_radius(InputIterator begin, InputIterator end)`. - * \sa `increase_scale(unsigned int iterations)`. - * \sa `reconstruct_surface(InputIterator begin, InputIterator end, unsigned int iterations)`. - */ - template < class InputIterator > - void construct_scale_space( InputIterator begin, InputIterator end, unsigned int iterations = 1, - typename boost::enable_if< CGAL::is_iterator >::type* = NULL ) { - clear(); - insert( begin, end ); - increase_scale( iterations ); - } - ///\endcond - -/// \} -private: - // constructs the scale-space from a triangulation. - void construct_scale_space( Triangulation& tr ) { - insert( tr.finite_vertices_begin(), tr.finite_vertices_end() ); - } - - // tries to perform a functor in parallel. - template< class F > void try_parallel( const F& func, std::size_t begin, std::size_t end, Sequential_tag ) const; - template< class F > void try_parallel( const F& func, std::size_t begin, std::size_t end, Parallel_tag ) const; - template< class F > void try_parallel( const F& func, std::size_t begin, std::size_t end ) const { try_parallel( func, begin, end, Ct() );} - -private: -/// \name Shape -/// \{ - /// constructs the shape of the points at a fixed scale. - /** The shape contains geometric and connectivity information - * of the scale space. - * - * If the neighborhood radius has not been set before, it is automatically - * estimated. - */ - void construct_shape() { - construct_shape( points_begin(), points_end() ); - } - - /// constructs the shape from an existing triangulation. - /** The shape contains geometric and connectivity information - * of the scale space. - * - * \param tr is the triangulation to construct the shape of. - * - * \note This does not set the current scale-space. - * To set this as well, use `construct_scale_space(Triangulation& tr)`. - * - * \note If the neighborhood radius has not been set before, it is automatically - * estimated. - */ - void construct_shape(Triangulation& tr ) { - deinit_shape(); - if( !has_neighborhood_squared_radius() ) - estimate_neighborhood_squared_radius(); - _shape = Shape_construction_3()( *tr, _squared_radius ); - } - - /// constructs the shape from a collection of points. - /** The shape contains geometric and connectivity information - * of the scale space. - * - * \tparam InputIterator is an iterator over the point sample. - * The value type of the iterator must be a `Point`. - * \param begin is an iterator to the first point of the collection. - * \param end is a past-the-end iterator for the point collection. - * - * \note This does not set the current scale-space. - * To set this as well, use `insert( InputIterator begin, InputIterator end )`. - * - * \note If the neighborhood radius has not been set before, it is automatically - * estimated. - * - * \sa `is_constructed()`. - */ - template < class InputIterator > -#ifdef DOXYGEN_RUNNING - void construct_shape( InputIterator begin, InputIterator end ); -#else // DOXYGEN_RUNNING - void construct_shape( InputIterator begin, InputIterator end, - typename boost::enable_if< CGAL::is_iterator >::type* = NULL ); -#endif // DOXYGEN_RUNNING - - // collects the surface mesh from the shape. - // If the sahep does not yet exist, it is constructed. - void collect_surface (bool separate_shells, bool force_manifold, FT border_angle ); - - // detects the non-manifold features of the shape - void find_two_other_vertices(const Facet& f, Vertex_handle v, - Vertex_handle& v1, Vertex_handle& v2); - void detect_bubbles( FT border_angle ); - void fix_nonmanifold_edges(); - void fix_nonmanifold_vertices(); - -/// \} - -public: -/// \name Surface Reconstruction -/// \{ - /// constructs a triangle mesh from the point set at a fixed scale. - /** The order of the points at the current scale is the same as the order - * at the original scale, meaning that the surface can interpolate the - * point set at the original scale by applying the indices of the surface - * triangles to the original point set. - * - * After construction, the triangles of the surface can be iterated using - * `surface_begin()` and `surface_end()`. - * - * If the neighborhood radius has not been set before, it is automatically - * estimated using `estimate_neighborhood_squared_radius()`. - * - * \param separate_shells determines whether to collect the surface per shell. - * \param force_manifold determines if the surface is forced to be 2-manifold. - * \param border_angle sets the maximal angle between two facets - * such that the edge is seen as a border. - * - * If the output is forced to be 2-manifold, some almost flat - * volume bubbles are detected. To do so, border edges must be - * estimated. - * - * An edge adjacent to 2 regular facets is considered as a border - * if it is also adjacent to a singular facet or if the angle - * between the two regular facets is lower than this parameter - * (set to 45° by default). - * - * \note This method processes the point set at the current scale. The - * points can be set with [insert(begin, end)](\ref insert). - * \note `border_angle` is not used if `force_manifold` is set to false. - * - * \sa `estimate_neighborhood_squared_radius()`. - * \sa `increase_scale(unsigned int iterations)`. - */ - void reconstruct_surface(bool separate_shells = true, - bool force_manifold = false, - FT border_angle = 45) - { - reconstruct_surface(0, separate_shells, force_manifold, border_angle); - } - - /// gives the number of triangles of the surface. - std::size_t number_of_triangles() const { return _surface.size(); } - - /// gives the number of points of the surface. - std::size_t number_of_points() const { return _tree.size(); } - - - /// gives the number of shells of the surface. - std::size_t number_of_shells() const { - return _shells.size(); - } - - /// \cond internal_doc - void reconstruct_surface( unsigned int iterations, bool separate_shells = true, - bool force_manifold = false, FT border_angle = 45); - /// \endcond - - /// \cond internal_doc - /// constructs a surface mesh from a collection of points at a fixed scale. - /** This method is equivalent to running - * `clear()` followed by - * [insert(begin, end)](\ref insert) and finally - * [reconstruct_surface(iterations)](\ref reconstruct_surface). - * - * If the neighborhood radius has not been set before, it is automatically - * estimated using `estimate_neighborhood_squared_radius()`. - * - * \tparam InputIterator is an iterator over the point collection. - * The value type of the iterator must be a `Point`. - * - * \param begin is an iterator to the first point of the collection. - * \param end is a past-the-end iterator for the point collection. - * \param iterations is the number of scale increase iterations to apply. - * If `iterations` is 0, the point set at the current scale is used. - * \param separate_shells determines whether to collect the surface per shell. - * \param force_manifold determines if the surface is forced to be 2-manifold. - * - * \sa `reconstruct_surface(unsigned int iterations)`. - * \sa `insert(InputIterator begin, InputIterator end)`. - * \sa `estimate_neighborhood_squared_radius(InputIterator begin, InputIterator end)`. - * \sa `construct_scale_space(InputIterator begin, InputIterator end, unsigned int iterations)`. - */ - template < class InputIterator > -#ifdef DOXYGEN_RUNNING - void reconstruct_surface( InputIterator begin, InputIterator end, - unsigned int iterations = 0, bool separate_shells = true, - bool force_manifold = false, FT border_angle = 45); -#else // DOXYGEN_RUNNING - void reconstruct_surface( InputIterator begin, InputIterator end, unsigned int iterations = 0, - bool separate_shells = true, bool force_manifold = false, FT border_angle = 45, - typename boost::enable_if< CGAL::is_iterator >::type* = NULL ); -#endif // DOXYGEN_RUNNING - /// \endcond -/// \} - -public: -/// \name Iterators -/// \{ - /// gives an iterator to the first point at the current scale. - Point_const_iterator points_begin() const { return _points.begin(); } - /// gives an iterator to the first point at the current scale. - /** \warning Changes to the scale-space do not cause an automatic update to - * the surface. - */ - Point_iterator points_begin() { return _points.begin(); } - - /// gives a past-the-end iterator of the points at the current scale. - Point_const_iterator points_end() const { return _points.end(); } - /// gives a past-the-end iterator of the points at the current scale. - /** \warning Changes to the scale-space do not cause an automatic update to - * the surface. - */ - Point_iterator points_end() { return _points.end(); } - - /// gives an iterator to the first triple in the surface. - Triple_const_iterator surface_begin() const { return _surface.begin(); } - /// gives an iterator to the first triple in the surface. - /** \warning Changes to the surface may change its topology. - */ - Triple_iterator surface_begin() { return _surface.begin(); } - - /// gives a past-the-end iterator of the triples in the surface. - Triple_const_iterator surface_end() const { return _surface.end(); } - /// gives a past-the-end iterator of the triples in the surface. - /** \warning Changes to the surface may change its topology. - */ - Triple_iterator surface_end() { return _surface.end(); } - - /// gives an iterator to the first triple in a given shell. - /** \param shell is the index of the shell to access. - * - * \pre `shell` is in the range [ 0, `number_of_shells()` ). - */ - Triple_const_iterator shell_begin( std::size_t shell ) const; - /// gives an iterator to the first triple in a given shell. - /** \param shell is the index of the shell to access. - * - * \pre `shell` is in the range [ 0, `number_of_shells()` ). - * - * \warning Changes to a shell may invalidate the topology of the surface. - */ - Triple_iterator shell_begin( std::size_t shell ); - - /// gives a past-the-end iterator of the triples in a given shell. - /** \param shell is the index of the shell to access. - * - * \pre `shell` is in the range [ 0, `number_of_shells()` ). - */ - Triple_const_iterator shell_end( std::size_t shell ) const; - - /// gives a past-the-end iterator of the triples in a given shell. - /** \param shell is the index of the shell to access. - * - * \pre `shell` is in the range [ 0, `number_of_shells()` ). - * - * \warning Changes to a shell may invalidate the topology of the surface. - */ - Triple_iterator shell_end( std::size_t shell ); - - /// gives an iterator to the first triple of the garbage facets - /// that may be discarded if 2-manifold output is required. - Triple_const_iterator garbage_begin() const { return _garbage.begin(); } - /// gives an iterator to the first triple of the garbage facets - /// that may be discarded if 2-manifold output is required. - Triple_iterator garbage_begin() { return _garbage.begin(); } - - /// gives a past-the-end iterator of the triples of the garbage facets - /// that may be discarded if 2-manifold output is required. - Triple_const_iterator garbage_end() const { return _garbage.end(); } - /// gives a past-the-end iterator of the triples of the garbage facets - /// that may be discarded if 2-manifold output is required. - Triple_iterator garbage_end() { return _garbage.end(); } - -/// \} -}; // class Scale_space_surface_reconstruction_3 +template +std::ostream& operator<< (std::ostream& os, const CGAL::Scale_space_surface_reconstruction_3& ss) +{ + typedef CGAL::Scale_space_surface_reconstruction_3 Reconstruction; + + os << "OFF" << std::endl + << ss.number_of_points() << " " << ss.number_of_facets() << " 0" << std::endl; + for (typename Reconstruction::Point_const_iterator it = ss.points_begin(); + it != ss.points_end(); ++ it) + os << *it << std::endl; + for (typename Reconstruction::Facet_const_iterator it = ss.facets_begin(); + it != ss.facets_end(); ++ it) + os << "3 " << (*it)[0] << " " << (*it)[1] << " " << (*it)[2] << std::endl; + return os; +} } // namespace CGAL -template< typename T > -std::ostream& -operator<<( std::ostream& os, const CGAL::cpp11::array< T, 3 >& t ) { - return os << t[0] << " " << t[1] << " " << t[2]; -} - -template< typename T > -std::istream& -operator>>( std::istream& is, CGAL::cpp11::array< T, 3 >& t ) { - return is >> get<0>(t) >> get<1>(t) >> get<2>(t); -} - -#include - #endif // CGAL_SCALE_SPACE_SURFACE_RECONSTRUCTION_3_H