diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md
index ca749a987028..9768c8694cbc 100644
--- a/Installation/CHANGES.md
+++ b/Installation/CHANGES.md
@@ -72,6 +72,15 @@
 
 -   Add a the non-zero rule, as well as functions to compute the conservative inner and outer hull of similar polygons.
 
+
+### [Shape Detection](https://doc.cgal.org/6.1/Manual/packages.html#PkgShapeDetection) 
+
+- Added the region type [`CGAL::Shape_detection::Polygon_mesh::Plane_face_region`](https://doc.cgal.org/6.1/Shape_detection/class_c_g_a_l_1_1_shape__detection_1_1_polygon__mesh_1_1_plane__face__region.html) that extends the support plane of the seed face without refitting the plane to the region
+- Added the sorting [`CGAL::Shape_detection::Polygon_mesh::Face_area_sorting`](https://doc.cgal.org/6.1/Shape_detection/class_c_g_a_l_1_1_shape__detection_1_1_polygon__mesh_1_1_face__area__sorting.html) that provides a sorting of faces in descending order of area
+- Added the named parameter `face_normal_map` to [`CGAL::Shape_detection::Polygon_mesh::Linear_least_squares_fit_region`](https://doc.cgal.org/6.1/Shape_detection/class_c_g_a_l_1_1_shape__detection_1_1_polygon__mesh_1_1_least__squares__plane__fit__region.html) ,[`CGAL::Shape_detection::Polygon_mesh::Plane_face_region`](https://doc.cgal.org/6.1/Shape_detection/class_c_g_a_l_1_1_shape__detection_1_1_polygon__mesh_1_1_plane__face__region.html) and [`CGAL::Polygon_mesh_processing::region_growing_of_planes_on_faces`](https://doc.cgal.org/6.1/Polygon_mesh_processing/group___pkg_polygon_mesh_processing_ref.html#ga50dcd2f6295f584d2e378b57290ae2af) to allow the use of face normal property maps instead of calculating the face normals.
+- The [`CGAL::Shape_detection::Polygon_mesh::Plane_face_region`](https://doc.cgal.org/6.1/Shape_detection/class_c_g_a_l_1_1_shape__detection_1_1_polygon__mesh_1_1_plane__face__region.html) and [`CGAL::Shape_detection::Polygon_mesh::Face_area_sorting`](https://doc.cgal.org/6.1/Shape_detection/class_c_g_a_l_1_1_shape__detection_1_1_polygon__mesh_1_1_face__area__sorting.html) are now used as the default in [`CGAL::Polygon_mesh_processing::region_growing_of_planes_on_faces`](https://doc.cgal.org/6.1/Polygon_mesh_processing/group___pkg_polygon_mesh_processing_ref.html#ga50dcd2f6295f584d2e378b57290ae2af)
+
+
 ### Triangulations
 -   All triangulations now offer the functions `point(Vertex_handle)` and `point(Simplex, int)`, which enables users to access the geometric position of a vertex and of the i-th vertex of a simplex of a triangulation.
 
diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h
index b71470a62625..7d875d069db6 100644
--- a/Kernel_23/include/CGAL/Kernel/function_objects.h
+++ b/Kernel_23/include/CGAL/Kernel/function_objects.h
@@ -275,6 +275,40 @@ namespace CommonKernelFunctors {
                                 * sq_length(ba)*sq_length(bc));
       }
     }
+
+    Comparison_result
+    operator()(const Point_3& a1, const Point_3& b1, const Point_3& c1,
+               const Point_3& a2, const Point_3& b2, const Point_3& c2,
+               const FT& cosine) const
+    {
+      typename K::Compute_scalar_product_3 scalar_product = K().compute_scalar_product_3_object();
+      typename K::Compute_squared_length_3 sq_length = K().compute_squared_length_3_object();
+      typename K::Construct_normal_3 normal = K().construct_normal_3_object();
+
+      Vector_3 n1 = normal(a1,b1,c1);
+      Vector_3 n2 = normal(a2,b2,c2);
+
+      typename K::FT sc_prod = scalar_product(n1, n2);
+
+      if (sc_prod >= 0)
+      {
+        if (cosine >= 0)
+          return CGAL::compare(CGAL::square(cosine)
+                                * sq_length(n1)*sq_length(n2),
+                               CGAL::square(sc_prod));
+        else
+          return SMALLER;
+      }
+      else
+      {
+        if (cosine >= 0)
+          return LARGER;
+        else
+          return CGAL::compare(CGAL::square(sc_prod),
+                               CGAL::square(cosine)
+                                * sq_length(n1)*sq_length(n2));
+      }
+    }
   };
 
   template <typename K>
@@ -899,6 +933,10 @@ namespace CommonKernelFunctors {
   {
     typedef typename K::Comparison_result  Comparison_result;
     typedef typename K::FT                 FT;
+    typedef typename K::Point_3            Point_3;
+    typedef typename K::Plane_3            Plane_3;
+    typename K::Construct_plane_3 construct_plane;
+  public:
 
   public:
     template <class T1, class T2>
@@ -914,6 +952,14 @@ namespace CommonKernelFunctors {
     {
       return internal::compare_squared_distance(p, q, K(), internal::squared_distance(r, s, K()));
     }
+
+
+    Needs_FT<Comparison_result>
+    operator()(const Point_3& p, const Point_3& q, const Point_3& r, const Point_3& query, const FT& sqd) const
+    {
+      Plane_3 plane = construct_plane(p,q,r);
+      return this->operator()(plane, query, sqd);
+    }
   };
 
  template <typename K>
diff --git a/Shape_detection/doc/Shape_detection/PackageDescription.txt b/Shape_detection/doc/Shape_detection/PackageDescription.txt
index e1efc07b3dd8..acc7d9177a98 100644
--- a/Shape_detection/doc/Shape_detection/PackageDescription.txt
+++ b/Shape_detection/doc/Shape_detection/PackageDescription.txt
@@ -96,7 +96,9 @@ to simplify the definition and creation of the different classes. The naming con
 - `CGAL::Shape_detection::Polygon_mesh::One_ring_neighbor_query<PolygonMesh>`
 - `CGAL::Shape_detection::Polygon_mesh::Polyline_graph<GeomTraits, PolygonMesh, VertexPointMap>`
 - `CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_region<GeomTraits, PolygonMesh, VertexToPointMap>`
+- `CGAL::Shape_detection::Polygon_mesh::Plane_face_region<GeomTraits, PolygonMesh, VertexToPointMap>`
 - `CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_sorting<GeomTraits, PolygonMesh, %NeighborQuery, VertexToPointMap>`
+- `CGAL::Shape_detection::Polygon_mesh::Face_area_sorting<GeomTraits, PolygonMesh, VertexToPointMap>`
 
 
 \defgroup PkgShapeDetectionRGConcepts Concepts
diff --git a/Shape_detection/doc/Shape_detection/Shape_detection.txt b/Shape_detection/doc/Shape_detection/Shape_detection.txt
index d58d86eeaccb..e1b4e3af8903 100644
--- a/Shape_detection/doc/Shape_detection/Shape_detection.txt
+++ b/Shape_detection/doc/Shape_detection/Shape_detection.txt
@@ -261,7 +261,7 @@ their own propagation and seeding conditions (see \ref Shape_detection_RegionGro
 
 \subsubsection Shape_detection_RegionGrowingFramework_examples Examples
 
-This toy example shows how to define a custom models of`NeighborQuery` and
+This toy example shows how to define a custom models of `NeighborQuery` and
 `RegionType` concept, which are used to parameterize the
 `CGAL::Shape_detection::Region_growing`.
 
@@ -483,6 +483,7 @@ In particular, it provides
 - `CGAL::Shape_detection::Polygon_mesh::One_ring_neighbor_query` class that retrieves all edge-adjacent faces of a face, and
 - `CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_region` class that fits a 3D plane to the vertices of
 all mesh faces, which have been added to the region so far, and controls the quality of this fit.
+- `CGAL::Shape_detection::Polygon_mesh::Plane_face_region` is another `RegionType` that uses the 3D plane of the seed face and does not refit the plane to the region.
 
 This component accepts any model of the concept `FaceListGraph` as a polygon mesh. \cgalFigureRef{Region_growing_on_surface_mesh}
 gives an \ref Shape_detection_RegionGrowingMesh_examples "example" of the region growing algorithm
@@ -495,10 +496,11 @@ A surface mesh depicted with one color per detected plane.
 The quality of region growing on a polygon mesh can be improved by slightly increasing the running time. To achieve
 this, one can sort the input faces with respect to some quality criteria. These quality criteria
 can be included in region growing through the seed range (see \ref Shape_detection_RegionGrowingFramework for more details).
-We provide such a quality sorting:
+We provide two sortings:
 
 - `CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_sorting` - the input faces are sorted with respect
 to the quality of the least squares plane fit applied to the neighbors of each face.
+- `CGAL::Shape_detection::Polygon_mesh::Face_area_sorting` - the input faces are sorted with respect to the area. The largest face are chosen as the first seed.
 
 
 \subsubsection Shape_detection_RegionGrowingMesh_parameters Parameters
@@ -520,7 +522,7 @@ as the one explained in Section \ref Shape_detection_RegionGrowingPoints_paramet
 In the example below, we show how to use region growing to detect planes on a polygon mesh that can be either stored as `CGAL::Surface_mesh` or `CGAL::Polyhedron_3`.
 
 We can improve the quality of region growing by providing a different seeding order (analogously to \ref Shape_detection_RegionGrowingPoints "Point Sets") that is why in this example we also sort the input faces using the `CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_sorting`
-and only then detect regions.
+and only then detect regions. The example also shows the use of `CGAL::Shape_detection::Polygon_mesh::Plane_face_region` which expands a region using the support plane of a face without refitting and `CGAL::Shape_detection::Polygon_mesh::Face_area_sorting` which sorts the faces in decreasing area.
 
 \cgalExample{Shape_detection/region_growing_planes_on_polygon_mesh.cpp}
 
@@ -528,7 +530,7 @@ and only then detect regions.
 \subsubsection Shape_detection_RegionGrowingMesh_performance Performance
 
 Since accessing neighbors of a face in a polygon mesh is fast, performance of the region growing on a polygon mesh
-mostly depends on the `RegionType` model's implementation, which is usually fast, too.
+mostly depends on the `RegionType` model's implementation, which is usually fast, too. `CGAL::Shape_detection::Polygon_mesh::Plane_face_region` is faster than `CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_region` as it does not perform a refitting of the plane to the region, whereas the latter may provide a tighter fit of the plane to the region and thus a lower number of larger regions.
 
 
 \subsection Shape_detection_RegionGrowingSegments Segment Set
diff --git a/Shape_detection/examples/Shape_detection/region_growing_planes_on_polygon_mesh.cpp b/Shape_detection/examples/Shape_detection/region_growing_planes_on_polygon_mesh.cpp
index 43d7b5b291ae..b7faf914dd86 100644
--- a/Shape_detection/examples/Shape_detection/region_growing_planes_on_polygon_mesh.cpp
+++ b/Shape_detection/examples/Shape_detection/region_growing_planes_on_polygon_mesh.cpp
@@ -21,58 +21,40 @@ using Polygon_mesh   = CGAL::Surface_mesh<Point_3>;
 #endif
 
 using Neighbor_query = CGAL::Shape_detection::Polygon_mesh::One_ring_neighbor_query<Polygon_mesh>;
-using Region_type    = CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_region<Kernel, Polygon_mesh>;
-using Sorting        = CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_sorting<Kernel, Polygon_mesh, Neighbor_query>;
-using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, Region_type>;
-
-int main(int argc, char *argv[]) {
-
-  // Load data either from a local folder or a user-provided file.
-  const bool is_default_input = argc > 1 ? false : true;
-  const std::string filename = is_default_input ? CGAL::data_file_path("meshes/building.off") : argv[1];
-  std::ifstream in(filename);
-  CGAL::IO::set_ascii_mode(in);
-
-  Polygon_mesh polygon_mesh;
-  if (!CGAL::IO::read_polygon_mesh(filename, polygon_mesh)) {
-    std::cerr << "ERROR: cannot read the input file!" << std::endl;
-    return EXIT_FAILURE;
-  }
-  const auto& face_range = faces(polygon_mesh);
-  std::cout << "* number of input faces: " << face_range.size() << std::endl;
-  assert(!is_default_input || face_range.size() == 32245);
-
-  // Default parameter values for the data file building.off.
-  const FT          max_distance    = FT(1);
-  const FT          max_angle       = FT(45);
-  const std::size_t min_region_size = 5;
+using LS_region_type    = CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_region<Kernel, Polygon_mesh>;
+using LS_sorting        = CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_sorting<Kernel, Polygon_mesh, Neighbor_query>;
+using Face_region_type  = CGAL::Shape_detection::Polygon_mesh::Plane_face_region<Kernel, Polygon_mesh>;
+using Face_area_Sorting = CGAL::Shape_detection::Polygon_mesh::Face_area_sorting<Kernel, Polygon_mesh>;
 
+template<class RegionType, class Sorting_type>
+void detect(Polygon_mesh &mesh, FT max_distance, FT max_angle, std::size_t min_region_size, const std::string &out_filename) {
   // Create instances of the classes Neighbor_query and Region_type.
-  Neighbor_query neighbor_query(polygon_mesh);
+  Neighbor_query neighbor_query(mesh);
 
-  Region_type region_type(
-    polygon_mesh,
+  RegionType region_type(
+    mesh,
     CGAL::parameters::
     maximum_distance(max_distance).
     maximum_angle(max_angle).
     minimum_region_size(min_region_size));
 
   // Sort face indices.
-  Sorting sorting(
-    polygon_mesh, neighbor_query);
+  Sorting_type sorting(
+    mesh, neighbor_query);
   sorting.sort();
 
+  using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, RegionType>;
+
   // Create an instance of the region growing class.
   Region_growing region_growing(
-    face_range, sorting.ordered(), neighbor_query, region_type);
+    faces(mesh), sorting.ordered(), neighbor_query, region_type);
 
   // Run the algorithm.
   std::vector<typename Region_growing::Primitive_and_region> regions;
   region_growing.detect(std::back_inserter(regions));
-  std::cout << "* number of found planes: " << regions.size() << std::endl;
-  assert(!is_default_input || regions.size() == 365);
+  std::cout << regions.size() << " regions found" << std::endl;
 
-  const Region_growing::Region_map& map = region_growing.region_map();
+  const typename Region_growing::Region_map& map = region_growing.region_map();
 
   for (std::size_t i = 0; i < regions.size(); i++)
     for (auto& item : regions[i].second) {
@@ -82,7 +64,7 @@ int main(int argc, char *argv[]) {
     }
 
   std::vector<typename Region_growing::Item> unassigned;
-  region_growing.unassigned_items(face_range, std::back_inserter(unassigned));
+  region_growing.unassigned_items(faces(mesh), std::back_inserter(unassigned));
 
   for (auto& item : unassigned) {
     if (std::size_t(-1) != get(map, item)) {
@@ -91,8 +73,33 @@ int main(int argc, char *argv[]) {
   }
 
   // Save regions to a file.
-  const std::string fullpath = (argc > 2 ? argv[2] : "planes_polygon_mesh.ply");
-  utils::save_polygon_mesh_regions(polygon_mesh, regions, fullpath);
+  utils::save_polygon_mesh_regions(mesh, regions, out_filename);
+}
+
+int main(int argc, char *argv[]) {
+
+  // Load data either from a local folder or a user-provided file.
+  const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/building.off") : argv[1];
+  std::ifstream in(filename);
+  CGAL::IO::set_ascii_mode(in);
+
+  Polygon_mesh polygon_mesh;
+  if (!CGAL::IO::read_polygon_mesh(filename, polygon_mesh)) {
+    std::cerr << "ERROR: cannot read the input file!" << std::endl;
+    return EXIT_FAILURE;
+  }
+  const auto& face_range = faces(polygon_mesh);
+  std::cout << "* number of input faces: " << face_range.size() << std::endl;
+
+  // Default parameter values for the data file building.off.
+  const FT          max_distance    = FT(1);
+  const FT          max_angle       = FT(45);
+  const std::size_t min_region_size = 5;
+
+  std::cout << "Region growing with Least_squares_plane_fit_region: ";
+  detect<LS_region_type, LS_sorting>(polygon_mesh, max_distance, max_angle, min_region_size, "least_squares_planes_polygon_mesh.ply");
+  std::cout << "Region growing with Plane_face_region: ";
+  detect<Face_region_type, Face_area_Sorting>(polygon_mesh, max_distance, max_angle, min_region_size, "face_planes_polygon_mesh.ply");
 
   return EXIT_SUCCESS;
 }
diff --git a/Shape_detection/include/CGAL/Polygon_mesh_processing/region_growing.h b/Shape_detection/include/CGAL/Polygon_mesh_processing/region_growing.h
index d6f78085c448..8c936623e2d0 100644
--- a/Shape_detection/include/CGAL/Polygon_mesh_processing/region_growing.h
+++ b/Shape_detection/include/CGAL/Polygon_mesh_processing/region_growing.h
@@ -156,6 +156,13 @@ class One_ring_neighbor_query_with_constraints
       \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
       \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.}
     \cgalParamNEnd
+      \cgalParamNBegin{face_normal_map}
+        \cgalParamDescription{a property map associating normal vectors to the faces of `pmesh`.}
+        \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<PolygonMesh>::%face_descriptor`
+                       as key type and `GT::Vector_3` as value type, `GT` being the type of the parameter `geom_traits`.}
+        \cgalParamDefault{If this parameter is omitted, face normals will be estimated using crossproducts of vectors created
+                          from consecutive vertices of the face.}
+      \cgalParamNEnd
     \cgalParamNBegin{region_primitive_map}
       \cgalParamDescription{a property map filled by this function and that will contain for each region
                             the plane (or only its orthognonal vector) estimated that approximates it.}
@@ -192,13 +199,13 @@ region_growing_of_planes_on_faces(const PolygonMesh& mesh,
                              Static_boolean_property_map<edge_descriptor, false>());
 
   using Neighbor_query = internal::One_ring_neighbor_query_with_constraints<PolygonMesh, ECM>;
-  using Region_type = RG_PM::Least_squares_plane_fit_region<Traits, PolygonMesh, VPM>;
-  using Sorting = RG_PM::Least_squares_plane_fit_sorting<Traits, PolygonMesh, Neighbor_query, VPM>;
+  using Region_type = RG_PM::Plane_face_region<Traits, PolygonMesh, VPM>;
+  using Sorting = RG_PM::Face_area_sorting<Traits, PolygonMesh, VPM>;
   using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, Region_type, RegionMap>;
 
   Neighbor_query neighbor_query(mesh, ecm);
   Region_type region_type(mesh, np);
-  Sorting sorting(mesh, neighbor_query, np);
+  Sorting sorting(mesh, np);
   sorting.sort();
 
   std::vector<typename Region_growing::Primitive_and_region> regions;
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_circle_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_circle_fit_region.h
index 15c5d39e13c7..8eae00029974 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_circle_fit_region.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_circle_fit_region.h
@@ -220,7 +220,7 @@ namespace Point_set {
     /*!
       \brief implements `RegionType::region_index_map()`.
 
-      This function creates an empty property map that maps iterators on the input range `Item` to std::size_t
+      This function creates an empty property map that maps iterators on the input range `Item` to `std::size_t`.
     */
     Region_index_map region_index_map() {
       return Region_index_map(m_region_map);
@@ -247,8 +247,7 @@ namespace Point_set {
       its normal and the circle radius is within the `maximum_angle`. If both conditions
       are satisfied, it returns `true`, otherwise `false`.
 
-      \param query
-      item of the query point
+      \param query item of the query point
 
       \param region
       inlier items of the region
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_cylinder_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_cylinder_fit_region.h
index 12543ef500f8..b23d32b33b88 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_cylinder_fit_region.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_cylinder_fit_region.h
@@ -219,7 +219,7 @@ namespace Point_set {
     /*!
       \brief implements `RegionType::region_index_map()`.
 
-      This function creates an empty property map that maps iterators on the input range `Item` to std::size_t
+      This function creates an empty property map that maps iterators on the input range `Item` to `std::size_t`.
     */
     Region_index_map region_index_map() {
       return Region_index_map(m_region_map);
@@ -246,8 +246,7 @@ namespace Point_set {
       its normal and the cylinder radius is within the `maximum_angle`. If both conditions
       are satisfied, it returns `true`, otherwise `false`.
 
-      \param query
-      item of the query point
+      \param query item of the query point
 
       \param region
       inlier items of the region
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_line_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_line_fit_region.h
index d61ecad6b1d6..c4f7e6db8d89 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_line_fit_region.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_line_fit_region.h
@@ -214,8 +214,7 @@ namespace Point_set {
       its normal and the line's normal is within the `maximum_angle`. If both conditions
       are satisfied, it returns `true`, otherwise `false`.
 
-      \param query
-      item of the query point
+      \param query item of the query point
 
       The last parameter is not used in this implementation.
 
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_plane_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_plane_fit_region.h
index 912114e693ae..6c113283019f 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_plane_fit_region.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_plane_fit_region.h
@@ -188,7 +188,7 @@ namespace Point_set {
     /*!
       \brief implements `RegionType::region_index_map()`.
 
-      This function creates an empty property map that maps iterators on the input range `Item` to std::size_t.
+      This function creates an empty property map that maps iterators on the input range `Item` to `std::size_t`.
     */
     Region_index_map region_index_map() {
       return Region_index_map(m_region_map);
@@ -215,8 +215,7 @@ namespace Point_set {
       its normal and the plane's normal is within the `maximum_angle`. If both conditions
       are satisfied, it returns `true`, otherwise `false`.
 
-      \param query
-      `Item` of the query point
+      \param query item of the query point
 
       The last parameter is not used in this implementation.
 
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_sphere_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_sphere_fit_region.h
index 55b935c357ed..d00c27de6a5d 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_sphere_fit_region.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Point_set/Least_squares_sphere_fit_region.h
@@ -209,7 +209,7 @@ namespace Point_set {
     /*!
       \brief implements `RegionType::region_index_map()`.
 
-      This function creates an empty property map that maps iterators on the input range `Item` to std::size_t
+      This function creates an empty property map that maps iterators on the input range `Item` to `std::size_t`.
     */
 
     Region_index_map region_index_map() {
@@ -238,8 +238,7 @@ namespace Point_set {
       its normal and the sphere radius is within the `maximum_angle`. If both conditions
       are satisfied, it returns `true`, otherwise `false`.
 
-      \param query
-      item of the query point
+      \param query item of the query point
 
       \param region
       inlier items of the region
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh.h
index 73443a13b8f5..6b571861b82c 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh.h
@@ -29,5 +29,7 @@
 
 #include <CGAL/Shape_detection/Region_growing/Polygon_mesh/Least_squares_plane_fit_region.h>
 #include <CGAL/Shape_detection/Region_growing/Polygon_mesh/Least_squares_plane_fit_sorting.h>
+#include <CGAL/Shape_detection/Region_growing/Polygon_mesh/Plane_face_region.h>
+#include <CGAL/Shape_detection/Region_growing/Polygon_mesh/Face_area_sorting.h>
 
 #endif // CGAL_SHAPE_DETECTION_REGION_GROWING_POLYGON_MESH_H
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh/Face_area_sorting.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh/Face_area_sorting.h
new file mode 100644
index 000000000000..c63f3e9be4a6
--- /dev/null
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh/Face_area_sorting.h
@@ -0,0 +1,245 @@
+// Copyright (c) 2024-2025 GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Sébastien Loriot
+//
+
+
+#ifndef CGAL_SHAPE_DETECTION_REGION_GROWING_POLYGON_MESH_FACE_AREA_SORTING_H
+#define CGAL_SHAPE_DETECTION_REGION_GROWING_POLYGON_MESH_FACE_AREA_SORTING_H
+
+#include <CGAL/license/Shape_detection.h>
+#include <CGAL/Shape_detection/Region_growing/internal/utils.h>
+
+namespace CGAL {
+namespace Shape_detection {
+namespace Polygon_mesh {
+
+  /*!
+    \ingroup PkgShapeDetectionRGOnMesh
+
+    \brief Sorting of polygon mesh faces with respect to their area.
+
+    `Items` of faces in a polygon mesh are sorted in decreasing area.
+
+    \tparam GeomTraits
+    a model of `Kernel`
+
+    \tparam PolygonMesh
+    a model of `FaceListGraph`
+
+    \tparam VertexToPointMap
+    a model of `ReadablePropertyMap` whose key type is the vertex type of a polygon mesh and
+    value type is `Kernel::Point_3`
+  */
+  template<
+  typename GeomTraits,
+  typename PolygonMesh,
+  typename VertexToPointMap = typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::const_type>
+  class Face_area_sorting
+  {
+    using vertex_descriptor = typename boost::graph_traits<PolygonMesh>::vertex_descriptor;
+    using face_descriptor = typename boost::graph_traits<PolygonMesh>::face_descriptor;
+    using halfedge_descriptor = typename boost::graph_traits<PolygonMesh>::halfedge_descriptor;
+    using FT = typename GeomTraits::FT;
+
+  public:
+    /// \name Types
+    /// @{
+
+    /// Item type.
+    using Item = typename boost::graph_traits<PolygonMesh>::face_descriptor;
+
+    /// Seed range.
+    using Seed_range = std::vector<Item>;
+
+    /// @}
+
+    /// \name Initialization
+    /// @{
+
+    /*!
+      \brief initializes all internal data structures.
+
+      \tparam NamedParameters
+      a sequence of \ref bgl_namedparameters "Named Parameters"
+
+      \param pmesh
+      an instance of `PolygonMesh` that represents a polygon mesh
+
+      \param np
+      a sequence of \ref bgl_namedparameters "Named Parameters"
+      among the ones listed below
+
+      \cgalNamedParamsBegin
+        \cgalParamNBegin{vertex_point_map}
+          \cgalParamDescription{an instance of `VertexToPointMap` that maps a polygon mesh
+          vertex to `Kernel::Point_3`}
+          \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`}
+        \cgalParamNEnd
+        \cgalParamNBegin{geom_traits}
+          \cgalParamDescription{an instance of `GeomTraits`}
+          \cgalParamDefault{`GeomTraits()`}
+        \cgalParamNEnd
+      \cgalNamedParamsEnd
+
+      \pre `faces(pmesh).size() > 0`
+    */
+    template<typename CGAL_NP_TEMPLATE_PARAMETERS>
+    Face_area_sorting(
+      const PolygonMesh& pmesh,
+      const CGAL_NP_CLASS& np = parameters::default_values())
+        : m_face_graph(pmesh)
+        , m_vpm(parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
+                                                             get_const_property_map(CGAL::vertex_point, pmesh)))
+        , m_traits(parameters::choose_parameter<GeomTraits>(parameters::get_parameter(np, internal_np::geom_traits)))
+    {
+      CGAL_precondition(faces(pmesh).size() > 0);
+
+      m_ordered.resize(faces(pmesh).size());
+
+      std::size_t index = 0;
+      for (Item item : faces(pmesh))
+        m_ordered[index++] = item;
+      m_scores.resize(m_ordered.size(), 0.);
+    }
+
+    /*!
+      \brief initializes all internal data structures.
+      3 Parameter constructor with dummy parameter provided for compatibility with other sorting types.
+
+      \tparam Dummy
+      Dummy parameter, not used.
+
+      \tparam NamedParameters
+      a sequence of \ref bgl_namedparameters "Named Parameters"
+
+      \param pmesh
+      an instance of `PolygonMesh` that represents a polygon mesh
+
+      \param np
+      a sequence of \ref bgl_namedparameters "Named Parameters"
+      among the ones listed below
+
+      \cgalNamedParamsBegin
+        \cgalParamNBegin{vertex_point_map}
+          \cgalParamDescription{an instance of `VertexToPointMap` that maps a polygon mesh
+          vertex to `Kernel::Point_3`}
+          \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`}
+        \cgalParamNEnd
+        \cgalParamNBegin{geom_traits}
+          \cgalParamDescription{an instance of `GeomTraits`}
+          \cgalParamDefault{`GeomTraits()`}
+        \cgalParamNEnd
+      \cgalNamedParamsEnd
+
+      \pre `faces(pmesh).size() > 0`
+    */
+    template<typename Dummy, typename CGAL_NP_TEMPLATE_PARAMETERS>
+      Face_area_sorting(
+        const PolygonMesh& pmesh,
+        const Dummy&,
+        const CGAL_NP_CLASS& np = parameters::default_values())
+      : m_face_graph(pmesh)
+      , m_vpm(parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
+        get_const_property_map(CGAL::vertex_point, pmesh)))
+      , m_traits(parameters::choose_parameter<GeomTraits>(parameters::get_parameter(np, internal_np::geom_traits)))
+    {
+      CGAL_precondition(faces(pmesh).size() > 0);
+
+      m_ordered.resize(faces(pmesh).size());
+
+      std::size_t index = 0;
+      for (Item item : faces(pmesh))
+        m_ordered[index++] = item;
+      m_scores.resize(m_ordered.size(), 0.);
+    }
+
+    /// @}
+
+    /// \name Sorting
+    /// @{
+
+    /*!
+      \brief sorts `Items` of input faces.
+    */
+    void sort() {
+      compute_scores();
+      CGAL_precondition(m_scores.size() > 0);
+
+      auto cmp = [this](const std::size_t i, const std::size_t j)
+      {
+        CGAL_precondition(i < m_scores.size());
+        CGAL_precondition(j < m_scores.size());
+        return m_scores[i] > m_scores[j];
+      };
+      std::vector<std::size_t> order(m_ordered.size());
+      std::iota(order.begin(), order.end(), 0);
+      std::sort(order.begin(), order.end(), cmp);
+
+      std::vector<Item> tmp(m_ordered.size());
+      for (std::size_t i = 0; i < m_ordered.size(); i++)
+        tmp[i] = m_ordered[order[i]];
+
+      m_ordered.swap(tmp);
+    }
+    /// @}
+
+    /// \name Access
+    /// @{
+
+    /*!
+      \brief returns an instance of `Seed_range` to access the ordered `Items`
+      of input faces.
+    */
+    const Seed_range &ordered() {
+      return m_ordered;
+    }
+    /// @}
+
+  private:
+    const PolygonMesh& m_face_graph;
+    VertexToPointMap m_vpm;
+    const GeomTraits m_traits;
+    Seed_range m_ordered;
+    std::vector<FT> m_scores;
+
+    void compute_scores()
+    {
+      auto squared_area = m_traits.compute_squared_area_3_object();
+      std::size_t idx = 0;
+      for (Item item : m_ordered)
+      {
+        halfedge_descriptor hd = halfedge(item, m_face_graph);
+        std::vector<typename GeomTraits::Point_3> pts;
+
+        for (vertex_descriptor v : vertices_around_face(hd,m_face_graph))
+          pts.push_back( get(m_vpm, v) );
+
+        if (pts.size()==3)
+          m_scores[idx++] = approximate_sqrt(squared_area(pts[0], pts[1], pts[2]));
+        else
+        {
+          std::vector<typename GeomTraits::Triangle_3> triangles;
+          internal::triangulate_face<GeomTraits>(pts, triangles);
+          for (const typename GeomTraits::Triangle_3& tr : triangles)
+            m_scores[idx] += approximate_sqrt(squared_area(tr[0], tr[1], tr[2]));
+          ++idx;
+        }
+      }
+    }
+  };
+
+} // namespace Polygon_mesh
+} // namespace Shape_detection
+} // namespace CGAL
+
+#endif // CGAL_SHAPE_DETECTION_REGION_GROWING_POLYGON_MESH_FACE_AREA_SORTING_H
+
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh/Least_squares_plane_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh/Least_squares_plane_fit_region.h
index fdc710665819..4375cbada0f2 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh/Least_squares_plane_fit_region.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh/Least_squares_plane_fit_region.h
@@ -137,6 +137,13 @@ namespace Polygon_mesh {
           vertex to `Kernel::Point_3`}
           \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`}
         \cgalParamNEnd
+        \cgalParamNBegin{face_normal_map}
+          \cgalParamDescription{a property map associating normal vectors to the faces of `pmesh`.}
+          \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<PolygonMesh>::%face_descriptor`
+                         as key type and `GT::Vector_3` as value type, `GT` being the type of the parameter `geom_traits`.}
+          \cgalParamDefault{If this parameter is omitted, face normals will be estimated using crossproducts of vectors created
+                            from consecutive vertices of the face.}
+        \cgalParamNEnd
         \cgalParamNBegin{geom_traits}
           \cgalParamDescription{an instance of `GeomTraits`}
           \cgalParamDefault{`GeomTraits()`}
