Shape_detection: Add a way to pass input normal to faces

Installation/CHANGES.md

@@ -46,6 +46,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.
 

Kernel_23/include/CGAL/Kernel/function_objects.h

@@ -284,6 +284,40 @@ namespace CommonKernelFunctors {
                                 * sq_length(ba)*sq_length(bc));
       }
     }
+
+    result_type
+    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>
@@ -930,6 +964,9 @@ namespace CommonKernelFunctors {
   class Compare_squared_distance_3
   {
     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:
     typedef typename K::Comparison_result  result_type;
 
@@ -946,6 +983,14 @@ namespace CommonKernelFunctors {
     {
       return internal::compare_squared_distance(p, q, K(), internal::squared_distance(r, s, K()));
     }
+
+
+    Needs_FT<result_type>
+    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>

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

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. Figure \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
@@ -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

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;

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

Shape_detection/include/CGAL/Shape_detection/Region_growing/Polygon_mesh/Face_area_sorting.h

@@ -0,0 +1,193 @@
+// 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 descending order of 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.);
+    }
+
+    /// @}
+
+    /// \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]));
+        }
+      }
+    }
+  };
+
+} // namespace Polygon_mesh
+} // namespace Shape_detection
+} // namespace CGAL
+
+#endif // CGAL_SHAPE_DETECTION_REGION_GROWING_POLYGON_MESH_FACE_AREA_SORTING_H
+