diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
index aef2bd7c852..88701b9e196 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
@@ -43,13 +43,7 @@ find_package(Eigen3 3.2.0) #(requires 3.2.0 or greater)
 # Creating entries for all .cpp/.C files with "main" routine
 # ##########################################################
 
-find_package( TBB )
 create_single_source_cgal_program( "hausdorff_distance_remeshing_example.cpp")
-if( TBB_FOUND )
-  CGAL_target_use_TBB(hausdorff_distance_remeshing_example)
-else()
-  message( STATUS "NOTICE: Intel TBB was not found. hausdorff_distance_example will use sequential code." )
-endif()
 
 if (EIGEN3_FOUND)
   # Executables that require Eigen 3.2
@@ -112,6 +106,13 @@ create_single_source_cgal_program( "triangulate_faces_example_OM.cpp")
 target_link_libraries( triangulate_faces_example_OM PRIVATE ${OPENMESH_LIBRARIES} )
 endif(OpenMesh_FOUND)
 
+find_package( TBB )
+if( TBB_FOUND )
+  CGAL_target_use_TBB(hausdorff_distance_remeshing_example)
+else()
+  message( STATUS "NOTICE: Intel TBB was not found. Sequential code will be used." )
+endif()
+
 find_package(Ceres QUIET)
 if(TARGET ceres)
   target_compile_definitions( mesh_smoothing_example PRIVATE CGAL_PMP_USE_CERES_SOLVER )
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h
index 931b42d1d8f..0c80c8c1872 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h
@@ -1,4 +1,4 @@
-// Copyright (c) 2018 GeometryFactory (France).
+// Copyright (c) 2018, 2019 GeometryFactory (France).
 // All rights reserved.
 //
 // This file is part of CGAL (www.cgal.org).
@@ -28,8 +28,6 @@
 #include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
 #include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
 #include <CGAL/Polygon_mesh_processing/border.h>
-#include <CGAL/Polygon_mesh_processing/connected_components.h>
-#include <CGAL/Polygon_mesh_processing/orientation.h>
 
 #include <CGAL/AABB_tree.h>
 #include <CGAL/AABB_traits.h>
@@ -48,6 +46,11 @@
 #include <boost/multi_index/ordered_index.hpp>
 #include <boost/multi_index/member.hpp>
 
+#ifdef CGAL_LINKED_WITH_TBB
+# include <tbb/concurrent_vector.h>
+# include <tbb/parallel_for.h>
+#endif
+
 #include <functional>
 #include <iostream>
 #include <iterator>
@@ -410,7 +413,7 @@ std::size_t snap_vertex_range_onto_vertex_range(const SourceHalfedgeRange& sourc
     CGAL_assertion(sp.sq_dist >= prev);
     CGAL_assertion_code(prev = sp.sq_dist;)
 
-#ifdef CGAL_PMP_SNAP_DEBUG
+#ifdef CGAL_PMP_SNAP_DEBUG_PP
     std::cout << "Snapping " /*<< vs*/ << " (" << get(svpm, vs) << ") "
               << " to " /*<< vt*/ << " (" << get(tvpm, vt) << ") at dist: " << sp.sq_dist << std::endl;
 #endif
@@ -585,9 +588,11 @@ class Projection_traits
   typedef CGAL::AABB_node<AABBTraits>                           Node;
 
 public:
-  explicit Projection_traits(const AABBTraits& traits)
+  explicit Projection_traits(const AABBTraits& traits,
+                             const FT sq_tolerance)
     : m_traits(traits),
-      m_closest_point_initialized(false)
+      m_closest_point_initialized(false),
+      m_sq_dist(sq_tolerance)
   {}
 
   bool go_further() const { return true; }
@@ -595,8 +600,8 @@ public:
   void intersection(const Point& query, const Primitive& primitive)
   {
     // skip the primitive if one of its endpoints is the query
-    typename Primitive::Datum s = primitive.datum(m_traits.shared_data());
-    if(s[0] == query || s[1] == query)
+    const typename Primitive::Datum& s = primitive.datum(m_traits.shared_data());
+    if(m_traits.equal_3_object()(s[0], query) || m_traits.equal_3_object()(s[1], query))
       return;
 
     if(!m_closest_point_initialized)
@@ -604,31 +609,34 @@ public:
       m_closest_point_initialized = true;
       m_closest_primitive = primitive.id();
       m_closest_point = primitive.reference_point(m_traits.shared_data());
+      m_sq_dist = m_traits.squared_distance_object()(query, m_closest_point);
     }
 
     Point new_closest_point = m_traits.closest_point_object()(query, primitive, m_closest_point);
-    if(new_closest_point != m_closest_point)
+    if(!m_traits.equal_3_object()(new_closest_point, m_closest_point))
     {
       m_closest_primitive = primitive.id();
-      m_closest_point = new_closest_point; // this effectively shrinks the sphere
+      m_closest_point = new_closest_point;
+      m_sq_dist = m_traits.squared_distance_object()(query, m_closest_point);
     }
   }
 
   bool do_intersect(const Point& query, const Node& node) const
   {
-    return !m_closest_point_initialized ||
-           m_traits.compare_distance_object()(query, node.bbox(), m_closest_point) == CGAL::SMALLER;
+    return m_traits.compare_distance_object()(query, node.bbox(), m_sq_dist) == CGAL::SMALLER;
   }
 
