Merge pull request #2853 from MaelRL/Generator-Fix_missing_domain_check-GF

Generator: Fix `random_points_in_triangle_mesh_2`

# Conflicts:
#	Generator/include/CGAL/point_generators_2.h

Generator/examples/Generator/random_points_on_triangle_mesh_2.cpp

@@ -1,8 +1,11 @@
-#include <CGAL/Polygon_2.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
 #include <CGAL/Constrained_Delaunay_triangulation_2.h>
-#include <CGAL/Triangulation_face_base_with_info_2.h>
 #include <CGAL/Delaunay_mesh_face_base_2.h>
+#include <CGAL/Delaunay_mesh_size_criteria_2.h>
+#include <CGAL/Delaunay_mesher_2.h>
+#include <CGAL/Triangulation_face_base_with_info_2.h>
+#include <CGAL/Polygon_2.h>
 #include <CGAL/point_generators_2.h>
 
 #include <iostream>
@@ -13,15 +16,17 @@ typedef CGAL::Triangulation_vertex_base_2<K>                      Vb;
 typedef CGAL::Delaunay_mesh_face_base_2<K>                        Fb;
 typedef CGAL::Triangulation_data_structure_2<Vb, Fb>              Tds;
 typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds>        CDT;
+
 typedef CDT::Point                                                Point;
 typedef CGAL::Polygon_2<K>                                        Polygon_2;
