From ec8a00dce7ce5b5e70e8cff9d4cd64d7749093d2 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Fri, 3 Mar 2017 11:15:44 +0100 Subject: [PATCH] Parallelize local eigen analysis with TBB --- .../Classification/Local_eigen_analysis.h | 135 ++++++++++++++---- 1 file changed, 105 insertions(+), 30 deletions(-) diff --git a/Classification/include/CGAL/Classification/Local_eigen_analysis.h b/Classification/include/CGAL/Classification/Local_eigen_analysis.h index ef0597afedf..b804644fe0d 100644 --- a/Classification/include/CGAL/Classification/Local_eigen_analysis.h +++ b/Classification/include/CGAL/Classification/Local_eigen_analysis.h @@ -28,6 +28,13 @@ #include #include +#ifdef CGAL_LINKED_WITH_TBB +#include +#include +#include +#include +#endif // CGAL_LINKED_WITH_TBB + namespace CGAL { namespace Classification { @@ -49,7 +56,14 @@ namespace Classification { for matrix diagonalization. */ template > + typename DiagonalizeTraits = CGAL::Default_diagonalize_traits, +#if defined(DOXYGEN_RUNNING) + typename ConcurrencyTag> +#elif defined(CGAL_LINKED_WITH_TBB) + typename ConcurrencyTag = CGAL::Parallel_tag> +#else + typename ConcurrencyTag = CGAL::Sequential_tag> +#endif class Local_eigen_analysis { public: @@ -61,6 +75,53 @@ public: private: + +#ifdef CGAL_LINKED_WITH_TBB + template + class Compute_eigen_values + { + Local_eigen_analysis& m_eigen; + const PointRange& m_input; + PointMap m_point_map; + const NeighborQuery& m_neighbor_query; + double& m_mean_range; + tbb::mutex& m_mutex; + + public: + + Compute_eigen_values (Local_eigen_analysis& eigen, + const PointRange& input, + PointMap point_map, + const NeighborQuery& neighbor_query, + double& mean_range, + tbb::mutex& mutex) + : m_eigen (eigen), m_input (input), m_point_map (point_map), + m_neighbor_query (neighbor_query), m_mean_range (mean_range), m_mutex (mutex) + { } + + void operator()(const tbb::blocked_range& r) const + { + for (std::size_t i = r.begin(); i != r.end(); ++ i) + { + std::vector neighbors; + m_neighbor_query (get(m_point_map, *(m_input.begin()+i)), std::back_inserter (neighbors)); + + std::vector neighbor_points; + for (std::size_t j = 0; j < neighbors.size(); ++ j) + neighbor_points.push_back (get(m_point_map, *(m_input.begin()+neighbors[j]))); + + m_mutex.lock(); + m_mean_range += CGAL::sqrt (CGAL::squared_distance (get(m_point_map, *(m_input.begin() + i)), + get(m_point_map, *(m_input.begin() + neighbors.back())))); + m_mutex.unlock(); + + m_eigen.compute (i, get(m_point_map, *(m_input.begin()+i)), neighbor_points); + } + } + + }; +#endif + std::vector m_eigenvalues; std::vector m_sum_eigenvalues; std::vector m_centroids; @@ -91,29 +152,43 @@ public: PointMap point_map, const NeighborQuery& neighbor_query) { - m_eigenvalues.reserve (input.size()); - m_centroids.reserve (input.size()); - m_smallest_eigenvectors.reserve (input.size()); - m_middle_eigenvectors.reserve (input.size()); - m_largest_eigenvectors.reserve (input.size()); + m_eigenvalues.resize (input.size()); + m_sum_eigenvalues.resize (input.size()); + m_centroids.resize (input.size()); + m_smallest_eigenvectors.resize (input.size()); + m_middle_eigenvectors.resize (input.size()); + m_largest_eigenvectors.resize (input.size()); m_mean_range = 0.; - for (std::size_t i = 0; i < input.size(); i++) +#ifndef CGAL_LINKED_WITH_TBB + CGAL_static_assertion_msg (!(boost::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else + if (boost::is_convertible::value) { - std::vector neighbors; - neighbor_query (get(point_map, *(input.begin()+i)), std::back_inserter (neighbors)); - - std::vector neighbor_points; - for (std::size_t j = 0; j < neighbors.size(); ++ j) - neighbor_points.push_back (get(point_map, *(input.begin()+neighbors[j]))); - - m_mean_range += CGAL::sqrt (CGAL::squared_distance (get(point_map, *(input.begin() + i)), - get(point_map, *(input.begin() + neighbors.back())))); - - compute (get(point_map, *(input.begin()+i)), neighbor_points); + tbb::mutex mutex; + Compute_eigen_values f(*this, input, point_map, neighbor_query, m_mean_range, mutex); + tbb::parallel_for(tbb::blocked_range(0, input.size ()), f); } + else +#endif + { + for (std::size_t i = 0; i < input.size(); i++) + { + std::vector neighbors; + neighbor_query (get(point_map, *(input.begin()+i)), std::back_inserter (neighbors)); + std::vector neighbor_points; + for (std::size_t j = 0; j < neighbors.size(); ++ j) + neighbor_points.push_back (get(point_map, *(input.begin()+neighbors[j]))); + + m_mean_range += CGAL::sqrt (CGAL::squared_distance (get(point_map, *(input.begin() + i)), + get(point_map, *(input.begin() + neighbors.back())))); + + compute (i, get(point_map, *(input.begin()+i)), neighbor_points); + } + } m_mean_range /= input.size(); } @@ -143,21 +218,21 @@ public: private: - void compute (const Point& query, std::vector& neighbor_points) + void compute (std::size_t index, const Point& query, std::vector& neighbor_points) { if (neighbor_points.size() == 0) { Eigenvalues v = {{ 0., 0., 0. }}; - m_eigenvalues.push_back (v); - m_centroids.push_back (query); - m_smallest_eigenvectors.push_back (Vector (0., 0., 1.)); - m_middle_eigenvectors.push_back (Vector (0., 1., 0.)); - m_largest_eigenvectors.push_back (Vector (1., 0., 0.)); + m_eigenvalues[index] = v; + m_centroids[index] = query; + m_smallest_eigenvectors[index] = Vector (0., 0., 1.); + m_middle_eigenvectors[index] = Vector (0., 1., 0.); + m_largest_eigenvectors[index] = Vector (1., 0., 0.); return; } Point centroid = CGAL::centroid (neighbor_points.begin(), neighbor_points.end()); - m_centroids.push_back (centroid); + m_centroids[index] = centroid; CGAL::cpp11::array covariance = {{ 0., 0., 0., 0., 0., 0. }}; @@ -185,11 +260,11 @@ private: if (sum > 0.) for (std::size_t i = 0; i < 3; ++ i) evalues[i] = evalues[i] / sum; - m_sum_eigenvalues.push_back(sum); - m_eigenvalues.push_back (evalues); - m_smallest_eigenvectors.push_back (Vector (evectors[0], evectors[1], evectors[2])); - m_middle_eigenvectors.push_back (Vector (evectors[3], evectors[4], evectors[5])); - m_largest_eigenvectors.push_back (Vector (evectors[6], evectors[7], evectors[8])); + m_sum_eigenvalues[index] = sum; + m_eigenvalues[index] = evalues; + m_smallest_eigenvectors[index] = Vector (evectors[0], evectors[1], evectors[2]); + m_middle_eigenvectors[index] = Vector (evectors[3], evectors[4], evectors[5]); + m_largest_eigenvectors[index] = Vector (evectors[6], evectors[7], evectors[8]); }