Merge pull request #1271 from lrineau/Mesh_3-distance_criterion_sizing_field-GF

Allow a sizing field for the distance criterion

Mesh_3/include/CGAL/Mesh_3/mesh_standard_facet_criteria.h

@@ -154,6 +154,11 @@ private:
 template <typename Tr, typename Visitor_>
 class Curvature_size_criterion :
   public Mesh_3::Abstract_criterion<Tr, Visitor_>
+{};
+
+template <typename Tr, typename Visitor_>
+class Uniform_curvature_size_criterion :
+  public Curvature_size_criterion<Tr, Visitor_>
 {
 private:
   typedef typename Tr::Facet Facet;
@@ -163,12 +168,12 @@ private:
   typedef typename Base::Quality Quality;
   typedef typename Base::Badness Badness;
 
-  typedef Curvature_size_criterion<Tr,Visitor_> Self;
+  typedef Uniform_curvature_size_criterion<Tr,Visitor_> Self;
 
 public:
   // Nb: the default bound of the criterion is such that the criterion
   // is always fulfilled
-  Curvature_size_criterion(const FT b = 0) : B_(b * b) {}
+  Uniform_curvature_size_criterion(const FT b = 0) : B_(b * b) {}
 
 protected:
   virtual void do_accept(Visitor_& v) const
@@ -218,8 +223,84 @@ protected:
 private:
   FT B_;
 
-};  // end Curvature_size_criterion
+};  // end Uniform_curvature_size_criterion
 
+// Variable size Criterion class
+template <typename Tr, typename Visitor_, typename SizingField>
+class Variable_curvature_size_criterion :
+  public Curvature_size_criterion<Tr, Visitor_>
+{
+private:
+  typedef typename Tr::Facet            Facet;
+  typedef typename Tr::Geom_traits::FT  FT;
+  typedef typename Tr::Vertex::Index    Index;
+
+  typedef Mesh_3::Abstract_criterion<Tr,Visitor_> Base;
+  typedef typename Base::Quality Quality;
+  typedef typename Base::Badness Badness;
+
+  typedef Variable_curvature_size_criterion<Tr,Visitor_,SizingField> Self;
+  typedef SizingField Sizing_field;
+
+public:
+  // Nb: the default bound of the criterion is such that the criterion
+  // is always fulfilled
+  Variable_curvature_size_criterion(const Sizing_field& s) : size_(s) {}
+
+protected:
+  virtual void do_accept(Visitor_& v) const
+  {
+    v.visit(*this);
+  }
+
+  virtual Self* do_clone() const
+  {
+    // Call copy ctor on this
+    return new Self(*this);
+  }
+
+  virtual Badness do_is_bad (const Facet& f) const
+  {
+    CGAL_assertion (f.first->is_facet_on_surface(f.second));
+
+    typedef typename Tr::Geom_traits Gt;
+    typedef typename Tr::Point Point_3;
+
+    typename Gt::Compute_squared_distance_3 distance =
+        Gt().compute_squared_distance_3_object();
+    typename Gt::Construct_weighted_circumcenter_3 circumcenter =
+        Gt().construct_weighted_circumcenter_3_object();
+
+    const Point_3& p1 = f.first->vertex((f.second+1)&3)->point();
+    const Point_3& p2 = f.first->vertex((f.second+2)&3)->point();
+    const Point_3& p3 = f.first->vertex((f.second+3)&3)->point();
+
+    const Point_3 c = circumcenter(p1,p2,p3);
+    const Point_3& ball_center = f.first->get_facet_surface_center(f.second);
+
+    const FT sq_dist = distance(c, ball_center);
+
+    const Index& index = f.first->get_facet_surface_center_index(f.second);
+
+    const FT sq_bound = CGAL::square(size_(ball_center, 2, index));
+    CGAL_assertion(sq_bound > FT(0));
+
+    if ( sq_dist > sq_bound )
+    {
+#ifdef CGAL_MESH_3_DEBUG_FACET_CRITERIA
+      std::cerr << "Bad facet (curvature size): sq_dist[" << sq_dist
+      << "] bound[" << sq_bound << "]\n";
+#endif
+      return Badness(Quality(sq_bound/sq_dist));
+    }
+    else
+      return Badness();
+  }
+
+private:
+  Sizing_field size_;
+
+};  // end Variable_curvature_size_criterion
 
 // Size Criterion base class
 template < typename Tr, typename Visitor_ >

Mesh_3/include/CGAL/Mesh_facet_criteria_3.h

@@ -65,10 +65,10 @@ public:
       init_aspect(angle_bound);
     
     if ( FT(0) != radius_bound )
-      init_radius(radius_bound);
+      init_radius_bound(radius_bound);
     
     if ( FT(0) != distance_bound )
-      init_distance(distance_bound);
+      init_distance_bound(distance_bound);
     
     init_topo(topology);
   }
@@ -86,15 +86,57 @@ public:
     if ( FT(0) != angle_bound )
       init_aspect(angle_bound);
     
-    init_radius(radius_bound);
+    init_radius_field(radius_bound);
     
     if ( FT(0) != distance_bound )
-      init_distance(distance_bound);
+      init_distance_bound(distance_bound);
     
     init_topo(topology);  
   }
   
