parallelize function "compute_update_sample_point"

Point_set_processing_3/include/CGAL/wlop_simplify_and_regularize_point_set.h

@@ -61,7 +61,8 @@ namespace simplify_and_regularize_internal{
 /// @tparam Kernel Geometric traits class.
 ///
 /// @return average term vector
-template <typename Kernel>
+template <typename Concurrency_tag,
+          typename Kernel>
 typename Kernel::Point_3
 compute_update_sample_point(
   const typename Kernel::Point_3& query, ///< 3D point to project
@@ -74,6 +75,7 @@ compute_update_sample_point(
   const std::vector<typename Kernel::FT>& sample_densities ///< 
 )
 {
+
   CGAL_point_set_processing_precondition(radius > 0);
   CGAL_point_set_processing_precondition(neighbor_original_points.size() >= 1);
   bool is_original_densities_empty = original_densities.empty();
@@ -86,6 +88,8 @@ compute_update_sample_point(
   typedef rich_grid_internal::Rich_point<Kernel> Rich_point;
 
   FT radius2 = radius * radius;
+  size_t nb_neighbor_sample_points = neighbor_sample_points.size();
+  size_t nb_neighbor_original_points = neighbor_original_points.size();
 
   //Compute average term
   FT weight = (FT)0.0, average_weight_sum = (FT)0.0;
@@ -94,21 +98,52 @@ compute_update_sample_point(
 
   std::vector<Rich_point>::const_iterator iter;
   iter = neighbor_original_points.begin();
-  for (; iter != neighbor_original_points.end(); ++iter)
+  
+#ifdef CGAL_LINKED_WITH_TBB
+  if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
   {
-    const Point& np = iter->pt;
-    unsigned int idx_of_original = iter->index;
-
-    FT dist2 = CGAL::squared_distance(query, np);
-    weight = exp(dist2 * iradius16);
-
-    if(!is_original_densities_empty)
+    tbb::parallel_for(
+      tbb::blocked_range<size_t>(0,nb_neighbor_original_points),
+      [&](const tbb::blocked_range<size_t>&r)
     {
-      weight *= original_densities[idx_of_original];
-    }
+      for (size_t i = r.begin(); i != r.end(); i++)
+      {
+        const Point& np = neighbor_original_points[i].pt;
+        unsigned int idx_of_original = neighbor_original_points[i].index;
 
-    average_weight_sum += weight;
-    average = average + (np - CGAL::ORIGIN) * weight;
+        FT dist2 = CGAL::squared_distance(query, np);
+        weight = exp(dist2 * iradius16);
+
+        if(!is_original_densities_empty)
+        {
+          weight *= original_densities[idx_of_original];
+        }
+
+        average_weight_sum += weight;
+        average = average + (np - CGAL::ORIGIN) * weight;
+      }
+    }
+    );
+  }
+  else
+#endif
+  {
+    for (; iter != neighbor_original_points.end(); ++iter)
+    {
+      const Point& np = iter->pt;
+      unsigned int idx_of_original = iter->index;
+
+      FT dist2 = CGAL::squared_distance(query, np);
+      weight = exp(dist2 * iradius16);
+
+      if(!is_original_densities_empty)
+      {
+        weight *= original_densities[idx_of_original];
+      }
+
+      average_weight_sum += weight;
+      average = average + (np - CGAL::ORIGIN) * weight;
+    }
   }
 
   // Finishing compute average term
@@ -126,25 +161,59 @@ compute_update_sample_point(
   FT repulsion_weight_sum = (FT)0.0;
 
   Vector repulsion = CGAL::NULL_VECTOR; 
-  iter = neighbor_sample_points.begin();
-  for (; iter != neighbor_sample_points.end(); ++iter)
+
+#ifdef CGAL_LINKED_WITH_TBB
+  if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
   {
-    const Point& np = iter->pt;
-    unsigned int idx_of_sample = iter->index;
-
-    Vector diff = query - np;
-
-    FT dist2 = CGAL::squared_distance(query, np);
-    FT dist = std::sqrt(dist2);
-
-    weight = std::exp(dist2 * iradius16) * std::pow(FT(1.0)/dist, 2);
-    if(!is_sample_densities_empty)
+    tbb::parallel_for(
+      tbb::blocked_range<size_t>(0,nb_neighbor_sample_points),
+      [&](tbb::blocked_range<size_t>&r)
     {
-      weight *= sample_densities[idx_of_sample];
-    }
+      for (size_t i = r.begin(); i != r.end(); i++)
+      {
+        const Point& np = neighbor_sample_points[i].pt;
+        unsigned int idx_of_sample = neighbor_sample_points[i].index;
 
-    repulsion_weight_sum += weight;
-    repulsion = repulsion + diff * weight;
+        Vector diff = query - np;
+
+        FT dist2 = CGAL::squared_distance(query, np);
+        FT dist = std::sqrt(dist2);
+
+        weight = std::exp(dist2 * iradius16) * std::pow(FT(1.0)/dist, 2);
+        if(!is_sample_densities_empty)
+        {
+          weight *= sample_densities[idx_of_sample];
+        }
+
+        repulsion_weight_sum += weight;
+        repulsion = repulsion + diff * weight;
+      }
+    }
+    );
+  }
+  else
+#endif
+  {
+    iter = neighbor_sample_points.begin();
+    for (; iter != neighbor_sample_points.end(); ++iter)
+    {
+      const Point& np = iter->pt;
+      unsigned int idx_of_sample = iter->index;
+
+      Vector diff = query - np;
+
+      FT dist2 = CGAL::squared_distance(query, np);
+      FT dist = std::sqrt(dist2);
+
+      weight = std::exp(dist2 * iradius16) * std::pow(FT(1.0)/dist, 2);
+      if(!is_sample_densities_empty)
+      {
+        weight *= sample_densities[idx_of_sample];
+      }
+
+      repulsion_weight_sum += weight;
+      repulsion = repulsion + diff * weight;
+    }
   }
 
   // Finishing compute repulsion term
@@ -465,9 +534,9 @@ wlop_simplify_and_regularize_point_set(
   std::vector<Point>::iterator update_sample_iter;
   for (unsigned int iteration = 0; iteration < max_iter_number; iteration++)
   {
-  #ifdef CGAL_DEBUG_MODE
-    task_timer.reset();
-  #endif
+    #ifdef CGAL_DEBUG_MODE
+      task_timer.reset();
+    #endif
     // Build rich-grid for sample-sample neighborhood
     rich_grid_internal::compute_ball_neighbors_one_self(sample_rich_points,
                                                         bbox, 
@@ -512,7 +581,7 @@ wlop_simplify_and_regularize_point_set(
         tbb::blocked_range<size_t>(0,nb_points_sample),
         [&](const tbb::blocked_range<size_t>&r)
       {
-        //parallel function "for" begin
+        //parallel section "for" begin
         for (size_t i = r.begin(); i != r.end(); i++)
         {
           Point& p = sample_points[i];
@@ -541,16 +610,16 @@ wlop_simplify_and_regularize_point_set(
             Rich_point rp(sample_points[idx_of_sample], idx_of_sample);
             rich_sample_neighbors.push_back(rp);
           }
-
+          
           update_sample_points[i] = simplify_and_regularize_internal::
-            compute_update_sample_point<Kernel>
+            compute_update_sample_point<CGAL::Parallel_tag, Kernel>
             (sample_points[i],
              rich_original_neighbors,
              rich_sample_neighbors,
              neighbor_radius,
              original_densities,
              sample_densities);
-        }//parallel function "for" end
+          }//parallel section "for" end
         }
       );
     }
@@ -591,13 +660,12 @@ wlop_simplify_and_regularize_point_set(
         *update_sample_iter = simplify_and_regularize_internal::
           compute_update_sample_point<Kernel>
           (*sample_iter,
-            rich_original_neighbors,
-            rich_sample_neighbors,
-            neighbor_radius,
-            original_densities,
-            sample_densities);
+           rich_original_neighbors,
+           rich_sample_neighbors,
+           neighbor_radius,
+           original_densities,
+           sample_densities);
       }
-
     }//parallel judge "else" end
 
     sample_iter = sample_points.begin();