diff --git a/BGL/include/CGAL/boost/graph/parameters_interface.h b/BGL/include/CGAL/boost/graph/parameters_interface.h
index 5f9d6865af7..30a8e83cf93 100644
--- a/BGL/include/CGAL/boost/graph/parameters_interface.h
+++ b/BGL/include/CGAL/boost/graph/parameters_interface.h
@@ -212,4 +212,5 @@ CGAL_add_named_parameter(maximum_normal_deviation_t, maximum_normal_deviation, m
 CGAL_add_named_parameter(distance_threshold_t, distance_threshold, distance_threshold)
 CGAL_add_named_parameter(angle_deg_threshold_t, angle_deg_threshold, angle_deg_threshold)
 CGAL_add_named_parameter(cos_value_threshold_t, cos_value_threshold, cos_value_threshold)
-CGAL_add_named_parameter(min_region_size_t, min_region_size, min_region_size)
\ No newline at end of file
+CGAL_add_named_parameter(min_region_size_t, min_region_size, min_region_size)
+CGAL_add_named_parameter(sort_regions_t, sort_regions, sort_regions)
\ No newline at end of file
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/K_neighbor_query.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/K_neighbor_query.h
index abc6c0e5ef2..4d2aea14b07 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/K_neighbor_query.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/K_neighbor_query.h
@@ -56,7 +56,7 @@ namespace Point_set {
     a model of `ConstRange` whose iterator type is `RandomAccessIterator`
 
     \tparam PointMap
-    a model of `ReadPropertyMap` whose key type is the value type of the input
+    a model of `LValuePropertyMap` whose key type is the value type of the input
     range and value type is `Kernel::Point_2` or `Kernel::Point_3`
 
     \cgalModels `NeighborQuery`
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_line_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_line_fit_region.h
index f0cdced3a64..19dd6001f4a 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_line_fit_region.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_line_fit_region.h
@@ -16,21 +16,10 @@
 
 #include <CGAL/license/Shape_detection.h>
 
-// STL includes.
-#include <vector>
-
 // Boost includes.
 #include <CGAL/boost/graph/named_params_helper.h>
 #include <CGAL/boost/graph/Named_function_parameters.h>
 
-// CGAL includes.
-#include <CGAL/assertions.h>
-#include <CGAL/number_utils.h>
-#include <CGAL/Cartesian_converter.h>
-#include <CGAL/Eigen_diagonalize_traits.h>
-#include <CGAL/linear_least_squares_fitting_2.h>
-#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-
 // Internal includes.
 #include <CGAL/Shape_detection/Region_growing/internal/utils.h>
 
@@ -56,11 +45,11 @@ namespace Point_set {
     a model of `ConstRange` whose iterator type is `RandomAccessIterator`
 
     \tparam PointMap
-    a model of `ReadPropertyMap` whose key type is the value type of the input
+    a model of `LValuePropertyMap` whose key type is the value type of the input
     range and value type is `Kernel::Point_2`.
 
     \tparam NormalMap
-    a model of `ReadPropertyMap` whose key type is the value type of the input
+    a model of `LValuePropertyMap` whose key type is the value type of the input
     range and value type is `Kernel::Vector_2`.
 
     \cgalModels `RegionType`
@@ -93,19 +82,10 @@ namespace Point_set {
     using Vector_2 = typename Traits::Vector_2;
     using Line_2 = typename Traits::Line_2;
 
-    using ITraits = Exact_predicates_inexact_constructions_kernel;
-    using IFT = typename ITraits::FT;
-    using IPoint_2 = typename ITraits::Point_2;
-    using ILine_2 = typename ITraits::Line_2;
-    using IConverter = Cartesian_converter<Traits, ITraits>;
-
     using Squared_length_2 = typename Traits::Compute_squared_length_2;
     using Squared_distance_2 = typename Traits::Compute_squared_distance_2;
     using Scalar_product_2 = typename Traits::Compute_scalar_product_2;
 
-    using Get_sqrt = internal::Get_sqrt<Traits>;
-    using Sqrt = typename Get_sqrt::Sqrt;
-
   public:
     /// \name Initialization
     /// @{
@@ -155,9 +135,7 @@ namespace Point_set {
     m_normal_map(normal_map),
     m_squared_length_2(traits.compute_squared_length_2_object()),
     m_squared_distance_2(traits.compute_squared_distance_2_object()),
-    m_scalar_product_2(traits.compute_scalar_product_2_object()),
-    m_sqrt(Get_sqrt::sqrt_object(traits)),
-    m_iconverter() {
+    m_scalar_product_2(traits.compute_scalar_product_2_object()) {
 
       CGAL_precondition(input_range.size() > 0);
       m_distance_threshold = parameters::choose_parameter(
@@ -177,6 +155,9 @@ namespace Point_set {
       m_cos_value_threshold = parameters::choose_parameter(
         parameters::get_parameter(np, internal_np::cos_value_threshold), cos_value_threshold);
       CGAL_precondition(m_cos_value_threshold >= FT(0) && m_cos_value_threshold <= FT(1));
+
+      m_sort_regions = parameters::choose_parameter(
+        parameters::get_parameter(np, internal_np::sort_regions), false);
     }
 
     /// @}
@@ -272,31 +253,10 @@ namespace Point_set {
 
       } else { // update reference line and normal
 
-        std::vector<IPoint_2> points;
-        points.reserve(region.size());
-        for (const std::size_t point_index : region) {
-          CGAL_precondition(point_index < m_input_range.size());
-          const auto& key = *(m_input_range.begin() + point_index);
-          const Point_2& point = get(m_point_map, key);
-          points.push_back(m_iconverter(point));
-        }
-        CGAL_postcondition(points.size() == region.size());
-
-        ILine_2 fitted_line;
-        IPoint_2 fitted_centroid;
-
         // The best fit line will be a line fitted to all region points with
         // its normal being perpendicular to the line.
-        CGAL::linear_least_squares_fitting_2(
-          points.begin(), points.end(),
-          fitted_line, fitted_centroid,
-          CGAL::Dimension_tag<0>(), ITraits(),
-          CGAL::Eigen_diagonalize_traits<IFT, 2>());
-
-        m_line_of_best_fit = Line_2(
-          static_cast<FT>(fitted_line.a()),
-          static_cast<FT>(fitted_line.b()),
-          static_cast<FT>(fitted_line.c()));
+        internal::create_line_from_points_2(
+          m_input_range, m_point_map, region, m_line_of_best_fit);
         m_normal_of_best_fit = m_line_of_best_fit.perpendicular(
           m_line_of_best_fit.point(0)).to_vector();
       }
@@ -310,6 +270,7 @@ namespace Point_set {
     FT m_distance_threshold;
     FT m_cos_value_threshold;
     std::size_t m_min_region_size;
+    bool m_sort_regions;
 
     const Point_map m_point_map;
     const Normal_map m_normal_map;
@@ -317,9 +278,6 @@ namespace Point_set {
     const Squared_length_2 m_squared_length_2;
     const Squared_distance_2 m_squared_distance_2;
     const Scalar_product_2 m_scalar_product_2;
-    const Sqrt m_sqrt;
-
-    const IConverter m_iconverter;
 
     Line_2 m_line_of_best_fit;
     Vector_2 m_normal_of_best_fit;
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_line_fit_sorting.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_line_fit_sorting.h
index b0340cc3d9c..f8a3bda74ed 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_line_fit_sorting.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_line_fit_sorting.h
@@ -16,17 +16,6 @@
 
 #include <CGAL/license/Shape_detection.h>
 
-// STL includes.
-#include <vector>
-#include <algorithm>
-
-// CGAL includes.
-#include <CGAL/assertions.h>
-#include <CGAL/Cartesian_converter.h>
-#include <CGAL/Eigen_diagonalize_traits.h>
-#include <CGAL/linear_least_squares_fitting_2.h>
-#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-
 // Internal includes.
 #include <CGAL/Shape_detection/Region_growing/internal/utils.h>
 #include <CGAL/Shape_detection/Region_growing/internal/property_map.h>
@@ -53,7 +42,7 @@ namespace Point_set {
     a model of `NeighborQuery`
 
     \tparam PointMap
-    a model of `ReadPropertyMap` whose key type is the value type of the input
+    a model of `LValuePropertyMap` whose key type is the value type of the input
     range and value type is `Kernel::Point_2`
   */
   template<
@@ -77,7 +66,7 @@ namespace Point_set {
 
     #ifdef DOXYGEN_RUNNING
       /*!
-        a model of `ReadPropertyMap` whose key and value type is `std::size_t`.
+        a model of `ReadablePropertyMap` whose key and value type is `std::size_t`.
         This map provides an access to the ordered indices of input points.
       */
       typedef unspecified_type Seed_map;
@@ -86,12 +75,9 @@ namespace Point_set {
     /// @}
 
   private:
-    using ITraits = Exact_predicates_inexact_constructions_kernel;
-    using IFT = typename ITraits::FT;
-    using IPoint_2 = typename ITraits::Point_2;
-    using ILine_2 = typename ITraits::Line_2;
-    using IConverter = Cartesian_converter<Traits, ITraits>;
-    using Compare_scores = internal::Compare_scores<IFT>;
+    using FT = typename Traits::FT;
+    using Line_2 = typename Traits::Line_2;
+    using Compare_scores = internal::Compare_scores<FT>;
 
   public:
     /// \name Initialization
@@ -119,8 +105,7 @@ namespace Point_set {
       const PointMap point_map = PointMap()) :
     m_input_range(input_range),
     m_neighbor_query(neighbor_query),
-    m_point_map(point_map),
-    m_iconverter() {
+    m_point_map(point_map) {
 
       CGAL_precondition(input_range.size() > 0);
       m_order.resize(m_input_range.size());
@@ -164,35 +149,18 @@ namespace Point_set {
     Neighbor_query& m_neighbor_query;
     const Point_map m_point_map;
     std::vector<std::size_t> m_order;
-    std::vector<IFT> m_scores;
-    const IConverter m_iconverter;
+    std::vector<FT> m_scores;
 
     void compute_scores() {
 
-      std::vector<IPoint_2> points;
+      Line_2 line;
       std::vector<std::size_t> neighbors;
       for (std::size_t i = 0; i < m_input_range.size(); ++i) {
-
         neighbors.clear();
         m_neighbor_query(i, neighbors);
         neighbors.push_back(i);
-
-        points.clear();
-        for (const std::size_t point_index : neighbors) {
-          CGAL_precondition(point_index < m_input_range.size());
-          const auto& key = *(m_input_range.begin() + point_index);
-          const auto& point = get(m_point_map, key);
-          points.push_back(m_iconverter(point));
-        }
-        CGAL_postcondition(points.size() == neighbors.size());
-
-        ILine_2  fitted_line;
-        IPoint_2 fitted_centroid;
-        m_scores[i] = CGAL::linear_least_squares_fitting_2(
-          points.begin(), points.end(),
-          fitted_line, fitted_centroid,
-          CGAL::Dimension_tag<0>(), ITraits(),
-          CGAL::Eigen_diagonalize_traits<IFT, 2>());
+        m_scores[i] = internal::create_line_from_points_2(
+          m_input_range, m_point_map, neighbors, line);
       }
     }
   };
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_plane_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_plane_fit_region.h
index 3aedc2754c9..d3e6756dd37 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_plane_fit_region.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_plane_fit_region.h
@@ -56,11 +56,11 @@ namespace Point_set {
     a model of `ConstRange` whose iterator type is `RandomAccessIterator`
 
     \tparam PointMap
-    a model of `ReadPropertyMap` whose key type is the value type of the input
+    a model of `LValuePropertyMap` whose key type is the value type of the input
     range and value type is `Kernel::Point_3`
 
     \tparam NormalMap
-    a model of `ReadPropertyMap` whose key type is the value type of the input
+    a model of `LValuePropertyMap` whose key type is the value type of the input
     range and value type is `Kernel::Vector_3`
 
     \cgalModels `RegionType`
@@ -177,6 +177,9 @@ namespace Point_set {
       m_cos_value_threshold = parameters::choose_parameter(
         parameters::get_parameter(np, internal_np::cos_value_threshold), cos_value_threshold);
       CGAL_precondition(m_cos_value_threshold >= FT(0) && m_cos_value_threshold <= FT(1));
+
+      m_sort_regions = parameters::choose_parameter(
+        parameters::get_parameter(np, internal_np::sort_regions), false);
     }
 
     /// @}
@@ -310,6 +313,7 @@ namespace Point_set {
     FT m_distance_threshold;
     FT m_cos_value_threshold;
     std::size_t m_min_region_size;
+    bool m_sort_regions;
 
     const Point_map m_point_map;
     const Normal_map m_normal_map;
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_plane_fit_sorting.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_plane_fit_sorting.h
index 6a4b7c5ad5a..10117f533fa 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_plane_fit_sorting.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_plane_fit_sorting.h
@@ -53,7 +53,7 @@ namespace Point_set {
     a model of `NeighborQuery`
 
     \tparam PointMap
-    a model of `ReadPropertyMap` whose key type is the value type of the input
+    a model of `LValuePropertyMap` whose key type is the value type of the input
     range and value type is `Kernel::Point_3`
   */
   template<
@@ -77,7 +77,7 @@ namespace Point_set {
 
     #ifdef DOXYGEN_RUNNING
       /*!
-        a model of `ReadPropertyMap` whose key and value type is `std::size_t`.
+        a model of `ReadablePropertyMap` whose key and value type is `std::size_t`.
         This map provides an access to the ordered indices of input points.
       */
       typedef unspecified_type Seed_map;
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Sphere_neighbor_query.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Sphere_neighbor_query.h
index f419cda4ffb..ceaa31ab38c 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Sphere_neighbor_query.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Sphere_neighbor_query.h
@@ -57,7 +57,7 @@ namespace Point_set {
     a model of `ConstRange` whose iterator type is `RandomAccessIterator`
 
     \tparam PointMap
-    a model of `ReadPropertyMap` whose key type is the value type of the input
+    a model of `LValuePropertyMap` whose key type is the value type of the input
     range and value type is `Kernel::Point_2` or `Kernel::Point_3`
 
     \cgalModels `NeighborQuery`
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_polygon_mesh/Least_squares_plane_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_polygon_mesh/Least_squares_plane_fit_region.h
index 2e6709da7f9..5b248cce1c2 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_polygon_mesh/Least_squares_plane_fit_region.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_polygon_mesh/Least_squares_plane_fit_region.h
@@ -69,7 +69,7 @@ namespace Polygon_mesh {
     value type is the face type of a polygon mesh
 
     \tparam VertexToPointMap
-    a model of `ReadPropertyMap` whose key type is the vertex type of a polygon mesh and
+    a model of `LValuePropertyMap` whose key type is the vertex type of a polygon mesh and
     value type is `Kernel::Point_3`
 
     \cgalModels `RegionType`
@@ -182,6 +182,9 @@ namespace Polygon_mesh {
       m_cos_value_threshold = parameters::choose_parameter(
         parameters::get_parameter(np, internal_np::cos_value_threshold), cos_value_threshold);
       CGAL_precondition(m_cos_value_threshold >= FT(0) && m_cos_value_threshold <= FT(1));
+
+      m_sort_regions = parameters::choose_parameter(
+        parameters::get_parameter(np, internal_np::sort_regions), false);
     }
 
     /// @}
@@ -354,6 +357,7 @@ namespace Polygon_mesh {
     FT m_distance_threshold;
     FT m_cos_value_threshold;
     std::size_t m_min_region_size;
+    bool m_sort_regions;
 
     const Vertex_to_point_map m_vertex_to_point_map;
 
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_polygon_mesh/Least_squares_plane_fit_sorting.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_polygon_mesh/Least_squares_plane_fit_sorting.h
index 1ee3286367c..39512dbca8b 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_polygon_mesh/Least_squares_plane_fit_sorting.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_polygon_mesh/Least_squares_plane_fit_sorting.h
@@ -65,7 +65,7 @@ namespace Polygon_mesh {
     value type is the face type of a polygon mesh
 
     \tparam VertexToPointMap
-    a model of `ReadPropertyMap` whose key type is the vertex type of a polygon mesh and
+    a model of `LValuePropertyMap` whose key type is the vertex type of a polygon mesh and
     value type is `Kernel::Point_3`
   */
   template<
@@ -91,7 +91,7 @@ namespace Polygon_mesh {
 
     #ifdef DOXYGEN_RUNNING
       /*!
-        a model of `ReadPropertyMap` whose key and value type is `std::size_t`.
+        a model of `ReadablePropertyMap` whose key and value type is `std::size_t`.
         This map provides an access to the ordered indices of polygon mesh faces.
       */
       typedef unspecified_type Seed_map;
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/property_map.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/property_map.h
index d8c590a139b..4cdd5d9b160 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/property_map.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/property_map.h
@@ -52,15 +52,12 @@ namespace internal {
 
     reference operator[](const key_type item_index) const {
       CGAL_precondition(item_index < m_item_range.size());
-
       const auto& key = *(m_item_range.begin() + item_index);
       return get(m_property_map, key);
     }
 
     friend inline reference get(
-      const Item_property_map& item_map,
-      const key_type key) {
-
+      const Item_property_map& item_map, const key_type key) {
       return item_map[key];
     }
 
@@ -80,7 +77,7 @@ namespace internal {
 
     using value_type = std::size_t;
     using key_type = Item;
-    using category = boost::lvalue_property_map_tag;
+    using category = boost::readable_property_map_tag;
 
     using Item_map = std::map<key_type, value_type>;
 
@@ -103,9 +100,7 @@ namespace internal {
     }
 
     friend inline value_type get(
-      const Item_to_index_property_map& item_to_index_map,
-      const key_type& key) {
-
+      const Item_to_index_property_map& item_to_index_map, const key_type& key) {
       return item_to_index_map[key];
     }
 
@@ -119,7 +114,7 @@ namespace internal {
   public:
     using key_type = std::size_t;
     using value_type = std::size_t;
-    using category = boost::lvalue_property_map_tag;
+    using category = boost::readable_property_map_tag;
 
     Seed_property_map(
       const std::vector<std::size_t>& seeds) :
@@ -131,9 +126,7 @@ namespace internal {
     }
 
     friend value_type get(
-      const Seed_property_map& seed_map,
-      const key_type key) {
-
+      const Seed_property_map& seed_map, const key_type key) {
       return seed_map[key];
     }
 
@@ -141,11 +134,7 @@ namespace internal {
     const std::vector<std::size_t>& m_seeds;
   };
 
-} // namespace internal
-
-namespace RG {
-
-  class Point_to_shape_index_map {
+  class Point_to_region_index_map {
 
   public:
     using key_type = std::size_t;
@@ -153,24 +142,24 @@ namespace RG {
     using reference = value_type;
     using category = boost::readable_property_map_tag;
 
-    Point_to_shape_index_map() { }
+    Point_to_region_index_map() { }
 
     template<typename PointRange>
-    Point_to_shape_index_map(
+    Point_to_region_index_map(
       const PointRange& points,
       const std::vector< std::vector<std::size_t> >& regions) :
     m_indices(new std::vector<int>(points.size(), -1)) {
 
       for (std::size_t i = 0; i < regions.size(); ++i)
-        for (const std::size_t idx : regions[i])
-          (*m_indices)[idx] = static_cast<int>(i);
+        for (const std::size_t index : regions[i])
+          (*m_indices)[index] = static_cast<int>(i);
     }
 
     inline friend value_type get(
-      const Point_to_shape_index_map& point_to_shape_index_map,
+      const Point_to_region_index_map& point_to_region_index_map,
       const key_type key) {
 
-      const auto& indices = *(point_to_shape_index_map.m_indices);
+      const auto& indices = *(point_to_region_index_map.m_indices);
       return indices[key];
     }
 
@@ -178,7 +167,7 @@ namespace RG {
     std::shared_ptr< std::vector<int> > m_indices;
   };
 
-} // namespace RG
+} // namespace internal
 } // namespace Shape_detection
 } // namespace CGAL
 
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h
index 69ff36de4c3..426a7f11ed1 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h
@@ -18,6 +18,7 @@
 
 // STL includes.
 #include <vector>
+#include <algorithm>
 
 // Boost headers.
 #include <boost/mpl/has_xxx.hpp>
@@ -27,6 +28,7 @@
 #include <CGAL/number_utils.h>
 #include <CGAL/Cartesian_converter.h>
 #include <CGAL/Eigen_diagonalize_traits.h>
+#include <CGAL/linear_least_squares_fitting_2.h>
 #include <CGAL/linear_least_squares_fitting_3.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
@@ -42,8 +44,7 @@ namespace internal {
     using FT = typename Traits::FT;
 
   public:
-    FT operator()(const FT value) const {
-
+    const FT operator()(const FT value) const {
       CGAL_precondition(value >= FT(0));
       return static_cast<FT>(CGAL::sqrt(CGAL::to_double(value)));
     }
@@ -80,29 +81,74 @@ namespace internal {
 
   template<typename FT>
   struct Compare_scores {
-    const std::vector<FT>& m_scores;
 
+    const std::vector<FT>& m_scores;
     Compare_scores(const std::vector<FT>& scores) :
     m_scores(scores)
     { }
 
     bool operator()(const std::size_t i, const std::size_t j) const {
-
       CGAL_precondition(i < m_scores.size());
       CGAL_precondition(j < m_scores.size());
       return m_scores[i] > m_scores[j];
     }
   };
 
+  template<
+  typename InputRange,
+  typename PointMap,
+  typename Line_2>
+  const typename Kernel_traits<Line_2>::Kernel::FT
+  create_line_from_points_2(
+    const InputRange& input_range, const PointMap point_map,
+    const std::vector<std::size_t>& region, Line_2& line) {
+
+    using Traits = typename Kernel_traits<Line_2>::Kernel;
+    using FT = typename Traits::FT;
+
+    using ITraits = CGAL::Exact_predicates_inexact_constructions_kernel;
+    using IConverter = Cartesian_converter<Traits, ITraits>;
+
+    using IFT = typename ITraits::FT;
+    using IPoint_2 = typename ITraits::Point_2;
+    using ILine_2 = typename ITraits::Line_2;
+
+    std::vector<IPoint_2> points;
+    CGAL_precondition(region.size() > 0);
+    points.reserve(region.size());
+    const IConverter iconverter = IConverter();
+
+    for (const std::size_t point_index : region) {
+      CGAL_precondition(point_index < input_range.size());
+      const auto& key = *(input_range.begin() + point_index);
+      const auto& point = get(point_map, key);
+      points.push_back(iconverter(point));
+    }
+    CGAL_postcondition(points.size() == region.size());
+
+    ILine_2 fitted_line;
+    IPoint_2 fitted_centroid;
+    const IFT score = CGAL::linear_least_squares_fitting_2(
+      points.begin(), points.end(),
+      fitted_line, fitted_centroid,
+      CGAL::Dimension_tag<0>(), ITraits(),
+      CGAL::Eigen_diagonalize_traits<IFT, 2>());
+
+    line = Line_2(
+      static_cast<FT>(fitted_line.a()),
+      static_cast<FT>(fitted_line.b()),
+      static_cast<FT>(fitted_line.c()));
+    return static_cast<FT>(score);
+  }
+
   template<
   typename InputRange,
   typename PointMap,
   typename Plane_3>
-  void create_planes_from_points(
-    const InputRange& input_range,
-    const PointMap point_map,
-    std::vector< std::vector<std::size_t> >& regions,
-    std::vector<Plane_3>& planes) {
+  const typename Kernel_traits<Plane_3>::Kernel::FT
+  create_plane_from_points(
+    const InputRange& input_range, const PointMap point_map,
+    const std::vector<std::size_t>& region, Plane_3& plane) {
 
     using Traits = typename Kernel_traits<Plane_3>::Kernel;
     using FT = typename Traits::FT;
@@ -114,98 +160,56 @@ namespace internal {
     using IPoint_3 = typename ITraits::Point_3;
     using IPlane_3 = typename ITraits::Plane_3;
 
+    std::vector<IPoint_3> points;
+    CGAL_precondition(region.size() > 0);
+    points.reserve(region.size());
+    const IConverter iconverter = IConverter();
+
+    for (const std::size_t point_index : region) {
+      CGAL_precondition(point_index < input_range.size());
+      const auto& key = *(input_range.begin() + point_index);
+      const auto& point = get(point_map, key);
+      points.push_back(iconverter(point));
+    }
+    CGAL_postcondition(points.size() == region.size());
+
+    IPlane_3 fitted_plane;
+    IPoint_3 fitted_centroid;
+    const IFT score = CGAL::linear_least_squares_fitting_3(
+      points.begin(), points.end(),
+      fitted_plane, fitted_centroid,
+      CGAL::Dimension_tag<0>(), ITraits(),
+      CGAL::Eigen_diagonalize_traits<IFT, 3>());
+
+    plane = Plane_3(
+      static_cast<FT>(fitted_plane.a()),
+      static_cast<FT>(fitted_plane.b()),
+      static_cast<FT>(fitted_plane.c()),
+      static_cast<FT>(fitted_plane.d()));
+    return static_cast<FT>(score);
+  }
+
+  template<
+  typename InputRange,
+  typename PointMap,
+  typename Plane_3>
+  void create_planes_from_points(
+    const InputRange& input_range, const PointMap point_map,
+    std::vector< std::vector<std::size_t> >& regions,
+    std::vector<Plane_3>& planes) {
+
     planes.clear();
     planes.reserve(regions.size());
 
-    std::vector<IPoint_3> points;
-    const IConverter iconverter = IConverter();
+    Plane_3 plane;
     for (const auto& region : regions) {
-      CGAL_assertion(region.size() > 0);
-
-      points.clear();
-      for (const std::size_t point_index : region) {
-        CGAL_precondition(point_index < input_range.size());
-        const auto& key = *(input_range.begin() + point_index);
-        const auto& point = get(point_map, key);
-        points.push_back(iconverter(point));
-      }
-      CGAL_postcondition(points.size() == region.size());
-
-      IPlane_3 fitted_plane;
-      IPoint_3 fitted_centroid;
-      CGAL::linear_least_squares_fitting_3(
-        points.begin(), points.end(),
-        fitted_plane, fitted_centroid,
-        CGAL::Dimension_tag<0>(), ITraits(),
-        CGAL::Eigen_diagonalize_traits<IFT, 3>());
-
-      const Plane_3 plane = Plane_3(
-        static_cast<FT>(fitted_plane.a()),
-        static_cast<FT>(fitted_plane.b()),
-        static_cast<FT>(fitted_plane.c()),
-        static_cast<FT>(fitted_plane.d()));
+      create_plane_from_points(input_range, point_map, region, plane);
       planes.push_back(plane);
     }
     CGAL_postcondition(planes.size() == regions.size());
   }
 
 } // namespace internal
-
-namespace RG {
-
-  template<typename GeomTraits>
-  class Plane {
-
-  public:
-    using Traits = GeomTraits;
-    using FT = typename Traits::FT;
-    using Point_2 = typename Traits::Point_2;
-    using Point_3 = typename Traits::Point_3;
-    using Plane_3 = typename Traits::Plane_3;
-    using Vector_3 = typename Traits::Vector_3;
-
-    template<
-    typename Input_range,
-    typename Point_map>
-    Plane(
-      const Input_range& input_range,
-      const Point_map point_map,
-      const std::vector<std::size_t>& region,
-      const Plane_3& plane) {
-
-      FT x = FT(0), y = FT(0), z = FT(0);
-      for (const std::size_t idx : region) {
-        const auto& p = get(point_map, *(input_range.begin() + idx));
-        x += p.x();
-        y += p.y();
-        z += p.z();
-      }
-      x /= static_cast<FT>(region.size());
-      y /= static_cast<FT>(region.size());
-      z /= static_cast<FT>(region.size());
-
-      m_centroid = Point_3(x, y, z);
-      m_base1 = plane.base1() / static_cast<FT>(CGAL::sqrt(
-        CGAL::to_double(plane.base1() * plane.base1())));
-      m_base2 = plane.base2() / static_cast<FT>(CGAL::sqrt(
-        CGAL::to_double(plane.base2() * plane.base2())));
-    }
-
-    Point_2 to_2d(const Point_3& query) const {
-      const Vector_3 v(m_centroid, query);
-      return Point_2(v * m_base1, v * m_base2);
-    }
-
-    Point_3 to_3d(const Point_2& query) const {
-      return m_centroid + query.x() * m_base1 + query.y() * m_base2;
-    }
-
-  private:
-    Point_3 m_centroid;
-    Vector_3 m_base1, m_base2;
-  };
-
-} // namespace RG
 } // namespace Shape_detection
 } // namespace CGAL
 
diff --git a/Shape_detection/test/Shape_detection/test_region_growing_on_degenerated_mesh.cpp b/Shape_detection/test/Shape_detection/test_region_growing_on_degenerated_mesh.cpp
index 204dcb84117..a6084e99d50 100644
--- a/Shape_detection/test/Shape_detection/test_region_growing_on_degenerated_mesh.cpp
+++ b/Shape_detection/test/Shape_detection/test_region_growing_on_degenerated_mesh.cpp
@@ -83,8 +83,8 @@ bool test_region_growing_on_degenerated_mesh(int argc, char *argv[]) {
   region_growing.detect(std::back_inserter(regions));
 
   // Test data.
-  assert(regions.size() >= 265 && regions.size() <= 269);
-  if (regions.size() < 265 || regions.size() > 269)
+  assert(regions.size() >= 260 && regions.size() <= 264);
+  if (regions.size() < 260 || regions.size() > 264)
     return false;
 
   for (const auto& region : regions)
@@ -94,8 +94,8 @@ bool test_region_growing_on_degenerated_mesh(int argc, char *argv[]) {
   std::vector<std::size_t> unassigned_faces;
   region_growing.unassigned_items(std::back_inserter(unassigned_faces));
 
-  assert(unassigned_faces.size() >= 501 && unassigned_faces.size() <= 521);
-  if (unassigned_faces.size() < 501 || unassigned_faces.size() > 521)
+  assert(unassigned_faces.size() >= 493 && unassigned_faces.size() <= 513);
+  if (unassigned_faces.size() < 493 || unassigned_faces.size() > 513)
     return false;
 
   return true;
diff --git a/Shape_detection/test/Shape_detection/test_region_growing_on_polygon_mesh.cpp b/Shape_detection/test/Shape_detection/test_region_growing_on_polygon_mesh.cpp
index 21708bd6b84..d012ea4a4a7 100644
--- a/Shape_detection/test/Shape_detection/test_region_growing_on_polygon_mesh.cpp
+++ b/Shape_detection/test/Shape_detection/test_region_growing_on_polygon_mesh.cpp
@@ -82,8 +82,8 @@ bool test_region_growing_on_polygon_mesh(int argc, char *argv[]) {
   region_growing.detect(std::back_inserter(regions));
 
   // Test data.
-  assert(regions.size() >= 324 && regions.size() <= 328);
-  if (regions.size() < 324 || regions.size() > 328)
+  assert(regions.size() >= 325 && regions.size() <= 334);
+  if (regions.size() < 325 || regions.size() > 334)
     return false;
 
   for (const auto& region : regions)
@@ -93,8 +93,8 @@ bool test_region_growing_on_polygon_mesh(int argc, char *argv[]) {
   std::vector<std::size_t> unassigned_faces;
   region_growing.unassigned_items(std::back_inserter(unassigned_faces));
 
-  assert(unassigned_faces.size() >= 904 && unassigned_faces.size() <= 924);
-  if (unassigned_faces.size() < 904 || unassigned_faces.size() > 924)
+  assert(unassigned_faces.size() >= 908 && unassigned_faces.size() <= 928);
+  if (unassigned_faces.size() < 908 || unassigned_faces.size() > 928)
     return false;
 
   return true;
diff --git a/Shape_detection/todo.md b/Shape_detection/todo.md
index 97060ca79c8..353a8cca150 100644
--- a/Shape_detection/todo.md
+++ b/Shape_detection/todo.md
@@ -7,4 +7,4 @@
 - cleanup the code
 - finish new tests
 - update the docs
-- remove sqrt
\ No newline at end of file
+- do not use Plane_3 from RG namespace in the polyhedron demo
\ No newline at end of file