@@ -161,10 +168,13 @@ namespace Polygon_mesh {
     m_squared_distance_3(m_traits.compute_squared_distance_3_object()),
     m_scalar_product_3(m_traits.compute_scalar_product_3_object()),
     m_cross_product_3(m_traits.construct_cross_product_vector_3_object()),
-    m_face_normals( get(CGAL::dynamic_face_property_t<Vector_3>(), pmesh) ),
+    m_face_normals(get(CGAL::dynamic_face_property_t<Vector_3>(), pmesh)),
     m_face_triangulations( get(CGAL::dynamic_face_property_t<std::vector<Triangle_3>>(), pmesh) )
     {
 
+      static constexpr bool use_input_face_normal =
+        !parameters::is_default_parameter<CGAL_NP_CLASS, internal_np::face_normal_t>::value;
+
 #ifdef CGAL_SD_RG_USE_PMP
     auto get_face_normal = [this](Item face, const PolygonMesh& pmesh)
     {
@@ -193,8 +203,19 @@ namespace Polygon_mesh {
     };
 #endif
 
+      if constexpr (!use_input_face_normal)
+      {
+        for (const Item &i : faces(pmesh))
+          put(m_face_normals, i, get_face_normal(i, pmesh));
+      }
+      else
+      {
+        auto fnm = parameters::get_parameter(np, internal_np::face_normal);
+        for (const Item &i : faces(pmesh))
+          put(m_face_normals, i, get(fnm, i));
+      }
+
       for (const Item &i : faces(pmesh)) {
-        put(m_face_normals, i, get_face_normal(i, pmesh));
         std::vector<Point_3> pts;
         auto h = halfedge(i, pmesh);
         auto s = h;
@@ -239,7 +260,7 @@ namespace Polygon_mesh {
     /*!
       \brief implements `RegionType::region_index_map()`.
 
-      This function creates an empty property map that maps each face to a std::size_t
+      This function creates an empty property map that maps each face to a `std::size_t`.
     */
     Region_index_map region_index_map() {
       return get(CGAL::dynamic_face_property_t<std::size_t>(), m_face_graph);
@@ -267,8 +288,7 @@ namespace Polygon_mesh {
       its normal and the plane's normal is within the `maximum_angle`. If both conditions
       are satisfied, it returns `true`, otherwise `false`.
 
-      \param query
-      `Item` of the query face
+      \param query item of the query face
 
       The last parameter is not used in this implementation.
 
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh/Plane_face_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh/Plane_face_region.h
new file mode 100644
index 000000000000..76e02e79c291
--- /dev/null
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh/Plane_face_region.h
@@ -0,0 +1,461 @@
+// Copyright (c) 2024-2025 GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Sébastien Loriot
+//
+
+#ifndef CGAL_SHAPE_DETECTION_REGION_GROWING_POLYGON_MESH_PLANE_FACE_REGION_H
+#define CGAL_SHAPE_DETECTION_REGION_GROWING_POLYGON_MESH_PLANE_FACE_REGION_H
+
+#include <CGAL/license/Shape_detection.h>
+
+// Internal includes.
+#include <CGAL/Shape_detection/Region_growing/internal/property_map.h>
+#include <CGAL/Shape_detection/Region_growing/internal/utils.h>
+#ifdef CGAL_SD_RG_USE_PMP
+#include <CGAL/Polygon_mesh_processing/compute_normal.h>
+#endif
+
+namespace CGAL {
+namespace Shape_detection {
+namespace Polygon_mesh {
+
+  /*!
+    \ingroup PkgShapeDetectionRGOnMesh
+
+    \brief Region type based on the plane of the first face selected.
+
+    This class uses the supporting plane of the first face picked for the region
+    and expands it for all faces with a normal close to that of the first face
+    (close being defined by the input angle threshold) and such that vertices are
+    not far from that supporting plane (far being defined by the input distance threshold).
+
+    \tparam GeomTraits
+    a model of `Kernel`
+
+    \tparam PolygonMesh
+    a model of `FaceListGraph`
+
+    \tparam VertexToPointMap
+    a model of `ReadablePropertyMap` whose key type is the vertex type of a polygon mesh and
+    value type is `Kernel::Point_3`
+
+    \cgalModels{RegionType}
+  */
+  template<
+  typename GeomTraits,
+  typename PolygonMesh,
+  typename VertexToPointMap = typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::const_type>
+  class Plane_face_region {
+
+  public:
+    /// \name Types
+    /// @{
+
+    /// \cond SKIP_IN_MANUAL
+    using Face_graph = PolygonMesh;
+    using Vertex_to_point_map = VertexToPointMap;
+
+    using face_descriptor = typename boost::graph_traits<PolygonMesh>::face_descriptor;
+    using halfedge_descriptor = typename boost::graph_traits<PolygonMesh>::halfedge_descriptor;
+    using vertex_descriptor = typename boost::graph_traits<PolygonMesh>::vertex_descriptor;
+    /// \endcond
+
+    /// Number type.
+    typedef typename GeomTraits::FT FT;
+
+    /// Item type.
+    using Item = face_descriptor;
+    using Region = std::vector<Item>;
+
+    /// Primitive
+    using Primitive = typename GeomTraits::Plane_3;
+
+    /// Region map
+    using Region_index_map = typename boost::property_map<Face_graph, CGAL::dynamic_face_property_t<std::size_t> >::const_type;
+
+    /// @}
+
+  private:
+    using Point_3 = typename GeomTraits::Point_3;
+    using Vector_3 = typename GeomTraits::Vector_3;
+    using Plane_3 = typename GeomTraits::Plane_3;
+    using Triangle_3 = typename GeomTraits::Triangle_3;
+
+    using Squared_length_3 = typename GeomTraits::Compute_squared_length_3;
+    using Squared_distance_3 = typename GeomTraits::Compute_squared_distance_3;
+    using Scalar_product_3 = typename GeomTraits::Compute_scalar_product_3;
+    using Cross_product_3 = typename GeomTraits::Construct_cross_product_vector_3;
+
+  public:
+    /// \name Initialization
+    /// @{
+
+    /*!
+      \brief initializes all internal data structures.
+
+      \tparam NamedParameters
+      a sequence of \ref bgl_namedparameters "Named Parameters"
+
+      \param pmesh
+      an instance of `PolygonMesh` that represents a polygon mesh
+
+      \param np
+      a sequence of \ref bgl_namedparameters "Named Parameters"
+      among the ones listed below
+
+      \cgalNamedParamsBegin
+        \cgalParamNBegin{maximum_distance}
+          \cgalParamDescription{the maximum distance from the furthest vertex of a face to a plane}
+          \cgalParamType{`GeomTraits::FT`}
+          \cgalParamDefault{1}
+        \cgalParamNEnd
+        \cgalParamNBegin{maximum_angle}
+          \cgalParamDescription{the maximum angle in degrees between
+          the normal of a face and the normal of a plane}
+          \cgalParamType{`GeomTraits::FT`}
+          \cgalParamDefault{25 degrees}
+        \cgalParamNEnd
+        \cgalParamNBegin{cosine_of_maximum_angle}
+          \cgalParamDescription{the cosine value `cos(maximum_angle * PI / 180)` to be used instead of the parameter `maximum_angle()`}
+          \cgalParamType{`GeomTraits::FT`}
+          \cgalParamDefault{`cos(25 * PI / 180)`}
+        \cgalParamNEnd
+        \cgalParamNBegin{minimum_region_size}
+          \cgalParamDescription{the minimum number of faces a region must have}
+          \cgalParamType{`std::size_t`}
+          \cgalParamDefault{1}
+        \cgalParamNEnd
+        \cgalParamNBegin{vertex_point_map}
+          \cgalParamDescription{an instance of `VertexToPointMap` that maps a polygon mesh
+          vertex to `Kernel::Point_3`}
+          \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`}
+        \cgalParamNEnd
+        \cgalParamNBegin{face_normal_map}
+          \cgalParamDescription{a property map associating normal vectors to the faces of `pmesh`.}
+          \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<PolygonMesh>::%face_descriptor`
+                         as key type and `GT::Vector_3` as value type, `GT` being the type of the parameter `geom_traits`.}
+          \cgalParamDefault{If this parameter is omitted, face normals will be estimated using crossproducts of vectors created
+                            from consecutive vertices of the face.}
+        \cgalParamNEnd
+        \cgalParamNBegin{geom_traits}
+          \cgalParamDescription{an instance of `GeomTraits`}
+          \cgalParamDefault{`GeomTraits()`}
+        \cgalParamNEnd
+      \cgalNamedParamsEnd
+
+      \pre `faces(tmesh).size() > 0`
+      \pre `maximum_distance >= 0`
+      \pre `maximum_angle >= 0 && maximum_angle <= 90`
+      \pre `cosine_of_maximum_angle >= 0 && cosine_of_maximum_angle <= 1`
+      \pre `minimum_region_size > 0`
+    */
+    template<typename CGAL_NP_TEMPLATE_PARAMETERS>
+    Plane_face_region(
+      const PolygonMesh& pmesh,
+      const CGAL_NP_CLASS& np = parameters::default_values()) :
+    m_face_graph(pmesh),
+    m_vpm(parameters::choose_parameter(parameters::get_parameter(
+      np, internal_np::vertex_point), get_const_property_map(CGAL::vertex_point, pmesh))),
+    m_traits(parameters::choose_parameter<GeomTraits>(parameters::get_parameter(np, internal_np::geom_traits))),
+    m_squared_length_3(m_traits.compute_squared_length_3_object()),
+    m_squared_distance_3(m_traits.compute_squared_distance_3_object()),
+    m_scalar_product_3(m_traits.compute_scalar_product_3_object()),
+    m_cross_product_3(m_traits.construct_cross_product_vector_3_object()),
+    m_face_normals(get(CGAL::dynamic_face_property_t<Vector_3>(), pmesh)),
+    m_face_triangulations( get(CGAL::dynamic_face_property_t<std::vector<Triangle_3>>(), pmesh) )
+    {
+      static constexpr bool use_input_face_normal =
+        !parameters::is_default_parameter<CGAL_NP_CLASS, internal_np::face_normal_t>::value;
+
+#ifdef CGAL_SD_RG_USE_PMP
+    auto get_face_normal = [this](Item face, const PolygonMesh& pmesh)
+    {
+      return Polygon_mesh_processing::compute_face_normal(face, pmesh, parameters::vertex_point_map(m_vpm));
+    };
+#else
+    auto get_face_normal = [this](Item face, const PolygonMesh& pmesh) -> Vector_3
+    {
+      const auto hedge = halfedge(face, pmesh);
+      const auto vertices = vertices_around_face(hedge, pmesh);
+      CGAL_precondition(vertices.size() >= 3);
+
+      auto vertex = vertices.begin();
+      const Point_3& p1 = get(m_vpm, *vertex); ++vertex;
+      const Point_3& p2 = get(m_vpm, *vertex); ++vertex;
+      Point_3 p3 = get(m_vpm, *vertex);
+      while(collinear(p1, p2, p3))
+      {
+        if (++vertex == vertices.end()) return NULL_VECTOR;
+        p3 = get(m_vpm, *vertex);
+      }
+
+      const Vector_3 u = p2 - p1;
+      const Vector_3 v = p3 - p1;
+      return m_cross_product_3(u, v);
+    };
+#endif
+
+      if constexpr (!use_input_face_normal)
+      {
+        for (const Item &i : faces(pmesh))
+          put(m_face_normals, i, get_face_normal(i, pmesh));
+      }
+      else
+      {
+        auto fnm = parameters::get_parameter(np, internal_np::face_normal);
+        for (const Item &i : faces(pmesh))
+          put(m_face_normals, i, get(fnm, i));
+      }
+
+      for (const Item &i : faces(pmesh)) {
+        std::vector<Point_3> pts;
+        auto h = halfedge(i, pmesh);
+        auto s = h;
+
+        do {
+          pts.push_back(get(m_vpm, target(h, pmesh)));
+          h = next(h, pmesh);
+        } while (h != s);
+
+        std::vector<Triangle_3> face_triangulation;
+        internal::triangulate_face<GeomTraits>(pts, face_triangulation);
+        put(m_face_triangulations, i, face_triangulation);
+      }
+
+      CGAL_precondition(faces(m_face_graph).size() > 0);
+      const FT max_distance = parameters::choose_parameter(
+        parameters::get_parameter(np, internal_np::maximum_distance), FT(1));
+      CGAL_precondition(max_distance >= FT(0));
+      m_distance_threshold = max_distance;
+
+      const FT max_angle = parameters::choose_parameter(
+        parameters::get_parameter(np, internal_np::maximum_angle), FT(25));
+      CGAL_precondition(max_angle >= FT(0) && max_angle <= FT(90));
+
+      m_min_region_size = parameters::choose_parameter(
+        parameters::get_parameter(np, internal_np::minimum_region_size), 1);
+      CGAL_precondition(m_min_region_size > 0);
+
+      const FT default_cos_value = static_cast<FT>(std::cos(CGAL::to_double(
+        (max_angle * static_cast<FT>(CGAL_PI)) / FT(180))));
+      const FT cos_value = parameters::choose_parameter(
+        parameters::get_parameter(np, internal_np::cosine_of_maximum_angle), default_cos_value);
+      CGAL_precondition(cos_value >= FT(0) && cos_value <= FT(1));
+      m_cos_value_threshold = cos_value;
+    }
+
+    /// @}
+
+    /// \name Access
+    /// @{
+
+    /*!
+      \brief implements `RegionType::region_index_map()`.
+
+      This function creates an empty property map that maps each face to a `std::size_t`.
+    */
+    Region_index_map region_index_map() {
+      return get(CGAL::dynamic_face_property_t<std::size_t>(), m_face_graph);
+    }
+
+    /*!
+      \brief implements `RegionType::primitive()`.
+
+      This function provides the last primitive that has been fitted with the region.
+
+      \return Primitive parameters that fits the region.
+
+      \pre `successful fitted primitive via successful call of update(region) with a sufficient large region`
+    */
+
+    // TODO: we probably want to return the m_seed_face instead
+    Primitive primitive() const {
+      return m_plane;
+    }
+
+    /*!
+      \brief implements `RegionType::is_part_of_region()`.
+
+      This function controls if the face `query` is within
+      the `maximum_distance` from the corresponding plane and if the angle between
+      its normal and the plane's normal is within the `maximum_angle`. If both conditions
+      are satisfied, it returns `true`, otherwise `false`.
+
+      \param query item of the query face
+
+      The last parameter is not used in this implementation.
+
+      \return Boolean `true` or `false`
+
+      \pre `query` is a valid const_iterator of `input_range`
+    */
+    bool is_part_of_region(
+      const Item query,
+      const Region&) const
+    {
+      halfedge_descriptor h = halfedge(m_seed_face, m_face_graph);
+
+      //TODO: store me!
+      const typename GeomTraits::Point_3& p=get(m_vpm,source(h, m_face_graph));
+      const typename GeomTraits::Point_3& q=get(m_vpm,target(h, m_face_graph));
+      typename GeomTraits::Point_3 r;
+      //TODO: add safety checks for degenerate faces
+      halfedge_descriptor guard = prev(h, m_face_graph);
+
+      do{
+        h=next(h, m_face_graph);
+        if (h == guard) return false;
+        r=get(m_vpm,target(h, m_face_graph));
+      }
+      while(collinear(p,q,r));
+
+
+      if (m_cos_value_threshold==1 || m_distance_threshold == 0)
+      {
+        h = halfedge(query, m_face_graph);
+        for (vertex_descriptor v : vertices_around_face(h, m_face_graph))
+        {
+          if (!coplanar(p,q,r, get(m_vpm, v)))
+            return false;
+        }
+        return true;
+      }
+      else
+      {
+        // test on distance of points to the plane of the seed face
+        const FT squared_distance_threshold = m_distance_threshold * m_distance_threshold;
+        h = halfedge(query, m_face_graph);
+        for (vertex_descriptor v : vertices_around_face(h, m_face_graph))
+        {
+          //TODO: that's a bit dummy that we retest points that are already in the region...
+          //      not sure caching in a vpm does worth it (need reset for each region)
+          if (typename GeomTraits::Compare_squared_distance_3()(p,q,r,get(m_vpm, v), squared_distance_threshold) != SMALLER)
+            return false;
+        }
+
+        const typename GeomTraits::Point_3& p2=get(m_vpm,source(h, m_face_graph));
+        const typename GeomTraits::Point_3& q2=get(m_vpm,target(h, m_face_graph));
+        typename GeomTraits::Point_3 r2;
+
+        halfedge_descriptor guard = prev(h, m_face_graph);
+        do{
+          h=next(h, m_face_graph);
+          if (h == guard) return true;
+          r2=get(m_vpm,target(h, m_face_graph));
+        }
+        while(collinear(p2,q2,r2));
+
+        // test on the normal of the query face to the normal of the seed face
+        return typename GeomTraits::Compare_angle_3()(p,q,r,
+                                                      p2,q2,r2,
+                                                      m_cos_value_threshold) == SMALLER;
+      }
+    }
+
+    /*!
+      \brief implements `RegionType::is_valid_region()`.
+
+      This function controls if the `region` contains at least `minimum_region_size` faces.
+
+      \param region
+      Faces of the region represented as `Items`.
+
+      \return Boolean `true` or `false`
+    */
+    inline bool is_valid_region(const Region& region) const {
+      return (region.size() >= m_min_region_size);
+    }
+
+    /*!
+      \brief implements `RegionType::update()`.
+
+      This function fits the least squares plane to all vertices of the faces
+      from the `region`.
+
+      \param region
+      Faces of the region represented as `Items`.
+
+      \return Boolean `true` if the plane fitting succeeded and `false` otherwise
+
+      \pre `region.size() > 0`
+    */
+    bool update(const Region& region) {
+
+      CGAL_precondition(region.size() > 0);
+      if (region.size() == 1) { // init reference plane and normal
+        m_seed_face = region[0];
+
+        // The best fit plane will be a plane through this face centroid with
+        // its normal being the face's normal.
+        const Point_3 face_centroid = get_face_centroid(m_seed_face);
+        const Vector_3 face_normal = get(m_face_normals, m_seed_face);
+        if (face_normal == CGAL::NULL_VECTOR) return false;
+
+        CGAL_precondition(face_normal != CGAL::NULL_VECTOR);
+        m_plane = Plane_3(face_centroid, face_normal);
+        m_normal = face_normal;
+      }
+      //TODO: shall we try to find a better seed face in the region?
+      return true;
+    }
+
+    /// @}
+
+  private:
+    const Face_graph& m_face_graph;
+    const Vertex_to_point_map m_vpm;
+    const GeomTraits m_traits;
+
+    FT m_distance_threshold;
+    FT m_cos_value_threshold;
+    std::size_t m_min_region_size;
+
+    const Squared_length_3 m_squared_length_3;
+    const Squared_distance_3 m_squared_distance_3;
+    const Scalar_product_3 m_scalar_product_3;
+    const Cross_product_3 m_cross_product_3;
+
+    typename boost::property_map<Face_graph, CGAL::dynamic_face_property_t<Vector_3> >::const_type m_face_normals;
+    typename boost::property_map<Face_graph, CGAL::dynamic_face_property_t<std::vector<Triangle_3>> >::const_type m_face_triangulations;
+
+    Plane_3 m_plane;
+    Vector_3 m_normal;
+    face_descriptor m_seed_face;
+
+    // Compute centroid of the face.
+    template<typename Face>
+    Point_3 get_face_centroid(const Face& face) const {
+
+      const auto hedge = halfedge(face, m_face_graph);
+      const auto vertices = vertices_around_face(hedge, m_face_graph);
+      CGAL_precondition(vertices.size() > 0);
+
+      FT sum = FT(0), x = FT(0), y = FT(0), z = FT(0);
+      for (const auto vertex : vertices) {
+        const Point_3& point = get(m_vpm, vertex);
+        x += point.x();
+        y += point.y();
+        z += point.z();
+        sum += FT(1);
+      }
+      CGAL_precondition(sum > FT(0));
+      x /= sum;
+      y /= sum;
+      z /= sum;
+      return Point_3(x, y, z);
+    }
+  };
+
+} // namespace Polygon_mesh
+} // namespace Shape_detection
+} // namespace CGAL
+
+#endif // CGAL_SHAPE_DETECTION_REGION_GROWING_POLYGON_MESH_PLANE_FACE_REGION_H
diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Segment_set/Least_squares_line_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Segment_set/Least_squares_line_fit_region.h
index 702cbaa0fb95..552ff6dd72fa 100644
--- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Segment_set/Least_squares_line_fit_region.h
+++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Segment_set/Least_squares_line_fit_region.h
@@ -193,7 +193,7 @@ namespace Segment_set {
     /*!
       \brief implements `RegionType::region_index_map()`.
 
-      This function creates an empty property map that maps iterators on the input range `Item` to std::size_t.
+      This function creates an empty property map that maps iterators on the input range `Item` to `std::size_t`.
     */
     Region_index_map region_index_map() {
       return Region_index_map(m_region_map);
@@ -220,8 +220,7 @@ namespace Segment_set {
       direction of this segment and the line's direction is within the `maximum_angle`.
       If both conditions are satisfied, it returns `true`, otherwise `false`.
 
-      \param query
-      `Item` of the query segment
+      \param query item of the query segment
 
       The last parameter is not used in this implementation.