-
+typedef CGAL::Delaunay_mesh_size_criteria_2<CDT>                  Mesh_2_criteria;
 
 using namespace CGAL;
 int main()
 {
  // Generated points are in that vector
   std::vector<Point> points;
+
   //Construct two non-intersecting nested polygons
   ::Polygon_2 polygon1;
   polygon1.push_back(Point(0,0));
@@ -38,9 +43,11 @@ int main()
   cdt.insert_constraint(polygon1.vertices_begin(), polygon1.vertices_end(), true);
   cdt.insert_constraint(polygon2.vertices_begin(), polygon2.vertices_end(), true);
 
+  // Refine the triangulation (and mark the faces as inside/outside)
+  CGAL::refine_Delaunay_mesh_2(cdt, Mesh_2_criteria(0.125, 0.5));
+
   // Create the generator, input is the Triangulation_2 cdt
-  Random_points_in_triangle_mesh_2<Point, CDT>
+  Random_points_in_triangle_mesh_2<Point, CDT> g(cdt);
-    g(cdt);
 
   // Get 100 random points in cdt
   CGAL::cpp11::copy_n(g, 100, std::back_inserter(points));

Generator/include/CGAL/internal/Generic_random_point_generator.h

@@ -22,12 +22,15 @@
 #ifndef CGAL_INTERNAL_GENERIC_RANDOM_POINT_GENERATOR_H
 #define CGAL_INTERNAL_GENERIC_RANDOM_POINT_GENERATOR_H
 
+#include <CGAL/assertions.h>
+#include <CGAL/Iterator_range.h>
 #include <CGAL/generators.h>
 #include <CGAL/Random.h>
 #include <CGAL/property_map.h>
-#include <vector>
+
 #include <boost/foreach.hpp>
-#include <CGAL/Iterator_range.h>
+
+#include <vector>
 
 namespace CGAL {
 
@@ -57,6 +60,8 @@ public:
     , random(rnd)
   {
     std::size_t input_size = input.size();
+    CGAL_precondition(input_size > 0);
+
     ids.reserve(input_size);
     weights.reserve(input_size);
 
@@ -71,9 +76,11 @@ public:
       total_weight += to_double( compute_weight(object) );
       weights.push_back(total_weight);
     }
+
     //generate the first point
     generate_point();
   }
+
   This& operator++()
   {
     generate_point();
@@ -85,6 +92,7 @@ public:
     ++(*this);
     return tmp;
   }
+
   double sum_of_weights() const
   {
     if (weights.empty())

Generator/include/CGAL/point_generators_2.h

@@ -33,11 +33,12 @@
 #include <CGAL/disable_warnings.h>
 
 #include <CGAL/generators.h>
-#include <iterator>
 #include <CGAL/number_type_basic.h>
 #include <CGAL/internal/Generic_random_point_generator.h>
 #include <CGAL/iterator.h>
 
+#include <iterator>
+
 namespace CGAL {
 
 template < class P, class Creator = 
@@ -551,6 +552,7 @@ void Random_points_in_triangle_2<P, Creator>::generate_point() {
 }
 
 namespace internal {
+
 //Functor returning Triangle_2 from Triangulation_2 Faces
 template <class T>
 class Triangle_from_face_2
@@ -560,22 +562,55 @@ public:
   typedef Triangle result_type;
   Triangle_from_face_2() {}
 
-  Triangle operator()(typename T::Face_handle face)const {
+  Triangle operator()(typename T::Finite_faces_iterator face) const {
     return Triangle(face->vertex(0)->point(), face->vertex(1)->point(), face->vertex(2)->point());
   }
 };
+
+struct Is_not_in_domain
+{
+  typedef bool                                   result_type;
+
+  template <class FH>
+  result_type operator()(const FH fh) const {
+    return (!fh->is_in_domain());
+  }
+};
+
+template <class T>
+class In_domain_finite_faces_iterator
+  : public Filter_iterator<typename T::Finite_faces_iterator, Is_not_in_domain>
+{
+  typedef CGAL::Filter_iterator<typename T::Finite_faces_iterator, Is_not_in_domain>     Base;
+  typedef In_domain_finite_faces_iterator<T>                                             Self;
+
+  typedef typename T::Face_handle                                                        Face_handle;
+  typedef typename T::Finite_faces_iterator                                              Finite_faces_iterator;
+
+public:
+  In_domain_finite_faces_iterator() : Base() {}
+  In_domain_finite_faces_iterator(const Base &b) : Base(b) {}
+  Self & operator++() { Base::operator++(); return *this; }
+  Self & operator--() { Base::operator--(); return *this; }
+  Self operator++(int) { Self tmp(*this); ++(*this); return tmp; }
+  Self operator--(int) { Self tmp(*this); --(*this); return tmp; }
+  operator Finite_faces_iterator() const { return Base::base(); }
+  operator Face_handle() const { return Face_handle(Base::base()); }
+};
+
 }//end namespace internal
+
 template <class P,
           class T,
-          class Creator =
+          class Creator = Creator_uniform_2<typename Kernel_traits<P>::Kernel::RT, P> >
-            Creator_uniform_2<typename Kernel_traits<P>::Kernel::RT,P>
+class Random_points_in_triangle_mesh_2
->
+  : public Generic_random_point_generator<internal::In_domain_finite_faces_iterator<T>,
-class Random_points_in_triangle_mesh_2 : public Generic_random_point_generator<
-  typename T::Face_handle ,
                                           internal::Triangle_from_face_2<T>,
-    Random_points_in_triangle_2<P> , P> {
+                                          Random_points_in_triangle_2<P, Creator>,
+                                          P>
+{
 public:
-  typedef Generic_random_point_generator<typename T::Face_handle,
+  typedef Generic_random_point_generator<internal::In_domain_finite_faces_iterator<T>,
                                          internal::Triangle_from_face_2<T>,
                                          Random_points_in_triangle_2<P, Creator>,
                                          P>                           Base;
@@ -583,15 +618,19 @@ public:
   typedef P                                                           result_type;
   typedef Random_points_in_triangle_mesh_2<P, T, Creator>             This;
 
-
   Random_points_in_triangle_mesh_2(const T& triangulation, Random& rnd = get_default_random())
-    : Base( CGAL::make_prevent_deref_range(triangulation.finite_faces_begin(),
+    : Base(CGAL::make_prevent_deref_range(
-                                           triangulation.finite_faces_end()),
+             CGAL::filter_iterator(triangulation.finite_faces_end(),
+                                   internal::Is_not_in_domain(),
+                                   triangulation.finite_faces_begin()),
+             CGAL::filter_iterator(triangulation.finite_faces_end(),
+                                   internal::Is_not_in_domain())),
            internal::Triangle_from_face_2<T>(),
            typename Kernel_traits<P>::Kernel::Compute_area_2(),
            rnd)
   {
   }
+
   This& operator++() {
     Base::generate_point();
     return *this;

Generator/test/Generator/generic_random_test.cpp

@@ -1,26 +1,30 @@
-#include <CGAL/Mesh_triangulation_3.h>
-#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
-#include <CGAL/Mesh_criteria_3.h>
-
-#include <CGAL/Polyhedral_mesh_domain_3.h>
-#include <CGAL/make_mesh_3.h>
-#include <CGAL/refine_mesh_3.h>
-
-
-
-#include <CGAL/Polygon_2.h>
-
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Simple_cartesian.h>
 
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+#include <CGAL/Mesh_polyhedron_3.h>
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
+#include <CGAL/make_mesh_3.h>
+#include <CGAL/refine_mesh_3.h>
+
+#include <CGAL/Polygon_2.h>
+
 #include <CGAL/Constrained_Delaunay_triangulation_2.h>
-#include <CGAL/Triangulation_face_base_with_info_2.h>
 #include <CGAL/Delaunay_mesh_face_base_2.h>
+#include <CGAL/Delaunay_mesh_size_criteria_2.h>
+#include <CGAL/Delaunay_mesher_2.h>
+#include <CGAL/Triangulation_data_structure_2.h>
+#include <CGAL/Triangulation_vertex_base_2.h>
+
 #include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
+#include <CGAL/boost/graph/helpers.h>
 #include <CGAL/point_generators_3.h>
 
-#include <CGAL/boost/graph/helpers.h>
 #include <iostream>
+#include <iterator>
+#include <vector>
 
 using namespace CGAL::parameters;
 
@@ -37,28 +41,27 @@ typedef CGAL::Triangulation_vertex_base_2<K>                      Vb;
 typedef CGAL::Delaunay_mesh_face_base_2<K>                        Fb;
 typedef CGAL::Triangulation_data_structure_2<Vb, Fb>              Tds;
 typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds>        CDT;
+typedef CGAL::Delaunay_mesh_size_criteria_2<CDT>                  Mesh_2_criteria;
 
 typedef CGAL::Polygon_2<K>                                        Polygon_2;
 
-typedef CGAL::Polyhedron_3<K>                                     Polyhedron;
+typedef CGAL::Mesh_polyhedron_3<K>::type                          Polyhedron;
 
 // Domain
-typedef CGAL::Polyhedral_mesh_domain_3<Polyhedron, K>             Mesh_domain;
+typedef CGAL::Polyhedral_mesh_domain_with_features_3<K>           Mesh_domain;
+typedef Mesh_domain::Corner_index                                 Corner_index;
+typedef Mesh_domain::Curve_segment_index                          Curve_segment_index;
 
-// Triangulation
+typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type             Tr;
-#ifdef CGAL_CONCURRENT_MESH_3
+typedef CGAL::Mesh_complex_3_in_triangulation_3<
-typedef CGAL::Parallel_tag                                        Concurrency_tag;
+          Tr, Corner_index, Curve_segment_index>                  C3t3;
-#else
-typedef CGAL::Sequential_tag                                      Concurrency_tag;
-#endif
-
-typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
-typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr>               C3t3;
 typedef CGAL::Mesh_criteria_3<Tr>                                 Mesh_criteria;
 typedef C3t3::Point                                               Point;
 
-int test_triangles_2()
+int test_triangles_2(const FT eps)
 {
+  std::cout << "test_triangles_2 (tolerance: " << eps << ")" << std::endl;
+
   // Generated points are in that vector
   std::vector<Point_2> points;
 
@@ -78,12 +81,11 @@ int test_triangles_2()
 
   BOOST_FOREACH(Point_2 p, points)
   {
-    bool on_quad = p.x() > -0.01 && p.x() < 0.51
+    bool on_quad = p.x() > -eps && p.x() < 0.5 + eps &&
-        && p.y() > -0.01 && p.y() < 0.51;
+                   p.y() > -eps && p.y() < 0.5 + eps;
     if(!on_quad)
     {
-      std::cerr<<p<<std::endl;
+      std::cerr << "ERROR : Generated point (" << p << ") is not on the square." << std::endl;
-      std::cerr<<"ERROR : Generated point is not on the triangle range."<<std::endl;
       return 0;
     }
   }
@@ -91,8 +93,10 @@ int test_triangles_2()
   return 1;
 }
 
-int test_triangles_3()
+int test_triangles_3(const FT eps)
 {
+  std::cout << "test_triangles_3 (tolerance: " << eps << ")" << std::endl;
+
   // Generated points are in that vector
   std::vector<Point_3> points;
 
@@ -114,18 +118,16 @@ int test_triangles_3()
 
   BOOST_FOREACH(Point_3 p, points)
   {
-    bool on_front = p.z() < 0.01 && p.z()> -0.01
+    bool on_front = p.z() <  eps && p.z() > -eps &&
-        && p.x() > -0.01 && p.x() < 0.51
+                    p.x() > -eps && p.x() < 0.5 + eps &&
-        && p.y() > -0.01 && p.y() < 0.51;
+                    p.y() > -eps && p.y() < 0.5 + eps;
-
+    bool on_right = p.x() < 0.5 + eps && p.x() > 0.5 - eps &&
-    bool on_right = p.x() < 0.51 && p.x()> 0.49
+                    p.z() > -eps && p.z() < 0.5 + eps &&
-        && p.z() > -0.01 && p.z() < 0.51
+                    p.y() > -eps && p.y() < 0.5 + eps;
-        && p.y() > -0.01 && p.y() < 0.51;
 
     if(!on_front && !on_right)
     {
-      std::cerr<<p<<std::endl;
+      std::cerr << "ERROR : Generated point (" << p << ") is not on a triangle of the range." << std::endl;
-      std::cerr<<"ERROR : Generated point is not on the triangle range."<<std::endl;
       return 0;
     }
   }
@@ -133,11 +135,12 @@ int test_triangles_3()
   return 1;
 }
 
-int test_T2()
+int test_T2(const FT eps)
 {
+  std::cout << "test_T2 (tolerance: " << eps << ")" << std::endl;
+
   std::vector<Point_2> points;
 
-  //construct two non-intersecting nested polygons
   Polygon_2 polygon1;
 
   polygon1.push_back(Point_2(0,0));
@@ -145,10 +148,19 @@ int test_T2()
   polygon1.push_back(Point_2(2,2));
   polygon1.push_back(Point_2(0,2));
 
-  //Insert the polygons into a constrained triangulation
+  // Insert the constraint...
   CDT cdt;
   cdt.insert_constraint(polygon1.vertices_begin(), polygon1.vertices_end(), true);
 
+  // ... and four points outside the box to create faces outside the domain
+  cdt.insert(Point_2(-1,-1));
+  cdt.insert(Point_2(3,-1));
+  cdt.insert(Point_2(3,3));
+  cdt.insert(Point_2(-1,3));
+
+  // Refine the triangulation
+  CGAL::refine_Delaunay_mesh_2(cdt, Mesh_2_criteria(0.125, 0.5));
+
   CGAL::Random_points_in_triangle_mesh_2<Point_2, CDT> g(cdt);
   CGAL::cpp11::copy_n(g, 300, std::back_inserter(points));
   for(std::size_t i=0; i<points.size(); ++i)
@@ -157,9 +169,9 @@ int test_T2()
     for(int j=0; j<2; ++j)
     {
       double coords[2] = {p.x(), p.y()};
-      if(coords[j]>2.05 || coords[j]<-0.05)
+      if(coords[j] > 2.0 + eps || coords[j] < -eps)
       {
-        std::cerr<<"ERROR : Generated point is not on the cube."<<std::endl;
+        std::cerr << "ERROR : Generated point (" << p << ") is not on a face of the domain." << std::endl;
         return 0;
       }
     }
@@ -168,17 +180,20 @@ int test_T2()
   return 1;
 }
 
-bool on_face(int face, double coord[3])
+bool on_face(int face, double coord[3], const FT eps)
 {
-    if(CGAL::abs(CGAL::abs(coord[face]) - 0.5) < 0.05
+  if(CGAL::abs(CGAL::abs(coord[face]) - 0.5) < eps &&
-       && CGAL::abs(coord[(face+1)%3]) - 0.5 < 0.05
+     CGAL::abs(coord[(face+1)%3]) - 0.5 < eps &&
-       && CGAL::abs(coord[(face+2)%3]) - 0.5 < 0.05)
+     CGAL::abs(coord[(face+2)%3]) - 0.5 < eps)
     return true;
+
   return false;
 }
 
-int test_volume_mesh(Polyhedron& polyhedron)
+int test_volume_mesh(Polyhedron& polyhedron, const FT eps)
 {
+  std::cout << "test_volume_mesh (tolerance: " << eps << ")" << std::endl;
+
   std::vector<Point_3> points;
   CGAL::Random_points_in_triangle_mesh_3<Polyhedron> g(polyhedron);
   CGAL::cpp11::copy_n(g, 300, std::back_inserter(points));
@@ -186,51 +201,65 @@ int test_volume_mesh(Polyhedron& polyhedron)
   {
     Point_3 p = points[i];
     double coords[3] = {p.x(), p.y(), p.z()};
-    if(!(on_face(0, coords) || on_face(1, coords) || on_face(2,coords)))
+    if(!(on_face(0, coords, eps) || on_face(1, coords, eps) || on_face(2, coords, eps)))
     {
-      std::cerr<<"ERROR : Generated point is not on the cube."<<std::endl;
+      std::cerr << "ERROR : Generated point (" << p << ") is not on a face." << std::endl;
       return 0;
     }
   }
+
   return 1;
 }
 
-int test_on_c3t3(const Polyhedron& polyhedron)
+int test_on_c3t3(const Polyhedron& polyhedron, const FT eps)
 {
+  std::cout << "test_on_c3t3 (tolerance: " << eps << ")" << std::endl;
+
   std::vector<Point_3> points;
 
   // Create domain
   Mesh_domain domain(polyhedron);
+  domain.detect_features();
 
   // Mesh criteria (no cell_size set)
-  Mesh_criteria criteria(facet_angle=25, facet_size=0.15, facet_distance=0.008,
+  Mesh_criteria criteria(facet_angle=25,
+                         facet_size=0.05,
+                         facet_distance=0.008,
                          cell_radius_edge_ratio=3);
+
   // Mesh generation
   C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb(), no_exude());
+
   CGAL::Random_points_in_tetrahedral_mesh_boundary_3<C3t3> g(c3t3);
   CGAL::cpp11::copy_n(g, 300, std::back_inserter(points));
   for(std::size_t i=0; i<points.size(); ++i)
   {
     Point_3 p = points[i];
     double coords[3] = {p.x(), p.y(), p.z()};
-    if(!(on_face(0, coords) || on_face(1, coords) || on_face(2,coords)))
+    if(!(on_face(0, coords, eps) || on_face(1, coords, eps) || on_face(2, coords, eps)))
     {
-      std::cerr<<"ERROR : Generated point is not on the cube."<<std::endl;
+      std::cerr << "ERROR : Generated point (" << p << ") is not on a face." << std::endl;
       return 0;
     }
   }
+
   return 1;
 }
 
-int test_in_c3t3(const Polyhedron& polyhedron)
+int test_in_c3t3(const Polyhedron& polyhedron, const FT eps)
 {
+  std::cout << "test_in_c3t3 (tolerance: " << eps << ")" << std::endl;
+
   std::vector<Point_3> points;
 
   // Create domain
   Mesh_domain domain(polyhedron);
+  domain.detect_features();
 
   // Mesh criteria (no cell_size set)
-  Mesh_criteria criteria(facet_angle=25, facet_size=0.15, facet_distance=0.008,
+  Mesh_criteria criteria(facet_angle=25,
+                         facet_size=0.05,
+                         facet_distance=0.008,
                          cell_radius_edge_ratio=3);
 
   // Mesh generation
@@ -243,12 +272,15 @@ int test_in_c3t3(const Polyhedron& polyhedron)
     Point_3 p = points[i];
     double coords[3] = {p.x(), p.y(), p.z()};
     for(int j=0; j<3; ++j)
-      if(CGAL::abs(coords[j]) > 0.501)
     {
-        std::cerr<<"ERROR : Generated point is not in the cube."<<std::endl;
+      if(CGAL::abs(coords[j]) > 0.5 + eps)
+      {
+        std::cerr << "ERROR : Generated point (" << p << ") is not on the cube." << std::endl;
         return 0;
       }
     }
+  }
+
   return 1;
 }
 
@@ -256,9 +288,11 @@ int main()
 {
   Polyhedron polyhedron;
 
+  // A cube
   make_hexahedron(Point_3(-0.5,-0.5,-0.5), Point_3(0.5,-0.5,-0.5), Point_3(0.5,0.5,-0.5), Point_3(-0.5,0.5,-0.5),
                   Point_3(-0.5,0.5,0.5), Point_3(-0.5,-0.5,0.5), Point_3(0.5,-0.5,0.5), Point_3(0.5,0.5,0.5),
                   polyhedron);
+
   boost::graph_traits<Polyhedron>::halfedge_descriptor facets[6];
   int i = 0;
   BOOST_FOREACH(boost::graph_traits<Polyhedron>::face_descriptor fd, faces(polyhedron))
@@ -267,14 +301,17 @@ int main()
   for(int i=0; i<6; ++i)
     CGAL::Euler::split_face(facets[i],next(next(facets[i], polyhedron), polyhedron), polyhedron);
 
+  const FT eps = 1e-10;
+
   int validity =
-      test_triangles_2()
+      test_triangles_2(eps)
-      *test_triangles_3()
+      *test_triangles_3(eps)
-      *test_volume_mesh(polyhedron)
+      *test_volume_mesh(polyhedron, eps)
-      *test_T2()
+      *test_T2(eps)
-      *test_on_c3t3(polyhedron)
+      *test_on_c3t3(polyhedron, eps)
-      *test_in_c3t3(polyhedron)
+      *test_in_c3t3(polyhedron, eps)
       ;
   assert(validity == 1);
+
   return 0;
 }