Merge remote-tracking branch 'cgal/master' into CGAL-move_semantic_for_triangulations-GF

Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_fct.cpp

@@ -60,6 +60,8 @@ int main(int argc, char* argv[])
 {
   std::ifstream in((argc>1)?argv[1]:"data/half.xyz");
   double per = (argc>2)?boost::lexical_cast<double>(argv[2]):0;
+  double radius_ratio_bound = (argc>3)?boost::lexical_cast<double>(argv[3]):5.0;
+
   std::vector<Point_3> points;
   std::vector<Facet> facets;
 
@@ -71,7 +73,8 @@ int main(int argc, char* argv[])
   CGAL::advancing_front_surface_reconstruction(points.begin(),
                                                points.end(),
                                                std::back_inserter(facets),
-                                               perimeter);
+                                               perimeter,
+                                               radius_ratio_bound);
 
   std::cout << "OFF\n" << points.size() << " " << facets.size() << " 0\n";
   std::copy(points.begin(),

BGL/examples/BGL_LCC/copy_lcc.cpp

@@ -55,6 +55,7 @@ int main(int argc, char* argv[])
     CGAL::copy_face_graph(T1, S, CGAL::parameters::vertex_to_vertex_output_iterator(std::inserter(v2v, v2v.end()))
                                                   .halfedge_to_halfedge_output_iterator(std::inserter(h2h, h2h.end())));
     std::ofstream out("reverse.off");
+    out.precision(17);
     out << S;
   }
   return 0;

BGL/examples/BGL_polyhedron_3/copy_polyhedron.cpp

@@ -43,6 +43,7 @@ int main(int argc, char* argv[])
   {
     CGAL::copy_face_graph(S, T1);
     std::ofstream out("sm.off");
+    out.precision(17);
     out << T1;
   }
 
@@ -81,6 +82,7 @@ int main(int argc, char* argv[])
     boost::unordered_map<source_face_descriptor, tm_face_descriptor> f2f;
     CGAL::copy_face_graph(T1, S, std::inserter(v2v, v2v.end()), std::inserter(h2h, h2h.end()));
     std::ofstream out("reverse.off");
+    out.precision(17);
     out << S;
     CGAL::copy_face_graph(T1, S, CGAL::parameters::vertex_to_vertex_map(boost::make_assoc_property_map(v2v))
                           .halfedge_to_halfedge_output_iterator(std::inserter(h2h, h2h.end()))

BGL/examples/BGL_surface_mesh/write_inp.cpp

@@ -20,6 +20,7 @@ int main()
   CGAL::make_quad(Point(0,0,1), Point(1,0,1),Point(1,1,1),Point(0,1,1), sm);
 
   std::ofstream out("out.inp");
+  out.precision(17);
   CGAL::write_inp(out, sm, "out.inp", "S4R");
   return 0;
 }

Barycentric_coordinates_2/benchmark/Barycentric_coordinates_2/Discrete_harmonic_coordinates/DH_speed_test_100_vertices.cpp

@@ -2,6 +2,7 @@
 // We test speed of discrete harmonic coordinates on a set of automatically generated
 // points inside a regular polygon with 100 vertices. We use inexact kernel.
 
+#include <CGAL/number_type_config.h>
 #include <CGAL/Real_timer.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Barycentric_coordinates_2/Discrete_harmonic_2.h>
