diff --git a/Point_set_processing_3/include/CGAL/grid_simplify_point_set.h b/Point_set_processing_3/include/CGAL/grid_simplify_point_set.h
index fb488710b73..851f48cf6fc 100644
--- a/Point_set_processing_3/include/CGAL/grid_simplify_point_set.h
+++ b/Point_set_processing_3/include/CGAL/grid_simplify_point_set.h
@@ -30,6 +30,7 @@
 #include <CGAL/point_set_processing_assertions.h>
 #include <CGAL/unordered.h>
 #include <CGAL/Iterator_range.h>
+#include <CGAL/function.h>
 #include <boost/functional/hash.hpp>
 
 #include <CGAL/boost/graph/named_function_params.h>
@@ -204,6 +205,8 @@ grid_simplify_point_set(
   
   typedef typename Point_set_processing_3::GetPointMap<PointRange, NamedParameters>::const_type PointMap;
   PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
+  const cpp11::function<bool(double)>& callback = choose_param(get_param(np, internal_np::callback),
+                                                               cpp11::function<bool(double)>());
 
   // actual type of input points
   typedef typename std::iterator_traits<typename PointRange::iterator>::value_type Enriched_point;
@@ -214,12 +217,15 @@ grid_simplify_point_set(
   // points_to_keep[] will contain 1 point per cell; the others will be in points_to_remove[].
   Epsilon_point_set_3<Enriched_point, PointMap> points_to_keep(epsilon, point_map);
   std::deque<Enriched_point> points_to_remove;
-  for (typename PointRange::iterator it = points.begin(); it != points.end(); it++)
+  std::size_t nb = 0, nb_points = points.size();
+  for (typename PointRange::iterator it = points.begin(); it != points.end(); it++, ++ nb)
   {
     std::pair<typename Epsilon_point_set_3<Enriched_point, PointMap>::iterator,bool> result;
     result = points_to_keep.insert(*it);
     if (!result.second) // if not inserted
       points_to_remove.push_back(*it);
+    if (callback && !callback ((nb+1) / double(nb_points)))
+      break;
   }
 
   // Replaces `[first, beyond)` range by the content of points_to_keep, then points_to_remove.
diff --git a/Point_set_processing_3/include/CGAL/hierarchy_simplify_point_set.h b/Point_set_processing_3/include/CGAL/hierarchy_simplify_point_set.h
index bbb94e0d81e..18f0cf25da3 100644
--- a/Point_set_processing_3/include/CGAL/hierarchy_simplify_point_set.h
+++ b/Point_set_processing_3/include/CGAL/hierarchy_simplify_point_set.h
@@ -39,6 +39,7 @@
 #include <CGAL/PCA_util.h>
 #include <CGAL/squared_distance_3.h>
 #include <CGAL/Iterator_range.h>
+#include <CGAL/function.h>
 
 #include <CGAL/boost/graph/named_function_params.h>
 #include <CGAL/boost/graph/named_params_helper.h>
@@ -167,6 +168,8 @@ namespace CGAL {
     PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
     unsigned int size = choose_param(get_param(np, internal_np::size), 10);
     double var_max = choose_param(get_param(np, internal_np::maximum_variation), 1./3.);
+    const cpp11::function<bool(double)>& callback = choose_param(get_param(np, internal_np::callback),
+                                                                 cpp11::function<bool(double)>());
 
     typedef typename std::iterator_traits<typename PointRange::iterator>::value_type Input_type;
 
@@ -191,6 +194,8 @@ namespace CGAL {
 
     std::list<Input_type> points_to_keep;
     std::list<Input_type> points_to_remove;
+
+    std::size_t nb_done = 0;
     
     while (!(clusters_stack.empty ()))
       {
@@ -203,6 +208,7 @@ namespace CGAL {
 	    points_to_keep.splice (points_to_keep.end (), current_cluster->first,
 				   current_cluster->first.begin ());
 	    clusters_stack.pop_front ();
+            ++ nb_done;
 	    continue;
 	  }
 
@@ -268,6 +274,7 @@ namespace CGAL {
 		cluster_iterator nonempty = (current_cluster->first.empty ()
 					     ? negative_side : current_cluster);
 
+                nb_done += nonempty->first.size();
 		// Compute the centroid
 		nonempty->second = internal::hsps_centroid (nonempty->first.begin (),
 							    nonempty->first.end (),
@@ -279,6 +286,7 @@ namespace CGAL {
 						 point_map,
 						 nonempty->second,
 						 Kernel ());
+                
 		clusters_stack.pop_front ();
 		clusters_stack.pop_front ();
 	      }
@@ -311,6 +319,7 @@ namespace CGAL {
 	// and output point
 	else
 	  {
+            nb_done += current_cluster->first.size();
 	    internal::hsc_terminate_cluster (current_cluster->first,
 					     points_to_keep,
 					     points_to_remove,
@@ -319,7 +328,14 @@ namespace CGAL {
 					     Kernel ());
 	    clusters_stack.pop_front ();
 	  }
+
+        if (callback && !callback (nb_done / double (points.size())))
+          break;
       }
+
+    if (callback)
+      callback (1.);
+    
     typename PointRange::iterator first_point_to_remove =
       std::copy (points_to_keep.begin(), points_to_keep.end(), points.begin());
     std::copy (points_to_remove.begin(), points_to_remove.end(), first_point_to_remove);
diff --git a/Point_set_processing_3/include/CGAL/jet_estimate_normals.h b/Point_set_processing_3/include/CGAL/jet_estimate_normals.h
index edc1cc6ee6a..bbd97d98e3e 100644
--- a/Point_set_processing_3/include/CGAL/jet_estimate_normals.h
+++ b/Point_set_processing_3/include/CGAL/jet_estimate_normals.h
@@ -40,9 +40,11 @@
 #include <list>
 
 #ifdef CGAL_LINKED_WITH_TBB
+#include <CGAL/internal/Parallel_callback.h>
 #include <tbb/parallel_for.h>
 #include <tbb/blocked_range.h>
 #include <tbb/scalable_allocator.h>  
+#include <tbb/atomic.h>
 #endif // CGAL_LINKED_WITH_TBB
 
 namespace CGAL {
@@ -125,17 +127,28 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
     const unsigned int degree_fitting;
     const std::vector<Point>& input;
     std::vector<Vector>& output;
+    tbb::atomic<std::size_t>& advancement;
+    tbb::atomic<bool>& interrupted;
 
   public:
     Jet_estimate_normals(Tree& tree, unsigned int k, std::vector<Point>& points,
-		     unsigned int degree_fitting, std::vector<Vector>& output)
+                         unsigned int degree_fitting, std::vector<Vector>& output,
+                         tbb::atomic<std::size_t>& advancement,
+                         tbb::atomic<bool>& interrupted)
       : tree(tree), k (k), degree_fitting (degree_fitting), input (points), output (output)
+      , advancement (advancement)
+      , interrupted (interrupted)
     { }
     
     void operator()(const tbb::blocked_range<std::size_t>& r) const
     {
       for( std::size_t i = r.begin(); i != r.end(); ++i)
+      {
+        if (interrupted)
+          break;
 	output[i] = CGAL::internal::jet_estimate_normal<Kernel,SvdTraits>(input[i], tree, k, degree_fitting);
+        ++ advancement;
+      }
     }
 
   };
@@ -212,6 +225,8 @@ jet_estimate_normals(
   PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
   NormalMap normal_map = choose_param(get_param(np, internal_np::normal_map), NormalMap());
   unsigned int degree_fitting = choose_param(get_param(np, internal_np::degree_fitting), 2);
+  const cpp11::function<bool(double)>& callback = choose_param(get_param(np, internal_np::callback),
+                                                               cpp11::function<bool(double)>());
 
   typedef typename Kernel::Point_3 Point;
 
@@ -254,28 +269,35 @@ jet_estimate_normals(
 #else
    if (boost::is_convertible<ConcurrencyTag,Parallel_tag>::value)
    {
+     internal::Point_set_processing_3::Parallel_callback
+       parallel_callback (callback, kd_tree_points.size());
+     
      std::vector<Vector> normals (kd_tree_points.size ());
      CGAL::internal::Jet_estimate_normals<Kernel, SvdTraits, Tree>
-       f (tree, k, kd_tree_points, degree_fitting, normals);
+       f (tree, k, kd_tree_points, degree_fitting, normals,
+          parallel_callback.advancement(),
+          parallel_callback.interrupted());
      tbb::parallel_for(tbb::blocked_range<size_t>(0, kd_tree_points.size ()), f);
-     unsigned int i = 0;
+     std::size_t i = 0;
      for(it = points.begin(); it != points.end(); ++ it, ++ i)
-       {
-	 put (normal_map, *it, normals[i]);
-       }
+       put (normal_map, *it, normals[i]);
+
+     parallel_callback.join();
    }
    else
 #endif
      {
-       for(it = points.begin(); it != points.end(); it++)
+       std::size_t nb = 0;
+       for(it = points.begin(); it != points.end(); it++, ++ nb)
 	 {
 	   Vector normal = internal::jet_estimate_normal<Kernel,SvdTraits,Tree>(
 										get(point_map,*it), 
 										tree, k, degree_fitting);
 
 	   put(normal_map, *it, normal); // normal_map[it] = normal
-    
-	 }
+           if (callback && !callback ((nb+1) / double(kd_tree_points.size())))
+             break;
+    	 }
      }
 
 
diff --git a/Point_set_processing_3/include/CGAL/jet_smooth_point_set.h b/Point_set_processing_3/include/CGAL/jet_smooth_point_set.h
index 83edb6dd1c0..c03062c64d9 100644
--- a/Point_set_processing_3/include/CGAL/jet_smooth_point_set.h
+++ b/Point_set_processing_3/include/CGAL/jet_smooth_point_set.h
@@ -39,9 +39,11 @@
 #include <list>
 
 #ifdef CGAL_LINKED_WITH_TBB
+#include <CGAL/internal/Parallel_callback.h>
 #include <tbb/parallel_for.h>
 #include <tbb/blocked_range.h>
 #include <tbb/scalable_allocator.h>  
+#include <tbb/atomic.h>
 #endif // CGAL_LINKED_WITH_TBB
 
 namespace CGAL {
@@ -125,20 +127,31 @@ jet_smooth_point(
     unsigned int degree_monge;
     const std::vector<Point>& input;
     std::vector<Point>& output;
+    tbb::atomic<std::size_t>& advancement;
+    tbb::atomic<bool>& interrupted;
 
   public:
     Jet_smooth_pwns (Tree& tree, unsigned int k, std::vector<Point>& points,
-		     unsigned int degree_fitting, unsigned int degree_monge, std::vector<Point>& output)
+		     unsigned int degree_fitting, unsigned int degree_monge, std::vector<Point>& output,
+                     tbb::atomic<std::size_t>& advancement,
+                     tbb::atomic<bool>& interrupted)
       : tree(tree), k (k), degree_fitting (degree_fitting),
 	degree_monge (degree_monge), input (points), output (output)
+      , advancement (advancement)
+      , interrupted (interrupted)
     { }
     
     void operator()(const tbb::blocked_range<std::size_t>& r) const
     {
       for( std::size_t i = r.begin(); i != r.end(); ++i)
+      {
+        if (interrupted)
+          break;
 	output[i] = CGAL::internal::jet_smooth_point<Kernel, SvdTraits>(input[i], tree, k,
 									degree_fitting,
 									degree_monge);
+        ++ advancement;
+      }
     }
 
   };
@@ -210,6 +223,8 @@ jet_smooth_point_set(
   PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
   unsigned int degree_fitting = choose_param(get_param(np, internal_np::degree_fitting), 2);
   unsigned int degree_monge = choose_param(get_param(np, internal_np::degree_monge), 2);
+  const cpp11::function<bool(double)>& callback = choose_param(get_param(np, internal_np::callback),
+                                                               cpp11::function<bool(double)>());
 
   typedef typename Kernel::Point_3 Point;
 
@@ -244,27 +259,35 @@ jet_smooth_point_set(
 #else
    if (boost::is_convertible<ConcurrencyTag,Parallel_tag>::value)
    {
+     internal::Point_set_processing_3::Parallel_callback
+       parallel_callback (callback, kd_tree_points.size());
+     
      std::vector<Point> mutated_points (kd_tree_points.size ());
      CGAL::internal::Jet_smooth_pwns<Kernel, SvdTraits, Tree>
        f (tree, k, kd_tree_points, degree_fitting, degree_monge,
-	  mutated_points);
+	  mutated_points,
+          parallel_callback.advancement(),
+          parallel_callback.interrupted());
      tbb::parallel_for(tbb::blocked_range<size_t>(0, kd_tree_points.size ()), f);
      unsigned int i = 0;
      for(it = points.begin(); it != points.end(); ++ it, ++ i)
-       {
-	 put(point_map, *it, mutated_points[i]);
+       put(point_map, *it, mutated_points[i]);
+
+     parallel_callback.join();
 
-       }
    }
    else
 #endif
      {
-       for(it = points.begin(); it != points.end(); it++)
+       std::size_t nb = 0;
+       for(it = points.begin(); it != points.end(); it++, ++ nb)
 	 {
 	   const typename boost::property_traits<PointMap>::reference p = get(point_map, *it);
 	   put(point_map, *it ,
 	       internal::jet_smooth_point<Kernel, SvdTraits>(
 							     p,tree,k,degree_fitting,degree_monge) );
+           if (callback && !callback ((nb+1) / double(kd_tree_points.size())))
+             break;
 	 }
      }
 }
diff --git a/Point_set_processing_3/include/CGAL/pca_estimate_normals.h b/Point_set_processing_3/include/CGAL/pca_estimate_normals.h
index 6011e1428cb..022837fc5d8 100644
--- a/Point_set_processing_3/include/CGAL/pca_estimate_normals.h
+++ b/Point_set_processing_3/include/CGAL/pca_estimate_normals.h
@@ -41,9 +41,11 @@
 #include <list>
 
 #ifdef CGAL_LINKED_WITH_TBB
+#include <CGAL/internal/Parallel_callback.h>
 #include <tbb/parallel_for.h>
 #include <tbb/blocked_range.h>
 #include <tbb/scalable_allocator.h>  
+#include <tbb/atomic.h>
 #endif // CGAL_LINKED_WITH_TBB
 
 namespace CGAL {
@@ -117,17 +119,28 @@ pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute
     const unsigned int k;
     const std::vector<Point>& input;
     std::vector<Vector>& output;
+    tbb::atomic<std::size_t>& advancement;
+    tbb::atomic<bool>& interrupted;
 
   public:
     PCA_estimate_normals(Tree& tree, unsigned int k, std::vector<Point>& points,
-			 std::vector<Vector>& output)
+			 std::vector<Vector>& output,
+                     tbb::atomic<std::size_t>& advancement,
+                     tbb::atomic<bool>& interrupted)
       : tree(tree), k (k), input (points), output (output)
+      , advancement (advancement)
+      , interrupted (interrupted)
     { }
     
     void operator()(const tbb::blocked_range<std::size_t>& r) const
     {
       for( std::size_t i = r.begin(); i != r.end(); ++i)
+      {
+        if (interrupted)
+          break;
 	output[i] = CGAL::internal::pca_estimate_normal<Kernel,Tree>(input[i], tree, k);
+        ++ advancement;
+      }
     }
 
   };
@@ -193,6 +206,8 @@ pca_estimate_normals(
 
   PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
   NormalMap normal_map = choose_param(get_param(np, internal_np::normal_map), NormalMap());
+  const cpp11::function<bool(double)>& callback = choose_param(get_param(np, internal_np::callback),
+                                                               cpp11::function<bool(double)>());
 
   typedef typename Kernel::Point_3 Point;
 
@@ -235,20 +250,26 @@ pca_estimate_normals(
 #else
   if (boost::is_convertible<ConcurrencyTag,Parallel_tag>::value)
     {
+      internal::Point_set_processing_3::Parallel_callback
+        parallel_callback (callback, kd_tree_points.size());
+     
       std::vector<Vector> normals (kd_tree_points.size ());
       CGAL::internal::PCA_estimate_normals<Kernel, Tree>
-	f (tree, k, kd_tree_points, normals);
+	f (tree, k, kd_tree_points, normals,
+           parallel_callback.advancement(),
+           parallel_callback.interrupted());
       tbb::parallel_for(tbb::blocked_range<size_t>(0, kd_tree_points.size ()), f);
       unsigned int i = 0;
       for(it = points.begin(); it != points.end(); ++ it, ++ i)
-	{
-	  put (normal_map, *it, normals[i]);
-	}
+        put (normal_map, *it, normals[i]);
+
+      parallel_callback.join();
     }
   else
 #endif
     {
-      for(it = points.begin(); it != points.end(); it++)
+      std::size_t nb = 0;
+      for(it = points.begin(); it != points.end(); it++, ++ nb)
 	{
 	  Vector normal = internal::pca_estimate_normal<Kernel,Tree>(      
 								     get(point_map,*it),
@@ -256,6 +277,8 @@ pca_estimate_normals(
 								     k);
 
 	  put(normal_map, *it, normal); // normal_map[it] = normal
+          if (callback && !callback ((nb+1) / double(kd_tree_points.size())))
+            break;
 	}
     }
    
diff --git a/Point_set_processing_3/include/CGAL/remove_outliers.h b/Point_set_processing_3/include/CGAL/remove_outliers.h
index d268a26adf4..9980ec462b2 100644
--- a/Point_set_processing_3/include/CGAL/remove_outliers.h
+++ b/Point_set_processing_3/include/CGAL/remove_outliers.h
@@ -29,6 +29,7 @@
 #include <CGAL/Orthogonal_k_neighbor_search.h>
 #include <CGAL/property_map.h>
 #include <CGAL/point_set_processing_assertions.h>
+#include <CGAL/function.h>
 
 #include <CGAL/boost/graph/named_function_params.h>
 #include <CGAL/boost/graph/named_params_helper.h>
@@ -163,6 +164,8 @@ remove_outliers(
   PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap());
   double threshold_percent = choose_param(get_param(np, internal_np::threshold_percent), 10.);
   double threshold_distance = choose_param(get_param(np, internal_np::threshold_distance), 0.);
+  const cpp11::function<bool(double)>& callback = choose_param(get_param(np, internal_np::callback),
+                                                               cpp11::function<bool(double)>());
   
   typedef typename Kernel::FT FT;
   
@@ -198,12 +201,15 @@ remove_outliers(
 
   // iterate over input points and add them to multimap sorted by distance to k
   std::multimap<FT,Enriched_point> sorted_points;
-  for(it = points.begin(); it != points.end(); it++)
+  std::size_t nb = 0;
+  for(it = points.begin(); it != points.end(); it++, ++ nb)
   {
     FT sq_distance = internal::compute_avg_knn_sq_distance_3<Kernel>(
       get(point_map,*it),
       tree, k);
     sorted_points.insert( std::make_pair(sq_distance, *it) );
+    if (callback && !callback ((nb+1) / double(kd_tree_points.size())))
+      break;
   }
 
   // Replaces [points.begin(), points.end()) range by the multimap content.
diff --git a/Point_set_processing_3/include/CGAL/wlop_simplify_and_regularize_point_set.h b/Point_set_processing_3/include/CGAL/wlop_simplify_and_regularize_point_set.h
index c47bb214eff..d9ecad56d7e 100644
--- a/Point_set_processing_3/include/CGAL/wlop_simplify_and_regularize_point_set.h
+++ b/Point_set_processing_3/include/CGAL/wlop_simplify_and_regularize_point_set.h
@@ -42,8 +42,11 @@
 #include <ctime>
 
 #ifdef CGAL_LINKED_WITH_TBB
+#include <CGAL/internal/Parallel_callback.h>
 #include <tbb/parallel_for.h>
 #include <tbb/blocked_range.h>
+#include <tbb/scalable_allocator.h>  
+#include <tbb/atomic.h>
 #endif // CGAL_LINKED_WITH_TBB
 
 #include <CGAL/Simple_cartesian.h>
@@ -354,6 +357,8 @@ class Sample_point_updater
   const typename Kernel::FT radius;  
   const std::vector<typename Kernel::FT> &original_densities;
   const std::vector<typename Kernel::FT> &sample_densities; 
+  tbb::atomic<std::size_t>& advancement;
+  tbb::atomic<bool>& interrupted;
 
 public:
   Sample_point_updater(
@@ -363,20 +368,25 @@ public:
     const Tree &_sample_kd_tree,              
     const typename Kernel::FT _radius,
     const std::vector<typename Kernel::FT> &_original_densities,
-    const std::vector<typename Kernel::FT> &_sample_densities): 
+    const std::vector<typename Kernel::FT> &_sample_densities,
+    tbb::atomic<std::size_t>& advancement,
+    tbb::atomic<bool>& interrupted):
   update_sample_points(out), 
     sample_points(in),
     original_kd_tree(_original_kd_tree),
     sample_kd_tree(_sample_kd_tree),
     radius(_radius),
     original_densities(_original_densities),
-    sample_densities(_sample_densities){} 
-
+    sample_densities(_sample_densities),
+    advancement (advancement),
+    interrupted (interrupted) {} 
 
   void operator() ( const tbb::blocked_range<size_t>& r ) const 
   { 
     for (size_t i = r.begin(); i != r.end(); ++i) 
     {
+      if (interrupted)
+        break;
       update_sample_points[i] = simplify_and_regularize_internal::
         compute_update_sample_point<Kernel, Tree, RandomAccessIterator>(
         sample_points[i], 
@@ -385,6 +395,7 @@ public:
         radius, 
         original_densities,
         sample_densities);
+      ++ advancement;
     }
   }
 };
@@ -466,6 +477,8 @@ wlop_simplify_and_regularize_point_set(
   double radius = choose_param(get_param(np, internal_np::neighbor_radius), -1);
   unsigned int iter_number = choose_param(get_param(np, internal_np::number_of_iterations), 35);
   bool require_uniform_sampling = choose_param(get_param(np, internal_np::require_uniform_sampling), false);
+  const cpp11::function<bool(double)>& callback = choose_param(get_param(np, internal_np::callback),
+                                                               cpp11::function<bool(double)>());
 
   typedef typename Kernel::Point_3   Point;
   typedef typename Kernel::FT        FT;
@@ -589,6 +602,9 @@ wlop_simplify_and_regularize_point_set(
     //parallel
     if (boost::is_convertible<ConcurrencyTag, Parallel_tag>::value)
     {
+      internal::Point_set_processing_3::Parallel_callback
+        parallel_callback (callback, iter_number * number_of_sample, iter_n * number_of_sample);
+     
       tbb::blocked_range<size_t> block(0, number_of_sample);
       Sample_point_updater<Kernel, Kd_Tree, typename PointRange::iterator> sample_updater(
                            update_sample_points,
@@ -597,15 +613,28 @@ wlop_simplify_and_regularize_point_set(
                            sample_kd_tree,
                            radius2,
                            original_density_weights,
-                           sample_density_weights);
+                           sample_density_weights,
+                           parallel_callback.advancement(),
+                           parallel_callback.interrupted());
 
        tbb::parallel_for(block, sample_updater);
+
+       bool interrupted = parallel_callback.interrupted();
+  
+       // We interrupt by hand as counter only goes halfway and won't terminate by itself
+       parallel_callback.interrupted() = true;
+       parallel_callback.join();       
+
+       // If interrupted during this step, nothing is computed, we return NaN
+       if (interrupted)
+         return output;
     }else
 #endif
     {
       //sequential
+      std::size_t nb = iter_n * number_of_sample;
       for (sample_iter = sample_points.begin();
-        sample_iter != sample_points.end(); ++sample_iter, ++update_iter)
+           sample_iter != sample_points.end(); ++sample_iter, ++update_iter, ++ nb)
       {
         *update_iter = simplify_and_regularize_internal::
           compute_update_sample_point<Kernel,
@@ -617,6 +646,8 @@ wlop_simplify_and_regularize_point_set(
                                        radius2,
                                        original_density_weights,
                                        sample_density_weights);
+        if (callback && !callback ((nb+1) / double(iter_number * number_of_sample)))
+          return output;
       }
     }