-  /// Destructor
+  // Nb: SFINAE (dummy) to avoid wrong matches with built-in numerical types
+  // as int.
+  template < typename Sizing_field >
+  Mesh_facet_criteria_3(const FT& angle_bound,
+                        const FT& radius_bound,
+                        const Sizing_field& distance_bound,
+                        const Mesh_facet_topology topology =
+                          FACET_VERTICES_ON_SURFACE,
+                        typename Sizing_field::FT /*dummy*/ = 0)
+  {
+    if ( FT(0) != angle_bound )
+      init_aspect(angle_bound);
+
+    if ( FT(0) != radius_bound )
+      init_radius_bound(radius_bound);
+
+    init_distance_field(distance_bound);
+
+    init_topo(topology);
+  }
+
+  // Nb: SFINAE (dummy) to avoid wrong matches with built-in numerical types
+  // as int.
+  template < typename Sizing_field, typename Sizing_field2 >
+  Mesh_facet_criteria_3(const FT& angle_bound,
+                        const Sizing_field & radius_bound,
+                        const Sizing_field2& distance_bound,
+                        const Mesh_facet_topology topology =
+                          FACET_VERTICES_ON_SURFACE,
+                        typename Sizing_field::FT /*dummy*/ = 0,
+                        typename Sizing_field2::FT /*dummy*/ = 0)
+  {
+    if ( FT(0) != angle_bound )
+      init_aspect(angle_bound);
+
+    init_radius_field(radius_bound);
+
+    init_distance_field(distance_bound);
+
+    init_topo(topology);
+  }
+
+/// Destructor
   ~Mesh_facet_criteria_3() { }
 
    /**
@@ -123,23 +165,32 @@ private:
     criteria_.add(new Aspect_criterion(angle_bound));
   }
   
-  void init_radius(const FT& radius_bound)
+  void init_radius_bound(const FT& radius_bound)
   {
     typedef Mesh_3::Uniform_size_criterion<Tr,Visitor> Uniform_size_criterion;
     criteria_.add(new Uniform_size_criterion(radius_bound));
   }
   
   template <typename Sizing_field>
-  void init_radius(const Sizing_field& radius_bound)
+  void init_radius_field(const Sizing_field& radius_bound)
   {
     typedef Mesh_3::Variable_size_criterion<Tr,Visitor,Sizing_field> Variable_size_criterion;
     criteria_.add(new Variable_size_criterion(radius_bound));
   }
   
-  void init_distance(const FT& distance_bound)
+  void init_distance_bound(const FT& distance_bound)
   {
-    typedef Mesh_3::Curvature_size_criterion<Tr,Visitor> Curvature_criterion;
-    criteria_.add(new Curvature_criterion(distance_bound));
+    typedef Mesh_3::Uniform_curvature_size_criterion<Tr,Visitor> Criterion;
+    criteria_.add(new Criterion(distance_bound));
+  }
+  
+  template <typename Sizing_field>
+  void init_distance_field(const Sizing_field& distance_bound)
+  {
+    typedef Mesh_3::Variable_curvature_size_criterion<Tr,
+                                                      Visitor,
+                                                      Sizing_field> Criterion;
+    criteria_.add(new Criterion(distance_bound));
   }
   
   void init_topo(const Mesh_facet_topology topology)

Mesh_3/test/Mesh_3/CMakeLists.txt

@@ -27,6 +27,7 @@ if ( CGAL_FOUND )
   create_single_source_cgal_program( "test_backward_compatibility.cpp" )
   create_single_source_cgal_program( "test_boost_has_xxx.cpp" )
   create_single_source_cgal_program( "test_c3t3.cpp" )
+  create_single_source_cgal_program( "test_mesh_capsule_var_distance_bound.cpp" )
   create_single_source_cgal_program( "test_implicit_multi_domain_to_labeling_function_wrapper.cpp" )
   create_single_source_cgal_program( "test_mesh_3_implicit_vector_to_labeled_function_wrapper.cpp" )
   create_single_source_cgal_program( "test_c3t3_io.cpp" )

Mesh_3/test/Mesh_3/test_mesh_capsule_var_distance_bound.cpp

@@ -0,0 +1,74 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+
+#include <CGAL/Implicit_mesh_domain_3.h>
+#include <CGAL/make_mesh_3.h>
+
+// Domain
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::FT FT;
+typedef K::Point_3 Point;
+typedef FT (Function)(const Point&);
+typedef CGAL::Implicit_mesh_domain_3<Function,K> Mesh_domain;
+
+#ifdef CGAL_CONCURRENT_MESH_3
+typedef CGAL::Parallel_tag Concurrency_tag;
+#else
+typedef CGAL::Sequential_tag Concurrency_tag;
+#endif
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3<Mesh_domain,K,Concurrency_tag>::type Tr;
+
+typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
+
+// Criteria
+typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+
+// To avoid verbose function and named parameters call
+using namespace CGAL::parameters;
+
+// Function, a capsule:
+//   a cylinder centered at x=y=0, with z in [-5, -5]
+//   with radius 1 and "round ends".
+FT capsule_function(const Point& p)
+{
+  const FT base = CGAL::square(p.x())+CGAL::square(p.y()) - 1;
+  const FT z = p.z();
+  if(z > FT(5)) return base+CGAL::square(z-5);
+  else if(z < FT(-5)) return base+CGAL::square(z+5);
+  else return base;
+}
+
+struct Field {
+  typedef ::FT FT;
+  
+  FT operator()(const Point& p, const int, const Mesh_domain::Index) const {
+    if(p.z() > 2) return 0.025;
+    if(p.z() < -3) return 0.01;
+    else return 1;
+  }
+} field;
+
+int main()
+{
+  Mesh_domain domain(capsule_function,
+                     K::Sphere_3(CGAL::ORIGIN, 49.));
+
+  // Mesh criteria
+  Mesh_criteria criteria(facet_angle=30, facet_size=0.5,
+			 facet_distance=field);
+  
+  // Mesh generation
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
+
+  // Output
+  std::ofstream medit_file("out.mesh");
+  c3t3.output_to_medit(medit_file);
+
+  return 0;
+}
+