@@ -28,12 +29,11 @@ void generate_regular_polygon(const int number_of_vertices, const double polygon
 {
     const int n = number_of_vertices;
     const double r =  polygon_radius;
-    const double number_pi = 3.14159;
 
     vertices.resize(n);
 
     for(int i = 0; i < n; ++i)    
-        vertices[i] = Point(Scalar(r*sin((number_pi / n) + ((i * 2.0 * number_pi) / n))), Scalar(-r*cos((number_pi / n) + ((i * 2.0 * number_pi) / n))));
+        vertices[i] = Point(Scalar(r*sin((CGAL_PI / n) + ((i * 2.0 * CGAL_PI) / n))), Scalar(-r*cos((CGAL_PI / n) + ((i * 2.0 * CGAL_PI) / n))));
 }
 
 int main()

Barycentric_coordinates_2/benchmark/Barycentric_coordinates_2/Mean_value_coordinates/MV_speed_test_100_vertices.cpp

@@ -2,6 +2,7 @@
 // We test speed of mean value coordinates on a set of automatically generated
 // points inside a regular polygon with 100 vertices. We use inexact kernel.
 
+#include <CGAL/number_type_config.h>
 #include <CGAL/Real_timer.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Barycentric_coordinates_2/Mean_value_2.h>
@@ -28,12 +29,11 @@ void generate_regular_polygon(const int number_of_vertices, const double polygon
 {
     const int n = number_of_vertices;
     const double r =  polygon_radius;
-    const double number_pi = 3.14159;
 
     vertices.resize(n);
 
     for(int i = 0; i < n; ++i)    
-        vertices[i] = Point(Scalar(r*sin((number_pi / n) + ((i * 2.0 * number_pi) / n))), Scalar(-r*cos((number_pi / n) + ((i * 2.0 * number_pi) / n))));
+        vertices[i] = Point(Scalar(r*sin((CGAL_PI / n) + ((i * 2.0 * CGAL_PI) / n))), Scalar(-r*cos((CGAL_PI / n) + ((i * 2.0 * CGAL_PI) / n))));
 }
 
 int main()

Barycentric_coordinates_2/benchmark/Barycentric_coordinates_2/Wachspress_coordinates/WP_speed_test_100_vertices.cpp

@@ -2,6 +2,7 @@
 // We test speed of Wachspress coordinates on a set of automatically generated
 // points inside a regular polygon with 100 vertices. We use inexact kernel.
 
+#include <CGAL/number_type_config.h>
 #include <CGAL/Real_timer.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Barycentric_coordinates_2/Wachspress_2.h>
@@ -28,12 +29,11 @@ void generate_regular_polygon(const int number_of_vertices, const double polygon
 {
     const int n = number_of_vertices;
     const double r =  polygon_radius;
-    const double number_pi = 3.14159;
 
     vertices.resize(n);
 
     for(int i = 0; i < n; ++i)    
-        vertices[i] = Point(Scalar(r*sin((number_pi / n) + ((i * 2.0 * number_pi) / n))), Scalar(-r*cos((number_pi / n) + ((i * 2.0 * number_pi) / n))));
+        vertices[i] = Point(Scalar(r*sin((CGAL_PI / n) + ((i * 2.0 * CGAL_PI) / n))), Scalar(-r*cos((CGAL_PI / n) + ((i * 2.0 * CGAL_PI) / n))));
 }
 
 int main()

Box_intersection_d/doc/Box_intersection_d/CGAL/box_intersection_d.h

@@ -242,6 +242,29 @@ void box_intersection_all_pairs_d(
   cutoff parameters are recommended. See also 
   Section \ref secboxintersperformance . 
 
+\cgalHeading{Concurrency}
+
+  The first template parameter of the function enables to choose whether
+  the algorithm is to be run in parallel, if `CGAL::Parallel_tag` is specified
+  and %CGAL has been linked with the Intel TBB library, or sequentially,
+  if `CGAL::Sequential_tag` - the default value - is specified.
+  The parallelization of the algorithm is based on a divide-and-conquer
+  approach: the two ranges are split in a number of smaller ranges, and
+  all combinations of subranges are treated in parallel. It is thus
+  recommended to use ranges of pointers to bounding boxes, to keep
+  these copies light.
+
+  \warning The parallel mode comes with a small overhead due to the
+  duplication and splitting of the input ranges. It is an improvement
+  for almost all inputs, but not all. A configuration where the two ranges
+  are small and entirely disjoint might result in a slightly worse
+  runtime when using the parallel version. Users should benchmark both
+  versions to verify that their data does not fall in this (small)
+  set of inputs.
+
+  \warning When using the parallel mode, the callback function must
+  be threadsafe.
+
 \cgalHeading{Example}
 
   The box implementation provided with 
@@ -274,7 +297,8 @@ void box_intersection_all_pairs_d(
 
 
 */
-template< class RandomAccessIterator1, 
+template< class ConcurrencyTag = CGAL::Sequential_tag,
+          class RandomAccessIterator1,
           class RandomAccessIterator2, 
           class Callback >
 void box_intersection_d(
@@ -291,7 +315,8 @@ void box_intersection_d(
   Invocation with custom box traits.
 
 */
-template< class RandomAccessIterator1,
+template< class ConcurrencyTag = CGAL::Sequential_tag,
+          class RandomAccessIterator1,
           class RandomAccessIterator2,
           class Callback, class BoxTraits >
 void box_intersection_d(
@@ -489,8 +514,11 @@ namespace CGAL {
 
 \cgalHeading{Implementation}
 
-  See the implementation section of the `box_intersection_d()` 
+  See the implementation section of the `box_intersection_d()` function.
-  function.
+
+\cgalHeading{Concurrency}
+
+  See the concurrency section of the `box_intersection_d()` function.
 
 \cgalHeading{Example}
 
@@ -520,7 +548,8 @@ namespace CGAL {
   `RandomAccessIterator`.
 
 */
-template< class RandomAccessIterator, class Callback >
+template< class ConcurrencyTag = CGAL::Sequential_tag,
+          class RandomAccessIterator, class Callback >
 void box_self_intersection_d(
   RandomAccessIterator begin, RandomAccessIterator end,
   Callback callback,
@@ -532,7 +561,8 @@ void box_self_intersection_d(
 
   Invocation with custom box traits.
 */
-template< class RandomAccessIterator,
+template< class ConcurrencyTag = CGAL::Sequential_tag
+          class RandomAccessIterator,
           class Callback, class BoxTraits >
 void box_self_intersection_d(
   RandomAccessIterator begin, RandomAccessIterator end,

Box_intersection_d/include/CGAL/Box_intersection_d/Box_with_info_d.h

@@ -16,34 +16,34 @@
 
 #include <CGAL/license/Box_intersection_d.h>
 
-
 #include <CGAL/basic.h>
 #include <CGAL/Box_intersection_d/Box_d.h>
 
 namespace CGAL {
-
 namespace Box_intersection_d {
 
-template<class NT_, int N, class Info_> 
+template<class NT_, int N, class Info_, class IdPolicy = ID_EXPLICIT>
-class Box_with_info_d : public Box_d< NT_, N, ID_FROM_BOX_ADDRESS> {
+class Box_with_info_d
+  : public Box_d< NT_, N, IdPolicy>
+{
 protected:
     Info_ m_info;
 public:
-    typedef Box_d< NT_, N, ID_FROM_BOX_ADDRESS> Base;
+    typedef Box_d< NT_, N, IdPolicy>    Base;
     typedef NT_                         NT;
     typedef Info_                       Info;
+
     Box_with_info_d() {}
     Box_with_info_d( Info h) : m_info(h) {}
     Box_with_info_d( bool complete, Info h): Base(complete), m_info(h) {}
     Box_with_info_d(NT l[N], NT h[N], Info n) : Base( l, h), m_info(n) {}
     Box_with_info_d( const Bbox_2& b, Info h) : Base( b), m_info(h) {}
     Box_with_info_d( const Bbox_3& b, Info h) : Base( b), m_info(h) {}
+
     Info info() const { return m_info; }
 };
 
-} // end namespace Box_intersection_d
+} // namespace Box_intersection_d
+} // namespace CGAL
 
-
+#endif // CGAL_BOX_INTERSECTION_D_BOX_WITH_INFO_D_H
-} //namespace CGAL
-
-#endif

Box_intersection_d/include/CGAL/Box_intersection_d/segment_tree.h

@@ -16,7 +16,6 @@
 
 #include <CGAL/license/Box_intersection_d.h>
 
-
 #include <CGAL/basic.h>
 #include <CGAL/Box_intersection_d/box_limits.h>
 
@@ -24,7 +23,6 @@
 #include <boost/random/uniform_int.hpp>
 #include <boost/random/variate_generator.hpp>
 
-
 #include <algorithm>
 #include <iterator>
 #include <functional>
@@ -93,6 +91,8 @@ void one_way_scan( RandomAccessIter1 p_begin, RandomAccessIter1 p_end,
                    bool in_order = true )
 {
     typedef typename Traits::Compare Compare;
+
+    // Putting a parallel sort here slows down the overall parallel algorithm
     std::sort( p_begin, p_end, Compare( 0 ) );
     std::sort( i_begin, i_end, Compare( 0 ) );
 
@@ -323,7 +323,7 @@ void segment_tree( RandomAccessIter1 p_begin, RandomAccessIter1 p_end,
                    RandomAccessIter2 i_begin, RandomAccessIter2 i_end,
                    T lo, T hi,
                    Callback callback, Predicate_traits traits,
-                   std::ptrdiff_t cutoff, int dim, bool in_order )
+                   std::ptrdiff_t cutoff, int dim, bool in_order)
 {
     typedef typename Predicate_traits::Spanning   Spanning;
     typedef typename Predicate_traits::Lo_less    Lo_less;

Box_intersection_d/include/CGAL/box_intersection_d.h

@@ -25,12 +25,191 @@
 #include <CGAL/Box_intersection_d/Box_traits_d.h>
 #include <CGAL/Box_intersection_d/box_limits.h>
 
+#include <CGAL/use.h>
+#include <CGAL/tags.h>
+
+#ifdef CGAL_LINKED_WITH_TBB
+#include <tbb/task_group.h>
+#endif
+
+#include <iterator>
 #include <vector>
 
+////////////////////////////////////////////////////////////////////////////////////////////////
+/// THE CALLBACK MUST BE THREADSAFE IF YOU ARE USING THE PARALLEL MODE
+////////////////////////////////////////////////////////////////////////////////////////////////
+
 namespace CGAL {
+namespace internal {
 
 // Generic call with custom predicate traits parameter.
-template< class RandomAccessIter1, class RandomAccessIter2,
+template< class ConcurrencyTag,
+          class RandomAccessIter1, class RandomAccessIter2,
+          class Callback, class Traits >
+void box_intersection_segment_tree_d(
+    RandomAccessIter1 begin1, RandomAccessIter1 end1,
+    RandomAccessIter2 begin2, RandomAccessIter2 end2,
+    Callback callback,
+    const Traits& traits,
+    const std::ptrdiff_t cutoff,
+    const bool in_order)
+{
+  typedef typename Traits::NT NT;
+
+  CGAL_assertion(Traits::dimension() > 0);
+  const int dim = Traits::dimension() - 1;
+
+  const NT inf = Box_intersection_d::box_limits<NT>::inf();
+  const NT sup = Box_intersection_d::box_limits<NT>::sup();
+
+#ifndef CGAL_LINKED_WITH_TBB
+  CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
+                             "Parallel_tag is enabled but TBB is unavailable.");
+#else // CGAL_LINKED_WITH_TBB
+  if(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value)
+  {
+    // Here is an illustration for n=2.
+    //
+    // Doing a R1-R2 intersection with a 2-split is doing 4 subpairs in parallel
+    // a-c and b-d
+    // r1[0][0] --a-- r1[0][1] --b-- r1[0][2]
+    // r2[0][0] --c-- r2[0][1] --d-- r2[0][2]
+    //
+    // a-d and b-c
+    // r1[1][0] --a-- r1[1][1] --b-- r1[1][2]
+    // r2[1][0] --c-- r2[1][1] --d-- r2[1][2]
+    //
+    // Ranges must be duplicates since sorting is performed
+
+    typedef typename std::iterator_traits<RandomAccessIter1>::value_type         val_t;
+    typedef typename std::iterator_traits<RandomAccessIter1>::difference_type    diff_size;
+
+    typedef std::vector<val_t>                                                   val_container;
+    typedef typename val_container::iterator                                     It;
+
+    static constexpr int n = 4;
+
+    const diff_size r1s = std::distance(begin1, end1);
+    const diff_size r2s = std::distance(begin2, end2);
+
+    val_container range_1_copies, range_2_copies;
+    range_1_copies.reserve(r1s * n);
+    range_2_copies.reserve(r2s * n);
+
+    const diff_size r1_step = r1s / n;
+    const diff_size r2_step = r2s / n;
+
+    for(int i=0; i<n; ++i)
+    {
+      range_1_copies.insert(range_1_copies.end(), begin1, end1);
+      range_2_copies.insert(range_2_copies.end(), begin2, end2);
+    }
+
+    // for example for n=2, there's 'begin', 'mid', and 'end' but we leave out 'end' for convenience
+    std::array<std::array<It, n>, n> range_1_iterators;
+    std::array<std::array<It, n>, n> range_2_iterators;
+
+    for(int i=0; i<n; ++i)
+    {
+      for(int j=0; j<n; ++j)
+      {
+        range_1_iterators[i][j] = range_1_copies.begin() + i * r1s + j * r1_step;
+        range_2_iterators[i][j] = range_2_copies.begin() + i * r2s + j * r2_step;
+      }
+    }
+
+    tbb::task_group g;
+
+    for(int i=0; i<n; ++i)
+    {
+      for(int j=0; j<n; ++j)
+      {
+        // 'j' vs 'j+i', meaning for each value of 'i' we're matching r1[i] with a shifted version of r2[i]
+        int r1_endi, r1_endj;
+        int r2_endi, r2_endj;
+
+        if(i == (n-1))
+        {
+          if(j == (n-1)) // both very last range and very last iterator
+          {
+            r1_endi = -1;
+            r1_endj = -1;
+          }
+          else
+          {
+            r1_endi = i;
+            r1_endj = j+1;
+          }
+        }
+        else // i != (n-1)
+        {
+          if(j == (n-1)) // end of that range, but not of the full container
+          {
+            r1_endi = i+1;
+            r1_endj = 0;
+          }
+          else
+          {
+            r1_endi = i;
+            r1_endj = j+1;
+          }
+        }
+
+        if(i == (n-1))
+        {
+          if(j+i == (n-1)) // both very last range and very last iterator
+          {
+            r2_endi = -1;
+            r2_endj = -1;
+          }
+          else
+          {
+            r2_endi = i;
+            r2_endj = (j+i+1)%n;
+          }
+        }
+        else // i != (n-1)
+        {
+          if(j+i == (n-1)) // end of that range, but not of the full container
+          {
+            r2_endi = i+1;
+            r2_endj = 0;
+          }
+          else
+          {
+            r2_endi = i;
+            r2_endj = (j+i+1)%n;
+          }
+        }
+
+        It r1_start = range_1_iterators[i][j];
+        It r1_end = (r1_endi == -1) ? range_1_copies.end() : range_1_iterators[r1_endi][r1_endj];
+        It r2_start = range_2_iterators[i][(j+i)%n];
+        It r2_end = (r2_endi == -1) ? range_2_copies.end() : range_2_iterators[r2_endi][r2_endj];
+        CGAL_assertion(range_1_copies.begin() <= r1_start && r1_start <= r1_end && r1_end <= range_1_copies.end());
+        CGAL_assertion(range_2_copies.begin() <= r2_start && r2_start <= r2_end && r2_end <= range_2_copies.end());
+
+        // Specify "copy by value" otherwise the values of iterators for next (i,j) iterations
+        // become shared with different lambdas being run in parallel, and things go wrong
+        g.run([=]{ Box_intersection_d::segment_tree( r1_start, r1_end, r2_start, r2_end,
+                                                     inf, sup, callback, traits, cutoff, dim, in_order); });
+      }
+    }
+
+    g.wait();
+  }
+  else
+#endif // CGAL_LINKED_WITH_TBB
+  {
+    Box_intersection_d::segment_tree(begin1, end1, begin2, end2, inf, sup, callback, traits, cutoff, dim, in_order);
+  }
+}
+
+} // namespace internal
+
+// Generic call with custom predicate traits parameter.
+template< class ConcurrencyTag = Sequential_tag,
+          class RandomAccessIter1, class RandomAccessIter2,
           class Callback, class BoxPredicateTraits >
 void box_intersection_custom_predicates_d(
     RandomAccessIter1 begin1, RandomAccessIter1 end1,
@@ -40,23 +219,15 @@ void box_intersection_custom_predicates_d(
     std::ptrdiff_t cutoff = 10,
     Box_intersection_d::Setting setting = Box_intersection_d::BIPARTITE)
 {
-    typedef BoxPredicateTraits Traits;
+  internal::box_intersection_segment_tree_d<ConcurrencyTag>(begin1, end1, begin2, end2, callback, traits, cutoff, true);
-    typedef typename Traits::NT NT;
-    CGAL_assertion( Traits::dimension() > 0 );
-    const int dim = Traits::dimension() - 1;
-    const NT inf = Box_intersection_d::box_limits<NT>::inf();
-    const NT sup = Box_intersection_d::box_limits<NT>::sup();
-    Box_intersection_d::segment_tree(begin1, end1, begin2, end2,
-                              inf, sup, callback, traits, cutoff, dim, true);
   if(setting == Box_intersection_d::BIPARTITE)
-        Box_intersection_d::segment_tree(begin2, end2, begin1, end1,
+    internal::box_intersection_segment_tree_d<ConcurrencyTag>(begin2, end2, begin1, end1, callback, traits, cutoff, false);
-                              inf, sup, callback, traits, cutoff, dim, false);
 }
 
-
 // Generic call with box traits parameter.
 // - make all default parameters explicit overloads (workaround)
-template< class RandomAccessIter1, class RandomAccessIter2,
+template< class ConcurrencyTag = Sequential_tag,
+          class RandomAccessIter1, class RandomAccessIter2,
           class Callback, class BoxTraits >
 void box_intersection_d(
     RandomAccessIter1 begin1, RandomAccessIter1 end1,
@@ -69,16 +240,17 @@ void box_intersection_d(
 {
   if (topology == Box_intersection_d::CLOSED) {
     typedef Box_intersection_d::Predicate_traits_d<BoxTraits,true> Traits;
-        box_intersection_custom_predicates_d(begin1, end1, begin2, end2,
+    box_intersection_custom_predicates_d<ConcurrencyTag>(begin1, end1, begin2, end2,
                                                          callback, Traits(), cutoff, setting);
   } else {
     typedef Box_intersection_d::Predicate_traits_d<BoxTraits,false> Traits;
-        box_intersection_custom_predicates_d(begin1, end1, begin2, end2,
+    box_intersection_custom_predicates_d<ConcurrencyTag>(begin1, end1, begin2, end2,
                                                          callback, Traits(), cutoff, setting);
   }
 }
 
-template< class RandomAccessIter1, class RandomAccessIter2,
+template< class ConcurrencyTag = Sequential_tag,
+          class RandomAccessIter1, class RandomAccessIter2,
           class Callback, class BoxTraits >
 void box_intersection_d(
     RandomAccessIter1 begin1, RandomAccessIter1 end1,
@@ -86,35 +258,39 @@ void box_intersection_d(
     Callback callback, BoxTraits box_traits, std::ptrdiff_t cutoff,
     Box_intersection_d::Topology topology)
 {
-    box_intersection_d( begin1, end1, begin2, end2, callback, box_traits, 
+    box_intersection_d<ConcurrencyTag>( begin1, end1, begin2, end2, callback, box_traits,
                                         cutoff, topology, Box_intersection_d::BIPARTITE);
 }
-template< class RandomAccessIter1, class RandomAccessIter2,
+template< class ConcurrencyTag = Sequential_tag,
+          class RandomAccessIter1, class RandomAccessIter2,
           class Callback, class BoxTraits >
 void box_intersection_d(
     RandomAccessIter1 begin1, RandomAccessIter1 end1,
     RandomAccessIter2 begin2, RandomAccessIter2 end2,
     Callback callback, BoxTraits box_traits, std::ptrdiff_t cutoff)
 {
-    box_intersection_d( begin1, end1, begin2, end2, callback, box_traits, 
+  box_intersection_d<ConcurrencyTag>( begin1, end1, begin2, end2, callback, box_traits,
                                       cutoff, Box_intersection_d::CLOSED,
                                       Box_intersection_d::BIPARTITE);
 }
-template< class RandomAccessIter1, class RandomAccessIter2,
+
+template< class ConcurrencyTag = Sequential_tag,
+          class RandomAccessIter1, class RandomAccessIter2,
           class Callback, class BoxTraits >
 void box_intersection_d(
     RandomAccessIter1 begin1, RandomAccessIter1 end1,
     RandomAccessIter2 begin2, RandomAccessIter2 end2,
     Callback callback, BoxTraits box_traits)
 {
-    box_intersection_d( begin1, end1, begin2, end2, callback, box_traits, 
+  box_intersection_d<ConcurrencyTag>( begin1, end1, begin2, end2, callback, box_traits,
                                       10, Box_intersection_d::CLOSED,
                                       Box_intersection_d::BIPARTITE);
 }
 
 // Specialized call with default box traits.
 // - make all default parameters explicit overloads (workaround)
-template< class RandomAccessIter1, class RandomAccessIter2, class Callback >
+template< class ConcurrencyTag = Sequential_tag,
+          class RandomAccessIter1, class RandomAccessIter2, class Callback >
 void box_intersection_d(
     RandomAccessIter1 begin1, RandomAccessIter1 end1,
     RandomAccessIter2 begin2, RandomAccessIter2 end2,
@@ -124,11 +300,13 @@ void box_intersection_d(
 {
   typedef typename std::iterator_traits<RandomAccessIter1>::value_type val_t;
   typedef Box_intersection_d::Box_traits_d< val_t>                     Box_traits;
-    box_intersection_d( begin1, end1, begin2, end2, callback, Box_traits(), 
+
+  box_intersection_d<ConcurrencyTag>( begin1, end1, begin2, end2, callback, Box_traits(),
                                       cutoff, topology, setting);
 }
 
-template< class RandomAccessIter1, class RandomAccessIter2, class Callback >
+template< class ConcurrencyTag = Sequential_tag,
+          class RandomAccessIter1, class RandomAccessIter2, class Callback >
 void box_intersection_d(
     RandomAccessIter1 begin1, RandomAccessIter1 end1,
     RandomAccessIter2 begin2, RandomAccessIter2 end2,
@@ -137,10 +315,12 @@ void box_intersection_d(
 {
   typedef typename std::iterator_traits<RandomAccessIter1>::value_type val_t;
   typedef Box_intersection_d::Box_traits_d< val_t>                     Box_traits;
-    box_intersection_d( begin1, end1, begin2, end2, callback, Box_traits(), 
+
+  box_intersection_d<ConcurrencyTag>( begin1, end1, begin2, end2, callback, Box_traits(),
                                       cutoff, topology, Box_intersection_d::BIPARTITE);
 }
-template< class RandomAccessIter1, class RandomAccessIter2, class Callback >
+template< class ConcurrencyTag = Sequential_tag,
+          class RandomAccessIter1, class RandomAccessIter2, class Callback >
 void box_intersection_d(
     RandomAccessIter1 begin1, RandomAccessIter1 end1,
     RandomAccessIter2 begin2, RandomAccessIter2 end2,
@@ -148,11 +328,12 @@ void box_intersection_d(
 {
   typedef typename std::iterator_traits<RandomAccessIter1>::value_type val_t;
   typedef Box_intersection_d::Box_traits_d< val_t>  Box_traits;
-    box_intersection_d( begin1, end1, begin2, end2, callback, Box_traits(), 
+  box_intersection_d<ConcurrencyTag>( begin1, end1, begin2, end2, callback, Box_traits(),
                                       cutoff, Box_intersection_d::CLOSED,
                                       Box_intersection_d::BIPARTITE);
 }
-template< class RandomAccessIter1, class RandomAccessIter2, class Callback >
+template< class ConcurrencyTag = Sequential_tag,
+          class RandomAccessIter1, class RandomAccessIter2, class Callback >
 void box_intersection_d(
     RandomAccessIter1 begin1, RandomAccessIter1 end1,
     RandomAccessIter2 begin2, RandomAccessIter2 end2,
@@ -160,15 +341,15 @@ void box_intersection_d(
 {
   typedef typename std::iterator_traits<RandomAccessIter1>::value_type val_t;
   typedef Box_intersection_d::Box_traits_d< val_t>  Box_traits;
-    box_intersection_d( begin1, end1, begin2, end2, callback, Box_traits(), 
+  box_intersection_d<ConcurrencyTag>( begin1, end1, begin2, end2, callback, Box_traits(),
                                       10, Box_intersection_d::CLOSED,
                                       Box_intersection_d::BIPARTITE);
 }
 
-
 // Generic call with box traits parameter, specialized for self-intersection.
 // - make all default parameters explicit overloads (workaround)
-template< class RandomAccessIter, class Callback, class BoxTraits >
+template< class ConcurrencyTag = Sequential_tag,
+          class RandomAccessIter, class Callback, class BoxTraits >
 void box_self_intersection_d(
     RandomAccessIter begin, RandomAccessIter end,
     Callback callback,
@@ -181,44 +362,47 @@ void box_self_intersection_d(
   // would be messed up by sorts on the second range otherwise.
   typedef typename std::iterator_traits<RandomAccessIter>::value_type val_t;
   std::vector< val_t> i( begin, end);
-    box_intersection_d( begin, end, i.begin(), i.end(),
+
+  box_intersection_d<ConcurrencyTag>( begin, end, i.begin(), i.end(),
                                       callback, box_traits, cutoff, topology,
                                       Box_intersection_d::COMPLETE);
 }
 
-template< class RandomAccessIter, class Callback, class BoxTraits >
+template< class ConcurrencyTag = Sequential_tag,
+          class RandomAccessIter, class Callback, class BoxTraits >
 void box_self_intersection_d(
     RandomAccessIter begin, RandomAccessIter end,
     Callback callback,
     BoxTraits box_traits,
     std::ptrdiff_t cutoff)
 {
-    return box_self_intersection_d(begin, end, callback, box_traits, cutoff,
+    return box_self_intersection_d<ConcurrencyTag>(begin, end, callback, box_traits, cutoff,
                                                    Box_intersection_d::CLOSED);
 }
 
-template< class RandomAccessIter, class Callback, class BoxTraits >
+template< class ConcurrencyTag = Sequential_tag,
+          class RandomAccessIter, class Callback, class BoxTraits >
 void box_self_intersection_d(
     RandomAccessIter begin, RandomAccessIter end,
     Callback callback,
     BoxTraits box_traits)
 {
-    return box_self_intersection_d(begin, end, callback, box_traits, 10);
+  return box_self_intersection_d<ConcurrencyTag>(begin, end, callback, box_traits, 10);
 }
 
 // Specialized call with default box traits, specialized for self-intersection.
 // - make all default parameters explicit overloads (workaround)
-template< class RandomAccessIter, class Callback >
+template< class ConcurrencyTag = Sequential_tag, class RandomAccessIter, class Callback >
 void box_self_intersection_d(
     RandomAccessIter begin, RandomAccessIter end,
     Callback callback)
 {
     typedef typename std::iterator_traits<RandomAccessIter>::value_type val_t;
     typedef Box_intersection_d::Box_traits_d< val_t>  Box_traits;
-    box_self_intersection_d(begin, end, callback, Box_traits());
+    box_self_intersection_d<ConcurrencyTag>(begin, end, callback, Box_traits());
 }
 
-template< class RandomAccessIter, class Callback >
+template< class ConcurrencyTag = Sequential_tag, class RandomAccessIter, class Callback >
 void box_self_intersection_d(
     RandomAccessIter begin, RandomAccessIter end,
     Callback callback,
@@ -226,10 +410,10 @@ void box_self_intersection_d(
 {
     typedef typename std::iterator_traits<RandomAccessIter>::value_type val_t;
     typedef Box_intersection_d::Box_traits_d< val_t>  Box_traits;
-    box_self_intersection_d(begin, end, callback, Box_traits(), cutoff);
+    box_self_intersection_d<ConcurrencyTag>(begin, end, callback, Box_traits(), cutoff);
 }
 
-template< class RandomAccessIter, class Callback >
+  template< class ConcurrencyTag = Sequential_tag, class RandomAccessIter, class Callback >
 void box_self_intersection_d(
     RandomAccessIter begin, RandomAccessIter end,
     Callback callback,
@@ -238,11 +422,10 @@ void box_self_intersection_d(
 {
     typedef typename std::iterator_traits<RandomAccessIter>::value_type val_t;
     typedef Box_intersection_d::Box_traits_d< val_t>  Box_traits;
-    box_self_intersection_d(begin, end, callback,
+    box_self_intersection_d<ConcurrencyTag>(begin, end, callback,
                                             Box_traits(), cutoff, topology );
 }
 
-
 // Generic call for trivial all-pairs algorithm with box traits parameter.
 // - make all default parameters explicit overloads (workaround)
 template< class ForwardIter1, class ForwardIter2,

Box_intersection_d/test/Box_intersection_d/CMakeLists.txt

@@ -5,20 +5,23 @@
 cmake_minimum_required(VERSION 3.1...3.15)
 project( Box_intersection_d_Tests )
 
+find_package( CGAL QUIET )
 
-find_package(CGAL QUIET)
+find_package( TBB )
 
 if ( CGAL_FOUND )
 
-  # create a target per cppfile
+  create_single_source_cgal_program( "automated_test.cpp" )
-  file(GLOB cppfiles RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp)
+  create_single_source_cgal_program( "benchmark_box_intersection.cpp" )
-  foreach(cppfile ${cppfiles})
+  create_single_source_cgal_program( "random_set_test.cpp" )
-    create_single_source_cgal_program( "${cppfile}" )
+  create_single_source_cgal_program( "test_box_grid.cpp" )
-  endforeach()
 
+  if( TBB_FOUND )
+    CGAL_target_use_TBB( test_box_grid )
+  else()
+    message( STATUS "NOTICE: Intel TBB was not found. Sequential code will be used." )
+  endif()
 else()
-  
     message(STATUS "This program requires the CGAL library, and will not be compiled.")
-  
 endif()
 

Box_intersection_d/test/Box_intersection_d/test_box_grid.cpp

@@ -1,33 +1,45 @@
 // file: test/Box_intersection_d/box_grid.C
 // similar to examples/Box_intersection_d/box_grid.C but stricter in checking
 // and more extensive in what is tested.
-#include <CGAL/box_intersection_d.h>
-#include <cassert>
-#include <vector>
-#include <algorithm>
-#include <iterator>
 
-typedef CGAL::Box_intersection_d::Box_d<int,2> Box;
+#include <CGAL/box_intersection_d.h>
+#include <CGAL/tags.h>
+
+#include <algorithm>
+#include <cassert>
+#include <iterator>
+#include <vector>
+
+#ifdef CGAL_LINKED_WITH_TBB
+#include <tbb/concurrent_vector.h>
+#endif
+
+typedef CGAL::Box_intersection_d::Box_d<int, 2> Box;
 
 // coordinates for 9 boxes of a grid
 int p[9*4]   = { 0,0,1,1,  1,0,2,1,  2,0,3,1, // lower
                  0,1,1,2,  1,1,2,2,  2,1,3,2, // middle
                  0,2,1,3,  1,2,2,3,  2,2,3,3};// upper
 // 9 boxes + 2 selected boxes as query; center and upper right
-Box init_boxes[11] = { Box( p,    p+ 2),  Box( p+ 4, p+ 6),  Box( p+ 8, p+10),
+const Box init_boxes[11] = { Box( p,    p+ 2),  Box( p+ 4, p+ 6),  Box( p+ 8, p+10),
                              Box( p+12, p+14),  Box( p+16, p+18),  Box( p+20, p+22),
                              Box( p+24, p+26),  Box( p+28, p+30),  Box( p+32, p+34),
                              Box( p+16, p+18),  Box( p+32, p+34)};
 Box  boxes[11];
 Box* query = boxes+9; 
 
-void init() {
+void init()
+{
   for ( int i = 0; i < 11; ++i)
     boxes[i] = init_boxes[i];
 }
 
-void check_result( const char* text, std::vector<std::size_t>& result, 
+template <typename Vector>
-                   const std::size_t* check, std::size_t size) {
+void check_result( const char* text,
+                   Vector& result,
+                   const std::size_t* check,
+                   std::size_t size)
+{
   // sort, show, and check result
   std::sort( result.begin(), result.end());
   std::cout << text << ": got " << result.size() << " elements, expected "
@@ -43,46 +55,112 @@ void check_result( const char* text, std::vector<std::size_t>& result,
 }
 
 // callback function object writing results to an output iterator
-template <class OutputIterator>
+template <class Container>
-struct Report {
+struct Report
-    OutputIterator it;
+{
-    Report( OutputIterator i) : it(i) {} // store iterator in object
+  Container& c_;
+
+  Report(Container& c) : c_(c) {} // store iterator in object
   // We encode both id-numbers in a single number, a.id() + 100 * b.id(),
   // and write that number to the output iterator.
-    void operator()( const Box& a, const Box& b) { *it++ = a.id()+100*b.id(); }
+  void operator()(const Box& a, const Box& b) { c_.push_back(a.id()+100*b.id()); }
+  void operator()(const Box* a, const Box* b) { c_.push_back(a->id()+100*b->id()); }
 };
-template <class Iter> // helper function to create the function object
+
-Report<Iter> report( Iter it) { return Report<Iter>(it); }
+template <class Container> // helper function to create the function object
+Report<Container> report(Container& c) { return Report<Container>(c); }
 
 // ---------------------------------------------------------------------
 // box_intersection_d
 // run the intersection algorithms and store results in a vector
 // ---------------------------------------------------------------------
-void test_box_intersection() {
+void test_box_intersection()
+{
+#ifdef CGAL_LINKED_WITH_TBB
+  tbb::concurrent_vector<std::size_t> result;
+#else
+  std::vector<std::size_t> result;
+#endif
+
+  // Some degenerate cases
+  init();
+  CGAL::box_intersection_d<CGAL::Parallel_if_available_tag>(boxes, boxes, boxes, boxes,
+                                                            report(result));
+  assert(result.empty());
+
+  init();
+  CGAL::box_intersection_d<CGAL::Parallel_if_available_tag>(boxes, boxes, query, query + 1,
+                                                            report(result));
+  assert(result.empty());
+
+  init();
+  CGAL::box_intersection_d<CGAL::Parallel_if_available_tag>(boxes, boxes + 3, query, query,
+                                                            report(result));
+  assert(result.empty());
+
+  // With pointers
+  init();
+  std::vector<const Box*> range_1 = {{ boxes, boxes+1, boxes+2, boxes+3, boxes+4, boxes+5,
+                                       boxes+6, boxes+7, boxes+8 }};
+  std::vector<const Box*> range_2 = {{ query, query+1 }};
+
+  CGAL::box_intersection_d<CGAL::Parallel_if_available_tag>(range_1.begin(), range_1.begin(),
+                                                            range_2.begin(), range_2.end(),
+                                                            report(result));
+  assert(result.empty());
+
+  CGAL::box_intersection_d<CGAL::Parallel_if_available_tag>(range_1.begin(), range_1.end(),
+                                                            range_2.begin(), range_2.begin(),
+                                                            report(result));
+  assert(result.empty());
+
+  CGAL::box_intersection_d<CGAL::Parallel_if_available_tag>(range_1.begin(), range_1.end(),
+                                                            range_2.begin(), range_2.end(),
+                                                            report(result));
+  std::size_t check0[13] = {900,901,902,903,904,905,906,907,908,
+                            1004,1005,1007,1008};
+  check_result( "Box inters. 3x3 (ptr), 2, closed", result, check0, 13);
+
   // intersect 3x3 with 2 query boxes, closed boxes
   init();
-    std::vector<std::size_t> result;
+  result.clear();
-    CGAL::box_intersection_d( boxes, boxes+9, query, query+2, 
+  CGAL::box_intersection_d<CGAL::Parallel_if_available_tag>(boxes, boxes+9, query, query+1, report(result));
-                              report( std::back_inserter( result)));
+  std::size_t check1[13] = {900,901,902,903,904,905,906,907,908};
-    std::size_t check1[13] = {900,901,902,903,904,905,906,907,908,
+  check_result( "Box inters. 3x3, 2, closed", result, check1, 9);
+
+  // intersect 3x3 with 2 query boxes, closed boxes
+  init();
+  result.clear();
+  CGAL::box_intersection_d<CGAL::Parallel_if_available_tag>(boxes, boxes+9, query, query+2, report(result));
+  std::size_t check1bis[13] = {900,901,902,903,904,905,906,907,908,
                             1004,1005,1007,1008};
-    check_result( "Box inters. 3x3, 2, closed", result, check1, 13);
+  check_result( "Box inters. 3x3, 2, closed", result, check1bis, 13);
 
   // intersect 3x3 with 2 query boxes, half-open boxes and changed cutoff
   init();
   result.clear();
   CGAL::box_intersection_d( boxes, boxes+9, query, query+2,
-                              report( std::back_inserter( result)),
+                            report(result),
                             std::ptrdiff_t(1),
                             CGAL::Box_intersection_d::HALF_OPEN);
   std::size_t check2[2]  = {904,1008};
   check_result( "Box inters. 3x3, 2, half-open", result, check2, 2);
 
+  // intersect 3x3 with 2 query boxes, half-open boxes and changed cutoff (reversed)
+  init();
+  result.clear();
+  CGAL::box_intersection_d<CGAL::Parallel_if_available_tag>( query, query+2, boxes, boxes+9,
+                                                             report(result),
+                                                             std::ptrdiff_t(1),
+                                                             CGAL::Box_intersection_d::HALF_OPEN);
+  std::size_t check2bis[2]  = {409,810};
+  check_result( "Box inters. 3x3 (reversed), 2, half-open", result, check2bis, 2);
+
   // self intersect 3x2, closed boxes
   init();
   result.clear();
   CGAL::box_self_intersection_d( boxes, boxes+6,
-                                   report( std::back_inserter( result)));
+                                 report(result));
   std::size_t check3[11] = {1,3,4,102,103,104,105,204,205,304,405};
   check_result( "Box self inters. 3x2, closed", result, check3, 11);
 
@@ -90,7 +168,7 @@ void test_box_intersection() {
   init();
   result.clear();
   CGAL::box_self_intersection_d( boxes, boxes+6,
-                                   report( std::back_inserter( result)),
+                                 report(result),
                                  std::ptrdiff_t(1),
                                  CGAL::Box_intersection_d::HALF_OPEN);
   std::size_t check4[1] = {9999};
@@ -100,7 +178,7 @@ void test_box_intersection() {
   init();
   result.clear();
   CGAL::box_self_intersection_d( boxes, boxes+11,
-                                   report( std::back_inserter( result)),
+                                 report(result),
                                  std::ptrdiff_t(1),
                                  CGAL::Box_intersection_d::HALF_OPEN);
   std::size_t check5[2] = {409,810};
@@ -112,7 +190,7 @@ void test_box_intersection() {
   result.clear();
   std::vector<Box> boxes2( boxes, boxes+11);
   CGAL::box_intersection_d( boxes, boxes+11, boxes2.begin(), boxes2.end(),
-                              report( std::back_inserter( result)),
+                            report(result),
                             std::ptrdiff_t(1),
                             CGAL::Box_intersection_d::HALF_OPEN,
                             CGAL::Box_intersection_d::COMPLETE);
@@ -125,7 +203,7 @@ void test_box_intersection() {
   result.clear();
   boxes2 = std::vector<Box>( boxes, boxes+11);
   CGAL::box_intersection_d( boxes, boxes+11, boxes2.begin(), boxes2.end(),
-                              report( std::back_inserter( result)),
+                            report(result),
                             std::ptrdiff_t(20),
                             CGAL::Box_intersection_d::HALF_OPEN,
                             CGAL::Box_intersection_d::BIPARTITE);
@@ -137,12 +215,13 @@ void test_box_intersection() {
 // box_intersection_all_pairs_d
 // run the intersection algorithms and store results in a vector
 // ---------------------------------------------------------------------
-void test_box_intersection_all_pairs() {
+void test_box_intersection_all_pairs()
+{
   // intersect 3x3 with 2 query boxes, closed boxes
   init();
   std::vector<std::size_t> result;
   CGAL::box_intersection_all_pairs_d( boxes, boxes+9, query, query+2,
-                                        report( std::back_inserter( result)));
+                                      report(result));
   std::size_t check1[13] = {900,901,902,903,904,905,906,907,908,
                             1004,1005,1007,1008};
   check_result( "All-pairs inters. 3x3, 2, closed", result, check1, 13);
@@ -151,7 +230,7 @@ void test_box_intersection_all_pairs() {
   init();
   result.clear();
   CGAL::box_intersection_all_pairs_d( boxes, boxes+9, query, query+2,
-                                        report( std::back_inserter( result)),
+                                      report(result),
                                       CGAL::Box_intersection_d::HALF_OPEN);
   std::size_t check2[2]  = {904,1008};
   check_result( "All-pairs inters. 3x3, 2, half-open", result, check2, 2);
@@ -160,7 +239,7 @@ void test_box_intersection_all_pairs() {
   init();
   result.clear();
   CGAL::box_self_intersection_all_pairs_d( boxes, boxes+6,
-                                   report( std::back_inserter( result)));
+                                           report(result));
   std::size_t check3[11] = {100,201,300,301,400,401,402,403,501,502,504};
   check_result( "All-pairs self inters. 3x2, closed", result, check3, 11);
 
@@ -168,7 +247,7 @@ void test_box_intersection_all_pairs() {
   init();
   result.clear();
   CGAL::box_self_intersection_all_pairs_d( boxes, boxes+6,
-                                   report( std::back_inserter( result)),
+                                           report(result),
                                            CGAL::Box_intersection_d::HALF_OPEN);
   std::size_t check4[1] = {9999};
   check_result( "All-pairs self inters. 3x2, half-open", result, check4, 0);
@@ -177,7 +256,7 @@ void test_box_intersection_all_pairs() {
   init();
   result.clear();
   CGAL::box_self_intersection_all_pairs_d( boxes, boxes+11,
-                                   report( std::back_inserter( result)),
+                                           report(result),
                                            CGAL::Box_intersection_d::HALF_OPEN);
   std::size_t check5[2] = {904,1008};
   check_result( "All-pairs self inters. 3x3+2, half-open", result, check5,2);
@@ -189,7 +268,7 @@ void test_box_intersection_all_pairs() {
   std::vector<Box> boxes2( boxes, boxes+11);
   CGAL::box_intersection_all_pairs_d( boxes, boxes+11,
                                       boxes2.begin(), boxes2.end(),
-                              report( std::back_inserter( result)),
+                                      report(result),
                                       CGAL::Box_intersection_d::HALF_OPEN,
                                       CGAL::Box_intersection_d::COMPLETE);
   check_result( "All-pairs self inters. 3x3+2, half-open", result, check5,2);
@@ -202,16 +281,17 @@ void test_box_intersection_all_pairs() {
   boxes2 = std::vector<Box>( boxes, boxes+11);
   CGAL::box_intersection_all_pairs_d( boxes, boxes+11,
                                       boxes2.begin(), boxes2.end(),
-                              report( std::back_inserter( result)),
+                                      report(result),
                                       CGAL::Box_intersection_d::HALF_OPEN,
                                       CGAL::Box_intersection_d::BIPARTITE);
   std::size_t check6[4] = {409,810,904,1008};
   check_result( "All-pairs inters. 3x3+2, half-open", result, check6, 4);
 }
 
-
+int main()
-int main() {
+{
   test_box_intersection();
   test_box_intersection_all_pairs();
-    return 0;
+
+  return EXIT_SUCCESS;
 }

Classification/examples/Classification/example_classification.cpp

@@ -11,16 +11,12 @@
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/Classification.h>
 #include <CGAL/bounding_box.h>
-
+#include <CGAL/tags.h>
 #include <CGAL/IO/read_ply_points.h>
 
 #include <CGAL/Real_timer.h>
 
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
 
 typedef CGAL::Simple_cartesian<double> Kernel;
 typedef Kernel::Point_3 Point;

Classification/examples/Classification/example_cluster_classification.cpp

@@ -18,11 +18,7 @@
 #include <CGAL/Shape_detection/Region_growing.h>
 #include <CGAL/Real_timer.h>
 
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag   Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
 
 typedef CGAL::Simple_cartesian<double> Kernel;
 typedef Kernel::Point_3                Point;

Classification/include/CGAL/Classification/ETHZ/Random_forest_classifier.h

@@ -127,11 +127,7 @@ public:
               std::size_t num_trees = 25,
               std::size_t max_depth = 20)
   {
-#ifdef CGAL_LINKED_WITH_TBB
+    train<CGAL::Parallel_if_available_tag>(ground_truth, reset_trees, num_trees, max_depth);
-    train<CGAL::Parallel_tag>(ground_truth, reset_trees, num_trees, max_depth);
-#else
-    train<CGAL::Sequential_tag>(ground_truth, reset_trees, num_trees, max_depth);
-#endif
   }
   /// \endcond
     
@@ -146,7 +142,7 @@ public:
     label.
 
     \tparam ConcurrencyTag enables sequential versus parallel
-    algorithm. Possible values are `Parallel_tag` (default value is
+    algorithm. Possible values are `Parallel_tag` (default value if
     %CGAL is linked with TBB) or `Sequential_tag` (default value
     otherwise).
 

Classification/include/CGAL/Classification/Local_eigen_analysis.h

@@ -226,7 +226,7 @@ public:
     is `CGAL::Point_3`.
     \tparam NeighborQuery model of `NeighborQuery`
     \tparam ConcurrencyTag enables sequential versus parallel
-    algorithm. Possible values are `Parallel_tag` (default value is %CGAL
+    algorithm. Possible values are `Parallel_tag` (default value if %CGAL
     is linked with TBB) or `Sequential_tag` (default value otherwise).
     \tparam DiagonalizeTraits model of `DiagonalizeTraits` used for
     matrix diagonalization. It can be omitted if Eigen 3 (or greater)
@@ -242,10 +242,8 @@ public:
             typename NeighborQuery,
 #if defined(DOXYGEN_RUNNING)
             typename ConcurrencyTag,
-#elif defined(CGAL_LINKED_WITH_TBB)
-            typename ConcurrencyTag = CGAL::Parallel_tag,
 #else
-            typename ConcurrencyTag = CGAL::Sequential_tag,
+            typename ConcurrencyTag = CGAL::Parallel_if_available_tag,
 #endif
 #if defined(DOXYGEN_RUNNING)
             typename DiagonalizeTraits>
@@ -310,7 +308,7 @@ public:
     \tparam FaceListGraph model of `FaceListGraph`. 
     \tparam NeighborQuery model of `NeighborQuery`
     \tparam ConcurrencyTag enables sequential versus parallel
-    algorithm. Possible values are `Parallel_tag` (default value is %CGAL
+    algorithm. Possible values are `Parallel_tag` (default value if %CGAL
     is linked with TBB) or `Sequential_tag` (default value otherwise).
     \tparam DiagonalizeTraits model of `DiagonalizeTraits` used for
     matrix diagonalization. It can be omitted: if Eigen 3 (or greater)
@@ -325,10 +323,8 @@ public:
             typename NeighborQuery,
 #if defined(DOXYGEN_RUNNING)
             typename ConcurrencyTag,
-#elif defined(CGAL_LINKED_WITH_TBB)
-            typename ConcurrencyTag = CGAL::Parallel_tag,
 #else
-            typename ConcurrencyTag = CGAL::Sequential_tag,
+            typename ConcurrencyTag = CGAL::Parallel_if_available_tag,
 #endif
 #if defined(DOXYGEN_RUNNING)
             typename DiagonalizeTraits>
@@ -396,7 +392,7 @@ public:
     `RandomAccessIterator` and its value type is the key type of
     `PointMap`.
     \tparam ConcurrencyTag enables sequential versus parallel
-    algorithm. Possible values are `Parallel_tag` (default value is %CGAL
+    algorithm. Possible values are `Parallel_tag` (default value if %CGAL
     is linked with TBB) or `Sequential_tag` (default value otherwise).
     \tparam DiagonalizeTraits model of `DiagonalizeTraits` used for
     matrix diagonalization. It can be omitted: if Eigen 3 (or greater)
@@ -409,10 +405,8 @@ public:
   template <typename ClusterRange,
 #if defined(DOXYGEN_RUNNING)
             typename ConcurrencyTag,
-#elif defined(CGAL_LINKED_WITH_TBB)
-            typename ConcurrencyTag = CGAL::Parallel_tag,
 #else
-            typename ConcurrencyTag = CGAL::Sequential_tag,
+            typename ConcurrencyTag = CGAL::Parallel_if_available_tag,
 #endif
 #if defined(DOXYGEN_RUNNING)
             typename DiagonalizeTraits>

Classification/include/CGAL/Classification/Mesh_feature_generator.h

@@ -79,7 +79,7 @@ namespace Classification {
   is `GeomTraits::Point_3`.
   \tparam ConcurrencyTag enables sequential versus parallel
   computation of `CGAL::Classification::Local_eigen_analysis`
-  objects. Possible values are `Parallel_tag` (default value is %CGAL
+  objects. Possible values are `Parallel_tag` (default value if %CGAL
   is linked with TBB) or `Sequential_tag` (default value otherwise).
   \tparam DiagonalizeTraits model of `DiagonalizeTraits` used for
   matrix diagonalization. It can be omitted: if Eigen 3 (or greater)
@@ -93,10 +93,8 @@ template <typename GeomTraits,
           typename PointMap,
 #if defined(DOXYGEN_RUNNING)
           typename ConcurrencyTag,
-#elif defined(CGAL_LINKED_WITH_TBB)
-          typename ConcurrencyTag = CGAL::Parallel_tag,
 #else
-          typename ConcurrencyTag = CGAL::Sequential_tag,
+          typename ConcurrencyTag = CGAL::Parallel_if_available_tag,
 #endif
           typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<float,3> >
 class Mesh_feature_generator

Classification/include/CGAL/Classification/Point_set_feature_generator.h

@@ -90,10 +90,8 @@ template <typename GeomTraits,
           typename PointMap,
 #if defined(DOXYGEN_RUNNING)
           typename ConcurrencyTag,
-#elif defined(CGAL_LINKED_WITH_TBB)
-          typename ConcurrencyTag = CGAL::Parallel_tag,
 #else
-          typename ConcurrencyTag = CGAL::Sequential_tag,
+          typename ConcurrencyTag = CGAL::Parallel_if_available_tag,
 #endif
 #if defined(DOXYGEN_RUNNING)
           typename DiagonalizeTraits>

Classification/include/CGAL/Classification/classify.h

@@ -334,7 +334,7 @@ namespace internal {
     suboptimal results.
 
     \tparam ConcurrencyTag enables sequential versus parallel
-    algorithm. Possible values are `Parallel_tag` or `Sequential_tag`.
+    algorithm. Possible values are `Parallel_if_available_tag`, `Parallel_tag` or `Sequential_tag`.
     
     \tparam ItemRange model of `ConstRange`. Its iterator type is
     `RandomAccessIterator`. Its value type depends on the data that is
@@ -424,7 +424,7 @@ namespace internal {
     efficiency and better quality results.
 
     \tparam ConcurrencyTag enables sequential versus parallel
-    algorithm. Possible values are `Parallel_tag` or `Sequential_tag`.
+    algorithm. Possible values are `Parallel_if_available_tag`, `Parallel_tag` or `Sequential_tag`.
     \tparam ItemRange model of `ConstRange`. Its iterator type is
     `RandomAccessIterator`.
     \tparam ItemMap model of `ReadablePropertyMap` whose key
@@ -502,7 +502,7 @@ namespace internal {
     results.
 
     \tparam ConcurrencyTag enables sequential versus parallel
-    algorithm. Possible values are `Parallel_tag` or `Sequential_tag`.
+    algorithm. Possible values are `Parallel_if_available_tag`, `Parallel_tag` or `Sequential_tag`.
     \tparam ItemRange model of `ConstRange`. Its iterator type is
     `RandomAccessIterator`.
     \tparam ItemMap model of `ReadablePropertyMap` whose key

Convex_hull_2/doc/Convex_hull_2/CGAL/ch_jarvis.h

@@ -12,12 +12,12 @@ cyclic sequence of extreme points is cut into a linear sequence.
 \pre The source range [`first`,`beyond`) does not contain `result`.
 
 The default traits class `Default_traits` is the kernel in which the 
-value type of `InputIterator` is defined. 
+value type of `ForwardIterator` is defined. 
 
 \cgalHeading{Requirements}
 
 <OL> 
-<LI>The value type of `InputIterator` and 
+<LI>The value type of `ForwardIterator` and 
 `OutputIterator` is equivalent to `Traits::Point_2`. 
 <LI>`Traits` defines the following subset of types from 
 the concept `ConvexHullTraits_2` and their corresponding member 
@@ -46,10 +46,10 @@ in the worst case for \f$ n\f$ input points with \f$ h\f$ extreme points.
 
 
 */
-template <class InputIterator, class OutputIterator, class Traits>
+template <class ForwardIterator, class OutputIterator, class Traits>
 OutputIterator
-ch_jarvis( InputIterator first,
+ch_jarvis( ForwardIterator first,
-InputIterator beyond,
+ForwardIterator beyond,
 OutputIterator result,
 const Traits & ch_traits = Default_traits);
 

Filtered_kernel/include/CGAL/internal/Static_filters/Coplanar_3.h

@@ -0,0 +1,54 @@
+// Copyright (c) 2011 GeometryFactory (France)
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Andreas Fabri, Laurent Rineau
+
+
+#ifndef CGAL_INTERNAL_STATIC_FILTERS_COPLANAR_3_H
+#define CGAL_INTERNAL_STATIC_FILTERS_COPLANAR_3_H
+
+namespace CGAL {
+
+namespace internal {
+
+namespace Static_filters_predicates {
+
+
+template < typename K_base, typename SFK >
+class Coplanar_3
+  : public K_base::Coplanar_3
+{
+  typedef typename K_base::Point_3 Point_3;
+  typedef typename K_base::Coplanar_3 Base;
+  typedef typename SFK::Orientation_3 Orientation_3;
+
+public:
+
+  typedef typename Base::result_type  result_type;
+
+
+
+  result_type 
+  operator()(const Point_3& p,const Point_3& q, const Point_3& r, const Point_3& s) const
+  {
+    return Orientation_3()(p,q,r,s) == COPLANAR;
+  }
+
+
+
+}; // end class Coplanar_3
+
+} // end namespace Static_filters_predicates
+
+} // end namespace internal
+
+} // end namespace CGAL
+
+#endif  // CGAL_INTERNAL_STATIC_FILTERS_COPLANAR_3_H

Filtered_kernel/include/CGAL/internal/Static_filters/Static_filters.h

@@ -41,7 +41,6 @@
 
 #endif // CGAL_DISABLE_STATIC_FILTERS_ADDED_2011
 
-
 #ifndef CGAL_NO_EQUAL_3_STATIC_FILTERS
 #  include <CGAL/internal/Static_filters/Equal_3.h>
 #  include <CGAL/internal/Static_filters/Equal_2.h>
@@ -65,6 +64,7 @@
 #  include <CGAL/internal/Static_filters/Do_intersect_2.h>
 #endif // NOT NOT CGAL_NO_DO_INTERSECT_STATIC_FILTERS
 
+#include <CGAL/internal/Static_filters/Coplanar_3.h>
 #include <CGAL/internal/Static_filters/Compare_y_at_x_2.h>
 #include <CGAL/internal/Static_filters/Side_of_oriented_circle_2.h>
 #include <CGAL/internal/Static_filters/Side_of_oriented_sphere_3.h>
@@ -121,6 +121,7 @@ public:
   typedef Static_filters_predicates::Side_of_oriented_circle_2<K_base>      Side_of_oriented_circle_2;
   typedef Static_filters_predicates::Side_of_oriented_sphere_3<K_base>      Side_of_oriented_sphere_3;
   typedef Static_filters_predicates::Compare_squared_radius_3<K_base>       Compare_squared_radius_3;
+  typedef Static_filters_predicates::Coplanar_3<K_base,Self>                Coplanar_3;
 
   typedef Static_filters_predicates::Compare_weighted_squared_radius_3<K_base>     Compare_weighted_squared_radius_3;
   typedef Static_filters_predicates::Power_side_of_oriented_power_sphere_3<K_base>                          Power_side_of_oriented_power_sphere_3;
@@ -152,7 +153,7 @@ public:
   compare_x_2_object() const
   { return Compare_x_2(); }
 
-Compare_y_2
+  Compare_y_2
   compare_y_2_object() const
   { return Compare_y_2(); }
 
@@ -182,7 +183,12 @@ Compare_y_2
   compare_squared_radius_3_object() const
   { return Compare_squared_radius_3(); }
 
-  Power_side_of_oriented_power_sphere_3 power_side_of_oriented_power_sphere_3_object() const
+  Coplanar_3
+  coplanar_3_object() const
+  { return Coplanar_3(); }
+
+  Power_side_of_oriented_power_sphere_3
+  power_side_of_oriented_power_sphere_3_object() const
   { return Power_side_of_oriented_power_sphere_3();}
 
   Compare_weighted_squared_radius_3

Generator/include/CGAL/point_generators_d.h

@@ -18,6 +18,8 @@
 
 #include <CGAL/generators.h>
 #include <CGAL/number_type_basic.h>
+#include <CGAL/number_type_config.h>
+
 #include <cmath>
 #include <vector>
 
@@ -56,11 +58,10 @@ generate_point() {
   for(int i=0; i<dimension; ++i) {
     // normal distribution 
     //( a product of normal distib is a normal distrib in higher dim)
-    const double pi = 3.141592653589793238462643;
     double a=this->_rnd.get_double();
     a = std::sqrt( -2* std::log(1-a) );
     double b=this->_rnd.get_double();
-    b = std::cos(2*pi*b);
+    b = std::cos(2*CGAL_PI*b);
     coord[i]= a*b;
     norm += coord[i]*coord[i];
   }
@@ -104,11 +105,10 @@ generate_point() {
   
   for(int i=0; i<dimension; ++i) {
     // normal distribution
-    const double pi = 3.141592653589793238462643;
     double a=this->_rnd.get_double();
     a = std::sqrt( -2* std::log(1-a) );
     double b=this->_rnd.get_double();
-    b = std::cos(2*pi*b);
+    b = std::cos(2*CGAL_PI*b);
     coord[i]= a*b;
     norm += coord[i]*coord[i];
   }

Hash_map/test/Hash_map/Hash.cpp

@@ -145,8 +145,9 @@ int main()
 
 #ifdef CGAL_LINKED_WITH_TBB
   Triangulation_3 T3;
-fct3(T3);
+  fct3(T3);
 #endif
+
   return 0;
 }
 

INSTALL.md

@@ -1,8 +1,7 @@
-NOTICE
+Building an Example or a Demo of CGAL
-======
+=====================================
 
-Since Version 5.0, CGAL is a header-only library it is not needed
+Since Version 5.0, CGAL is a header-only library, hence it is not needed to build it. Usage of CGAL should simply amount to:
-to build and install it. Usage of CGAL should thus simply amount to:
 
 ``` {.bash}
 git clone https://github.com/CGAL/cgal.git /path/to/cgal.git
@@ -13,43 +12,41 @@ cmake -DCMAKE_BUILD_TYPE=Debug -DCGAL_DIR=/path/to/cgal.git
 make
 ```
 
-in the case of the building of an example in debug mode.
+in the case of building some CGAL-provided examples in debug mode.
-
-For more information head over to the [CGAL manual](https://doc.cgal.org/latest/Manual/general_intro.html).
-Note that this page describes the setting of CGAL as a sources release and, as such,
-files are organized in a slightly different way, see the [Layout of the CGAL Git Repository](README.md).
 
+Note that although CGAL is a header-only library, some parts of it must link to several external libraries, such as GMP, MPFR, etc.
 
 Building a Program Using CGAL
 =============================
 
-To compile a program using CGAL, simply set `CGAL_DIR` to the location
+If you wish to build a program that is not provided with CGAL and does not already have a `CMakeLists.txt`,
-of the directory containing `CGALConfig.cmake` (for example the root 
+you can trigger the creation of a basic `CMakeLists.txt` by calling the script [`cgal_create_cmake_script`](Scripts/scripts/cgal_create_cmake_script)
-of the extracted source archive or the root of a git checkout).
+found in `/path/to/cgal.git/Scripts/scripts/` at the root of your program directory.
-
-Here is an example of how to build in debug the examples from the 3D Triangulations package:
 
 ``` {.bash}
- cd /path/to/cgal.git/Triangulation_3/examples/Triangulation_3
+ git clone https://github.com/CGAL/cgal.git /path/to/cgal.git
- mkdir -p build/debug
+ cd /path/to/your/program
- cd build/debug
- cmake -DCGAL_DIR:PATH=/path/to/cgal.git ../..
- make
-```
-
-If you are trying to build examples or tests that do not already have a `CMakeLists.txt`,
-you can trigger its creation by calling the script [`cgal_create_cmake_script`](Scripts/scripts/cgal_create_cmake_script)
-found in `/path/to/cgal.git/Scripts/scripts/` at the root of the example/test directory.
-Here is an example for the examples of the 2D Triangulation package:
-
-``` {.bash}
- cd /path/to/cgal.git/Triangulation_2/examples/Triangulation_2
  /path/to/cgal.git/Scripts/scripts/cgal_create_cmake_script
- cd /path/to/cgal.git/Triangulation_2/examples/Triangulation_2
  mkdir -p build/debug
  cd build/debug
- cmake -DCGAL_DIR:PATH=/path/to/cgal.git ../..
+ cmake -DCMAKE_BUILD_TYPE=Debug -DCGAL_DIR:PATH=/path/to/cgal.git ../..
- make
+ make your_program
 ```
 
-For more information head over to the [CGAL manual](https://doc.cgal.org/latest/Manual/general_intro.html).
+Since the basic `CMakeLists.txt` created by the script `cgal_create_cmake_script` cannot
+guess which part(s) of CGAL you are using, it does not link with any optional third party 
+dependency of CGAL. You should look at the documentation of the package(s) that you
+are using to learn which dependencies you must add. The `CMakeLists.txt`
+of the examples and demos provided with the package(s) that you are using can be used
+to complete your basic `CMakeLists.txt`.
+
+
+Repository Structure
+====================
+
+If you have downloaded a source release instead of cloning the Git repository, the files will be organized in a slightly different way, see the [Layout of the CGAL Git Repository](README.md).
+
+Documentation
+=============
+
+For more information see the [CGAL manual](https://doc.cgal.org/latest/Manual/general_intro.html).

Installation/CHANGES.md

@@ -18,6 +18,9 @@ Release date: June 2020
 -   **Breaking change**: the internal search tree is now lazily constructed. To disable it, one must call
     the new function `do_not_accelerate_distance_queries()` before the first distance query.
 
+### Intersecting Sequences of dD Iso-oriented Boxes 
+ -   Added parallel versions of the functions `CGAL::box_intersection_d()` and `CGAL::box_self_intersection_d()`.
+
 ### Polygon Mesh Processing
 
 -   Introduced a new function, `CGAL::Polygon_mesh_processing::remove_connected_components_of_negligible_size()`, 
@@ -28,6 +31,8 @@ Release date: June 2020
     components that would be removed with the specified threshold, but without actually removing them.
 -   The function `CGAL::Polygon_mesh_processing::stitch_borders()` now returns the number
     of halfedge pairs that were stitched.
+ -   Added parallel versions of the functions `CGAL::Polygon_mesh_processing::does_self_intersect()` 
+     and `CGAL::Polygon_mesh_processing::self_intersections()`.
 
 ### 2D Triangulations
 -   To fix an inconsistency between code and documentation and to clarify which types of intersections
@@ -65,6 +70,10 @@ Release date: June 2020
         k-NN search to interrupt some distance computations before its end,
         saving precious milliseconds, in particular in medium-to-high dimension.
 
+### Spatial Sorting
+ -   Added parallel versions of `hilbert_sort()` and `spatial_sort()` in 2D and 3D when the median policy is used.
+     The parallel versions use up to four threads in 2D, and up to eight threads in 3D.
+
 ### dD Geometry Kernel
 -   Epick\_d and Epeck\_d gain 2 new functors: `Power_side_of_bounded_power_sphere_d` and
     `Compute_squared_radius_smallest_orthogonal_sphere_d`. Those are
@@ -75,6 +84,10 @@ Release date: June 2020
 - The concept "EdgeProfile" has been removed. This concept was not actually in use as the CGAL-provided model `CGAL::Edge_profile`
   was imposed to the user. Other concepts have been clarified to reflect the fact that the API uses this particular class.
 
+### STL Extensions for CGAL
+ -   Added a new concurrency tag: `CGAL::Parallel_if_available_tag`. This tag is a convenience typedef to `CGAL::Parallel_tag`
+     if the third party library TBB has been found and linked with, and to `CGAL::Sequential_tag` otherwise.
+
 [Release 5.0](https://github.com/CGAL/cgal/releases/tag/releases%2FCGAL-5.0)
 -----------
 

Installation/CMakeLists.txt

@@ -568,10 +568,6 @@ message("== Detect external libraries ==")
 #  Coin is used in KDS, but no FindCoin or FindCOIN exists
 #  There exists FindIPE, FindMKL, but they are only used to support supporting libs
 list (INSERT CGAL_SUPPORTING_3RD_PARTY_LIBRARIES 0 GMP MPFR ZLIB OpenGL LEDA MPFI RS RS3 OpenNL Eigen3 ESBTL Coin3D NTL IPE)
-if (NOT WIN32)
-  # GMPXX is not supported on WIN32 machines
-  list (INSERT CGAL_SUPPORTING_3RD_PARTY_LIBRARIES 1 GMPXX)
-endif()
 
 # Where CMake is run several times, to avoid duplicates
 list(REMOVE_DUPLICATES CGAL_SUPPORTING_3RD_PARTY_LIBRARIES)
@@ -611,8 +607,6 @@ if(CGAL_DISABLE_GMP)
   unset(WITH_GMP CACHE)
   unset(CGAL_USE_GMP)
   unset(CGAL_USE_GMP CACHE)
-  unset(WITH_GMPXX)
-  unset(WITH_GMPXX CACHE)
   unset(WITH_MPFR)
   unset(WITH_MPFR CACHE)
 
@@ -621,9 +615,6 @@ if(CGAL_DISABLE_GMP)
     "#error GMP is disabled by the CMake option CGAL_DISABLE_GMP")
 else()
   list(APPEND      CGAL_ESSENTIAL_3RD_PARTY_LIBRARIES GMP MPFR)
-  if(WITH_GMPXX)
-    list(APPEND      CGAL_ESSENTIAL_3RD_PARTY_LIBRARIES GMPXX)
-  endif()
 
   # When CMake is run several times, to avoid duplicates
   list (REMOVE_DUPLICATES CGAL_ESSENTIAL_3RD_PARTY_LIBRARIES)

Installation/cmake/modules/CGAL_SetupGMP.cmake

@@ -23,12 +23,12 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CGAL_MODULES_DIR})
 
 find_package(GMP REQUIRED)
 find_package(MPFR REQUIRED)
+find_package(GMPXX QUIET)
 
-if(NOT DEFINED WITH_GMPXX)
+if(NOT GMPXX_FOUND)
   option(CGAL_WITH_GMPXX "Use CGAL with GMPXX: use C++ classes of GNU MP instead of CGAL wrappers" OFF)
-endif()
+else()
-if(WITH_GMPXX OR CGAL_WITH_GMPXX)
+  option(CGAL_WITH_GMPXX "Use CGAL with GMPXX: use C++ classes of GNU MP instead of CGAL wrappers" ON)
-  find_package(GMPXX REQUIRED)
 endif()
 
 #.rst:

Intersections_3/test/Intersections_3/call_test.cpp

@@ -52,11 +52,13 @@ void call_do_intersect_with_kernel(const A& a, const B& b, const K&) {
 }
 
 
-int main()
+int main(int argc, char**)
 {
   CGAL::Interval_nt_advanced::Protector p;
   CGAL_USE(p);
-  try {
+  //we only want to check compilation
+  if(argc > 666 )
+  {
     call_intersection_global(S(), S());
     call_intersection_global(S(), L());
     call_intersection_global(S(), Pl());
@@ -264,10 +266,6 @@ int main()
     call_do_intersect_with_kernel(Bbox_3(), S(), K());
     call_do_intersect_with_kernel(Bbox_3(), Sph(), K());
     call_do_intersect_with_kernel(Bbox_3(), Tr(), K());
-    // There is no kernel to call
-    // call_do_intersect_with_kernel(Bbox_3(), Bbox_3(), K());
-  } catch(...) {
-    // as long as this test compiles, it is fine
   }
   return EXIT_SUCCESS;
 }

Kernel_23/include/CGAL/Kernel/interface_macros.h

@@ -488,7 +488,7 @@ CGAL_Kernel_pred_RT(Coplanar_orientation_3,
                     coplanar_orientation_3_object)
 CGAL_Kernel_pred_RT(Coplanar_side_of_bounded_circle_3,
                     coplanar_side_of_bounded_circle_3_object)
-CGAL_Kernel_pred(Coplanar_3,
+CGAL_Kernel_pred_RT(Coplanar_3,
 		 coplanar_3_object)
 CGAL_Kernel_pred(Counterclockwise_in_between_2,
 		 counterclockwise_in_between_2_object)

Mesh_3/doc/Mesh_3/CGAL/Mesh_triangulation_3.h

@@ -12,8 +12,8 @@ type to be used for the 3D triangulation embedding the mesh.
 and defaults to `Kernel_traits<MD>::%Kernel`.
 
 \tparam Concurrency_tag enables sequential versus parallel meshing and optimization algorithms.
-                        Possible values are `Sequential_tag` (the default) and
+                        Possible values are `Sequential_tag` (the default), `Parallel_tag`,
-                        `Parallel_tag`.
+                        and `Parallel_if_available_tag`.
 
 \tparam Vertex_base must be a model of `MeshVertexBase_3` or `Default`
 and defaults to `Mesh_vertex_base_3<Gt, MD>`.

Mesh_3/test/Mesh_3/CMakeLists.txt

@@ -65,15 +65,15 @@ if ( CGAL_FOUND )
       test_meshing_verbose
       test_meshing_polyhedron_with_features
       test_meshing_utilities.h
-      test_mesh_implicit_domains
       test_meshing_implicit_function
       test_meshing_3D_image
       test_meshing_3D_gray_image
       test_meshing_unit_tetrahedron
-      test_backward_compatibility
       test_meshing_polyhedron
       test_meshing_polyhedral_complex
+      test_mesh_capsule_var_distance_bound
       test_mesh_3_issue_1554
+      test_mesh_polyhedral_domain_with_features_deprecated
       )
     if(TBB_FOUND AND TARGET ${target})
       CGAL_target_use_TBB(${target})

Mesh_3/test/Mesh_3/test_mesh_3_issue_1554.cpp

@@ -18,11 +18,7 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
 typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;
 
-#ifdef CGAL_CONCURRENT_MESH_3
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
 
 // Triangulation
 typedef CGAL::Mesh_triangulation_3<Mesh_domain,

Mesh_3/test/Mesh_3/test_mesh_capsule_var_distance_bound.cpp

@@ -16,11 +16,7 @@ typedef K::Point_3 Point;
 typedef FT (Function)(const Point&);
 typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
 
-#ifdef CGAL_CONCURRENT_MESH_3
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-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;

Mesh_3/test/Mesh_3/test_mesh_polyhedral_domain_with_features_deprecated.cpp

@@ -14,12 +14,7 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
 typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;
 
-
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-#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,CGAL::Default,Concurrency_tag>::type Tr;

Number_types/include/CGAL/FPU.h

@@ -408,7 +408,7 @@ typedef unsigned int FPU_CW_t;
 
 #elif defined __mips__
 #define CGAL_IA_SETFPCW(CW) asm volatile ("ctc1 %0,$31" : :"r" (CW))
-#define CGAL_IA_GETFPCW(CW) asm volatile ("cfc1 %0,$31" : "=r" (CW))
+#define CGAL_IA_GETFPCW(CW) asm volatile ("cfc1 %0,$31" : "=r" (CW)); CW &= 3
 typedef unsigned int FPU_CW_t;
 #define CGAL_FE_TONEAREST    (0x0)
 #define CGAL_FE_TOWARDZERO   (0x1)

Optimal_transportation_reconstruction_2/demo/Optimal_transportation_reconstruction_2/scene.h

@@ -18,9 +18,11 @@
 #define cimg_display 0 // To avoid X11 or Windows-GDI dependency
 #include <CImg.h>
 #endif
-#include <CGAL/Random.h>
 #include <utility>      // std::pair
 #include <vector>
+
+#include <CGAL/number_type_config.h>
+#include <CGAL/Random.h>
 #include <CGAL/property_map.h>
 #include <CGAL/value_type_traits.h>
 #include <CGAL/compute_average_spacing.h>
@@ -953,7 +955,7 @@ public:
 
   void append_star(const int nb_branches, const int density) {
     std::cerr << "append star...";
-    const double deg_in_rad = 3.1415926535897932384626 / 180.0;
+    const double deg_in_rad = CGAL_PI / 180.0;
     const double incr = 180.0 / nb_branches;
     double angle = 0.0;
     const Point center(0.5, 0.5);
@@ -970,7 +972,7 @@ public:
 
   void append_predefined_increasingly_sharp_angles(const int density,
     const double min_angle) {
-    const double deg_in_rad = 3.1415926535897932384626 / 180.0;
+    const double deg_in_rad = CGAL_PI / 180.0;
     double prev_angle = 0.0;
     double curr_angle = min_angle;
     double incr = min_angle;

Periodic_3_triangulation_3/test/Periodic_3_triangulation_3/include/CGAL/_test_cls_periodic_3_tds_cell.h

@@ -88,11 +88,6 @@ _test_cell_tds_3(const Tds &)
    assert(c2->index(v1)==0);
    assert(c2->index(v2)==1);
    assert(c2->index(v3)==2);
-   //   c2->set_vertices();
-   //   assert(c2->vertex(0)==NULL);
-   //   assert(c2->vertex(1)==NULL);
-   //   assert(c2->vertex(2)==NULL);
-   //   assert(c2->vertex(3)==NULL);
    c2->set_vertices(v0, v1, v2, v3);
    assert(c2->index(v0)==0);
    assert(c2->index(v1)==1);
@@ -107,11 +102,6 @@ _test_cell_tds_3(const Tds &)
    assert(c2->index(n1)==0);
    assert(c2->index(n2)==1);
    assert(c2->index(n3)==2);
-   //   c2->set_neighbors();
-   //   assert(c2->neighbor(0)==NULL);
-   //   assert(c2->neighbor(1)==NULL);
-   //   assert(c2->neighbor(2)==NULL);
-   //   assert(c2->neighbor(3)==NULL);
    c2->set_neighbors(n0, n1, n2, n3);
    assert(c2->index(n0)==0);
    assert(c2->index(n1)==1);

Point_set_processing_3/examples/Point_set_processing_3/average_spacing_example.cpp

@@ -16,12 +16,7 @@ typedef Kernel::Point_3 Point;
 typedef boost::tuple<int, Point, int, int, int> IndexedPointWithColorTuple;
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
-
 
 int main(int argc, char*argv[])
 {

Point_set_processing_3/examples/Point_set_processing_3/bilateral_smooth_point_set_example.cpp

@@ -17,12 +17,7 @@ typedef Kernel::Vector_3 Vector;
 typedef std::pair<Point, Vector> PointVectorPair;
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
-
 
 int main(int argc, char*argv[])
 {

Point_set_processing_3/examples/Point_set_processing_3/callback_example.cpp

@@ -18,12 +18,7 @@ typedef Kernel::Point_3 Point;
 typedef CGAL::Random_points_on_sphere_3<Point> Generator;
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
-
 
 // instance of std::function<bool(double)>
 struct Progress_to_std_cerr_callback

Point_set_processing_3/examples/Point_set_processing_3/edge_aware_upsample_point_set_example.cpp

@@ -16,12 +16,7 @@ typedef Kernel::Vector_3 Vector;
 typedef std::pair<Point, Vector> PointVectorPair;
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
-
 
 int main(int argc, char* argv[])
 {

Point_set_processing_3/examples/Point_set_processing_3/jet_smoothing_example.cpp

@@ -8,12 +8,7 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 typedef Kernel::Point_3 Point;
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
-
 
 int main(void)
 {

Point_set_processing_3/examples/Point_set_processing_3/normal_estimation.cpp

@@ -48,12 +48,7 @@ typedef std::pair<Point, Vector> PointVectorPair;
 typedef std::vector<PointVectorPair> PointList;
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
-
 
 // ----------------------------------------------------------------------------
 // Private functions

Point_set_processing_3/examples/Point_set_processing_3/normals_example.cpp

@@ -18,11 +18,7 @@ typedef Kernel::Vector_3 Vector;
 typedef std::pair<Point, Vector> PointVectorPair;
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
 
 int main(int argc, char*argv[])
 {

Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp

@@ -12,11 +12,7 @@
 #include <fstream>
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
 
 // Types
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;

Point_set_processing_3/examples/Point_set_processing_3/wlop_simplify_and_regularize_point_set_example.cpp

@@ -12,11 +12,7 @@ typedef CGAL::Simple_cartesian<double> Kernel;
 typedef Kernel::Point_3 Point;
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
 
 int main(int argc, char** argv)
 {

Point_set_processing_3/include/CGAL/bilateral_smooth_point_set.h

@@ -16,6 +16,7 @@
 
 #include <CGAL/disable_warnings.h>
 
+#include <CGAL/number_type_config.h>
 #include <CGAL/Search_traits_3.h>
 #include <CGAL/Orthogonal_k_neighbor_search.h>
 #include <CGAL/Point_set_processing_3/internal/neighbor_query.h>
@@ -144,7 +145,7 @@ compute_denoise_projection(
   FT project_weight_sum = FT(0.0);
   Vector normal_sum = CGAL::NULL_VECTOR; 
 
-  FT cos_sigma = cos(sharpness_angle / 180.0 * 3.1415926);
+  FT cos_sigma = cos(sharpness_angle * CGAL_PI / 180.0);
   FT sharpness_bandwidth = std::pow((CGAL::max)(1e-8, 1 - cos_sigma), 2);
 
   typename std::vector<Pwn,CGAL_PSP3_DEFAULT_ALLOCATOR<Pwn> >::const_iterator 
@@ -350,9 +351,8 @@ public:
    \pre Normals must be unit vectors
    \pre k >= 2
 
-   \tparam ConcurrencyTag enables sequential versus parallel algorithm.
+   \tparam ConcurrencyTag enables sequential versus parallel algorithm. Possible values are `Sequential_tag`,
-   Possible values are `Sequential_tag`
+                          `Parallel_tag`, and `Parallel_if_available_tag`.
-   And `Parallel_tag`.
    \tparam PointRange is a model of `Range`. The value type of
    its iterator is the key type of the named parameter `point_map`.
 

Point_set_processing_3/include/CGAL/compute_average_spacing.h

@@ -150,9 +150,8 @@ compute_average_spacing(const typename Kernel::Point_3& query, ///< 3D point who
 
    \pre `k >= 2.`
 
-   \tparam ConcurrencyTag enables sequential versus parallel algorithm.
+   \tparam ConcurrencyTag enables sequential versus parallel algorithm. Possible values are `Sequential_tag`,
-   Possible values are `Sequential_tag`
+                          `Parallel_tag`, and `Parallel_if_available_tag`.
-   and `Parallel_tag`.
    \tparam PointRange is a model of `ConstRange`. The value type of
    its iterator is the key type of the named parameter `point_map`.
 

Point_set_processing_3/include/CGAL/edge_aware_upsample_point_set.h

@@ -281,7 +281,7 @@ update_new_point(
  
    \tparam ConcurrencyTag enables sequential versus parallel versions
    of `compute_average_spacing()` (called internally). Possible
-   values are `Sequential_tag` and `Parallel_tag`.
+   values are `Sequential_tag`, `Parallel_tag`, and `Parallel_if_available_tag`.
    \tparam PointRange is a model of `ConstRange`. The value type of
    its iterator is the key type of the named parameter `point_map`.
    \tparam OutputIterator Type of the output iterator. 

Point_set_processing_3/include/CGAL/jet_estimate_normals.h

@@ -152,9 +152,8 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
 
    \pre `k >= 2`
 
-   \tparam ConcurrencyTag enables sequential versus parallel algorithm.
+   \tparam ConcurrencyTag enables sequential versus parallel algorithm. Possible values are `Sequential_tag`,
-   Possible values are `Sequential_tag`
+                          `Parallel_tag`, and `Parallel_if_available_tag`.
-   and `Parallel_tag`.
    \tparam PointRange is a model of `Range`. The value type of
    its iterator is the key type of the named parameter `point_map`.
 

Point_set_processing_3/include/CGAL/jet_smooth_point_set.h

@@ -156,9 +156,8 @@ jet_smooth_point(
 
    \pre `k >= 2`
 
-   \tparam ConcurrencyTag enables sequential versus parallel algorithm.
+   \tparam ConcurrencyTag enables sequential versus parallel algorithm. Possible values are `Sequential_tag`,
-   Possible values are `Sequential_tag`
+                          `Parallel_tag`, and `Parallel_if_available_tag`.
-   and `Parallel_tag`.
    \tparam PointRange is a model of `Range`. The value type of
    its iterator is the key type of the named parameter `point_map`.
 

Point_set_processing_3/include/CGAL/pca_estimate_normals.h

@@ -142,9 +142,8 @@ pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
 
    \pre `k >= 2`
 
-   \tparam ConcurrencyTag enables sequential versus parallel algorithm.
+   \tparam ConcurrencyTag enables sequential versus parallel algorithm. Possible values are `Sequential_tag`,
-   Possible values are `Sequential_tag`
+                          `Parallel_tag`, and `Parallel_if_available_tag`.
-   and `Parallel_tag`.
    \tparam PointRange is a model of `Range`. The value type of
    its iterator is the key type of the named parameter `point_map`.
 

Point_set_processing_3/include/CGAL/wlop_simplify_and_regularize_point_set.h

@@ -411,9 +411,8 @@ public:
    See the <a href="https://www.threadingbuildingblocks.org/documentation">TBB documentation</a> 
    for more details.
 
-   \tparam ConcurrencyTag enables sequential versus parallel algorithm.
+   \tparam ConcurrencyTag enables sequential versus parallel algorithm. Possible values are `Sequential_tag`,
-   Possible values are `Sequential_tag`
+                          `Parallel_tag`, and `Parallel_if_available_tag`.
-   and `Parallel_tag`.
    \tparam PointRange is a model of `Range`. The value type of
    its iterator is the key type of the named parameter `point_map`.
    \tparam OutputIterator Type of the output iterator. 

Point_set_processing_3/test/Point_set_processing_3/analysis_test.cpp

@@ -33,12 +33,7 @@ typedef Kernel::FT FT;
 typedef Kernel::Point_3 Point;
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
-
 
 // ----------------------------------------------------------------------------
 // Tests

Point_set_processing_3/test/Point_set_processing_3/edge_aware_upsample_test.cpp

@@ -36,12 +36,7 @@ typedef Kernel::Vector_3 Vector;
 typedef std::pair<Point, Vector> PointVectorPair;
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
-
 
 // ----------------------------------------------------------------------------
 // Tests

Point_set_processing_3/test/Point_set_processing_3/normal_estimation_test.cpp

@@ -45,12 +45,7 @@ typedef CGAL::Point_with_normal_3<Kernel> Point_with_normal; // position + norma
 typedef std::vector<Point_with_normal> PointList;
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
-
 
 // ----------------------------------------------------------------------------
 // Tests

Point_set_processing_3/test/Point_set_processing_3/smoothing_test.cpp

@@ -36,11 +36,7 @@ typedef Kernel::Point_3 Point;
 typedef Kernel::Vector_3 Vector;
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
 
 // ----------------------------------------------------------------------------
 // Tests

Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt

@@ -412,10 +412,13 @@ one incident non-null face) after each hole filling step.
 Holes are filled one after the other, and the process stops when there is no border edge left.
 
 This process is illustrated by the example below, where holes are
-iteratively filled, refined and faired to get a faired mesh with no hole.
+iteratively filled, refined and faired. Optionally, only holes
+not exceeding a certain diameter or number of edges can be filled.
+This example assumes that the mesh is stored in a `Surface_mesh`
+datastructure. Analogous examples when using the `Polyhedron_3`
+class and a few others are part of the code base.
 
-
+\cgalExample{Polygon_mesh_processing/hole_filling_example_SM.cpp}
-\cgalExample{Polygon_mesh_processing/hole_filling_example.cpp}
 
  \cgalFigureBegin{Triangulated_fork, fork.jpg}
  Holes in the fork model are filled with triangle patches.

Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt

@@ -2,6 +2,7 @@
 
 \example Polygon_mesh_processing/self_intersections_example.cpp
 \example Polygon_mesh_processing/hole_filling_example.cpp
+\example Polygon_mesh_processing/hole_filling_example_SM.cpp
 \example Polygon_mesh_processing/stitch_borders_example.cpp
 \example Polygon_mesh_processing/compute_normals_example.cpp
 \example Polygon_mesh_processing/point_inside_example.cpp

Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt

@@ -109,6 +109,7 @@ endif(OpenMesh_FOUND)
 
 find_package( TBB )
 if( TBB_FOUND )
+  CGAL_target_use_TBB(self_intersections_example)
   CGAL_target_use_TBB(hausdorff_distance_remeshing_example)
 else()
   message( STATUS "NOTICE: Intel TBB was not found. Sequential code will be used." )

Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_LCC.cpp

@@ -34,6 +34,7 @@ int main(int argc, char* argv[])
             << num_vertices(mesh2) << "\n";
 
   std::ofstream output("mesh1_refined.off");
+  output.precision(17);
   CGAL::write_off(output, mesh1);
   output.close();
   output.open("mesh2_refined.off");

Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_SM.cpp

@@ -41,6 +41,7 @@ int main(int argc, char* argv[])
             << num_vertices(mesh2) << "\n";
 
   std::ofstream output("mesh1_refined.off");
+  output.precision(17);
   output << mesh1;
   output.close();
   output.open("mesh2_refined.off");

Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_consecutive_bool_op.cpp

@@ -136,6 +136,7 @@ int main(int argc, char* argv[])
     {
       std::cout << "Intersection and union were successfully computed\n";
       std::ofstream output("inter_union.off");
+      output.precision(17);
       output << mesh2;
       return 0;
     }

Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_difference_remeshed.cpp

@@ -67,6 +67,7 @@ int main(int argc, char* argv[])
   {
     std::cout << "Difference was successfully computed\n";
     std::ofstream output("difference.off");
+    output.precision(17);
     output << mesh1;
   }
   else{
@@ -111,6 +112,7 @@ int main(int argc, char* argv[])
     params::edge_is_constrained_map(is_constrained_map) );
 
   std::ofstream output("difference_remeshed.off");
+  output.precision(17);
   output << mesh1;
 
   return 0;

Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union.cpp

@@ -37,6 +37,7 @@ int main(int argc, char* argv[])
   {
     std::cout << "Union was successfully computed\n";
     std::ofstream output("union.off");
+    output.precision(17);
     output << out;
     return 0;
   }

Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union_and_intersection.cpp

@@ -55,6 +55,7 @@ int main(int argc, char* argv[])
   {
     std::cout << "Union was successfully computed\n";
     std::ofstream output("union.off");
+    output.precision(17);
     output << out_union;
   }
   else
@@ -64,6 +65,7 @@ int main(int argc, char* argv[])
   {
     std::cout << "Intersection was successfully computed\n";
     std::ofstream output("intersection.off");
+    output.precision(17);
     output << out_intersection;
   }
   else

Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union_with_attributes.cpp

@@ -94,6 +94,7 @@ int main(int argc, char* argv[])
   {
     std::cout << "Union was successfully computed\n";
     std::ofstream output("union.off");
+    output.precision(17);
     output << out;
     return 0;
   }

Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_polyhedron_union.cpp

@@ -39,6 +39,7 @@ int main(int argc, char* argv[])
   {
     std::cout << "Union was successfully computed\n";
     std::ofstream output("union.off");
+    output.precision(17);
     output << out;
     return 0;
   }

Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_distance_remeshing_example.cpp

@@ -3,11 +3,7 @@
 #include <CGAL/Polygon_mesh_processing/distance.h>
 #include <CGAL/Polygon_mesh_processing/remesh.h>
 
-#if defined(CGAL_LINKED_WITH_TBB)
+#define TAG CGAL::Parallel_if_available_tag
-#define TAG CGAL::Parallel_tag
-#else
-#define TAG CGAL::Sequential_tag
-#endif
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 typedef K::Point_3                     Point;

Polygon_mesh_processing/examples/Polygon_mesh_processing/hole_filling_example_SM.cpp

@@ -2,10 +2,14 @@
 #include <CGAL/Surface_mesh.h>
 
 #include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
+#include <CGAL/Polygon_mesh_processing/border.h>
+
+#include <boost/lexical_cast.hpp>
 
 #include <iostream>
 #include <fstream>
 #include <vector>
+#include <set>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 typedef Kernel::Point_3 Point;
@@ -15,23 +19,58 @@ typedef boost::graph_traits<Mesh>::halfedge_descriptor    halfedge_descriptor;
 typedef boost::graph_traits<Mesh>::face_descriptor       face_descriptor;
 typedef boost::graph_traits<Mesh>::vertex_descriptor     vertex_descriptor;
 
+bool is_small_hole(halfedge_descriptor h, Mesh & mesh,
+                   double max_hole_diam, int max_num_hole_edges)
+{
+  int num_hole_edges = 0;
+  CGAL::Bbox_3 hole_bbox;
+  for (halfedge_descriptor hc : CGAL::halfedges_around_face(h, mesh))
+  {
+    const Point& p = mesh.point(target(hc, mesh));
+
+    hole_bbox += p.bbox();
+    ++num_hole_edges;
+
+    // Exit early, to avoid unnecessary traversal of large holes
+    if (num_hole_edges > max_num_hole_edges) return false;
+    if (hole_bbox.xmax() - hole_bbox.xmin() > max_hole_diam) return false;
+    if (hole_bbox.ymax() - hole_bbox.ymin() > max_hole_diam) return false;
+    if (hole_bbox.zmax() - hole_bbox.zmin() > max_hole_diam) return false;
+  }
+
+  return true;
+}
+
+// Incrementally fill the holes that are no larger than given diameter
+// and with no more than a given number of edges (if specified).
+
 int main(int argc, char* argv[])
 {
   const char* filename = (argc > 1) ? argv[1] : "data/mech-holes-shark.off";
-  std::ifstream input(filename);
 
+  // Both of these must be positive in order to be considered
+  double max_hole_diam   = (argc > 2) ? boost::lexical_cast<double>(argv[2]): -1.0;
+  int max_num_hole_edges = (argc > 3) ? boost::lexical_cast<int>(argv[3]) : -1;
+
+  std::ifstream input(filename);
   Mesh mesh;
   if ( !input || !(input >> mesh) ) {
     std::cerr << "Not a valid off file." << std::endl;
     return 1;
   }
 
-  // Incrementally fill the holes
   unsigned int nb_holes = 0;
-  for(halfedge_descriptor h : halfedges(mesh))
+  std::vector<halfedge_descriptor> border_cycles;
-  {
+
-    if(is_border(h,mesh))
+  // collect one halfedge per boundary cycle
+  CGAL::Polygon_mesh_processing::extract_boundary_cycles(mesh, std::back_inserter(border_cycles));
+
+  for(halfedge_descriptor h : border_cycles)
   {
+    if(max_hole_diam > 0 && max_num_hole_edges > 0 &&
+       !is_small_hole(h, mesh, max_hole_diam, max_num_hole_edges))
+      continue;
+
     std::vector<face_descriptor>  patch_facets;
     std::vector<vertex_descriptor> patch_vertices;
     bool success = std::get<0>(
@@ -39,21 +78,20 @@ int main(int argc, char* argv[])
                 mesh,
                 h,
                 std::back_inserter(patch_facets),
-                  std::back_inserter(patch_vertices),
+                std::back_inserter(patch_vertices)) );
-     CGAL::Polygon_mesh_processing::parameters::vertex_point_map(get(CGAL::vertex_point, mesh)).
-                  geom_traits(Kernel())) );
 
     std::cout << "* Number of facets in constructed patch: " << patch_facets.size() << std::endl;
     std::cout << "  Number of vertices in constructed patch: " << patch_vertices.size() << std::endl;
     std::cout << "  Is fairing successful: " << success << std::endl;
-      nb_holes++;
+    ++nb_holes;
-    }
   }
 
   std::cout << std::endl;
   std::cout << nb_holes << " holes have been filled" << std::endl;
 
-  std::ofstream out("filled_SM.off");
+  std::string outfile = "filled_SM.off";
+  std::ofstream out(outfile.c_str());
+  std::cout << "Mesh written to: " << outfile << std::endl;
   out.precision(17);
   out << mesh << std::endl;
   return 0;

Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_of_patch_example.cpp

@@ -74,6 +74,7 @@ int main(int argc, char* argv[])
                            );
 
   std::ofstream out("out.off");
+  out.precision(17);
   out << mesh;
   std::cout << "Remeshing done." << std::endl;
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_smoothing_example.cpp

@@ -48,6 +48,7 @@ int main(int argc, char** argv)
                                          .edge_is_constrained_map(eif));
 
   std::ofstream output("mesh_smoothed.off");
+  output.precision(17);
   output << mesh;
 
   std::cout << "Done!" << std::endl;

Polygon_mesh_processing/examples/Polygon_mesh_processing/orient_polygon_soup_example.cpp

@@ -44,11 +44,13 @@ int main(int argc, char* argv[])
     CGAL::Polygon_mesh_processing::orient_to_bound_a_volume(mesh);
 
   std::ofstream out("tet-oriented1.off");
+  out.precision(17);
   out << mesh;
   out.close();
 
   CGAL::Polygon_mesh_processing::reverse_face_orientations(mesh);
   std::ofstream out2("tet-oriented2.off");
+  out2.precision(17);
   out2 << mesh;
   out2.close();
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/random_perturbation_SM_example.cpp

@@ -33,6 +33,7 @@ int main(int argc, char* argv[])
     PMP::parameters::vertex_point_map(mesh.points()).geom_traits(K()));
 
   std::ofstream out("data/eight_perturbed.off");
+  out.precision(17);
   out << mesh;
   out.close();
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/refine_fair_example.cpp

@@ -62,6 +62,7 @@ int main(int argc, char* argv[])
                   CGAL::Polygon_mesh_processing::parameters::density_control_factor(2.));
 
   std::ofstream refined_off("refined.off");
+  refined_off.precision(17);
   refined_off << poly;
   refined_off.close();
   std::cout << "Refinement added " << new_vertices.size() << " vertices." << std::endl;
@@ -75,6 +76,7 @@ int main(int argc, char* argv[])
   std::cout << "Fairing : " << (success ? "succeeded" : "failed") << std::endl;
 
   std::ofstream faired_off("faired.off");
+  faired_off.precision(17);
   faired_off << poly;
   faired_off.close();
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/self_intersections_example.cpp

@@ -2,7 +2,10 @@
 #include <CGAL/Surface_mesh.h>
 
 #include <CGAL/Polygon_mesh_processing/self_intersections.h>
+#include <CGAL/Real_timer.h>
+#include <CGAL/tags.h>
 
+#include <iostream>
 #include <fstream>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
@@ -20,20 +23,25 @@ int main(int argc, char* argv[])
   if (!input || !(input >> mesh) || !CGAL::is_triangle_mesh(mesh))
   {
     std::cerr << "Not a valid input file." << std::endl;
-    return 1;
+    return EXIT_FAILURE;
   }
 
-  bool intersecting = PMP::does_self_intersect(mesh,
+  std::cout << "Using parallel mode? " << std::is_same<CGAL::Parallel_if_available_tag, CGAL::Parallel_tag>::value << std::endl;
-      PMP::parameters::vertex_point_map(get(CGAL::vertex_point, mesh)));
 
-  std::cout
+  CGAL::Real_timer timer;
-    << (intersecting ? "There are self-intersections." : "There is no self-intersection.")
+  timer.start();
-    << std::endl;
+
+  bool intersecting = PMP::does_self_intersect<CGAL::Parallel_if_available_tag>(mesh, CGAL::parameters::vertex_point_map(get(CGAL::vertex_point, mesh)));
+  std::cout << (intersecting ? "There are self-intersections." : "There is no self-intersection.") << std::endl;
+  std::cout << "Elapsed time (does self intersect): " << timer.time() << std::endl;
+
+  timer.reset();
 
   std::vector<std::pair<face_descriptor, face_descriptor> > intersected_tris;
-  PMP::self_intersections(mesh, std::back_inserter(intersected_tris));
+  PMP::self_intersections<CGAL::Parallel_if_available_tag>(faces(mesh), mesh, std::back_inserter(intersected_tris));
-
   std::cout << intersected_tris.size() << " pairs of triangles intersect." << std::endl;
 
-  return 0;
+  std::cout << "Elapsed time (self intersections): " << timer.time() << std::endl;
+
+  return EXIT_SUCCESS;
 }

Polygon_mesh_processing/examples/Polygon_mesh_processing/shape_smoothing_example.cpp

@@ -41,6 +41,7 @@ int main(int argc, char* argv[])
                                                 .vertex_is_constrained_map(vcmap));
 
   std::ofstream output("mesh_shape_smoothed.off");
+  output.precision(17);
   output << mesh;
 
   std::cout << "Done!" << std::endl;

Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example.cpp

@@ -34,6 +34,7 @@ int main(int argc, char* argv[])
   std::cout << "\t Number of facets    :\t" << mesh.size_of_facets() << std::endl;
 
   std::ofstream output("mesh_stitched.off");
+  output.precision(17);
   output << std::setprecision(17) << mesh;
 
   return 0;

Polygon_mesh_processing/examples/Polygon_mesh_processing/surface_mesh_intersection.cpp

@@ -43,6 +43,7 @@ int main(int argc, char* argv[])
 
   //dump polylines
   std::ofstream output("intersection_polylines.cgal");
+  output.precision(17);
   for(const std::vector<Point>& polyline : polylines)
   {
     output << polyline.size() << " ";

Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_example.cpp

@@ -32,6 +32,7 @@ int main(int argc, char* argv[])
       std::cerr << "Error: non-triangular face left in mesh." << std::endl;
 
   std::ofstream cube_off(outfilename);
+  cube_off.precision(17);
   cube_off << mesh;
 
   return 0;

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h

@@ -266,7 +266,7 @@ clip_to_bbox(const Plane_3& plane,
   * clips `tm` by keeping the part that is inside the volume \link coref_def_subsec bounded \endlink
   * by `clipper`.
   * If `tm` is closed, the clipped part can be closed too if the named parameter `clip_volume` is set to `true`.
-  * See subsection \ref coref_clip for more details.
+  * See Subsection \ref coref_clip for more details.
   * \attention With the current implementation, `clipper` will be modified (refined with the intersection with `tm`).
   *
   * \pre \link CGAL::Polygon_mesh_processing::does_self_intersect() `!CGAL::Polygon_mesh_processing::does_self_intersect(tm1)` \endlink
@@ -305,7 +305,7 @@ clip_to_bbox(const Plane_3& plane,
   *   \cgalParamEnd
   *   \cgalParamBegin{clip_volume} if `true` and `tm` is closed, the clipping will be done on
   *      the volume \link coref_def_subsec bounded \endlink by `tm` rather than on its surface
-  *      (i.e. `tm` will be kept closed).
+  *      (i.e., `tm` will be kept closed).
   *   \cgalParamEnd
   *   \cgalParamBegin{use_compact_clipper} if `false` and `clip_volume` is `false` and `tm` is open, the parts of `tm` coplanar with `clipper`
   *                                        will not be part of the output.
@@ -358,9 +358,10 @@ bool dispatch_clip_call(TriangleMesh& tm, TriangleMesh& clipper,
   * \ingroup PMP_corefinement_grp
   * clips `tm` by keeping the part that is on the negative side of `plane` (side opposite to its normal vector).
   * If `tm` is closed, the clipped part can be closed too if the named parameter `clip_volume` is set to `true`.
-  * See subsection \ref coref_clip for more details.
+  * See Subsection \ref coref_clip for more details.
   *
-  * \note In the current implementation it is not possible to set the vertex point map and the default will be used.
+  * \note In the current implementation it is not possible to set the vertex point map and the default will be used. `Plane_3` must be
+  * from the same %Kernel as the point of the vertex point map.
   * \pre \link CGAL::Polygon_mesh_processing::does_self_intersect() `!CGAL::Polygon_mesh_processing::does_self_intersect(tm)` \endlink
   *
   * @tparam TriangleMesh a model of `MutableFaceGraph`, `HalfedgeListGraph` and `FaceListGraph`.
@@ -386,7 +387,7 @@ bool dispatch_clip_call(TriangleMesh& tm, TriangleMesh& clipper,
   *   \cgalParamEnd
   *   \cgalParamBegin{clip_volume} if `true` and `tm` is closed, the clipping will be done on
   *      the volume \link coref_def_subsec bounded \endlink by `tm` rather than on its surface
-  *      (i.e. `tm` will be kept closed).
+  *      (i.e., `tm` will be kept closed).
   *   \cgalParamEnd
   *   \cgalParamBegin{use_compact_clipper} if `false` and `clip_volume` is `false` and `tm` is open, the parts of `tm` coplanar with `plane`
   *                                        will not be part of the output.
@@ -432,6 +433,70 @@ bool clip(      TriangleMesh& tm,
                                       np, CGAL::graph_has_property<TriangleMesh, CGAL::face_index_t>());
 }
 
+/**
+  * \ingroup PMP_corefinement_grp
+  * clips `tm` by keeping the part that is inside `iso_cuboid`.
+  * If `tm` is closed, the clipped part can be closed too if the named parameter `clip_volume` is set to `true`.
+  * See Subsection \ref coref_clip for more details.
+  *
+  * \note In the current implementation it is not possible to set the vertex point map and the default will be used. `Iso_cuboid_3` must be
+  * from the same %Kernel as the point of the vertex point map.
+  * \pre \link CGAL::Polygon_mesh_processing::does_self_intersect() `!CGAL::Polygon_mesh_processing::does_self_intersect(tm)` \endlink
+  *
+  * @tparam TriangleMesh a model of `MutableFaceGraph`, `HalfedgeListGraph` and `FaceListGraph`.
+  *                      If `TriangleMesh` has an internal property map for `CGAL::face_index_t`,
+  *                      as a named parameter, then it must be initialized.
+  *                      An internal property map for `CGAL::vertex_point_t` must be available.
+  *
+  * @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
+  *
+  * @param tm input triangulated surface mesh
+  * @param iso_cuboid iso-cuboid used to clip `tm`.
+  * @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
+  *
+  * \cgalNamedParamsBegin
+  *   \cgalParamBegin{visitor} a class model of `PMPCorefinementVisitor`
+  *                            that is used to track the creation of new faces.
+  *   \cgalParamEnd
+  *   \cgalParamBegin{throw_on_self_intersection} if `true`,
+  *      the set of triangles closed to the intersection of `tm` and `iso_cuboid` will be
+  *      checked for self-intersections and `CGAL::Polygon_mesh_processing::Corefinement::Self_intersection_exception`
+  *      will be thrown if at least one is found.
+  *   \cgalParamEnd
+  *   \cgalParamBegin{clip_volume} if `true` and `tm` is closed, the clipping will be done on
+  *      the volume \link coref_def_subsec bounded \endlink by `tm` rather than on its surface
+  *      (i.e., `tm` will be kept closed).
+  *   \cgalParamEnd
+  *   \cgalParamBegin{use_compact_clipper} if `false` and `clip_volume` is `false` and `tm` is open, the parts of `tm` coplanar with `is_cuboid`
+  *                                        will not be part of the output.
+  * \cgalNamedParamsEnd
+  *
+  * @return `true` if the output surface mesh is manifold.
+  *         If `false` is returned `tm` is only refined by the intersection with `iso_cuboid`.
+  */
+template <class TriangleMesh,
+          class NamedParameters>
+bool clip(      TriangleMesh& tm,
+          #ifdef DOXYGEN_RUNNING
+          const Iso_cuboid_3& iso_cuboid,
+          #else
+          const typename GetGeomTraits<TriangleMesh, NamedParameters>::type::Iso_cuboid_3& iso_cuboid,
+          #endif
+          const NamedParameters& np)
+{
+  if( boost::begin(faces(tm))==boost::end(faces(tm)) ) return true;
+  TriangleMesh clipper;
+
+  make_hexahedron(iso_cuboid[0], iso_cuboid[1], iso_cuboid[2], iso_cuboid[3],
+                  iso_cuboid[4], iso_cuboid[5], iso_cuboid[6], iso_cuboid[7],
+                  clipper);
+  triangulate_faces(clipper);
+
+  // dispatch is needed because face index map for tm and clipper have to be of the same time
+  return internal::dispatch_clip_call(tm, clipper,
+                                      np, CGAL::graph_has_property<TriangleMesh, CGAL::face_index_t>());
+}
+
 /// \cond SKIP_IN_MANUAL
 
 // convenience overloads
@@ -442,6 +507,14 @@ bool clip(      TriangleMesh& tm,
   return clip(tm, plane, parameters::all_default());
 }
 
+// convenience overloads
+template <class TriangleMesh>
+bool clip(      TriangleMesh& tm,
+          const typename GetGeomTraits<TriangleMesh>::type::Iso_cuboid_3& iso_cuboid)
+{
+  return clip(tm, iso_cuboid, parameters::all_default());
+}
+
 // convenience overload
 template <class TriangleMesh,
           class NamedParameters1>

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h

@@ -27,6 +27,7 @@
 #include <boost/graph/graph_traits.hpp>
 
 #include <iostream>
+#include <limits>
 #include <utility>
 #include <vector>
 
@@ -264,6 +265,8 @@ bool does_enclose_other_normals(const std::size_t i, const std::size_t j, const
       continue;
 
     const Vector_ref nl = get(face_normals, incident_faces[l]);
+    if(nl == CGAL::NULL_VECTOR)
+      continue;
 
     // this is a bound on how much the scalar product between (v1,v2) and (v1, v3) can change
     // when the angle changes theta_bound := 0.01°
@@ -271,6 +274,8 @@ bool does_enclose_other_normals(const std::size_t i, const std::size_t j, const
     // with theta_j - theta_i = theta_bound
     const FT sp_diff_bound = nbn * 0.00017453292431333;
     const FT sp_bl = sp(nb, nl);
+
+    // norm of nl is 1 by construction
     if(CGAL::abs(sp_bi - sp_bl) <= sp_diff_bound)
       continue;
 
@@ -308,18 +313,18 @@ typename GT::Vector_3 compute_normals_bisector(const typename GT::Vector_3& ni,
 
   Vector_3 nb = cv_3(CGAL::NULL_VECTOR);
 
-  if(almost_equal(ni, nj, traits))
+  if(almost_equal(ni, nj, traits) || nk == CGAL::NULL_VECTOR)
   {
     if(almost_equal(nj, nk, traits))
       nb = ni;
     else // ni == nj, but nij != nk
       nb = compute_normals_bisector(nj, nk, traits);
   }
-  else if(almost_equal(ni, nk, traits)) // ni != nj
+  else if(almost_equal(ni, nk, traits) || nj == CGAL::NULL_VECTOR) // ni != nj
   {
     nb = compute_normals_bisector(nj, nk, traits);
   }
-  else if(almost_equal(nj, nk, traits)) // ni != nj, ni != nk
+  else if(almost_equal(nj, nk, traits) || ni == CGAL::NULL_VECTOR) // ni != nj, ni != nk
   {
     nb = compute_normals_bisector(ni, nk, traits);
   }
@@ -342,7 +347,7 @@ typename GT::Vector_3 compute_normals_bisector(const typename GT::Vector_3& ni,
       return csv_3(csv_3(cslv_3(ni, third), cslv_3(nj, third)), cslv_3(nk, third));
     }
 
-    nb = cv_3(CGAL::ORIGIN, c); // note that this isn't normalized
+    nb = cv_3(CGAL::ORIGIN, c); // not normalized
   }
 
   return nb;
@@ -382,9 +387,8 @@ compute_most_visible_normal_2_points(std::vector<typename boost::graph_traits<Po
       if(traits.equal_3_object()(nb, CGAL::NULL_VECTOR))
         return CGAL::NULL_VECTOR;
 
-      const FT sp_bi = sp_3(nb, ni);
+      FT sp_bi = sp_3(nb, ni);
-
+      sp_bi = (std::max)(FT(0), sp_bi);
-      CGAL_assertion(sp_bi >= 0);
       if(sp_bi <= min_sp)
         continue;
 
@@ -428,6 +432,9 @@ compute_most_visible_normal_3_points(const std::vector<typename boost::graph_tra
         const Vector_ref nj = get(face_normals, incident_faces[j]);
         const Vector_ref nk = get(face_normals, incident_faces[k]);
 
+        if(ni == CGAL::NULL_VECTOR || nj == CGAL::NULL_VECTOR || nk == CGAL::NULL_VECTOR)
+          continue;
+
         Vector_3 nb = compute_normals_bisector(ni, nj, nk, traits);
         if(traits.equal_3_object()(nb, CGAL::NULL_VECTOR))
           return nb;
@@ -487,7 +494,12 @@ compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits<Polyg
   if(res != CGAL::NULL_VECTOR) // found a valid normal through 2 point min circle
     return res;
 
-  CGAL_assertion(incident_faces.size() > 2);
+  // The vertex has only two incident faces with opposite normals (fold)...
+  // @todo devise something based on the directions of the 2/3/4 incident edges?
+  if(incident_faces.size() == 2 && res == CGAL::NULL_VECTOR)
+    return res;
+
+  CGAL_assertion(incident_faces.size() >= 2);
 
   return compute_most_visible_normal_3_points<PolygonMesh>(incident_faces, face_normals, traits);
 }
@@ -531,10 +543,14 @@ compute_vertex_normal_as_sum_of_weighted_normals(typename boost::graph_traits<Po
         const Vector_3 v1 = cv_3(get(vpmap, v), get(vpmap, source(h, pmesh)));
         const Vector_3 v2 = cv_3(get(vpmap, v), get(vpmap, target(next(h, pmesh), pmesh)));
 
-        //v(i) and v(i+1) must me seen in ccw order, from v, so we reverse v1 and v2
+        //v(i) and v(i+1) must be seen in ccw order, from v, so we reverse v1 and v2
         Vector_3 n = traits.construct_cross_product_vector_3_object()(v2, v1);
-        n = traits.construct_scaled_vector_3_object()(n, CGAL::approximate_sqrt(FT(1)/(csl_3(v1) * csl_3(v2))));
+        const FT den = CGAL::approximate_sqrt(csl_3(v1) * csl_3(v2));
 
+        if(den == FT(0))
+          return compute_vertex_normal_as_sum_of_weighted_normals(v, NO_WEIGHT, face_normals, vpmap, pmesh, traits);
+
+        n = traits.construct_scaled_vector_3_object()(n, FT(1) / den);
         normal = traits.construct_sum_of_vectors_3_object()(normal, n);
       }
       else

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/connected_components.h

@@ -23,8 +23,6 @@
 #include <CGAL/boost/graph/Named_function_parameters.h>
 #include <CGAL/boost/graph/helpers.h>
 #include <boost/graph/graph_traits.hpp>
-#include <boost/graph/filtered_graph.hpp>
-#include <boost/graph/connected_components.hpp>
 #include <boost/property_map/vector_property_map.hpp>
 
 #include <CGAL/assertions.h>
@@ -206,6 +204,10 @@ connected_components(const PolygonMesh& pmesh,
   using parameters::choose_parameter;
   using parameters::get_parameter;
 
+  typedef boost::graph_traits<PolygonMesh> GT;
+  typedef typename GT::halfedge_descriptor halfedge_descriptor;
+  typedef typename GT::face_descriptor face_descriptor;
+
   typedef typename internal_np::Lookup_named_param_def <
     internal_np::edge_is_constrained_t,
     NamedParameters,
@@ -215,22 +217,40 @@ connected_components(const PolygonMesh& pmesh,
     = choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
                    internal::No_constraint<PolygonMesh>());
 
-  typedef Dual<PolygonMesh>                              Dual;
+  typedef typename GetFaceIndexMap<PolygonMesh, NamedParameters>::const_type  FaceIndexMap;
-  typedef boost::filtered_graph<Dual,
+  FaceIndexMap fimap = choose_parameter(get_parameter(np, internal_np::face_index),
-    internal::No_border<PolygonMesh,EdgeConstraintMap> > FiniteDual;
-  Dual dual(pmesh);
-
-  FiniteDual finite_dual(dual,
-    internal::No_border<PolygonMesh, EdgeConstraintMap>(pmesh, ecmap));
-
-  typename GetFaceIndexMap<PolygonMesh, NamedParameters>::const_type
-    fimap = choose_parameter(get_parameter(np, internal_np::face_index),
                                         get_const_property_map(boost::face_index, pmesh));
 
-  return boost::connected_components(finite_dual,
+  typename boost::property_traits<FaceComponentMap>::value_type i=0;
-    fcm,
+  std::vector<bool> handled(num_faces(pmesh), false);
-    boost::vertex_index_map(fimap)
+  for (face_descriptor f : faces(pmesh))
-  );
+  {
+    if (handled[get(fimap,f)]) continue;
+    std::vector<face_descriptor> queue;
+    queue.push_back(f);
+    while(!queue.empty())
+    {
+      face_descriptor fq = queue.back();
+      queue.pop_back();
+      typename boost::property_traits<FaceIndexMap>::value_type  fq_id = get(fimap,fq);
+      if ( handled[fq_id]) continue;
+      handled[fq_id]=true;
+      put(fcm, fq, i);
+      for (halfedge_descriptor h : halfedges_around_face(halfedge(fq, pmesh), pmesh))
+      {
+        if ( get(ecmap, edge(h, pmesh)) ) continue;
+        halfedge_descriptor opp = opposite(h, pmesh);
+        face_descriptor fqo = face(opp, pmesh);
+        if ( fqo != GT::null_face() )
+        {
+          if ( !handled[get(fimap,fqo)] )
+            queue.push_back(fqo);
+        }
+      }
+    }
+    ++i;
+  }
+  return i;
 }
 
 template <typename PolygonMesh, typename FaceComponentMap>

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h

@@ -589,8 +589,7 @@ double approximate_Hausdorff_distance(
  * for more details.
  *
  * @tparam Concurrency_tag enables sequential versus parallel algorithm.
- *                         Possible values are `Sequential_tag`
+ *                         Possible values are `Sequential_tag`, `Parallel_tag`, and `Parallel_if_available_tag`.
- *                         and `Parallel_tag`.
  * @tparam TriangleMesh a model of the concept `FaceListGraph`
  * @tparam NamedParameters1 a sequence of \ref pmp_namedparameters "Named Parameters" for `tm1`
  * @tparam NamedParameters2 a sequence of \ref pmp_namedparameters "Named Parameters" for `tm2`

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/AABB_traversal_traits_with_transformation.h

@@ -26,7 +26,6 @@
 #include <CGAL/Filtered_predicate.h>
 #include <CGAL/Simple_cartesian.h>
 
-#include <CGAL/Box_intersection_d/Box_with_info_d.h>
 #include <CGAL/Aff_transformation_3.h>
 #include <boost/mpl/if.hpp>
 

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_callbacks.h

@@ -15,7 +15,7 @@
 
 #include <CGAL/license/Polygon_mesh_processing/corefinement.h>
 
-
+#include <CGAL/Box_intersection_d/Box_with_info_d.h>
 #include <CGAL/property_map.h>
 #include <CGAL/enum.h>
 #include <CGAL/Polygon_mesh_processing/self_intersections.h>
@@ -38,7 +38,9 @@ protected:
   typedef boost::graph_traits<TriangleMesh> Graph_traits;
   typedef typename Graph_traits::face_descriptor face_descriptor;
   typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor;
-  typedef typename CGAL::Box_intersection_d::Box_with_info_d<double, 3, halfedge_descriptor> Box;
+
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, halfedge_descriptor, Box_policy> Box;
 
 public:
   Collect_face_bbox_per_edge_bbox(
@@ -80,9 +82,12 @@ protected:
   typedef boost::graph_traits<TriangleMesh> Graph_traits;
   typedef typename Graph_traits::face_descriptor face_descriptor;
   typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor;
-  typedef typename CGAL::Box_intersection_d::Box_with_info_d<double, 3, halfedge_descriptor> Box;
+
   typedef typename boost::property_traits<VertexPointMap>::reference Point;
 
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, halfedge_descriptor, Box_policy> Box;
+
 public:
   Collect_face_bbox_per_edge_bbox_with_coplanar_handling(
     const TriangleMesh& tm_faces,
@@ -162,7 +167,10 @@ protected:
   typedef typename Graph_traits::face_descriptor face_descriptor;
   typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor;
   typedef typename Graph_traits::vertex_descriptor vertex_descriptor;
-  typedef typename CGAL::Box_intersection_d::Box_with_info_d<double, 3, halfedge_descriptor> Box;
+
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, halfedge_descriptor, Box_policy> Box;
+
   typedef typename boost::property_traits<VertexPointMap>::reference Point;
 
   bool is_edge_target_incident_to_face(halfedge_descriptor hd,

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h

@@ -18,7 +18,6 @@
 
 #include <boost/graph/graph_traits.hpp>
 #include <CGAL/box_intersection_d.h>
-#include <CGAL/Box_intersection_d/Box_with_info_d.h>
 #include <CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_callbacks.h>
 #include <CGAL/Polygon_mesh_processing/internal/Corefinement/Intersection_type.h>
 #include <CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_of_coplanar_triangles_3.h>
@@ -156,7 +155,8 @@ class Intersection_of_triangle_meshes
   typedef typename graph_traits::halfedge_descriptor halfedge_descriptor;
   typedef typename graph_traits::vertex_descriptor vertex_descriptor;
 
-  typedef typename CGAL::Box_intersection_d::Box_with_info_d<double, 3, halfedge_descriptor> Box;
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, halfedge_descriptor, Box_policy> Box;
 
   typedef boost::unordered_set<face_descriptor> Face_set;
   typedef boost::unordered_map<edge_descriptor, Face_set> Edge_to_faces;
@@ -250,11 +250,12 @@ class Intersection_of_triangle_meshes
         if (callback_si.self_intersections_found())
          throw Self_intersection_exception();
     }
-    else
+    else {
       CGAL::box_intersection_d( face_boxes_ptr.begin(), face_boxes_ptr.end(),
                                 edge_boxes_ptr.begin(), edge_boxes_ptr.end(),
                                 callback, cutoff );
     }
+  }
 
   // for autorefinement
   void filter_intersections(const TriangleMesh& tm,

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/repair_extra.h

@@ -110,7 +110,8 @@ void collect_close_stitchable_boundary_edges(PM& pm,
   typedef boost::unordered_map<halfedge_descriptor, int> Halfedge_multiplicity;
   typedef std::vector<std::pair<halfedge_descriptor, halfedge_descriptor> > Halfedge_pairs;
 
-  typedef typename Box_intersection_d::Box_with_info_d<double, 3, edge_descriptor> Box;
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, edge_descriptor, Box_policy> Box;
 
   typedef Union_find<vertex_descriptor> UF_vertices;
   typedef std::map<vertex_descriptor, typename UF_vertices::handle> Handle_map;
@@ -141,7 +142,6 @@ void collect_close_stitchable_boundary_edges(PM& pm,
   for(Box& b : boxes)
     box_ptrs.push_back(&b);
 
-
   Halfedge_multiplicity multiplicity;
   Halfedge_pairs matching_hedges;
 

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h

@@ -17,20 +17,27 @@
 
 #include <CGAL/disable_warnings.h>
 
-#include <CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h>
-#include <boost/type_traits/is_same.hpp>
-#include <boost/utility/enable_if.hpp>
-#include <CGAL/Polygon_mesh_processing/bbox.h>
-#include <CGAL/boost/iterator/counting_iterator.hpp>
-#include <boost/mpl/if.hpp>
-#include <CGAL/Polygon_mesh_processing/bbox.h>
-#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
-#include <CGAL/Polygon_mesh_processing/connected_components.h>
 #include <CGAL/AABB_face_graph_triangle_primitive.h>
 #include <CGAL/AABB_traits.h>
 #include <CGAL/AABB_tree.h>
+#include <CGAL/boost/iterator/counting_iterator.hpp>
+#include <CGAL/box_intersection_d.h>
+#include <CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h>
+#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+#include <CGAL/Polygon_mesh_processing/connected_components.h>
 #include <CGAL/Side_of_triangle_mesh.h>
 
+#include <boost/function_output_iterator.hpp>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/utility/enable_if.hpp>
+#include <boost/mpl/if.hpp>
+
+#include <exception>
+#include <iterator>
+#include <utility>
+#include <vector>
+
 namespace CGAL {
 namespace Polygon_mesh_processing{
 namespace internal {
@@ -401,7 +408,9 @@ compute_face_face_intersection(const FaceRange& face_range1,
 
   typedef TriangleMesh TM;
   typedef typename boost::graph_traits<TM>::face_descriptor face_descriptor;
-  typedef typename CGAL::Box_intersection_d::Box_with_info_d<double, 3, face_descriptor> Box;
+
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, face_descriptor, Box_policy> Box;
 
   CGAL::Bbox_3 b1 = CGAL::Polygon_mesh_processing::bbox(tm1, np1),
                b2 = CGAL::Polygon_mesh_processing::bbox(tm2, np2);
@@ -414,12 +423,8 @@ compute_face_face_intersection(const FaceRange& face_range1,
   // make one box per facet
   std::vector<Box> boxes1;
   std::vector<Box> boxes2;
-  boxes1.reserve(
+  boxes1.reserve(std::distance(boost::begin(face_range1), boost::end(face_range1)));
-        std::distance( boost::begin(face_range1), boost::end(face_range1) )
+  boxes2.reserve(std::distance(boost::begin(face_range2), boost::end(face_range2)));
-        );
-  boxes2.reserve(
-        std::distance( boost::begin(face_range2), boost::end(face_range2) )
-        );
 
   typedef typename GetVertexPointMap<TM, NamedParameters1>::const_type VertexPointMap1;
   typedef typename GetVertexPointMap<TM, NamedParameters2>::const_type VertexPointMap2;
@@ -433,6 +438,7 @@ compute_face_face_intersection(const FaceRange& face_range1,
        typename boost::property_traits<VertexPointMap1>::value_type,
        typename boost::property_traits<VertexPointMap2>::value_type
        >::value) );
+
   for(face_descriptor f : face_range1)
   {
     boxes1.push_back(Box(Polygon_mesh_processing::face_bbox(f, tm1), f));
@@ -449,7 +455,6 @@ compute_face_face_intersection(const FaceRange& face_range1,
   std::vector<const Box*> box2_ptr(boost::make_counting_iterator<const Box*>(&boxes2[0]),
                                    boost::make_counting_iterator<const Box*>(&boxes2[0]+boxes2.size()));
 
-
   // compute intersections filtered out by boxes
   typedef typename GetGeomTraits<TM, NamedParameters1>::type GeomTraits;
   GeomTraits gt = choose_parameter(get_parameter(np1, internal_np::geom_traits), GeomTraits());
@@ -459,10 +464,7 @@ compute_face_face_intersection(const FaceRange& face_range1,
                             Box,
                             OutputIterator,
                             VertexPointMap1,
-                            VertexPointMap2> Intersect_faces(tm1, tm2,
+                            VertexPointMap2> Intersect_faces(tm1, tm2, out, vpmap1, vpmap2, gt);
-                                                             out,
-                                                             vpmap1, vpmap2,
-                                                             gt);
 
   std::ptrdiff_t cutoff = 2000;
   CGAL::box_intersection_d(box1_ptr.begin(), box1_ptr.end(),
@@ -541,21 +543,16 @@ compute_face_polyline_intersection( const FaceRange& face_range,
         typename boost::range_value<Polyline>::type>::value));
 
   std::vector<face_descriptor> faces;
-  faces.reserve(std::distance( boost::begin(face_range), boost::end(face_range) ));
+  faces.reserve(std::distance(boost::begin(face_range), boost::end(face_range)));
 
-  typedef typename CGAL::Box_intersection_d::Box_with_info_d<double, 3, std::size_t> Box;
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, std::size_t, Box_policy> Box;
 
   // make one box per facet
   std::vector<Box> boxes1;
   std::vector<Box> boxes2;
-  boxes1.reserve(
+  boxes1.reserve(std::distance(boost::begin(face_range), boost::end(face_range)));
-        std::distance( boost::begin(face_range), boost::end(face_range) )
+  boxes2.reserve(std::distance(boost::begin(polyline), boost::end(polyline)) - 1);
-        );
-
-  boxes2.reserve(
-        std::distance( boost::begin(polyline), boost::end(polyline) ) - 1
-        );
-
 
   for(face_descriptor f : face_range)
   {
@@ -571,7 +568,6 @@ compute_face_polyline_intersection( const FaceRange& face_range,
   }
 
   // generate box pointers
-
   std::vector<const Box*> box1_ptr(boost::make_counting_iterator<const Box*>(&boxes1[0]),
                                    boost::make_counting_iterator<const Box*>(&boxes1[0]+boxes1.size()));
   std::vector<const Box*> box2_ptr(boost::make_counting_iterator<const Box*>(&boxes2[0]),
@@ -587,12 +583,7 @@ compute_face_polyline_intersection( const FaceRange& face_range,
                                     OutputIterator,
                                     Polyline,
                                     VertexPointMap>
-                                    Intersect_face_polyline(tm,
+                                    Intersect_face_polyline(tm, faces, polyline, out, vpmap, gt);
-                                                            faces,
-                                                            polyline,
-                                                            out,
-                                                            vpmap,
-                                                            gt);
 
   std::ptrdiff_t cutoff = 2000;
   CGAL::box_intersection_d(box1_ptr.begin(), box1_ptr.end(),
@@ -672,34 +663,32 @@ compute_face_polylines_intersection(const FaceRange& face_range,
                                           get_const_property_map(boost::vertex_point, tm));
   typedef typename boost::property_traits<VertexPointMap>::value_type Point;
   typedef typename boost::range_value<PolylineRange>::type Polyline;
-  CGAL_static_assertion(
+  CGAL_static_assertion((boost::is_same<Point, typename boost::range_value<Polyline>::type>::value));
-        (boost::is_same<Point,
-        typename boost::range_value<Polyline>::type>::value));
 
   std::vector<face_descriptor> faces;
   faces.reserve(std::distance( boost::begin(face_range), boost::end(face_range) ));
 
-  typedef typename CGAL::Box_intersection_d::Box_with_info_d<double, 3, std::pair<std::size_t, std::size_t> > Box;
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, std::pair<std::size_t, std::size_t>, Box_policy> Box;
 
   // make one box per facet
   std::vector<Box> boxes1;
   std::vector<Box> boxes2;
-  boxes1.reserve(
+  boxes1.reserve(std::distance(boost::begin(face_range), boost::end(face_range)));
-        std::distance( boost::begin(face_range), boost::end(face_range) )
-        );
 
   std::size_t polylines_size = 0;
   for(Polyline poly : polyline_range)
   {
     polylines_size += std::distance( boost::begin(poly), boost::end(poly) ) -1;
   }
-  boxes2.reserve( polylines_size );
+  boxes2.reserve(polylines_size);
 
   for(face_descriptor f : face_range)
   {
     faces.push_back(f);
     boxes1.push_back(Box(Polygon_mesh_processing::face_bbox(f, tm), std::make_pair(0, faces.size()-1)));
   }
+
   std::size_t range_size = std::distance( boost::begin(polyline_range), boost::end(polyline_range) );
   for(std::size_t j = 0; j < range_size; ++j)
   {
@@ -730,12 +719,7 @@ compute_face_polylines_intersection(const FaceRange& face_range,
                                      PolylineRange,
                                      OutputIterator,
                                      VertexPointMap>
-                                     Intersect_face_polyline(tm,
+                                     Intersect_face_polyline(tm, faces, polyline_range, out, vpmap, gt);
-                                                             faces,
-                                                             polyline_range,
-                                                             out,
-                                                             vpmap,
-                                                             gt);
 
   std::ptrdiff_t cutoff = 2000;
   CGAL::box_intersection_d(box1_ptr.begin(), box1_ptr.end(),
@@ -773,18 +757,15 @@ compute_polyline_polyline_intersection(const Polyline& polyline1,
                                        OutputIterator out,
                                        const Kernel& K)
 {
-  typedef typename CGAL::Box_intersection_d::Box_with_info_d<double, 3, std::size_t> Box;
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, std::size_t, Box_policy> Box;
+
   typedef typename Kernel::Point_3 Point;
   // make one box per facet
   std::vector<Box> boxes1;
   std::vector<Box> boxes2;
-  boxes1.reserve(
+  boxes1.reserve(std::distance(boost::begin(polyline1), boost::end(polyline1)) - 1);
-        std::distance( boost::begin(polyline1), boost::end(polyline1) ) - 1
+  boxes2.reserve(std::distance(boost::begin(polyline2), boost::end(polyline2)) - 1);
-        );
-
-  boxes2.reserve(
-        std::distance( boost::begin(polyline2), boost::end(polyline2) ) - 1
-        );
 
   for(std::size_t i =0; i< polyline1.size()-1; ++i)
   {
@@ -813,10 +794,7 @@ compute_polyline_polyline_intersection(const Polyline& polyline1,
                                 Kernel,
                                 Box,
                                 OutputIterator>
-                                intersect_polylines(polyline1,
+                                intersect_polylines(polyline1, polyline2, out, K);
-                                                    polyline2,
-                                                    out,
-                                                    K);
 
   std::ptrdiff_t cutoff = 2000;
   CGAL::box_intersection_d(box1_ptr.begin(), box1_ptr.end(),
@@ -856,7 +834,9 @@ compute_polylines_polylines_intersection(const PolylineRange& polylines1,
                                          const Kernel& K)
 {
   //info.first is the index of the polyline in the range, info.second is the index of the point in the polyline
-  typedef typename CGAL::Box_intersection_d::Box_with_info_d<double, 3, std::pair<std::size_t, std::size_t> > Box;
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, std::pair<std::size_t, std::size_t>, Box_policy> Box;
+
   typedef typename Kernel::Point_3 Point;
   typedef typename boost::range_value<PolylineRange>::type Polyline;
 
@@ -871,6 +851,7 @@ compute_polylines_polylines_intersection(const PolylineRange& polylines1,
     b1 += CGAL::bbox_3(poly.begin(), poly.end());
   }
   boxes1.reserve( polylines_size );
+
   polylines_size = 0;
   for(Polyline poly : polylines2)
   {
@@ -881,6 +862,7 @@ compute_polylines_polylines_intersection(const PolylineRange& polylines1,
 
   if(!CGAL::do_overlap(b1,b2))
     return out;
+
   std::size_t range_size = std::distance( boost::begin(polylines1), boost::end(polylines1) );
   for(std::size_t j = 0; j < range_size; ++j)
   {
@@ -915,15 +897,11 @@ compute_polylines_polylines_intersection(const PolylineRange& polylines1,
 
 
   // compute intersections filtered out by boxes
-
   internal::Intersect_polyline_ranges<PolylineRange,
                                       Kernel,
                                       Box,
                                       OutputIterator>
-                                      intersect_polylines(polylines1,
+                                      intersect_polylines(polylines1, polylines2, out, K);
-                                                          polylines2,
-                                                          out,
-                                                          K);
 
   std::ptrdiff_t cutoff = 2000;
   CGAL::box_intersection_d(box1_ptr.begin(), box1_ptr.end(),
@@ -1595,7 +1573,9 @@ OutputIterator intersecting_meshes(const TriangleMeshRange& range,
 
   bool report_overlap =  choose_parameter(get_parameter(np, internal_np::overlap_test),false);
 
-  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, TriangleMeshIterator> Mesh_box;
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, TriangleMeshIterator, Box_policy> Mesh_box;
+
   std::vector<Mesh_box> boxes;
   boxes.reserve(std::distance(range.begin(), range.end()));
 
@@ -1604,8 +1584,7 @@ OutputIterator intersecting_meshes(const TriangleMeshRange& range,
     boxes.push_back( Mesh_box(Polygon_mesh_processing::bbox(*it), it) );
   }
 
-  std::vector<Mesh_box*> boxes_ptr(
+  std::vector<Mesh_box*> boxes_ptr(boost::make_counting_iterator(&boxes[0]),
-        boost::make_counting_iterator(&boxes[0]),
                                    boost::make_counting_iterator(&boxes[0]+boxes.size()));
 
   typedef typename boost::range_value<NamedParametersRange>::type NP_rng;
@@ -1613,12 +1592,10 @@ OutputIterator intersecting_meshes(const TriangleMeshRange& range,
   typedef typename GetGeomTraits<TriangleMesh, NamedParameters, NP_rng>::type GT;
   GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
 
-
   //get all the pairs of meshes intersecting (no strict inclusion test)
   std::ptrdiff_t cutoff = 2000;
   internal::Mesh_callback<TriangleMeshRange, GT, OutputIterator, NamedParametersRange> callback(range, out, report_overlap, gt, nps);
-  CGAL::box_self_intersection_d(boxes_ptr.begin(), boxes_ptr.end(),
+  CGAL::box_self_intersection_d(boxes_ptr.begin(), boxes_ptr.end(), callback, cutoff);
-                                callback, cutoff);
   return callback.m_iterator;
 }
 
@@ -1827,9 +1804,10 @@ surface_self_intersection(const TriangleMesh& tm,
     CGAL::Polygon_mesh_processing::parameters::all_default()
   );
 }
-} //end of namespace experimental
 
-} } //end of namespace CGAL::Polygon_mesh_processing
+} //end of namespace experimental
+} //end of namespace Polygon_mesh_processing
+} //end of namespace CGAL
 
 #include <CGAL/enable_warnings.h>
 

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h

@@ -21,215 +21,378 @@
 
 #include <CGAL/disable_warnings.h>
 
-#include <CGAL/box_intersection_d.h>
-#include <CGAL/intersections.h>
-#include <CGAL/Bbox_3.h>
-
-#include <CGAL/Kernel/global_functions_3.h>
-
-#include <vector>
-#include <exception>
-#include <boost/range.hpp>
-
-#include <boost/function_output_iterator.hpp>
-#include <boost/type_traits/is_const.hpp>
-#include <boost/graph/graph_traits.hpp>
-#include <CGAL/boost/graph/helpers.h>
-#include <CGAL/boost/graph/properties.h>
-
 #include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
 #include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
 
+#include <CGAL/algorithm.h>
+#include <CGAL/Bbox_3.h>
+#include <CGAL/boost/graph/helpers.h>
+#include <CGAL/boost/graph/properties.h>
+#include <CGAL/box_intersection_d.h>
+#include <CGAL/exceptions.h>
+#include <CGAL/intersections.h>
+#include <CGAL/iterator.h>
+#include <CGAL/Kernel/global_functions_3.h>
+#include <CGAL/Random.h>
+#include <CGAL/use.h>
+
+#ifdef CGAL_LINKED_WITH_TBB
+#include <tbb/parallel_for.h>
+#include <tbb/blocked_range.h>
+#include <tbb/concurrent_vector.h>
+#endif
+
+#include <boost/function_output_iterator.hpp>
+
+#include <exception>
+#include <sstream>
+#include <type_traits>
+#include <vector>
+
 #ifdef DOXYGEN_RUNNING
 #define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters
 #define CGAL_PMP_NP_CLASS NamedParameters
 #endif
 
 namespace CGAL {
+namespace Polygon_mesh_processing {
 namespace internal {
-template <class TM,//TriangleMesh
+
-          class Kernel,
+// Checks for 'real' intersections, i.e. not simply a shared vertex or edge
-          class Box,
+template <class GT, class TM, class VPM>
-          class OutputIterator,
+bool do_faces_intersect(typename boost::graph_traits<TM>::halfedge_descriptor h,
-          class VertexPointMap>
+                        typename boost::graph_traits<TM>::halfedge_descriptor g,
-struct Intersect_facets
+                        const TM& tmesh,
+                        const VPM vpmap,
+                        const typename GT::Construct_segment_3& construct_segment,
+                        const typename GT::Construct_triangle_3& construct_triangle,
+                        const typename GT::Do_intersect_3& do_intersect)
 {
-  // wrapper to check whether anything is inserted to output iterator
+  typedef typename boost::graph_traits<TM>::vertex_descriptor                 vertex_descriptor;
-  struct Output_iterator_with_bool
-  {
-    Output_iterator_with_bool(OutputIterator* out, bool* intersected)
-      : m_iterator(out), m_intersected(intersected) { }
-
-    template<class T>
-    void operator()(const T& t) {
-      *m_intersected = true;
-      *(*m_iterator)++ = t;
-    }
-
-    OutputIterator* m_iterator;
-    bool* m_intersected;
-  };
-// typedefs
-  typedef typename Kernel::Segment_3    Segment;
-  typedef typename Kernel::Triangle_3   Triangle;
   typedef typename boost::graph_traits<TM>::halfedge_descriptor               halfedge_descriptor;
-  typedef typename boost::property_map<TM, boost::vertex_point_t>::const_type Ppmap;
 
-// members
+  typedef typename GT::Segment_3                                              Segment;
-  const TM& m_tmesh;
+  typedef typename GT::Triangle_3                                             Triangle;
-  const VertexPointMap m_vpmap;
-  mutable OutputIterator  m_iterator;
-  mutable bool            m_intersected;
-  mutable boost::function_output_iterator<Output_iterator_with_bool> m_iterator_wrapper;
 
-  typename Kernel::Construct_segment_3  segment_functor;
+  CGAL_assertion(!is_border(h, tmesh));
-  typename Kernel::Construct_triangle_3 triangle_functor;
+  CGAL_assertion(!is_border(g, tmesh));
-  typename Kernel::Do_intersect_3       do_intersect_3_functor;
 
+  vertex_descriptor hv[3], gv[3];
+  hv[0] = target(h, tmesh);
+  hv[1] = target(next(h, tmesh), tmesh);
+  hv[2] = source(h, tmesh);
 
-  Intersect_facets(const TM& tmesh, OutputIterator it, VertexPointMap vpmap, const Kernel& kernel)
+  gv[0] = target(g, tmesh);
-    :
+  gv[1] = target(next(g, tmesh), tmesh);
-    m_tmesh(tmesh),
+  gv[2] = source(g, tmesh);
-    m_vpmap(vpmap),
-    m_iterator(it),
-    m_intersected(false),
-    m_iterator_wrapper(Output_iterator_with_bool(&m_iterator, &m_intersected)),
-    segment_functor(kernel.construct_segment_3_object()),
-    triangle_functor(kernel.construct_triangle_3_object()),
-    do_intersect_3_functor(kernel.do_intersect_3_object())
-  { }
 
-  void operator()(const Box* b, const Box* c) const
+  // check for shared edge
+  for(unsigned int i=0; i<3; ++i)
   {
-    halfedge_descriptor h = halfedge(b->info(), m_tmesh);
+    halfedge_descriptor opp_h = opposite(h, tmesh);
-    halfedge_descriptor opp_h;
+    if(face(opp_h, tmesh) == face(g, tmesh))
+    {
+      // there is an intersection if the four points are coplanar and the triangles overlap
+      if(CGAL::coplanar(get(vpmap, hv[i]),
+                        get(vpmap, hv[(i+1)%3]),
+                        get(vpmap, hv[(i+2)%3]),
+                        get(vpmap, target(next(opp_h, tmesh), tmesh))) &&
+         CGAL::coplanar_orientation(get(vpmap, hv[(i+2)%3]),
+                                    get(vpmap, hv[i]),
+                                    get(vpmap, hv[(i+1)%3]),
+                                    get(vpmap, target(next(opp_h, tmesh), tmesh)))
+           == CGAL::POSITIVE)
+      {
+        return true;
+      }
+      else
+      {
+        // there is a shared edge but no intersection
+        return false;
+      }
+    }
 
-    // check for shared egde
+    h = next(h, tmesh);
-    for(unsigned int i=0; i<3; ++i){
-      opp_h = opposite(h, m_tmesh);
-      if(face(opp_h, m_tmesh) == c->info()){
-        // there is an intersection if the four points are coplanar and
-        // the triangles overlap
-        if(CGAL::coplanar(get(m_vpmap, target(h, m_tmesh)),
-                          get(m_vpmap, target(next(h, m_tmesh), m_tmesh)),
-                          get(m_vpmap, source(h, m_tmesh)),
-                          get(m_vpmap, target(next(opp_h, m_tmesh), m_tmesh))) &&
-           CGAL::coplanar_orientation(get(m_vpmap, source(h, m_tmesh)),
-                                      get(m_vpmap, target(h, m_tmesh)),
-                                      get(m_vpmap, target(next(h, m_tmesh), m_tmesh)),
-                                      get(m_vpmap, target(next(opp_h, m_tmesh), m_tmesh)))
-             == CGAL::POSITIVE){
-          *m_iterator_wrapper++ = std::make_pair(b->info(), c->info());
-          return;
-        } else { // there is a shared edge but no intersection
-          return;
-        }
-      }
-      h = next(h, m_tmesh);
   }
 
   // check for shared vertex --> maybe intersection, maybe not
-    halfedge_descriptor g = halfedge(c->info(),m_tmesh);
+  int i(0), j(0);
-    halfedge_descriptor v;
+  bool shared = false;
-
+  for(; i<3 && (! shared); ++i)
-    if(target(h,m_tmesh) == target(g,m_tmesh))
+  {
-      v = g;
+    for(j=0; j<3 && (! shared); ++j)
-    if(target(h,m_tmesh) == target(next(g,m_tmesh),m_tmesh))
+    {
-      v = next(g,m_tmesh);
+      if(hv[i] == gv[j])
-    if(target(h,m_tmesh) == target(next(next(g,m_tmesh),m_tmesh),m_tmesh))
+      {
-      v = next(next(g,m_tmesh),m_tmesh);
+        shared = true;
-
+        break;
-    if(v == halfedge_descriptor()){
-      h = next(h,m_tmesh);
-      if(target(h,m_tmesh) == target(g,m_tmesh))
-        v = g;
-      if(target(h,m_tmesh) == target(next(g,m_tmesh),m_tmesh))
-        v = next(g,m_tmesh);
-      if(target(h,m_tmesh) == target(next(next(g,m_tmesh),m_tmesh),m_tmesh))
-        v = next(next(g,m_tmesh),m_tmesh);
-      if(v == halfedge_descriptor()){
-        h = next(h,m_tmesh);
-        if(target(h,m_tmesh) == target(g,m_tmesh))
-          v = g;
-        if(target(h,m_tmesh) == target(next(g,m_tmesh),m_tmesh))
-          v = next(g,m_tmesh);
-        if(target(h,m_tmesh) == target(next(next(g,m_tmesh),m_tmesh),m_tmesh))
-          v = next(next(g,m_tmesh),m_tmesh);
       }
     }
 
-    if(v != halfedge_descriptor()){
+    if(shared)
+      break;
+  }
+
+  if(shared)
+  {
     // found shared vertex:
-      CGAL_assertion(target(h,m_tmesh) == target(v,m_tmesh));
+    CGAL_assertion(hv[i] == gv[j]);
+
     // geometric check if the opposite segments intersect the triangles
-      Triangle t1 = triangle_functor( get(m_vpmap,target(h,m_tmesh)),
+    const Triangle t1 = construct_triangle(get(vpmap, hv[0]), get(vpmap, hv[1]), get(vpmap, hv[2]));
-                                      get(m_vpmap, target(next(h,m_tmesh),m_tmesh)),
+    const Triangle t2 = construct_triangle(get(vpmap, gv[0]), get(vpmap, gv[1]), get(vpmap, gv[2]));
-                                      get(m_vpmap, target(next(next(h,m_tmesh),m_tmesh),m_tmesh)));
-      Triangle t2 = triangle_functor( get(m_vpmap, target(v,m_tmesh)),
-                                      get(m_vpmap, target(next(v,m_tmesh),m_tmesh)),
-                                      get(m_vpmap, target(next(next(v,m_tmesh),m_tmesh),m_tmesh)));
 
-      Segment s1 = segment_functor( get(m_vpmap, target(next(h,m_tmesh),m_tmesh)),
+    const Segment s1 = construct_segment(get(vpmap, hv[(i+1)%3]), get(vpmap, hv[(i+2)%3]));
-                                    get(m_vpmap, target(next(next(h,m_tmesh),m_tmesh),m_tmesh)));
+    const Segment s2 = construct_segment(get(vpmap, gv[(j+1)%3]), get(vpmap, gv[(j+2)%3]));
-      Segment s2 = segment_functor( get(m_vpmap, target(next(v,m_tmesh),m_tmesh)),
-                                    get(m_vpmap, target(next(next(v,m_tmesh),m_tmesh),m_tmesh)));
 
-      if(do_intersect_3_functor(t1,s2)){
+    if(do_intersect(t1, s2))
-        *m_iterator_wrapper++ = std::make_pair(b->info(), c->info());
+      return true;
-      } else if(do_intersect_3_functor(t2,s1)){
+    else if(do_intersect(t2, s1))
-        *m_iterator_wrapper++ = std::make_pair(b->info(), c->info());
+      return true;
-      }
+
-      return;
+    return false;
   }
 
   // check for geometric intersection
-    Triangle t1 = triangle_functor( get(m_vpmap, target(h,m_tmesh)),
+  const Triangle th = construct_triangle(get(vpmap, hv[0]), get(vpmap, hv[1]), get(vpmap, hv[2]));
-                                    get(m_vpmap, target(next(h,m_tmesh),m_tmesh)),
+  const Triangle tg = construct_triangle(get(vpmap, gv[0]), get(vpmap, gv[1]), get(vpmap, gv[2]));
-                                    get(m_vpmap, target(next(next(h,m_tmesh),m_tmesh),m_tmesh)));
+  if(do_intersect(th, tg))
-    Triangle t2 = triangle_functor( get(m_vpmap, target(g,m_tmesh)),
+    return true;
-                                    get(m_vpmap, target(next(g,m_tmesh),m_tmesh)),
-                                    get(m_vpmap, target(next(next(g,m_tmesh),m_tmesh),m_tmesh)));
-    if(do_intersect_3_functor(t1, t2)){
-      *m_iterator_wrapper++ = std::make_pair(b->info(), c->info());
-    }
-  } // end operator ()
-}; // end struct Intersect_facets
 
-struct Throw_at_output {
+  return false;
-  class Throw_at_output_exception: public std::exception
+};
-  { };
 
-  template<class T>
+template <class Box, class TM, class VPM, class GT,
-  void operator()(const T& /* t */) const {
+          class OutputIterator>
-    throw Throw_at_output_exception();
+struct Strict_intersect_faces // "strict" as in "not sharing a subface"
+{
+  typedef typename boost::graph_traits<TM>::halfedge_descriptor halfedge_descriptor;
+
+  mutable OutputIterator m_iterator;
+  const TM& m_tmesh;
+  const VPM m_vpmap;
+  typename GT::Construct_segment_3 m_construct_segment;
+  typename GT::Construct_triangle_3 m_construct_triangle;
+  typename GT::Do_intersect_3 m_do_intersect;
+
+  Strict_intersect_faces(const TM& tmesh, VPM vpmap, const GT& gt, OutputIterator it)
+    :
+      m_iterator(it),
+      m_tmesh(tmesh),
+      m_vpmap(vpmap),
+      m_construct_segment(gt.construct_segment_3_object()),
+      m_construct_triangle(gt.construct_triangle_3_object()),
+      m_do_intersect(gt.do_intersect_3_object())
+  {}
+
+  void operator()(const Box* b, const Box* c) const
+  {
+    const halfedge_descriptor h = halfedge(b->info(), m_tmesh);
+    const halfedge_descriptor g = halfedge(c->info(), m_tmesh);
+
+    if(do_faces_intersect<GT>(h, g, m_tmesh, m_vpmap, m_construct_segment, m_construct_triangle, m_do_intersect))
+      *m_iterator++ = std::make_pair(b->info(), c->info());
   }
 };
 
-}// namespace internal
+template <class ConcurrencyTag,
-
+          class TriangleMesh,
-namespace Polygon_mesh_processing {
+          class FaceRange,
-
+          class FacePairOutputIterator,
-#ifndef DOXYGEN_RUNNING
+          class NamedParameters>
-template <class TriangleMesh
+FacePairOutputIterator
-        , class FaceRange
+self_intersections_impl(const FaceRange& face_range,
-        , class OutputIterator
-        , class NamedParameters
->
-OutputIterator
-self_intersections( const FaceRange& face_range,
                         const TriangleMesh& tmesh,
-                    OutputIterator out,
+                        FacePairOutputIterator out,
-                    const NamedParameters& np);
+                        const bool throw_on_SI,
+                        const NamedParameters& np)
+{
+  CGAL_precondition(CGAL::is_triangle_mesh(tmesh));
+
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+
+  typedef TriangleMesh                                                                   TM;
+  typedef typename boost::graph_traits<TM>::halfedge_descriptor                          halfedge_descriptor;
+  typedef typename boost::graph_traits<TM>::face_descriptor                              face_descriptor;
+
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS                                  Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, face_descriptor, Box_policy> Box;
+
+  typedef typename GetGeomTraits<TM, NamedParameters>::type                              GT;
+  GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
+
+  typedef typename GetVertexPointMap<TM, NamedParameters>::const_type                    VPM;
+  VPM vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
+                               get_const_property_map(boost::vertex_point, tmesh));
+
+  const unsigned int seed = choose_parameter(get_parameter(np, internal_np::random_seed), 0);
+  CGAL_USE(seed); // used in the random shuffle of the range, which is only done to balance tasks in parallel
+
+  const std::ptrdiff_t cutoff = 2000;
+
+  // make one box per face
+  std::vector<Box> boxes;
+  boxes.reserve(std::distance(std::begin(face_range), std::end(face_range)));
+
+  // This loop is very cheap, so there is hardly anything to gain from parallelizing it
+  for(face_descriptor f : face_range)
+  {
+    halfedge_descriptor h = halfedge(f, tmesh);
+    typename boost::property_traits<VPM>::reference
+      p = get(vpmap, target(h,tmesh)),
+      q = get(vpmap, target(next(h, tmesh), tmesh)),
+      r = get(vpmap, target(prev(h, tmesh), tmesh));
+
+    // tiny fixme: if f is degenerate, we might still have a real intersection between f
+    // and another face f', but right now we are not creating a box for f and thus not returning those
+    if(collinear(p, q, r))
+    {
+      if(throw_on_SI)
+        throw CGAL::internal::Throw_at_output_exception();
+      else
+        *out++= std::make_pair(f, f);
+    }
+    else
+    {
+      boxes.push_back(Box(p.bbox() + q.bbox() + r.bbox(), f));
+    }
+  }
+
+  // generate box pointers
+  std::vector<const Box*> box_ptr;
+  box_ptr.reserve(boxes.size());
+
+  for(Box& b : boxes)
+    box_ptr.push_back(&b);
+
+  // In case we are throwing, like in `does_self_intersect()`, we keep the geometric test to throw ASAP.
+  // This is obviously not optimal if there are no or few self-intersections: it would be a greater speed-up
+  // to do the same as for `self_intersections()`. However, doing like `self_intersections()` would
+  // be a major slow-down over sequential code if there are a lot of self-intersections...
+  typedef boost::function_output_iterator<CGAL::internal::Throw_at_output>               Throwing_output_iterator;
+  typedef internal::Strict_intersect_faces<Box, TM, VPM, GT, Throwing_output_iterator>   Throwing_filter;
+  Throwing_filter throwing_filter(tmesh, vpmap, gt, Throwing_output_iterator());
+
+#if !defined(CGAL_LINKED_WITH_TBB)
+  CGAL_static_assertion_msg (!(std::is_convertible<ConcurrencyTag, Parallel_tag>::value),
+                             "Parallel_tag is enabled but TBB is unavailable.");
+#else
+  if(std::is_convertible<ConcurrencyTag, Parallel_tag>::value)
+  {
+    // We are going to split the range into a number of smaller ranges. To handle
+    // smaller trees of roughly the same size, we first apply a random shuffle to the range
+    CGAL::Random rng(seed);
+    CGAL::cpp98::random_shuffle(box_ptr.begin(), box_ptr.end(), rng);
+
+    // Write in a concurrent vector all pairs that intersect
+    typedef tbb::concurrent_vector<std::pair<face_descriptor, face_descriptor> >         Face_pairs;
+    typedef std::back_insert_iterator<Face_pairs>                                        Face_pairs_back_inserter;
+    typedef internal::Strict_intersect_faces<Box, TM, VPM, GT, Face_pairs_back_inserter> Intersecting_faces_filter;
+
+    Face_pairs face_pairs;
+    Intersecting_faces_filter callback(tmesh, vpmap, gt, std::back_inserter(face_pairs));
+
+    if(throw_on_SI)
+      CGAL::box_self_intersection_d<ConcurrencyTag>(box_ptr.begin(), box_ptr.end(), throwing_filter, cutoff);
+    else
+      CGAL::box_self_intersection_d<ConcurrencyTag>(box_ptr.begin(), box_ptr.end(), callback, cutoff);
+
+    // Sequentially write into the output iterator
+    for(std::size_t i=0; i<face_pairs.size(); ++i)
+      *out ++= face_pairs[i];
+
+    return out;
+  }
 #endif
 
+  // Sequential version of the code
+  // Compute self-intersections filtered out by boxes
+  typedef internal::Strict_intersect_faces<Box, TM, VPM, GT, FacePairOutputIterator> Intersecting_faces_filter;
+  Intersecting_faces_filter intersect_faces(tmesh, vpmap, gt, out);
+
+  if(throw_on_SI)
+    CGAL::box_self_intersection_d<CGAL::Sequential_tag>(box_ptr.begin(), box_ptr.end(), throwing_filter, cutoff);
+  else
+    CGAL::box_self_intersection_d<CGAL::Sequential_tag>(box_ptr.begin(), box_ptr.end(), intersect_faces, cutoff);
+
+  return intersect_faces.m_iterator;
+}
+
+} // namespace internal
+
+/*!
+ * \ingroup PMP_intersection_grp
+ * collects intersections between a subset of faces of a triangulated surface mesh.
+ * Two faces are said to intersect if the corresponding triangles intersect
+ * and the intersection is not an edge nor a vertex incident to both faces.
+ *
+ * This function depends on the package \ref PkgBoxIntersectionD
+ *
+ * @pre `CGAL::is_triangle_mesh(tmesh)`
+ *
+ * @tparam ConcurrencyTag enables sequential versus parallel algorithm.
+ *                        Possible values are `Sequential_tag`, `Parallel_tag`, and `Parallel_if_available_tag`.
+ * @tparam FaceRange a model of `ConstRange` with value type `boost::graph_traits<TriangleMesh>::%face_descriptor`.
+ * @tparam TriangleMesh a model of `FaceListGraph`
+ * @tparam FacePairOutputIterator a model of `OutputIterator` holding objects of type
+ *   `std::pair<boost::graph_traits<TriangleMesh>::%face_descriptor, boost::graph_traits<TriangleMesh>::%face_descriptor>`
+ * @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
+ *
+ * @param face_range the range of faces to check for self-intersection.
+ * @param tmesh the triangulated surface mesh to be checked
+ * @param out output iterator to be filled with all pairs of non-adjacent faces that intersect
+ * @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
+ *
+ * \cgalNamedParamsBegin
+ *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`.
+ *   If this parameter is omitted, an internal property map for
+ *   `CGAL::vertex_point_t` must be available in `TriangleMesh`\cgalParamEnd
+ *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPSelfIntersectionTraits` \cgalParamEnd
+ * \cgalNamedParamsEnd
+ */
+template < class ConcurrencyTag = Sequential_tag,
+           class TriangleMesh,
+           class FaceRange,
+           class FacePairOutputIterator,
+           class NamedParameters>
+FacePairOutputIterator
+self_intersections(const FaceRange& face_range,
+                   const TriangleMesh& tmesh,
+                         FacePairOutputIterator out,
+                   const NamedParameters& np)
+{
+  return internal::self_intersections_impl<ConcurrencyTag>(face_range, tmesh, out, false /*don't throw*/, np);
+}
+
+/// \cond SKIP_IN_MANUAL
+template <class ConcurrencyTag = Sequential_tag,
+          class TriangleMesh,
+          class FaceRange,
+          class FacePairOutputIterator>
+FacePairOutputIterator
+self_intersections(const FaceRange& face_range,
+                   const TriangleMesh& tmesh,
+                         FacePairOutputIterator out)
+{
+  return self_intersections<ConcurrencyTag>(face_range, tmesh, out, CGAL::parameters::all_default());
+}
+/// \endcond
+
 /**
  * \ingroup PMP_intersection_grp
- * detects and records self-intersections of a triangulated surface mesh.
+ * collects intersections between all the faces of a triangulated surface mesh.
+ * Two faces are said to intersect if the corresponding triangles intersect
+ * and the intersection is not an edge nor a vertex incident to both faces.
+ *
  * This function depends on the package \ref PkgBoxIntersectionD
+ *
  * @pre `CGAL::is_triangle_mesh(tmesh)`
  *
+ * @tparam ConcurrencyTag enables sequential versus parallel algorithm.
+ *                         Possible values are `Sequential_tag`, `Parallel_tag`, and `Parallel_if_available_tag`.
  * @tparam TriangleMesh a model of `FaceListGraph`
- * @tparam OutputIterator a model of `OutputIterator` holding objects of type
+ * @tparam FacePairOutputIterator a model of `OutputIterator` holding objects of type
  *   `std::pair<boost::graph_traits<TriangleMesh>::%face_descriptor, boost::graph_traits<TriangleMesh>::%face_descriptor>`
  * @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
  *
@@ -249,181 +412,35 @@ self_intersections( const FaceRange& face_range,
  *
  * @return `out`
  */
-template <class TriangleMesh
+template <class ConcurrencyTag = Sequential_tag,
-        , class OutputIterator
+          class TriangleMesh,
-#ifdef DOXYGEN_RUNNING
+          class FacePairOutputIterator,
-        , class NamedParameters
+          class CGAL_PMP_NP_TEMPLATE_PARAMETERS>
-#else //avoid ambiguity with self_intersections(faces, tmesh, out)
+FacePairOutputIterator
-        , class P, class T, class R
+self_intersections(const TriangleMesh& tmesh,
-#endif
+                         FacePairOutputIterator out,
->
+                   const CGAL_PMP_NP_CLASS& np)
-OutputIterator
-self_intersections(const TriangleMesh& tmesh
-                 , OutputIterator out
-#ifdef DOXYGEN_RUNNING
-                 , const NamedParameters& np)
-#else
-                 , const Named_function_parameters<P,T,R>& np)
-#endif
 {
-  return self_intersections(faces(tmesh), tmesh, out, np);
+  return self_intersections<ConcurrencyTag>(faces(tmesh), tmesh, out, np);
 }
 
 /// \cond SKIP_IN_MANUAL
-template <class TriangleMesh, class OutputIterator>
+template <class ConcurrencyTag = Sequential_tag, class TriangleMesh, class FacePairOutputIterator>
-OutputIterator
+FacePairOutputIterator
-self_intersections(const TriangleMesh& tmesh, OutputIterator out)
+self_intersections(const TriangleMesh& tmesh, FacePairOutputIterator out)
 {
-  return self_intersections(tmesh, out,
+  return self_intersections<ConcurrencyTag>(faces(tmesh), tmesh, out, parameters::all_default());
-    CGAL::Polygon_mesh_processing::parameters::all_default());
 }
 /// \endcond
 
-/*!
- * \ingroup PMP_intersection_grp
- * Same as above but the self-intersections reported
- * are only limited to the faces in `face_range`.
- *
- * @pre `CGAL::is_triangle_mesh(tmesh)`
- *
- * @tparam FaceRange range of `boost::graph_traits<TriangleMesh>::%face_descriptor`,
- *  model of `Range`.
- * Its iterator type is `RandomAccessIterator`.
- * @tparam TriangleMesh a model of `FaceListGraph`
- * @tparam OutputIterator a model of `OutputIterator` holding objects of type
- *   `std::pair<boost::graph_traits<TriangleMesh>::%face_descriptor, boost::graph_traits<TriangleMesh>::%face_descriptor>`
- * @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
- *
- * @param face_range the range of faces to check for self-intersection.
- * @param tmesh the triangulated surface mesh to be checked
- * @param out output iterator to be filled with all pairs of non-adjacent faces that intersect
- * @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
- *
- * \cgalNamedParamsBegin
- *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`.
- *   If this parameter is omitted, an internal property map for
- *   `CGAL::vertex_point_t` must be available in `TriangleMesh`\cgalParamEnd
- *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPSelfIntersectionTraits` \cgalParamEnd
- * \cgalNamedParamsEnd
-
- */
-template <class TriangleMesh
-        , class FaceRange
-        , class OutputIterator
-        , class NamedParameters
->
-OutputIterator
-self_intersections( const FaceRange& face_range,
-                    const TriangleMesh& tmesh,
-                    OutputIterator out,
-                    const NamedParameters& np)
-{
-  CGAL_precondition(CGAL::is_triangle_mesh(tmesh));
-
-  typedef TriangleMesh TM;
-  typedef typename boost::graph_traits<TM>::face_descriptor face_descriptor;
-  typedef typename CGAL::Box_intersection_d::Box_with_info_d<double, 3, face_descriptor> Box;
-
-  // make one box per facet
-  std::vector<Box> boxes;
-  boxes.reserve(
-    std::distance( boost::begin(face_range), boost::end(face_range) )
-  );
-
-  typedef typename GetVertexPointMap<TM, NamedParameters>::const_type VertexPointMap;
-  VertexPointMap vpmap = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
-                                                      get_const_property_map(boost::vertex_point, tmesh));
-
-  for(face_descriptor f : face_range)
-  {
-    typename boost::property_traits<VertexPointMap>::reference
-      p = get(vpmap, target(halfedge(f,tmesh),tmesh)),
-      q = get(vpmap, target(next(halfedge(f, tmesh), tmesh), tmesh)),
-      r = get(vpmap, target(next(next(halfedge(f, tmesh), tmesh), tmesh), tmesh));
-
-    if ( collinear(p, q, r) )
-      *out++= std::make_pair(f,f);
-    else
-      boxes.push_back(Box(p.bbox() + q.bbox() + r.bbox(), f));
-  }
-  // generate box pointers
-  std::vector<const Box*> box_ptr;
-  box_ptr.reserve(num_faces(tmesh));
-
-  for(Box& b : boxes)
-    box_ptr.push_back(&b);
-
-  // compute self-intersections filtered out by boxes
-  typedef typename GetGeomTraits<TM, NamedParameters>::type GeomTraits;
-  CGAL::internal::Intersect_facets<TM,GeomTraits,Box,OutputIterator,VertexPointMap>
-    intersect_facets(tmesh, out, vpmap,
-      parameters::choose_parameter(parameters::get_parameter(np, internal_np::geom_traits), GeomTraits()));
-
-  std::ptrdiff_t cutoff = 2000;
-  CGAL::box_self_intersection_d(box_ptr.begin(), box_ptr.end(),intersect_facets,cutoff);
-  return intersect_facets.m_iterator;
-}
-
-/// \cond SKIP_IN_MANUAL
-template <class TriangleMesh
-        , class FaceRange
-        , class OutputIterator
->
-OutputIterator self_intersections(const FaceRange& face_range,
-                                  const TriangleMesh& tmesh,
-                                  OutputIterator out)
-{
-  return self_intersections(face_range, tmesh, out,
-    CGAL::Polygon_mesh_processing::parameters::all_default());
-}
-/// \endcond
-
-/**
- * \ingroup PMP_intersection_grp
- * tests if a triangulated surface mesh self-intersects.
- * This function depends on the package \ref PkgBoxIntersectionD
- * @pre `CGAL::is_triangle_mesh(tmesh)`
- *
- * @tparam TriangleMesh a model of `FaceListGraph`
- * @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
- *
- * @param tmesh the triangulated surface mesh to be tested
- * @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
- *
- * \cgalNamedParamsBegin
- *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tmesh`.
- *   If this parameter is omitted, an internal property map for
- *   `CGAL::vertex_point_t` must be available in `TriangleMesh`\cgalParamEnd
- *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPSelfIntersectionTraits` \cgalParamEnd
- * \cgalNamedParamsEnd
- *
- * @return true if `tmesh` self-intersects
- */
-template <class TriangleMesh
-        , class CGAL_PMP_NP_TEMPLATE_PARAMETERS
-          >
-bool does_self_intersect(const TriangleMesh& tmesh
-                        , const CGAL_PMP_NP_CLASS& np)
-{
-  CGAL_precondition(CGAL::is_triangle_mesh(tmesh));
-
-  try
-  {
-    typedef boost::function_output_iterator<CGAL::internal::Throw_at_output> OutputIterator;
-    self_intersections(tmesh, OutputIterator(), np);
-  }
-  catch( CGAL::internal::Throw_at_output::Throw_at_output_exception& )
-  { return true; }
-
-  return false;
-}
-
 /**
  * \ingroup PMP_intersection_grp
  * tests if a set of faces of a triangulated surface mesh self-intersects.
  * This function depends on the package \ref PkgBoxIntersectionD
  * @pre `CGAL::is_triangle_mesh(tmesh)`
  *
+ * @tparam ConcurrencyTag enables sequential versus parallel algorithm.
+ *                        Possible values are `Sequential_tag`, `Parallel_tag`, and `Parallel_if_available_tag`.
  * @tparam FaceRange a range of `face_descriptor`
  * @tparam TriangleMesh a model of `FaceListGraph`
  * @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
@@ -439,12 +456,12 @@ bool does_self_intersect(const TriangleMesh& tmesh
  *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `SelfIntersectionTraits` \cgalParamEnd
  * \cgalNamedParamsEnd
  *
- * @return true if the faces in `face_range` self-intersect
+ * @return `true` if the faces in `face_range` self-intersect
  */
-template <class FaceRange,
+template <class ConcurrencyTag = Sequential_tag,
+          class FaceRange,
           class TriangleMesh,
-          class NamedParameters
+          class NamedParameters>
-          >
 bool does_self_intersect(const FaceRange& face_range,
                          const TriangleMesh& tmesh,
                          const NamedParameters& np)
@@ -453,35 +470,65 @@ bool does_self_intersect(const FaceRange& face_range,
 
   try
   {
-    typedef boost::function_output_iterator<CGAL::internal::Throw_at_output> OutputIterator;
+    CGAL::Emptyset_iterator unused_out;
-    self_intersections(face_range, tmesh, OutputIterator(), np);
+    internal::self_intersections_impl<ConcurrencyTag>(face_range, tmesh, unused_out, true /*throw*/, np);
+  }
+  catch(CGAL::internal::Throw_at_output_exception&)
+  {
+    return true;
   }
-  catch( CGAL::internal::Throw_at_output::Throw_at_output_exception& )
-  { return true; }
 
   return false;
 }
 
-/// \cond SKIP_IN_MANUAL
+/**
-template <class TriangleMesh>
+ * \ingroup PMP_intersection_grp
-bool does_self_intersect(const TriangleMesh& tmesh)
+ * tests if a triangulated surface mesh self-intersects.
+ * This function depends on the package \ref PkgBoxIntersectionD
+ * @pre `CGAL::is_triangle_mesh(tmesh)`
+ *
+ * @tparam ConcurrencyTag enables sequential versus parallel algorithm.
+ *                        Possible values are `Sequential_tag`, `Parallel_tag`, and `Parallel_if_available_tag`.
+ * @tparam TriangleMesh a model of `FaceListGraph`
+ * @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
+ *
+ * @param tmesh the triangulated surface mesh to be tested
+ * @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
+ *
+ * \cgalNamedParamsBegin
+ *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tmesh`.
+ *   If this parameter is omitted, an internal property map for
+ *   `CGAL::vertex_point_t` must be available in `TriangleMesh`\cgalParamEnd
+ *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPSelfIntersectionTraits` \cgalParamEnd
+ * \cgalNamedParamsEnd
+ *
+ * @return `true` if `tmesh` self-intersects
+ */
+template <class ConcurrencyTag = Sequential_tag,
+          class TriangleMesh,
+          class CGAL_PMP_NP_TEMPLATE_PARAMETERS>
+bool does_self_intersect(const TriangleMesh& tmesh,
+                         const CGAL_PMP_NP_CLASS& np)
 {
-  return does_self_intersect(tmesh,
+  return does_self_intersect<ConcurrencyTag>(faces(tmesh), tmesh, np);
-    CGAL::Polygon_mesh_processing::parameters::all_default());
 }
 
-template <class FaceRange, class TriangleMesh>
+/// \cond SKIP_IN_MANUAL
+template <class ConcurrencyTag = Sequential_tag, class TriangleMesh>
+bool does_self_intersect(const TriangleMesh& tmesh)
+{
+  return does_self_intersect<ConcurrencyTag>(faces(tmesh), tmesh, CGAL::parameters::all_default());
+}
+
+template <class ConcurrencyTag = Sequential_tag, class FaceRange, class TriangleMesh>
 bool does_self_intersect(const FaceRange& face_range,
                          const TriangleMesh& tmesh)
 {
-  return does_self_intersect(face_range, tmesh,
+  return does_self_intersect<ConcurrencyTag>(face_range, tmesh, CGAL::parameters::all_default());
-    CGAL::Polygon_mesh_processing::parameters::all_default());
 }
-
 /// \endcond
 
-}// end namespace Polygon_mesh_processing
+}// namespace Polygon_mesh_processing
-
 }// namespace CGAL
 
 #include <CGAL/enable_warnings.h>

Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt

@@ -101,6 +101,7 @@ endif()
 
   if( TBB_FOUND )
     CGAL_target_use_TBB(test_pmp_distance)
+    CGAL_target_use_TBB(self_intersection_surface_mesh_test)
   else()
     message( STATUS "NOTICE: Intel TBB was not found. test_pmp_distance will use sequential code." )
   endif()

Polygon_mesh_processing/test/Polygon_mesh_processing/data_degeneracies/deg_on_border.off

@@ -0,0 +1,10 @@
+OFF
+5 3 0
+0 0 0
+1 0 0
+2 -1 0
+3 -1 0
+4 0 0
+3 0 1 4
+3 1 2 3
+3 1 3 4

Polygon_mesh_processing/test/Polygon_mesh_processing/data_degeneracies/existing_flip.off

@@ -0,0 +1,19 @@
+OFF
+7 10 0
+0.314482809221655 -0.012877345839067217 -0.24601469358849659
+-0.081905665816609866 0.11959178909305103 0.59582950965713732
+1 1 0
+0.5 0.5 1
+0.75 0.75 0.5
+0.67302179212919955 0.77012616600142381 -0.44072953863256931
+0.14300722468117555 0.12954038296069673 0.88891433254868446
+3 0 5 2
+3 0 3 6
+3 5 0 4
+3 2 4 3
+3 4 0 1
+3 4 1 6
+3 4 2 5
+3 3 0 2
+3 6 3 4
+3 6 1 0

Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp

@@ -1,4 +1,4 @@
-// #define CGAL_PMP_COMPUTE_NORMAL_DEBUG
+// #define CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Exact_predicates_exact_constructions_kernel_with_sqrt.h>
@@ -6,7 +6,10 @@
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/Polyhedron_3.h>
 
+#include <CGAL/centroid.h>
 #include <CGAL/Polygon_mesh_processing/compute_normal.h>
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
 
 #include <iostream>
 #include <fstream>
@@ -27,6 +30,7 @@ void test(const Mesh& mesh,
   typedef typename K::Vector_3                                            Vector;
 
   typedef typename boost::graph_traits<Mesh>::vertex_descriptor           vertex_descriptor;
+  typedef typename boost::graph_traits<Mesh>::halfedge_descriptor         halfedge_descriptor;
   typedef typename boost::graph_traits<Mesh>::face_descriptor             face_descriptor;
 
   typedef typename CGAL::GetVertexPointMap<Mesh>::const_type              VPMap;
@@ -63,12 +67,58 @@ void test(const Mesh& mesh,
   PMP::compute_normals(mesh, vnormals, fnormals, CGAL::parameters::vertex_point_map(vpmap)
                                                                   .geom_traits(K()));
 
-  for(vertex_descriptor v : vertices(mesh)) {
+#if 1//def CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
-    assert(get(vnormals, v) != CGAL::NULL_VECTOR);
+  std::ofstream vn_out("vertex_normals.cgal.polylines.txt");
+  std::ofstream fn_out("face_normals.cgal.polylines.txt");
+
+  const CGAL::Bbox_3 bb = PMP::bbox(mesh);
+  const auto bbox_diagonal = CGAL::sqrt(CGAL::square(bb.xmax() - bb.xmin()) +
+                                        CGAL::square(bb.ymax() - bb.ymin()) +
+                                        CGAL::square(bb.zmax() - bb.zmin()));
+
+  for(vertex_descriptor v : vertices(mesh))
+    vn_out << "2 " << get(vpmap, v) << " " << get(vpmap, v) + 0.1 * bbox_diagonal * get(vnormals, v) << "\n";
+  vn_out.close();
+
+  for(face_descriptor f : faces(mesh))
+  {
+    std::list<typename K::Point_3> vertices;
+    for(halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(f, mesh), mesh))
+      vertices.push_back(get(vpmap, target(h, mesh)));
+
+    const typename K::Point_3& c = CGAL::centroid(vertices.begin(), vertices.end());
+    fn_out << "2 " << c << " " << c + 0.1 * bbox_diagonal * get(fnormals, f) << "\n";
+  }
+  fn_out.close();
+#endif
+
+  // Check sanity of output
+  for(const face_descriptor f : faces(mesh))
+  {
+    // tests on non triangular meshes are @todo
+    if(CGAL::is_triangle(halfedge(f, mesh), mesh))
+    {
+      if(PMP::is_degenerate_triangle_face(f, mesh))
+        assert(get(fnormals, f) == CGAL::NULL_VECTOR);
+      else
+        assert(get(fnormals, f) != CGAL::NULL_VECTOR);
+    }
   }
 
-  for(face_descriptor f : faces(mesh)) {
+  for(const vertex_descriptor v : vertices(mesh))
-    assert(get(fnormals, f) != CGAL::NULL_VECTOR);
+  {
+    if(get(vnormals, v) == CGAL::NULL_VECTOR)
+    {
+      for(halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(v, mesh), mesh))
+      {
+        if(!is_border(h, mesh))
+        {
+          // There are other cases where a vertex normal can be null without the face normals being null,
+          // (only two incident faces with opposite normals, for example), but they are not tested for now.
+          assert(get(fnormals, face(h, mesh)) == CGAL::NULL_VECTOR);
+        }
+      }
+    }
   }
 }
 
@@ -161,6 +211,17 @@ int main()
   test("data/mannequin-devil.off");
   test("data/U.off");
 
+  test("data_degeneracies/deg_on_border.off");
+  test("data_degeneracies/degtri_edge.off");
+  test("data_degeneracies/degtri_three.off");
+  test("data_degeneracies/degtri_four.off");
+  test("data_degeneracies/degtri_nullface.off");
+  test("data_degeneracies/degtri_single.off");
+  test("data_degeneracies/existing_flip.off");
+  test("data_degeneracies/fused_vertices.off");
+  test("data_degeneracies/small_ccs.off");
+  test("data_degeneracies/trihole.off");
+
   std::cerr << "All done." << std::endl;
   return EXIT_SUCCESS;
 }

Polygon_mesh_processing/test/Polygon_mesh_processing/self_intersection_surface_mesh_test.cpp

@@ -1,22 +1,25 @@
-#include <cstdlib>
-#include <iostream>
-#include <fstream>
-#include <sstream>
-
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/Polygon_mesh_processing/self_intersections.h>
-
+#include <CGAL/tags.h>
 #include <CGAL/Timer.h>
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel     Epic;
+#include <cstdlib>
-typedef CGAL::Exact_predicates_exact_constructions_kernel       Epec;
+#include <iostream>
+#include <fstream>
+#include <sstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel     EPICK;
+typedef CGAL::Exact_predicates_exact_constructions_kernel       EPECK;
+
+namespace PMP = ::CGAL::Polygon_mesh_processing;
+namespace CP = ::CGAL::parameters;
 
 template <typename K>
-int
+int test_self_intersections(const char* filename,
-test_self_intersections(const char* filename, const bool expected)
+                            const bool expected)
 {
   typedef CGAL::Surface_mesh<typename K::Point_3>                Mesh;
   typedef typename boost::graph_traits<Mesh>::face_descriptor    face_descriptor;
@@ -35,18 +38,36 @@ test_self_intersections(const char* filename, const bool expected)
   timer.start();
 
   std::vector<std::pair<face_descriptor, face_descriptor> > intersected_tris;
-  CGAL::Polygon_mesh_processing::self_intersections(
+
-    m,
+  if(std::is_same<K, EPECK>::value) // EPECK isn't threadsafe
-    std::back_inserter(intersected_tris),
+  {
-    CGAL::Polygon_mesh_processing::parameters::vertex_index_map(get(CGAL::vertex_point, m)));
+    PMP::self_intersections<CGAL::Sequential_tag>(
+      m, std::back_inserter(intersected_tris), CP::vertex_index_map(get(CGAL::vertex_point, m)));
+  }
+  else
+  {
+    PMP::self_intersections<CGAL::Parallel_if_available_tag>(
+      m, std::back_inserter(intersected_tris), CP::vertex_index_map(get(CGAL::vertex_point, m)));
+  }
   bool intersecting_1 = !intersected_tris.empty();
 
   std::cout << "self_intersections test took " << timer.time() << " sec." << std::endl;
   std::cout << intersected_tris.size() << " pairs of triangles are intersecting." << std::endl;
 
   timer.reset();
-  bool intersecting_2 = CGAL::Polygon_mesh_processing::does_self_intersect(m,
+
-    CGAL::Polygon_mesh_processing::parameters::vertex_index_map(get(CGAL::vertex_point, m)));
+  bool intersecting_2;
+
+  if(std::is_same<K, EPECK>::value) // EPECK isn't threadsafe
+  {
+    intersecting_2 = PMP::does_self_intersect<CGAL::Sequential_tag>(
+                            m, CP::vertex_index_map(get(CGAL::vertex_point, m)));
+  }
+  else
+  {
+    intersecting_2 = PMP::does_self_intersect<CGAL::Parallel_if_available_tag>(
+                            m, CP::vertex_index_map(get(CGAL::vertex_point, m)));
+  }
 
   std::cout << "does_self_intersect test took " << timer.time() << " sec." << std::endl;
   std::cout << (intersecting_2 ? "There is a self-intersection." :
@@ -77,11 +98,11 @@ int main(int argc, char** argv)
     assert(!ss.fail()); // make sure that argv[2] is either 'true' or 'false'
   }
 
-  std::cout << "First test (Epic):" << std::endl;
+  std::cout << "First test (EPICK):" << std::endl;
-  int r = test_self_intersections<Epic>(filename, expected);
+  int r = test_self_intersections<EPICK>(filename, expected);
 
-  std::cout << "First test (Epec):" << std::endl;
+  std::cout << "First test (EPECK):" << std::endl;
-  r += test_self_intersections<Epec>(filename, expected);
+  r += test_self_intersections<EPECK>(filename, expected);
 
   // Second test ---------------------------------------------------------------
   expected = true;
@@ -93,11 +114,11 @@ int main(int argc, char** argv)
     assert(!ss.fail());
   }
 
-  std::cout << "Second test (Epic):" << std::endl;
+  std::cout << "Second test (EPICK):" << std::endl;
-  r += test_self_intersections<Epic>(filename, expected);
+  r += test_self_intersections<EPICK>(filename, expected);
 
-  std::cout << "Second test (Epec):" << std::endl;
+  std::cout << "Second test (EPECK):" << std::endl;
-  r += test_self_intersections<Epec>(filename, expected);
+  r += test_self_intersections<EPECK>(filename, expected);
 
   // Third test ----------------------------------------------------------------
   expected = true;
@@ -109,11 +130,11 @@ int main(int argc, char** argv)
     assert(!ss.fail());
   }
 
-  std::cout << "Third test (Epic):" << std::endl;
+  std::cout << "Third test (EPICK):" << std::endl;
-  r += test_self_intersections<Epic>(filename, expected);
+  r += test_self_intersections<EPICK>(filename, expected);
 
-  std::cout << "Third test (Epec):" << std::endl;
+  std::cout << "Third test (EPECK):" << std::endl;
-  r += test_self_intersections<Epec>(filename, expected);
+  r += test_self_intersections<EPECK>(filename, expected);
 
   // Fourth test ----------------------------------------------------------------
   expected = true;
@@ -125,11 +146,11 @@ int main(int argc, char** argv)
     assert(!ss.fail());
   }
 
-  std::cout << "Fourth test (Epic):" << std::endl;
+  std::cout << "Fourth test (EPICK):" << std::endl;
-  r += test_self_intersections<Epic>(filename, expected);
+  r += test_self_intersections<EPICK>(filename, expected);
 
-  std::cout << "Fourth test (Epec):" << std::endl;
+  std::cout << "Fourth test (EPECK):" << std::endl;
-  r += test_self_intersections<Epec>(filename, expected);
+  r += test_self_intersections<EPECK>(filename, expected);
 
   return r;
 }

Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_clip.cpp

@@ -51,6 +51,17 @@ void test()
   CGAL::clear(tm1);
   CGAL::clear(tm2);
 
+
+  // test with a iso-cuboid
+  input.open("data-coref/elephant.off");
+  input >> tm1;
+  input.close();
+  K::Iso_cuboid_3 iso_cuboid(K::Point_3(0,0,0), K::Point_3(0.4, 0.6, 0.4));
+
+  PMP::clip(tm1, iso_cuboid, params::clip_volume(true));
+  assert(CGAL::is_closed(tm1));
+  CGAL::clear(tm1);
+
   // test with a plane
   input.open("data-coref/cube.off");
   input >> tm1;

Polygonal_surface_reconstruction/include/CGAL/internal/compute_confidences.h

@@ -24,17 +24,10 @@
 
 
 // Concurrency
-#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
-typedef CGAL::Parallel_tag Concurrency_tag;
-#else
-typedef CGAL::Sequential_tag Concurrency_tag;
-#endif
-
 
 namespace CGAL {
-
+namespace internal {
-
-	namespace internal {
 
 		/**
 		*	Computes the confidences of the candidate faces.

Polyhedron/demo/Polyhedron/CMakeLists.txt

@@ -252,6 +252,9 @@ if(CGAL_Qt5_FOUND AND Qt5_FOUND)
   add_item(scene_image_item Scene_image_item.cpp)
   add_item(scene_surface_mesh_item Scene_surface_mesh_item.cpp)
   
+  if(TBB_FOUND)
+    CGAL_target_use_TBB(scene_surface_mesh_item)
+  endif()
   # special
 
   add_item(scene_item_decorator Scene_polyhedron_item_decorator.cpp )
@@ -262,6 +265,9 @@ if(CGAL_Qt5_FOUND AND Qt5_FOUND)
 
   add_item(scene_selection_item Scene_polyhedron_selection_item.cpp)
   target_link_libraries(scene_selection_item PUBLIC scene_item_decorator  scene_k_ring_selection)
+ if(TBB_FOUND)
+    CGAL_target_use_TBB(scene_selection_item)
+  endif()
 
   add_item(scene_shortest_path_item Plugins/Surface_mesh/Scene_polyhedron_shortest_path_item.cpp)
   target_link_libraries(scene_shortest_path_item PUBLIC scene_item_decorator scene_surface_mesh_item scene_polylines_item)