Delaunay refinement for Poisson is moved from Poisson_implicit_function class to new poisson_refinement_3() function.

poisson_refinement_3() is implemented with mesher levels. See Mesh_3 package.

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

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"
 				>
 			</File>
+			<File
+				RelativePath="..\..\..\include\CGAL\poisson_refinement_3.h"
+				>
+			</File>
 			<File
 				RelativePath="..\..\..\include\CGAL\estimate_normals_3.h"
 				>
@@ -400,6 +404,10 @@
 				RelativePath="..\..\..\include\CGAL\Implicit_fct_delaunay_triangulation_3.h"
 				>
 			</File>
+			<File
+				RelativePath="..\..\..\include\CGAL\Iso_cuboid_oracle_3.h"
+				>
+			</File>
 			<File
 				RelativePath="..\..\..\include\CGAL\k_nearest_neighbor.h"
 				>

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_average_spacing = 0.0;
-	m_tag = -1;
+    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_average_spacing = 0.0;
-	m_tag = -1;
+    m_tag = -1;
 	
   }
 
@@ -173,8 +174,8 @@ public:
     m_type = 0;
     m_constrained = false;
     m_index = 0;
-	m_average_spacing = 0.0;
+    m_average_spacing = 0.0;
-	m_tag = -1;
+    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;
+    m_average_spacing = 0.0;
-	tag = -1;
+    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 <class CellIt>
+  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()

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/Surface_mesh_traits_generator_3.h>
+
+CGAL_BEGIN_NAMESPACE
+
+
+/// Mesher level oracle which tests if a point belongs to a 3D cube.
+template <class GT>
+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 <typename Kernel>
+struct Surface_mesh_traits_generator_3<CGAL::Iso_cuboid_3<Kernel> >
+{
+  typedef Iso_cuboid_oracle_3<Kernel> Type;
+  typedef Type type; // for Boost compatiblity (meta-programming)
+};
+
+
+CGAL_END_NAMESPACE
+
+#endif  // CGAL_ISO_CUBOID_ORACLE_3_H

Surface_reconstruction_3/include/CGAL/Poisson_implicit_function.h

@@ -36,6 +36,7 @@
 #include <CGAL/surface_reconstruction_assertions.h>
 #include <CGAL/remove_outliers_wrt_avg_knn_sq_distance_3.h>
 #include <CGAL/Memory_sizer.h>
+#include <CGAL/poisson_refinement_3.h>
 
 CGAL_BEGIN_NAMESPACE
 
@@ -179,12 +180,6 @@ private:
   typedef typename CGAL::Point_vertex_handle_3<Vertex_handle> Point_vertex_handle_3;
   K_nearest_neighbor m_nn_search;
 
-  // delaunay refinement
-  typedef typename CGAL::Candidate<Vertex_handle,Point> Candidate;
-  typedef typename std::priority_queue< Candidate,
-                                        std::vector<Candidate>,
-                                        less<Candidate> > 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 radius_edge_ratio_bound = 2.5;
-    const FT quality = 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_vertices_added << " Steiner points, "
-    CGAL_TRACE_STREAM << "Delaunay refinement: " << "added " << nb_vertices2-nb_vertices << " 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 delaunay_refinement(FT radius_edge_ratio_bound, ///< radius edge ratio bound (ignored if zero)
-                                   unsigned int maximum,
+                                   FT cell_radius_bound, ///< cell radius bound (ignored if zero)
-                                   const FT enlarge_ratio,
+                                   unsigned int max_vertices, ///< number of vertices bound
-                                   unsigned int restart_each)
+                                   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<Cell_handle> cells;
-        m_dt.incident_cells(v,std::back_inserter(cells));
-        typename std::list<Cell_handle>::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);
+        Point circumcenter = m_dt.dual(cell);
-        Vertex_handle v = m_dt.insert(point, Triangulation::STEINER);
+        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,
+  /// Compute enlarged geometric bounding sphere of the embedded triangulation.
-                   const FT threshold,
+  Sphere enlarged_bounding_sphere(FT ratio) const
-                   Iso_cuboid& enlarged_bbox)
   {
-    // clear queue
+    Sphere bbox = bounding_sphere(); // triangulation's bounding sphere
-    while(!queue.empty())
+    return Sphere(bbox.center(), bbox.squared_radius() * ratio*ratio);
-      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));
-      }
-    }
   }
 
   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);

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 <CGAL/Mesher_level.h>
+#include <CGAL/Mesh_3/Refine_tets.h>
+#include <CGAL/Mesh_criteria_3.h>
+#include <CGAL/Iso_cuboid_oracle_3.h>
+#include <CGAL/surface_reconstruction_assertions.h>
+
+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 Tr,
+          class Criteria,
+          class Surface,
+          class Oracle,
+          class Container = Meshes::Double_map_container<typename Tr::Cell_handle,
+                                                         typename Criteria::Quality>
+>
+class Poisson_mesher_level_impl_base :
+    public Mesh_3::Refine_tets_with_oracle_base<Tr, Criteria, Surface, Oracle, Container>
+{
+protected:
+  typedef Mesh_3::Refine_tets_with_oracle_base<Tr, Criteria, Surface, Oracle, Container> 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<Cell_handle> 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 <typename Tr,
+          typename Criteria,
+          typename Surface,
+          typename Oracle = typename CGAL::Surface_mesh_traits_generator_3<Surface>::type
+ >
+class Poisson_mesher_level : 
+  public Poisson_mesher_level_impl_base<Tr, Criteria, Surface, Oracle>, 
+  public Mesher_level <
+    Tr,
+    Poisson_mesher_level<Tr, Criteria, Surface, Oracle>,
+    typename Tr::Cell_handle,
+    Null_mesher_level,
+    Triangulation_mesher_level_traits_3<Tr>
+  >
+{
+  typedef Poisson_mesher_level_impl_base<Tr, Criteria, Surface, Oracle> Base;
+
+public:
+  typedef Mesher_level <
+    Tr,
+    Poisson_mesher_level<Tr, Criteria, Surface, Oracle>,
+    typename Tr::Cell_handle,
+    Null_mesher_level,
+    Triangulation_mesher_level_traits_3<Tr>
+  > 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 <typename Tr,
+          typename Surface>
+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<Tr> Tets_criteria;
+  typedef Poisson_mesher_level<Tr, Tets_criteria, Surface> 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

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"
 			>
+			<File
+				RelativePath="..\..\..\include\CGAL\poisson_refinement_3.h"
+				>
+			</File>
 			<File
 				RelativePath="..\..\..\include\CGAL\Implicit_fct_delaunay_triangulation_3.h"
 				>