fix bugs + add a robust version of the inside region predicate

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;
 
@@ -947,6 +984,14 @@ namespace CommonKernelFunctors {
       return CGAL::compare(internal::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/include/CGAL/Shape_detection/Region_growing/Polygon_mesh/Plane_face_region.h

@@ -30,13 +30,12 @@ namespace Polygon_mesh {
   /*!
     \ingroup PkgShapeDetectionRGOnMesh
 
-    \brief Region type based on the quality of the least squares plane
-    fit applied to faces of a polygon mesh.
+    \brief Region type based on the plane of the first face selected.
 
-    This class fits a plane, using \ref PkgPrincipalComponentAnalysisDRef "PCA",
-    to chunks of faces in a polygon mesh and controls the quality of this fit.
-    If all quality conditions are satisfied, the chunk is accepted as a valid region,
-    otherwise rejected.
+    This class uses the supporting plane of the first face picked for the region
+    and expand it for all faces with a normal close that 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`
@@ -173,9 +172,8 @@ namespace Polygon_mesh {
     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;
+        !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)
@@ -205,7 +203,7 @@ namespace Polygon_mesh {
     };
 #endif
 
-      if constexpr (use_input_face_normal)
+      if constexpr (!use_input_face_normal)
       {
         for (const Item &i : faces(pmesh))
           put(m_face_normals, i, get_face_normal(i, pmesh));
@@ -280,7 +278,7 @@ namespace Polygon_mesh {
 
     // TODO: we probably want to return the m_seed_face instead
     Primitive primitive() const {
-      return m_plane_of_best_fit;
+      return m_plane;
     }
 
     /*!
@@ -300,6 +298,7 @@ namespace Polygon_mesh {
 
       \pre `query` is a valid const_iterator of `input_range`
     */
+#ifdef FAST // TODO: bench me
     bool is_part_of_region(
       const Item query,
       const Region&) const
@@ -358,6 +357,70 @@ namespace Polygon_mesh {
 
       return true;
     }
+#else
+    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;
+      }
+    }
+#endif
 
     /*!
       \brief implements `RegionType::is_valid_region()`.
@@ -389,7 +452,7 @@ namespace Polygon_mesh {
     bool update(const Region& region) {
 
       CGAL_precondition(region.size() > 0);
-      if (region.size() == 1) { // create new reference plane and normal
+      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
@@ -399,8 +462,8 @@ namespace Polygon_mesh {
         if (face_normal == CGAL::NULL_VECTOR) return false;
 
         CGAL_precondition(face_normal != CGAL::NULL_VECTOR);
-        m_plane_of_best_fit = Plane_3(face_centroid, face_normal);
-        m_normal_of_best_fit = face_normal;
+        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;
@@ -425,8 +488,8 @@ namespace Polygon_mesh {
     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_of_best_fit;
-    Vector_3 m_normal_of_best_fit;
+    Plane_3 m_plane;
+    Vector_3 m_normal;
     face_descriptor m_seed_face;
 
     // Compute centroid of the face.