diff --git a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/PoissonDoc.cpp b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/PoissonDoc.cpp index f01c237faed..656728d7c51 100644 --- a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/PoissonDoc.cpp +++ b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/PoissonDoc.cpp @@ -758,10 +758,13 @@ void CPoissonDoc::OnReconstructionDelaunayRefinement() status_message("Delaunay refinement..."); CGAL::Timer task_timer; task_timer.start(); - const double quality = 2.5; + const double radius_edge_ratio_bound = 2.5; const unsigned int max_vertices = (unsigned int)1e7; // max 10M vertices const double enlarge_ratio = 1.5; - unsigned int nb_vertices_added = m_poisson_function->delaunay_refinement(quality,max_vertices,enlarge_ratio,50000); + const double size = sqrt(m_poisson_function->bounding_sphere().squared_radius()); // get triangulation's radius + const double cell_radius_bound = size/5.; // large + unsigned int nb_vertices_added = m_poisson_function->delaunay_refinement(radius_edge_ratio_bound,cell_radius_bound,max_vertices,enlarge_ratio); + m_triangulation_refined = true; status_message("Delaunay refinement...done (%.2lf s, %d vertices inserted)", @@ -786,7 +789,6 @@ void CPoissonDoc::OnAlgorithmsRefineInShell() status_message("Delaunay refinement in surface shell..."); CGAL::Timer task_timer; task_timer.start(); - const double quality = 2.5; const unsigned int max_vertices = (unsigned int)1e7; // max 10M vertices const double enlarge_ratio = 1.5; unsigned int nb_vertices_added = m_poisson_function->delaunay_refinement_shell(m_dr_shell_size,m_dr_sizing,m_dr_max_vertices); @@ -921,7 +923,7 @@ void CPoissonDoc::OnReconstructionPoissonSurfaceMeshing() return; } - // Get implicit surface's size + // Get implicit surface's radius Sphere bounding_sphere = m_poisson_function->bounding_sphere(); FT size = sqrt(bounding_sphere.squared_radius()); @@ -1253,7 +1255,7 @@ void CPoissonDoc::OnReconstructionApssReconstruction() return; } - // Get implicit surface's size + // Get implicit surface's radius Sphere bounding_sphere = m_apss_function->bounding_sphere(); FT size = sqrt(bounding_sphere.squared_radius()); @@ -1396,7 +1398,7 @@ void CPoissonDoc::OnPointCloudSimplificationByClustering() status_message("Point cloud simplification by clustering..."); CGAL::Timer task_timer; task_timer.start(); - // Get point set's size + // Get point set's radius Sphere bounding_sphere = m_points.bounding_sphere(); FT size = sqrt(bounding_sphere.squared_radius()); diff --git a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/Poisson_vc80.vcproj b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/Poisson_vc80.vcproj index f1f954d75f5..ceb3642c01a 100644 --- a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/Poisson_vc80.vcproj +++ b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/Poisson_vc80.vcproj @@ -45,7 +45,7 @@ AdditionalOptions="/Zm900" Optimization="0" AdditionalIncludeDirectories="include;..\..\..\include;..\..\..\..\Spatial_searching\include;..\..\..\..\Jet_fitting_3\include;"$(CGALROOT)\include\CGAL\config\msvc";"$(CGALROOT)\include";"$(CGALROOT)\auxiliary\gmp\include";"$(CGALROOT)\auxiliary\taucs\include";"$(BOOSTROOT)"" - PreprocessorDefinitions="WIN32;_SCL_SECURE_NO_DEPRECATE;_CRT_SECURE_NO_DEPRECATE;_HAS_ITERATOR_DEBUGGING=0;CGAL_USE_GMP;CGAL_USE_TAUCS" + PreprocessorDefinitions="WIN32;_SCL_SECURE_NO_DEPRECATE;_CRT_SECURE_NO_DEPRECATE;_HAS_ITERATOR_DEBUGGING=0;CGAL_USE_GMP;CGAL_USE_TAUCS;DEBUG_TRACE" MinimalRebuild="true" BasicRuntimeChecks="3" RuntimeLibrary="3" @@ -380,6 +380,10 @@ RelativePath="..\..\..\include\CGAL\average_spacing_3.h" > + + @@ -400,6 +404,10 @@ RelativePath="..\..\..\include\CGAL\Implicit_fct_delaunay_triangulation_3.h" > + + diff --git a/Surface_reconstruction_3/include/CGAL/Implicit_fct_delaunay_triangulation_3.h b/Surface_reconstruction_3/include/CGAL/Implicit_fct_delaunay_triangulation_3.h index bb5f28f7ff2..dcad2a470ed 100644 --- a/Surface_reconstruction_3/include/CGAL/Implicit_fct_delaunay_triangulation_3.h +++ b/Surface_reconstruction_3/include/CGAL/Implicit_fct_delaunay_triangulation_3.h @@ -131,11 +131,12 @@ public: // data members private: + FT m_f; // value of the implicit function // PA: should we make a separate type instead? // (so that the user can decide to run in float or double mode) bool m_constrained; // is vertex constrained? - unsigned char m_type; // INPUT or STEINER + unsigned char m_type; // INPUT or STEINER. Default type is INPUT. unsigned int m_index; // index in matrix double m_average_spacing; // average spacing int m_tag; @@ -150,8 +151,8 @@ public: m_type = 0; m_constrained = false; m_index = 0; - m_average_spacing = 0.0; - m_tag = -1; + m_average_spacing = 0.0; + m_tag = -1; } Implicit_fct_delaunay_triangulation_vertex_base_3(const Point& p) @@ -161,8 +162,8 @@ public: m_type = 0; m_constrained = false; m_index = 0; - m_average_spacing = 0.0; - m_tag = -1; + m_average_spacing = 0.0; + m_tag = -1; } @@ -173,8 +174,8 @@ public: m_type = 0; m_constrained = false; m_index = 0; - m_average_spacing = 0.0; - m_tag = -1; + m_average_spacing = 0.0; + m_tag = -1; } Implicit_fct_delaunay_triangulation_vertex_base_3(Cell_handle c) @@ -184,8 +185,8 @@ public: m_type = 0; m_constrained = false; m_index = 0; - m_average_spacing = 0.0; - tag = -1; + m_average_spacing = 0.0; + tag = -1; } /// is vertex constrained? @@ -469,20 +470,22 @@ public: m_bounding_box_is_valid = false; } - /// Insert point to the triangulation. + /// Insert point in the triangulation. + /// Default type is INPUT. Vertex_handle insert(const Point& p, unsigned char type = INPUT /* INPUT or STEINER */, Cell_handle start = Cell_handle()) { Vertex_handle v = Base::insert(p, start); + v->type() = type; - invalidate_bounding_box(); return v; } - /// Insert points to the triangulation using a spatial sort. + /// Insert points in the triangulation using a spatial sort. + /// Default type is INPUT. /// /// Precondition: the value type of InputIterator must 'Point'. /// @@ -513,6 +516,22 @@ public: return number_of_vertices() - n; } + /// Delaunay refinement callback: + /// insert STEINER point in the triangulation. + template + Vertex_handle + insert_in_hole(const Point & p, CellIt cell_begin, CellIt cell_end, + Cell_handle begin, int i, + unsigned char type = STEINER /* INPUT or STEINER */) + { + Vertex_handle v = Base::insert_in_hole(p, cell_begin, cell_end, begin, i); + + v->type() = type; + invalidate_bounding_box(); + + return v; + } + /// Index all (finite) vertices following the order of Finite_vertices_iterator. /// @return the number of (finite) vertices. unsigned int index_vertices() diff --git a/Surface_reconstruction_3/include/CGAL/Iso_cuboid_oracle_3.h b/Surface_reconstruction_3/include/CGAL/Iso_cuboid_oracle_3.h new file mode 100644 index 00000000000..edcd7cfcd40 --- /dev/null +++ b/Surface_reconstruction_3/include/CGAL/Iso_cuboid_oracle_3.h @@ -0,0 +1,72 @@ +// Copyright (c) 2006-2008 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Laurent RINEAU, Laurent Saboret + +#ifndef CGAL_ISO_CUBOID_ORACLE_3_H +#define CGAL_ISO_CUBOID_ORACLE_3_H + +#include + +CGAL_BEGIN_NAMESPACE + + +/// Mesher level oracle which tests if a point belongs to a 3D cube. +template +class Iso_cuboid_oracle_3 +{ +public: + + // Public types + typedef GT Geom_traits; + typedef typename GT::Point_3 Point; + typedef typename GT::Iso_cuboid_3 Iso_cuboid_3; + + typedef Iso_cuboid_3 Surface_3; + + typedef Point Intersection_point; + +public: + + // Constructors + Iso_cuboid_oracle_3 () + { + } + + // Predicates + bool is_in_volume(const Surface_3& cube, const Point& p) const + { + typename GT::Has_on_bounded_side_3 on_bounded_side_of_cube = + GT().has_on_bounded_side_3_object(); + + return on_bounded_side_of_cube(cube, p); + } + +}; // end Iso_cuboid_oracle_3 + + +/// Specialization of Surface_mesh_traits_generator_3 for Iso_cuboid_3. +template +struct Surface_mesh_traits_generator_3 > +{ + typedef Iso_cuboid_oracle_3 Type; + typedef Type type; // for Boost compatiblity (meta-programming) +}; + + +CGAL_END_NAMESPACE + +#endif // CGAL_ISO_CUBOID_ORACLE_3_H diff --git a/Surface_reconstruction_3/include/CGAL/Poisson_implicit_function.h b/Surface_reconstruction_3/include/CGAL/Poisson_implicit_function.h index 382cf275823..e6765d3b2e7 100644 --- a/Surface_reconstruction_3/include/CGAL/Poisson_implicit_function.h +++ b/Surface_reconstruction_3/include/CGAL/Poisson_implicit_function.h @@ -36,6 +36,7 @@ #include #include #include +#include CGAL_BEGIN_NAMESPACE @@ -179,12 +180,6 @@ private: typedef typename CGAL::Point_vertex_handle_3 Point_vertex_handle_3; K_nearest_neighbor m_nn_search; - // delaunay refinement - typedef typename CGAL::Candidate Candidate; - typedef typename std::priority_queue< Candidate, - std::vector, - less > Refinement_pqueue; - // contouring and meshing Point m_sink; // Point with the minimum value of f() mutable Cell_handle m_hint; // last cell found = hint for next search @@ -294,16 +289,16 @@ public: CGAL_TRACE_STREAM << "Delaunay refinement...\n"; // Delaunay refinement - int nb_vertices = m_dt.number_of_vertices(); - const FT quality = 2.5; + const FT radius_edge_ratio_bound = 2.5; const unsigned int max_vertices = (unsigned int)1e7; // max 10M vertices const FT enlarge_ratio = 1.5; - delaunay_refinement(quality,max_vertices,enlarge_ratio,50000); + const FT size = sqrt(bounding_sphere().squared_radius()); // get triangulation's radius + const FT cell_radius_bound = size/5.; // large + unsigned int nb_vertices_added = delaunay_refinement(radius_edge_ratio_bound,cell_radius_bound,max_vertices,enlarge_ratio); // Print status long memory = CGAL::Memory_sizer().virtual_size(); - int nb_vertices2 = m_dt.number_of_vertices(); - CGAL_TRACE_STREAM << "Delaunay refinement: " << "added " << nb_vertices2-nb_vertices << " Steiner points, " + CGAL_TRACE_STREAM << "Delaunay refinement: " << "added " << nb_vertices_added << " Steiner points, " << task_timer.time() << " seconds, " << (memory>>20) << " Mb allocated" << std::endl; @@ -407,92 +402,33 @@ public: /// bad means badly shaped or too big). The normal of /// Steiner points is set to zero. /// Return the number of vertices inserted. - unsigned int delaunay_refinement(const FT threshold, - unsigned int maximum, - const FT enlarge_ratio, - unsigned int restart_each) + unsigned int delaunay_refinement(FT radius_edge_ratio_bound, ///< radius edge ratio bound (ignored if zero) + FT cell_radius_bound, ///< cell radius bound (ignored if zero) + unsigned int max_vertices, ///< number of vertices bound + FT enlarge_ratio) ///< bounding box enlarge ratio { - CGAL_TRACE("Call delaunay_refinement()\n"); + CGAL_TRACE("Call delaunay_refinement(radius_edge_ratio_bound=%lf, cell_radius_bound=%lf, max_vertices=%u, enlarge_ratio=%lf)\n", + radius_edge_ratio_bound, cell_radius_bound, max_vertices, enlarge_ratio); - // create enlarged bounding box +#define DELAUNAY_REFINEMENT_USE_BOUNDING_BOX +#ifdef DELAUNAY_REFINEMENT_USE_BOUNDING_BOX Iso_cuboid enlarged_bbox = enlarged_bounding_box(enlarge_ratio); + unsigned int nb_vertices_added = poisson_refinement_3(m_dt,radius_edge_ratio_bound,cell_radius_bound,max_vertices,enlarged_bbox); +#else + Sphere enlarged_bbox = enlarged_bounding_sphere(enlarge_ratio); + unsigned int nb_vertices_added = poisson_refinement_3(m_dt,radius_edge_ratio_bound,cell_radius_bound,max_vertices,enlarged_bbox); +#endif - long memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20); - CGAL_TRACE(" Create queue\n"); - - long max_memory = memory; // max memory allocated by this algorithm - - // push bad cells to the queue - Refinement_pqueue queue; - init_queue(queue,threshold,enlarged_bbox); - - unsigned int nb = 0; - while(!queue.empty()) - { - memory = CGAL::Memory_sizer().virtual_size(); max_memory = (std::max)(max_memory, memory); - - if(nb%restart_each == 0) - { - reset_queue(queue,threshold,enlarged_bbox); - if(queue.empty()) - break; - } - - Candidate candidate = queue.top(); - queue.pop(); - Vertex_handle v0 = candidate.v0(); - Vertex_handle v1 = candidate.v1(); - Vertex_handle v2 = candidate.v2(); - Vertex_handle v3 = candidate.v3(); - Cell_handle cell = NULL; - if(m_dt.is_cell(v0,v1,v2,v3,cell)) - { - Point point = m_dt.dual(cell); - Vertex_handle v = m_dt.insert(point,Triangulation::STEINER,cell); - - if(nb++ > maximum) - break; // premature ending - - // iterate over incident cells and feed queue - std::list cells; - m_dt.incident_cells(v,std::back_inserter(cells)); - typename std::list::iterator it; - for(it = cells.begin(); - it != cells.end(); - it++) - { - Cell_handle c = *it; - if(m_dt.is_infinite(c)) - continue; - - FT rer = radius_edge_ratio(c); - Point point = m_dt.dual(c); - bool inside = enlarged_bbox.has_on_bounded_side(point); - if(inside && rer > threshold) - { - Vertex_handle v0 = c->vertex(0); - Vertex_handle v1 = c->vertex(1); - Vertex_handle v2 = c->vertex(2); - Vertex_handle v3 = c->vertex(3); - float score = (float)max_edge_len(cell); - queue.push(Candidate(v0,v1,v2,v3,score)); - } - } - } - } m_dt.invalidate_bounding_box(); - - CGAL_TRACE(" Max allocation = %ld Mb\n", max_memory>>20); - - /*long*/ memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20); + CGAL_TRACE("End of delaunay_refinement()\n"); - return nb; + return nb_vertices_added; } unsigned int delaunay_refinement_shell(FT size_shell, FT sizing, - unsigned int maximum) + unsigned int max_vertices) { // make parameters relative to size Sphere bounding_sphere = m_dt.bounding_sphere(); @@ -539,10 +475,10 @@ public: Cell_handle cell = NULL; if(m_dt.is_cell(v0,v1,v2,v3,cell)) { - Point point = m_dt.dual(cell); - Vertex_handle v = m_dt.insert(point, Triangulation::STEINER); + Point circumcenter = m_dt.dual(cell); + Vertex_handle v = m_dt.insert(circumcenter, Triangulation::STEINER); - if(nb++ > maximum) + if(nb++ > max_vertices) return nb; // premature ending // iterate over incident cells and feed queue @@ -1068,31 +1004,6 @@ private: d = std::fabs(td.volume() / v); } - // radius-edge ratio (the ratio of the circumradius to - // the shortest edge length of tetrahedron) - // check template type - double radius_edge_ratio(Cell_handle c) const - { - double r = circumradius(c); - double l = len_shortest_edge(c); - if(l != 0.0) - return r/l; - else - return 1e38; - } - - FT len_shortest_edge(Cell_handle c) const - { - FT d1 = distance(c->vertex(0),c->vertex(1)); - FT d2 = distance(c->vertex(0),c->vertex(2)); - FT d3 = distance(c->vertex(0),c->vertex(3)); - FT d4 = distance(c->vertex(1),c->vertex(2)); - FT d5 = distance(c->vertex(1),c->vertex(3)); - FT d6 = distance(c->vertex(2),c->vertex(3)); - return (std::min)((std::min)((std::min)(d1,d2),d3), - (std::min)((std::min)(d4,d5),d6)); - } - FT distance(Vertex_handle v1, Vertex_handle v2) const { const Point& a = v1->point(); @@ -1551,11 +1462,11 @@ private: return Edge(cell,i1,i2); } - /// Compute enlarged geometric bounding box of the embedded triangulation + /// Compute enlarged geometric bounding box of the embedded triangulation. Iso_cuboid enlarged_bounding_box(FT ratio) const { // Get triangulation's bounding box - Iso_cuboid bbox = m_dt.bounding_box(); + Iso_cuboid bbox = bounding_box(); // Its center point is: FT mx = 0.5 * (bbox.xmax() + bbox.xmin()); @@ -1572,41 +1483,11 @@ private: return Iso_cuboid(p,q); } - void reset_queue(Refinement_pqueue& queue, - const FT threshold, - Iso_cuboid& enlarged_bbox) + /// Compute enlarged geometric bounding sphere of the embedded triangulation. + Sphere enlarged_bounding_sphere(FT ratio) const { - // clear queue - while(!queue.empty()) - queue.pop(); - - // push bad cells to the queue - init_queue(queue,threshold,enlarged_bbox); - } - - // push bad cells to the queue - void init_queue(Refinement_pqueue& queue, - const FT threshold, - Iso_cuboid& enlarged_bbox) - { - Finite_cells_iterator c; - for(c = m_dt.finite_cells_begin(); - c != m_dt.finite_cells_end(); - c++) - { - FT rer = radius_edge_ratio(c); - Point point = m_dt.dual(c); - bool inside = enlarged_bbox.has_on_bounded_side(point); - if(inside && rer > threshold) - { - Vertex_handle v0 = c->vertex(0); - Vertex_handle v1 = c->vertex(1); - Vertex_handle v2 = c->vertex(2); - Vertex_handle v3 = c->vertex(3); - float score = (float)max_edge_len(c); - queue.push(Candidate(v0,v1,v2,v3,score)); - } - } + Sphere bbox = bounding_sphere(); // triangulation's bounding sphere + return Sphere(bbox.center(), bbox.squared_radius() * ratio*ratio); } void init_nn_search_shell() @@ -1662,23 +1543,6 @@ private: return std::sqrt(CGAL::squared_distance(a,b)); } - FT max_edge_len(Cell_handle cell) const - { - const Point& a = cell->vertex(0)->point(); - const Point& b = cell->vertex(1)->point(); - const Point& c = cell->vertex(2)->point(); - const Point& d = cell->vertex(3)->point(); - FT ab = (a-b)*(a-b); - FT ac = (a-c)*(a-c); - FT bc = (c-b)*(c-b); - FT ad = (a-d)*(a-d); - FT bd = (d-b)*(d-b); - FT cd = (c-d)*(c-d); - FT sq_max = (std::max)((std::max)((std::max)(ab,ac), - (std::max)(bc,ad)), (std::max)(bd,cd)); - return std::sqrt(sq_max); - } - FT circumradius(Cell_handle c) const { Point center = m_dt.dual(c); diff --git a/Surface_reconstruction_3/include/CGAL/poisson_refinement_3.h b/Surface_reconstruction_3/include/CGAL/poisson_refinement_3.h new file mode 100644 index 00000000000..dd2ff57571b --- /dev/null +++ b/Surface_reconstruction_3/include/CGAL/poisson_refinement_3.h @@ -0,0 +1,243 @@ +// Copyright (c) 2007-2008 INRIA (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Laurent RINEAU, Laurent Saboret + +#ifndef CGAL_POISSON_REFINEMENT_3_H +#define CGAL_POISSON_REFINEMENT_3_H + +// CGAL +#include +#include +#include +#include +#include + +CGAL_BEGIN_NAMESPACE + + +/// Utility class for poisson_refinement_3(): +/// implements Delaunay refinement (break bad tetrahedra, where +/// bad means badly shaped or too big). +/// +/// This class must be derived to inherit from Mesher_level. +template +> +class Poisson_mesher_level_impl_base : + public Mesh_3::Refine_tets_with_oracle_base +{ +protected: + typedef Mesh_3::Refine_tets_with_oracle_base Base; + +public: + typedef typename Tr::Geom_traits Geom_traits; + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Cell_handle Cell_handle; + typedef typename Tr::Point Point; + typedef typename Base::Quality Quality; + + using Base::triangulation_ref_impl; + +public: + /** \name CONSTRUCTORS */ + + Poisson_mesher_level_impl_base(Tr& t, Criteria crit, unsigned int max_vert, Surface& surface, Oracle& oracle) + : Base(t, crit, surface, oracle), + max_vertices(max_vert), ///< number of vertices bound (ignored if zero) + max_memory(CGAL::Memory_sizer().virtual_size()) ///< max memory allocated by this algorithm + { + } + +protected: + /* --- protected functions --- */ + + bool test_if_cell_is_bad(const Cell_handle c) + { + Quality q; + if( is_in_domain(c) && should_be_refined(c, q) ) + { + this->add_bad_element(c, q); + return true; + } + return false; + } + + bool is_in_domain(const Cell_handle c) + { + return oracle.is_in_volume(surface, triangulation_ref_impl().dual(c)); + } + +public: + /* Overriden functions of this level: */ + /* we override all methods that call test_if_cell_is_bad() */ + + void scan_triangulation_impl() + { + for(typename Tr::Finite_cells_iterator cit = triangulation_ref_impl().finite_cells_begin(); + cit != triangulation_ref_impl().finite_cells_end(); + ++cit) + test_if_cell_is_bad(cit); + } + + void after_insertion_impl(const Vertex_handle& v) + { + update_star(v); + + // Update used memory + long memory = CGAL::Memory_sizer().virtual_size(); + max_memory = (std::max)(max_memory, memory); + } + + void update_star(const Vertex_handle& v) + { + // scan refiner + typedef std::vector Cells; + typedef typename Cells::iterator Cell_iterator; + Cells incident_cells; + + triangulation_ref_impl().incident_cells(v, std::back_inserter(incident_cells)); + + for(Cell_iterator cit = incident_cells.begin(); + cit != incident_cells.end(); + ++cit) + { + if( ! triangulation_ref_impl().is_infinite(*cit) ) + test_if_cell_is_bad(*cit); + } + } + + /// Tells if the algorithm is done. + bool no_longer_element_to_refine_impl() const + { + return Base::no_longer_element_to_refine_impl() || + (max_vertices > 0 && triangulation_ref_impl().number_of_vertices() >= max_vertices); + } + + /// Give max memory allocated by this algorithm. + long max_memory_allocated() const + { + return max_memory; + } + +private: + /* --- private datas --- */ + unsigned int max_vertices; ///< number of vertices bound (ignored if zero) + long max_memory; ///< max memory allocated by this algorithm + +}; // end Poisson_mesher_level_impl_base + + +/// Utility class for poisson_refinement_3(): +/// glue class that inherits from both Mesher_level +/// and Poisson_mesher_level_impl_base. +template ::type + > +class Poisson_mesher_level : + public Poisson_mesher_level_impl_base, + public Mesher_level < + Tr, + Poisson_mesher_level, + typename Tr::Cell_handle, + Null_mesher_level, + Triangulation_mesher_level_traits_3 + > +{ + typedef Poisson_mesher_level_impl_base Base; + +public: + typedef Mesher_level < + Tr, + Poisson_mesher_level, + typename Tr::Cell_handle, + Null_mesher_level, + Triangulation_mesher_level_traits_3 + > Mesher; + + Poisson_mesher_level(Tr& t, Criteria criteria, unsigned int max_vertices, Surface& surface, Oracle& oracle = Oracle()) + : Base(t, criteria, max_vertices, surface, oracle), + Mesher(Null_mesher_level()) + { + } + +}; // end class Poisson_mesher_level + + +/// Delaunay refinement (break bad tetrahedra, where +/// bad means badly shaped or too big). +/// @return the number of vertices inserted. +/// +/// Precondition: +/// convergence is guaranteed if radius_edge_ratio_bound >= 1.0. +/// +/// @heading Parameters: +/// @param Tr 3D Delaunay triangulation. +/// @param Surface Sphere_3 or Iso_cuboid_3. +template +unsigned int poisson_refinement_3(Tr& tr, + double radius_edge_ratio_bound, ///< radius edge ratio bound (>= 1.0) + double cell_radius_bound, ///< cell radius bound (ignored if zero) + unsigned int max_vertices, ///< number of vertices bound (ignored if zero) + Surface& enlarged_bbox) ///< new bounding sphere or box +{ + CGAL_TRACE("Call poisson_refinement_3()\n"); + + // Convergence is guaranteed if radius_edge_ratio_bound >= 1.0 + CGAL_precondition(radius_edge_ratio_bound >= 1.0); + + // Basic geometric types + typedef typename Tr::Geom_traits Gt; + typedef typename Gt::FT FT; + typedef typename Gt::Point_3 Point; + + // Mesher_level types + typedef Mesh_criteria_3 Tets_criteria; + typedef Poisson_mesher_level Refiner; + + long memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20); + CGAL_TRACE(" Create queue\n"); + + int nb_vertices = tr.number_of_vertices(); // get former #vertices + + // Delaunay refinement + Tets_criteria tets_criteria(radius_edge_ratio_bound*radius_edge_ratio_bound, cell_radius_bound); + Refiner refiner(tr, tets_criteria, max_vertices, enlarged_bbox); + refiner.scan_triangulation(); // Push bad cells to the queue + refiner.refine(Null_mesh_visitor()); // Refine triangulation until queue is empty + + int nb_vertices_added = tr.number_of_vertices() - nb_vertices; + + long max_memory = refiner.max_memory_allocated(); + CGAL_TRACE(" Max allocation = %ld Mb\n", max_memory>>20); + + /*long*/ memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20); + CGAL_TRACE("End of poisson_refinement_3()\n"); + + return (unsigned int) nb_vertices_added; +} + + +CGAL_END_NAMESPACE + +#endif // CGAL_POISSON_REFINEMENT_3_H diff --git a/Surface_reconstruction_3/test/Surface_reconstruction_3/VC/poisson_reconstruction_test_vc80.vcproj b/Surface_reconstruction_3/test/Surface_reconstruction_3/VC/poisson_reconstruction_test_vc80.vcproj index d2d1220eaa3..c62cce818d2 100644 --- a/Surface_reconstruction_3/test/Surface_reconstruction_3/VC/poisson_reconstruction_test_vc80.vcproj +++ b/Surface_reconstruction_3/test/Surface_reconstruction_3/VC/poisson_reconstruction_test_vc80.vcproj @@ -41,7 +41,7 @@ AdditionalOptions="/Zm900" Optimization="0" AdditionalIncludeDirectories="..\include;..\..\..\include;..\..\..\..\Spatial_searching\include;..\..\..\..\Jet_fitting_3\include;"$(CGALROOT)\include\CGAL\config\msvc";"$(CGALROOT)\include";"$(CGALROOT)\auxiliary\gmp\include";"$(CGALROOT)\auxiliary\taucs\include";"$(BOOSTROOT)"" - PreprocessorDefinitions="WIN32;_SCL_SECURE_NO_DEPRECATE;_CRT_SECURE_NO_DEPRECATE;_HAS_ITERATOR_DEBUGGING=0;CGAL_USE_GMP;CGAL_USE_TAUCS" + PreprocessorDefinitions="WIN32;_SCL_SECURE_NO_DEPRECATE;_CRT_SECURE_NO_DEPRECATE;_HAS_ITERATOR_DEBUGGING=0;CGAL_USE_GMP;CGAL_USE_TAUCS;DEBUG_TRACE" MinimalRebuild="true" BasicRuntimeChecks="3" RuntimeLibrary="3" @@ -205,6 +205,10 @@ Name="Surface_reconstruction_3 package" Filter="h;hpp;hxx;hm;inl;inc" > + +