-  Point closest_point() const { return m_closest_point; }
+  const Point& closest_point() const { return m_closest_point; }
   typename Primitive::Id closest_primitive_id() const { return m_closest_primitive; }
   bool closest_point_initialized() const { return m_closest_point_initialized; }
 
 private:
-  Point m_closest_point;
-  typename Primitive::Id m_closest_primitive;
   const AABBTraits& m_traits;
   bool m_closest_point_initialized;
+  Point m_closest_point;
+  FT m_sq_dist;
+
+  typename Primitive::Id m_closest_primitive;
 };
 
 template <typename PolygonMesh, typename GeomTraits, typename VPM>
@@ -650,6 +658,116 @@ struct Compare_points_along_edge
   Point m_src;
 };
 
+template <typename TriangleMesh, typename EdgeToSplitContainer, typename AABBTree, typename VPM, typename GT>
+void find_splittable_edge(const std::pair<typename boost::property_traits<VPM>::value_type,
+                                          typename GT::FT>& point_with_tolerance,
+                          EdgeToSplitContainer& edges_to_split,
+                          const AABBTree* aabb_tree_ptr,
+                          const TriangleMesh& /* pms */,
+                          const VPM& /* vpms */,
+                          const TriangleMesh& pmt,
+                          const VPM& vpmt,
+                          const GT& gt)
+{
+  typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor   halfedge_descriptor;
+  typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor       edge_descriptor;
+
+  typedef typename GT::FT                                                   FT;
+  typedef typename boost::property_traits<VPM>::value_type                  Point;
+
+  typedef typename AABBTree::AABB_traits                                    AABB_traits;
+
+  const Point& query = point_with_tolerance.first;
+  const FT sq_tolerance = CGAL::square(point_with_tolerance.second);
+
+  // use the current halfedge as hint
+  internal::Projection_traits<AABB_traits> traversal_traits(aabb_tree_ptr->traits(), sq_tolerance);
+  aabb_tree_ptr->traversal(query, traversal_traits);
+
+  if(!traversal_traits.closest_point_initialized())
+    return;
+
+  const Point& closest_p = traversal_traits.closest_point();
+
+  // The filtering in the AABB tree checks the dist query <-> node bbox, which might be smaller than
+  // the actual distance between the query <-> closest point
+  const FT sq_dist_to_closest = gt.compute_squared_distance_3_object()(query, closest_p);
+  bool is_close_enough = (sq_dist_to_closest <= sq_tolerance);
+
+#ifdef CGAL_PMP_SNAP_DEBUG_PP
+  std::cout << "  Closest point: (" << closest_p << ")" << std::endl
+            << "  at distance " << gt.compute_squared_distance_3_object()(query, closest_p)
+            << " with a tolerance of " << point_with_tolerance.second << " at vertex (squared tol: " << sq_tolerance
+            << " && close enough? " << is_close_enough << ")" << std::endl;
+#endif
+
+  if(!is_close_enough)
+    return;
+
+#ifdef CGAL_PMP_SNAP_DEBUG_PP
+  std::cout << "\t and is thus beneath tolerance" << std::endl;
+#endif
+
+  edge_descriptor closest = traversal_traits.closest_primitive_id();
+  CGAL_assertion(get(vpmt, source(closest, pmt)) != query &&
+                                                    get(vpmt, target(closest, pmt)) != query);
+
+  halfedge_descriptor clos_hd = halfedge(closest, pmt);
+  CGAL_assertion(is_border(edge(clos_hd, pmt), pmt));
+
+  if(!is_border(clos_hd, pmt))
+    clos_hd = opposite(clos_hd, pmt);
+
+#ifdef CGAL_PMP_SNAP_DEBUG_PP
+  std::cout << "splitting position: " << query << std::endl;
+#endif
+
+  edges_to_split.push_back(std::make_pair(query, clos_hd));
+}
+
+#ifdef CGAL_LINKED_WITH_TBB
+template <typename PointWithToleranceContainer,
+          typename TriangleMesh, typename EdgeToSplitContainer,
+          typename AABBTree, typename VPM, typename GT>
+struct Find_splittable_edge_for_parallel_for
+{
+  Find_splittable_edge_for_parallel_for(const PointWithToleranceContainer& points_with_tolerance,
+                                        EdgeToSplitContainer& edges_to_split,
+                                        const AABBTree* aabb_tree_ptr,
+                                        const TriangleMesh& pms,
+                                        const VPM& vpms,
+                                        const TriangleMesh& pmt,
+                                        const VPM& vpmt,
+                                        const GT& gt)
+    :
+      m_points_with_tolerance(points_with_tolerance),
+      m_edges_to_split(edges_to_split), m_aabb_tree_ptr(aabb_tree_ptr),
+      m_pms(pms), m_vpms(vpms), m_pmt(pmt), m_vpmt(vpmt), m_gt(gt)
+  { }
+
+  void operator()(const tbb::blocked_range<size_t>& r) const
+  {
+    for(std::size_t i=r.begin(); i!=r.end(); ++i)
+    {
+      find_splittable_edge(m_points_with_tolerance[i], m_edges_to_split,
+                           m_aabb_tree_ptr,
+                           m_pms, m_vpms, m_pmt, m_vpmt,
+                           m_gt);
+    }
+  }
+
+private:
+  const PointWithToleranceContainer& m_points_with_tolerance;
+  EdgeToSplitContainer& m_edges_to_split;
+  const AABBTree* m_aabb_tree_ptr;
+  const TriangleMesh& m_pms;
+  const VPM m_vpms;
+  const TriangleMesh& m_pmt;
+  const VPM m_vpmt;
+  const GT& m_gt;
+};
+#endif
+
 } // namespace internal
 
 namespace experimental {
@@ -668,13 +786,13 @@ std::size_t snap_vertex_range_onto_vertex_range_non_conforming(const HalfedgeRan
 {
   typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor     vertex_descriptor;
   typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor   halfedge_descriptor;
-  typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor       edge_descriptor;
 
   typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::type   VPM;
+  typedef typename boost::property_traits<VPM>::value_type                  Point;
+  typedef typename boost::property_traits<VPM>::reference                   Point_ref;
 
   typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type       GT;
   typedef typename GT::FT                                                   FT;
-  typedef typename GT::Point_3                                              Point;
   typedef typename GT::Vector_3                                             Vector;
 
   typedef internal::Compare_points_along_edge<TriangleMesh, GT, VPM>        Point_along_edge_comparer;
@@ -697,8 +815,6 @@ std::size_t snap_vertex_range_onto_vertex_range_non_conforming(const HalfedgeRan
 
   const GT gt = choose_parameter(get_parameter(nps, internal_np::geom_traits), GT());
 
-  const bool is_same_mesh = (&pms == &pmt); // @todo probably not optimal
-
   // start by snapping vertices together to simplify things
   snapped_n = snap_vertex_range_onto_vertex_range(source_hrange, pms, target_hrange, pmt, tol_pmap, nps, npt);
 
@@ -708,121 +824,100 @@ std::size_t snap_vertex_range_onto_vertex_range_non_conforming(const HalfedgeRan
             << std::distance(target_hrange.begin(), target_hrange.end()) << std::endl;
 #endif
 
-  typedef std::map<Point, std::set<vertex_descriptor /*target vd*/> >     Occurrence_map;
-  Occurrence_map occurrences_as_target;
-  for(halfedge_descriptor hd : target_hrange)
-  {
-    vertex_descriptor vd = target(hd, pmt);
-
-    std::set<vertex_descriptor> corresponding_vd;
-    corresponding_vd.insert(vd);
-
-    std::pair<typename Occurrence_map::iterator, bool> is_insert_successful =
-      occurrences_as_target.insert(std::make_pair(get(vpmt, vd), corresponding_vd));
-    if(!is_insert_successful.second) // point already existed in the map
-      is_insert_successful.first->second.insert(vd);
-  }
-
-  // Since we're inserting primitives one by one, we can't pass this shared data in the constructor of the tree
-  AABB_Traits aabb_traits;
-  aabb_traits.set_shared_data(pmt, vpmt);
-
-  // Fill the AABB-tree
-  AABB_tree aabb_tree(aabb_traits);
-  for(halfedge_descriptor hd : target_hrange)
-  {
-    CGAL_precondition(is_border(edge(hd, pmt), pmt));
-    aabb_tree.insert(Primitive(edge(hd, pmt), pmt, vpmt));
-  }
-
   // Collect border points that can be projected onto a border edge
-  std::vector<std::pair<Point, halfedge_descriptor> > edges_to_split;
-  std::unordered_set<vertex_descriptor> unique_vertices;
+  typedef std::map<Point, FT>                       Unique_vertices_with_tolerance;  // unordered map is empirically slower
+
+  // If the point appears multiple times on the source border, use it only once to snap target edges
+  // use the smallest tolerance from all source vertices
+
+#ifdef CGAL_PMP_SNAP_DEBUG
+  std::cout << "gather unique points in source range" << std::endl;
+#endif
+
+  Unique_vertices_with_tolerance unique_source_vertices;
   for(halfedge_descriptor hd : source_hrange)
   {
-    const vertex_descriptor vd = target(hd, pms);
-    if(!unique_vertices.insert(vd).second)
-    {
-      // if 'vd' appears multiple times on the source border, use it only once to snap target edges
-      continue;
-    }
-
-    const Point& query = get(vpms, vd);
-    const std::set<vertex_descriptor>& occurrences = occurrences_as_target[query];
-
-    // Skip points that are already attached to another border. Keeping it in two 'continue' for clarity.
-
-    // If we are working with a single mesh, the vertex is only blocked if another vertex has the same
-    // position (that is, if occurrences.size() > 1)
-    if(is_same_mesh && occurrences.size() > 1)
-      continue;
-
-    // If it's not the same mesh, then block as soon as a vertex in the target range has already that position
-    if(!is_same_mesh && !occurrences.empty())
-      continue;
-
     // Skip the source vertex if its two incident halfedges are geometrically identical (it means that
     // the two halfedges are already stitchable and we don't want this common vertex to be used
     // to split a halfedge somewhere else)
     if(get(vpms, source(hd, pms)) == get(vpms, target(next(hd, pms), pms)))
       continue;
 
-#ifdef CGAL_PMP_SNAP_DEBUG
-    std::cout << "Query: " << vd << " (" << query << ")" << std::endl;
+    const vertex_descriptor vd = target(hd, pms);
+    const Point_ref query = get(vpms, vd);
+    const FT tolerance = get(tol_pmap, vd);
+
+    std::pair<typename Unique_vertices_with_tolerance::iterator, bool> is_insert_successful =
+      unique_source_vertices.insert(std::make_pair(query, tolerance));
+
+    // Keep the lowest tolerance out of all the source vertices with the same position
+    is_insert_successful.first->second = (std::min)(is_insert_successful.first->second, tolerance);
+
+#ifdef CGAL_PMP_SNAP_DEBUG_PP
+    if(is_insert_successful.second)
+      std::cout << "Query: " << vd << " (" << query << ")" << std::endl;
 #endif
-
-    // use the current halfedge as hint
-    internal::Projection_traits<AABB_Traits> traversal_traits(aabb_tree.traits());
-    aabb_tree.traversal(query, traversal_traits);
-
-    if(!traversal_traits.closest_point_initialized())
-      continue;
-
-    const FT sq_tolerance = CGAL::square(get(tol_pmap, vd));
-    const Point& closest_p = traversal_traits.closest_point();
-    const FT sq_dist_to_closest = gt.compute_squared_distance_3_object()(query, closest_p);
-    bool is_close_enough = (sq_dist_to_closest <= sq_tolerance);
-
-#ifdef CGAL_PMP_SNAP_DEBUG
-    std::cout << "  Closest point: (" << closest_p << ")" << std::endl
-              << "  at distance " << gt.compute_squared_distance_3_object()(query, closest_p)
-              << " with a tolerance of " << get(tol_pmap, vd) << " at vertex (squared: " << sq_tolerance
-              << " && close enough? " << is_close_enough << ")" << std::endl;;
-#endif
-
-    if(is_close_enough)
-    {
-#ifdef CGAL_PMP_SNAP_DEBUG
-      std::cout << "\t and is thus beneath tolerance" << std::endl;
-#endif
-
-      edge_descriptor closest = traversal_traits.closest_primitive_id();
-      CGAL_assertion(get(vpmt, source(closest, pmt)) != query &&
-                     get(vpmt, target(closest, pmt)) != query);
-
-      halfedge_descriptor clos_hd = halfedge(closest, pmt);
-      CGAL_assertion(is_border(edge(clos_hd, pmt), pmt));
-
-      if(!is_border(clos_hd, pmt))
-        clos_hd = opposite(clos_hd, pmt);
-
-      Point new_pos = query;
-      put(vpms, vd, new_pos);
-
-#ifdef CGAL_PMP_SNAP_DEBUG
-      std::cout << "splitting position: " << new_pos << std::endl;
-#endif
-
-      edges_to_split.push_back(std::make_pair(new_pos, clos_hd));
-    }
   }
 
+  // Since we're inserting primitives one by one, we can't pass this shared data in the constructor of the tree
+  AABB_Traits aabb_traits;
+  aabb_traits.set_shared_data(pmt, vpmt);
+  AABB_tree aabb_tree(aabb_traits);
+
+  for(halfedge_descriptor hd : target_hrange)
+  {
+    CGAL_precondition(is_border(edge(hd, pmt), pmt));
+    aabb_tree.insert(Primitive(edge(hd, pmt), pmt, vpmt));
+  }
+
+  // Now, check which edges of the target range ought to be split by source vertices
+#ifdef CGAL_PMP_SNAP_DEBUG_PP
+  std::cout << "Collect edges to split w/ " << unique_source_vertices.size() << " unique vertices" << std::endl;
+#endif
+
+#ifdef CGAL_LINKED_WITH_TBB
+#ifdef CGAL_PMP_SNAP_DEBUG
+  std::cout << "Parallel find splittable edges!" << std::endl;
+#endif
+
+  typedef tbb::concurrent_vector<std::pair<Point, halfedge_descriptor> >        Concurrent_edge_to_split_container;
+
+  // Can't parallel for a map...
+  std::vector<std::pair<Point, FT> > unique_vertices_with_tolerance_vector(unique_source_vertices.begin(),
+                                                                           unique_source_vertices.end());
+
+  typedef internal::Find_splittable_edge_for_parallel_for<std::vector<std::pair<Point, FT> >,
+                                                          TriangleMesh,
+                                                          Concurrent_edge_to_split_container,
+                                                          AABB_tree, VPM, GT>    Functor;
+
+  Concurrent_edge_to_split_container edges_to_split;
+  Functor f(unique_vertices_with_tolerance_vector, edges_to_split, &aabb_tree, pms, vpms, pmt, vpmt, gt);
+  tbb::parallel_for(tbb::blocked_range<std::size_t>(0, unique_vertices_with_tolerance_vector.size()), f);
+#else // CGAL_LINKED_WITH_TBB
+#ifdef CGAL_PMP_SNAP_DEBUG
+  std::cout << "Sequential find splittable edges!" << std::endl;
+#endif
+  typedef std::pair<Point, FT>                      Point_with_tolerance;
+
+  std::vector<std::pair<Point, halfedge_descriptor> > edges_to_split;
+  for(const Point_with_tolerance& point_with_tolerance : unique_source_vertices)
+  {
+    internal::find_splittable_edge(point_with_tolerance, edges_to_split, &aabb_tree,
+                                   pms, vpms, pmt, vpmt, gt);
+  }
+#endif
+
   // Sort points falling on the same edge and split the edge
   std::map<halfedge_descriptor, std::vector<Point> > points_per_edge;
   typedef std::pair<Point, halfedge_descriptor> Pair_type;
   for(const Pair_type& p : edges_to_split)
     points_per_edge[p.second].push_back(p.first);
 
+#ifdef CGAL_PMP_SNAP_DEBUG
+  std::cout << "split " << points_per_edge.size() << " edges" << std::endl;
+#endif
+
   typedef std::pair<const halfedge_descriptor, std::vector<Point> > Map_type;
   for(Map_type& mt : points_per_edge)
   {
@@ -832,7 +927,7 @@ std::size_t snap_vertex_range_onto_vertex_range_non_conforming(const HalfedgeRan
 
     if(mt.second.size() > 1)
     {
-#ifdef CGAL_PMP_SNAP_DEBUG
+#ifdef CGAL_PMP_SNAP_DEBUG_PP
       std::cout << " MUST SORT ON BORDER!" << std::endl;
 #endif
       /// \todo Sorting the projected positions is too simple (for example, this doesn't work
@@ -847,7 +942,7 @@ std::size_t snap_vertex_range_onto_vertex_range_non_conforming(const HalfedgeRan
     halfedge_descriptor hd_to_split = mt_hd;
     for(const Point& p : mt.second)
     {
-#ifdef CGAL_PMP_SNAP_DEBUG
+#ifdef CGAL_PMP_SNAP_DEBUG_PP
       std::cout << "SPLIT " << edge(hd_to_split, pmt) << " |||| "
                 << " vs " << source(hd_to_split, pmt) << " (" << pmt.point(source(hd_to_split, pmt)) << ")"
                 << " --- vt " << target(hd_to_split, pmt) << " (" << pmt.point(target(hd_to_split, pmt)) << ")" << std::endl;
@@ -878,9 +973,9 @@ std::size_t snap_vertex_range_onto_vertex_range_non_conforming(const HalfedgeRan
        *       opp
        */
 
-      const Point left_pt = get(vpmt, source(res, pmt));
-      const Point right_pt = get(vpmt, target(hd_to_split, pmt));
-      const Point opp = get(vpmt, target(next(opposite(res, pmt), pmt), pmt));
+      const Point_ref left_pt = get(vpmt, source(res, pmt));
+      const Point_ref right_pt = get(vpmt, target(hd_to_split, pmt));
+      const Point_ref opp = get(vpmt, target(next(opposite(res, pmt), pmt), pmt));
 
       // Check if 'p' is "visible" from 'opp' (i.e. its projection on the plane 'Pl(left, opp, right)'
       // falls in the cone with apex 'opp' and sides given by 'left' and 'right')
@@ -894,7 +989,7 @@ std::size_t snap_vertex_range_onto_vertex_range_non_conforming(const HalfedgeRan
 
       const bool is_visible = (!left_of_left && !right_of_right);
 
-#ifdef CGAL_PMP_SNAP_DEBUG
+#ifdef CGAL_PMP_SNAP_DEBUG_PP
       std::cout << "Left/Right: " << left_of_left << " " << right_of_right << std::endl;
       std::cout << "visible from " << opp << " ? " << is_visible << std::endl;
 #endif
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h
index 0b3de74bd38..d9f9b0f42c0 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h
@@ -470,7 +470,7 @@ std::size_t remove_isolated_points_in_polygon_soup(PointRange& points,
     for(std::size_t i=0, polygon_size = polygon.size(); i<polygon_size; ++i)
     {
       polygon[i] = id_remapping[polygon[i]];
-      CGAL_postcondition(polygon[i] < points.size());
+      CGAL_postcondition(static_cast<std::size_t>(polygon[i]) < points.size());
     }
   }
 
diff --git a/Stream_support/include/CGAL/IO/PLY.h b/Stream_support/include/CGAL/IO/PLY.h
index 6f2b90c93dd..59e3beb067f 100644
--- a/Stream_support/include/CGAL/IO/PLY.h
+++ b/Stream_support/include/CGAL/IO/PLY.h
@@ -23,13 +23,16 @@
 
 #include <CGAL/Kernel_traits.h>
 #include <CGAL/IO/io.h>
-
 #include <CGAL/property_map.h>
 
-#include <tuple>
+#include <boost/cstdint.hpp>
+
 #include <iostream>
 #include <sstream>
 #include <string>
+#include <tuple>
+#include <utility>
+#include <vector>
 
 #define TRY_TO_GENERATE_PROPERTY(STD_TYPE, T_TYPE, TYPE)                \
   if (type == STD_TYPE  || type == T_TYPE)                              \