updated docs, region growing now properly handles the first time call for all regions

Shape_detection/doc/Shape_detection/Concepts/NeighborQuery.h

@@ -8,7 +8,7 @@ to access neighbors of an item.
 \cgalHasModel
 - `CGAL::Shape_detection::Point_set::K_neighbor_query`
 - `CGAL::Shape_detection::Point_set::Sphere_neighbor_query`
-- `CGAL::Shape_detection::Segment_set::Polyline_graph_neighbor_query`
+- `CGAL::Shape_detection::Polygon_mesh::Polyline_graph`
 - `CGAL::Shape_detection::Polygon_mesh::One_ring_neighbor_query`
 - `CGAL::Shape_detection::Polyline::One_ring_neighbor_query`
 */

Shape_detection/doc/Shape_detection/Concepts/RegionType.h

@@ -56,8 +56,12 @@ public:
     condition `indices.size() == 1`. This function is also called periodically
     when enlarging the region. This case can be identified by checking the
     condition `indices.size() > 1`.
+
+    This function also returns a Boolean at the first call when a new region
+    with one seed item is being created. When it is `true`, the new region is
+    further propagated, otherwise, it is rejected.
   */
-  void update(
+  bool update(
     const std::vector<std::size_t>& indices) {
 
   }

Shape_detection/doc/Shape_detection/Shape_detection.txt

@@ -237,8 +237,8 @@ Shapes are detected by growing regions from seed items, where each region is cre
 Together with the generic algorithm's implementation `CGAL::Shape_detection::Region_growing`, five particular instances of this algorithm are provided:
 
 - Line detection in a \ref Shape_detection_RegionGrowingPoints "2D point set";
-- Line detection in a \ref Shape_detection_RegionGrowingPoints "2D/3D segment set";
+- Line detection in a \ref Shape_detection_RegionGrowingSegments "2D/3D segment set";
-- Line detection on a \ref Shape_detection_RegionGrowingPoints "2D/3D polyline";
+- Line detection on a \ref Shape_detection_RegionGrowingPolyline "2D/3D polyline";
 - Plane detection in a \ref Shape_detection_RegionGrowingPoints "3D point set";
 - Plane detection on a \ref Shape_detection_RegionGrowingMesh "polygon mesh".
 
@@ -279,7 +279,7 @@ all necessary region requirements and can be added to a region. It is called per
 - `RegionType::is_valid_region()` This function checks if a region satisfies
 all necessary region requirements. It is called per region.
 
-- `RegionType::update()` This utility function  enables to update any information,
+- `RegionType::update()` This function enables to update any information,
 which is maintained with the region.
 
 

Shape_detection/examples/Shape_detection/region_growing_on_polyline.cpp

@@ -41,15 +41,13 @@ int main(int argc, char *argv[]) {
       return EXIT_FAILURE;
     }
   }
-  polyline_3.pop_back();
   in.close();
   std::cout << "* number of input vertices: " << polyline_3.size() << std::endl;
-  assert(polyline_3.size() == 248);
+  assert(polyline_3.size() == 249);
 
   // Default parameter values for the data file polyline_3.polylines.txt.
-  const FT          max_distance_to_line = FT(45) / FT(10);
+  const FT max_distance_to_line = FT(45) / FT(10);
-  const FT          max_accepted_angle   = FT(45);
+  const FT max_accepted_angle   = FT(45);
-  const std::size_t min_region_size      = 5;
 
   // Create instances of the classes Neighbor_query and Region_type.
   Neighbor_query neighbor_query(polyline_3);
@@ -58,8 +56,7 @@ int main(int argc, char *argv[]) {
     polyline_3,
     CGAL::parameters::
     distance_threshold(max_distance_to_line).
-    angle_deg_threshold(max_accepted_angle).
+    angle_deg_threshold(max_accepted_angle));
-    min_region_size(min_region_size));
 
   // Create an instance of the region growing class.
   Region_growing region_growing(
@@ -69,7 +66,7 @@ int main(int argc, char *argv[]) {
   std::vector< std::vector<std::size_t> > regions;
   region_growing.detect(std::back_inserter(regions));
   std::cout << "* number of found 3D regions: " << regions.size() << std::endl;
-  assert(regions.size() == 10);
+  assert(regions.size() == 12);
 
   // Save 3D regions to a file.
   std::string fullpath = (argc > 2 ? argv[2] : "regions_polyline_3.ply");
@@ -98,10 +95,9 @@ int main(int argc, char *argv[]) {
     polyline_2, std::back_inserter(regions),
     CGAL::parameters::
     distance_threshold(max_distance_to_line).
-    angle_deg_threshold(max_accepted_angle).
+    angle_deg_threshold(max_accepted_angle));
-    min_region_size(min_region_size));
   std::cout << "* number of found 2D regions: " << regions.size() << std::endl;
-  assert(regions.size() == 6);
+  assert(regions.size() == 5);
 
   // Save 2D regions to a file.
   fullpath = (argc > 2 ? argv[2] : "regions_polyline_2.ply");

Shape_detection/examples/Shape_detection/region_growing_with_custom_classes.cpp

@@ -71,8 +71,9 @@ namespace Custom {
       return m_is_valid;
     }
 
-    void update(const std::vector<std::size_t>&) {
+    bool update(const std::vector<std::size_t>&) {
       m_is_valid = true;
+      return m_is_valid;
     }
   };
 

Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing.h

@@ -147,13 +147,14 @@ namespace Shape_detection {
 
         // Try to grow a new region from the index of the seed item.
         if (!m_visited[seed_index]) {
-          propagate(seed_index, region);
+          const bool is_success = propagate(seed_index, region);
 
           // Check global conditions.
-          if (!m_region_type.is_valid_region(region))
+          if (!is_success || !m_region_type.is_valid_region(region)) {
             revert(region);
-          else
+          } else {
             *(regions++) = region;
+          }
         }
       }
       return regions;
@@ -217,7 +218,7 @@ namespace Shape_detection {
     const Seed_map m_seed_map;
     Visited_items m_visited;
 
-    void propagate(const std::size_t seed_index, Indices& region) {
+    bool propagate(const std::size_t seed_index, Indices& region) {
       region.clear();
 
       // Use two queues, while running on this queue, push to the other queue;
@@ -232,7 +233,8 @@ namespace Shape_detection {
       region.push_back(seed_index);
 
       // Update internal properties of the region.
-      m_region_type.update(region);
+      const bool is_well_created = m_region_type.update(region);
+      if (!is_well_created) return false;
 
       Indices neighbors;
       while (
@@ -274,6 +276,7 @@ namespace Shape_detection {
           depth_index = !depth_index;
         }
       }
+      return true;
     }
 
     void revert(const Indices& region) {

Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_line_fit_region.h

@@ -208,7 +208,12 @@ namespace Point_set {
       const auto& key = *(m_input_range.begin() + query_index);
       const Point_2& query_point = get(m_point_map, key);
       const Vector_2& query_normal = get(m_normal_map, key);
-      CGAL_precondition(query_normal != Vector_2());
+
+      const FT a = CGAL::abs(m_line_of_best_fit.a());
+      const FT b = CGAL::abs(m_line_of_best_fit.b());
+      const FT c = CGAL::abs(m_line_of_best_fit.c());
+      if (a == FT(0) && b == FT(0) && c == FT(0))
+        return false;
 
       const FT squared_distance_to_fitted_line =
         m_squared_distance_2(query_point, m_line_of_best_fit);
@@ -251,9 +256,11 @@ namespace Point_set {
       \param region
       indices of points included in the region
 
+      \return Boolean `true` if the line fitting succeeded and `false` otherwise
+
       \pre `region.size() > 0`
     */
-    void update(const std::vector<std::size_t>& region) {
+    bool update(const std::vector<std::size_t>& region) {
 
       CGAL_precondition(region.size() > 0);
       if (region.size() == 1) { // create new reference line and normal
@@ -265,8 +272,9 @@ namespace Point_set {
         const auto& key = *(m_input_range.begin() + point_index);
         const Point_2& point = get(m_point_map, key);
         const Vector_2& normal = get(m_normal_map, key);
-        CGAL_precondition(normal != Vector_2());
+        if (normal == CGAL::NULL_VECTOR) return false;
 
+        CGAL_precondition(normal != CGAL::NULL_VECTOR);
         m_line_of_best_fit = Line_2(point, normal).perpendicular(point);
         m_normal_of_best_fit = m_line_of_best_fit.perpendicular(
           m_line_of_best_fit.point(0)).to_vector();
@@ -276,6 +284,7 @@ namespace Point_set {
         std::tie(m_line_of_best_fit, m_normal_of_best_fit) =
           get_line_and_normal(region);
       }
+      return true;
     }
 
     /// @}
@@ -300,7 +309,6 @@ namespace Point_set {
         CGAL_precondition(normal_index < m_input_range.size());
         const auto& key = *(m_input_range.begin() + normal_index);
         const Vector_2& normal = get(m_normal_map, key);
-        CGAL_precondition(normal != Vector_2());
         const bool agrees =
           m_scalar_product_2(normal, unoriented_normal_of_best_fit) > FT(0);
         votes_to_keep_normal += (agrees ? 1 : -1);

Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_plane_fit_region.h

@@ -208,7 +208,6 @@ namespace Point_set {
       const auto& key = *(m_input_range.begin() + query_index);
       const Point_3& query_point = get(m_point_map, key);
       const Vector_3& query_normal = get(m_normal_map, key);
-      CGAL_precondition(query_normal != Vector_3());
 
       const FT a = CGAL::abs(m_plane_of_best_fit.a());
       const FT b = CGAL::abs(m_plane_of_best_fit.b());
@@ -258,9 +257,11 @@ namespace Point_set {
       \param region
       indices of points included in the region
 
+      \return Boolean `true` if the plane fitting succeeded and `false` otherwise
+
       \pre `region.size() > 0`
     */
-    void update(const std::vector<std::size_t>& region) {
+    bool update(const std::vector<std::size_t>& region) {
 
       CGAL_precondition(region.size() > 0);
       if (region.size() == 1) { // create new reference plane and normal
@@ -272,17 +273,19 @@ namespace Point_set {
         const auto& key = *(m_input_range.begin() + point_index);
         const Point_3& point = get(m_point_map, key);
         const Vector_3& normal = get(m_normal_map, key);
-        CGAL_precondition(normal != Vector_3());
+        if (normal == CGAL::NULL_VECTOR) return false;
 
+        CGAL_precondition(normal != CGAL::NULL_VECTOR);
         m_plane_of_best_fit = Plane_3(point, normal);
         m_normal_of_best_fit = m_plane_of_best_fit.orthogonal_vector();
 
       } else { // update reference plane and normal
-        if (region.size() <= 3) return;
+        if (region.size() < 3) return false;
         CGAL_precondition(region.size() >= 3);
         std::tie(m_plane_of_best_fit, m_normal_of_best_fit) =
           get_plane_and_normal(region);
       }
+      return true;
     }
 
     /// @}
@@ -306,7 +309,6 @@ namespace Point_set {
         CGAL_precondition(normal_index < m_input_range.size());
         const auto& key = *(m_input_range.begin() + normal_index);
         const Vector_3& normal = get(m_normal_map, key);
-        CGAL_precondition(normal != Vector_3());
         const bool agrees =
           m_scalar_product_3(normal, unoriented_normal_of_best_fit) > FT(0);
         votes_to_keep_normal += (agrees ? 1 : -1);

Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_polygon_mesh/Least_squares_plane_fit_region.h

@@ -245,9 +245,11 @@ namespace Polygon_mesh {
       \param region
       indices of faces included in the region
 
+      \return Boolean `true` if the plane fitting succeeded and `false` otherwise
+
       \pre `region.size() > 0`
     */
-    void update(const std::vector<std::size_t>& region) {
+    bool update(const std::vector<std::size_t>& region) {
 
       CGAL_precondition(region.size() > 0);
       if (region.size() == 1) { // create new reference plane and normal
@@ -259,7 +261,9 @@ namespace Polygon_mesh {
         const auto face = *(m_face_range.begin() + face_index);
         const Point_3 face_centroid = get_face_centroid(face);
         const Vector_3 face_normal = get_face_normal(face);
+        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 = m_plane_of_best_fit.orthogonal_vector();
 
@@ -268,6 +272,7 @@ namespace Polygon_mesh {
         std::tie(m_plane_of_best_fit, m_normal_of_best_fit) =
           get_plane_and_normal(region);
       }
+      return true;
     }
 
     /// @}
@@ -373,7 +378,6 @@ namespace Polygon_mesh {
       const Vector_3 u = point2 - point1;
       const Vector_3 v = point3 - point1;
       const Vector_3 face_normal = m_cross_product_3(u, v);
-      CGAL_postcondition(face_normal != Vector_3());
       return face_normal;
     }
 
@@ -399,7 +403,7 @@ namespace Polygon_mesh {
         const FT squared_distance = m_squared_distance_3(point, m_plane_of_best_fit);
         max_squared_distance = (CGAL::max)(squared_distance, max_squared_distance);
       }
-      CGAL_postcondition(max_squared_distance >= FT(0));
+      CGAL_precondition(max_squared_distance >= FT(0));
       return max_squared_distance;
     }
   };

Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_polygon_mesh/Least_squares_plane_fit_sorting.h

@@ -72,7 +72,7 @@ namespace Polygon_mesh {
     #ifdef DOXYGEN_RUNNING
       /*!
         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.
+        This map provides an access to the ordered indices of input faces.
       */
       typedef unspecified_type Seed_map;
     #endif
@@ -129,7 +129,7 @@ namespace Polygon_mesh {
     /// @{
 
     /*!
-      \brief sorts indices of polygon mesh faces.
+      \brief sorts indices of input faces.
     */
     void sort() {
 
@@ -146,7 +146,7 @@ namespace Polygon_mesh {
 
     /*!
       \brief returns an instance of `Seed_map` to access the ordered indices
-      of polygon mesh faces.
+      of input faces.
     */
     Seed_map seed_map() {
       return Seed_map(m_order);

Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_polyline/Least_squares_line_fit_region.h

@@ -210,21 +210,22 @@ namespace Polyline {
       const Point& input_point = get(m_point_map, key1);
       const Point& query_point = get(m_point_map, key2);
 
-      if (region.size() == 1) { // update new reference line and direction
+      // Update new reference line and direction.
+      if (m_direction_of_best_fit == CGAL::NULL_VECTOR) {
+        if (input_point == query_point) return true;
         CGAL_precondition(input_point != query_point);
         m_line_of_best_fit = Line(input_point, query_point);
         m_direction_of_best_fit = m_line_of_best_fit.to_vector();
         return true;
       }
+      CGAL_precondition(m_direction_of_best_fit != CGAL::NULL_VECTOR);
 
-      CGAL_precondition(region.size() >= 2);
+      // Add equal points to the previously defined region.
       if (input_point == query_point) return true;
       CGAL_precondition(input_point != query_point);
       const Vector query_direction(input_point, query_point);
 
-      CGAL_precondition(m_line_of_best_fit != Line());
+      // Check real conditions.
-      CGAL_precondition(m_direction_of_best_fit != Vector());
-
       const FT squared_distance_to_fitted_line =
         m_squared_distance(query_point, m_line_of_best_fit);
       const FT squared_distance_threshold =
@@ -255,6 +256,8 @@ namespace Polyline {
       \return Boolean `true` or `false`
     */
     inline bool is_valid_region(const std::vector<std::size_t>& region) const {
+      if (m_direction_of_best_fit == CGAL::NULL_VECTOR)
+        return false; // all points are equal
       return (region.size() >= m_min_region_size);
     }
 
@@ -266,19 +269,23 @@ namespace Polyline {
       \param region
       indices of vertices included in the region
 
+      \return Boolean `true` if the line fitting succeeded and `false` otherwise
+
       \pre `region.size() > 0`
     */
-    void update(const std::vector<std::size_t>& region) {
+    bool update(const std::vector<std::size_t>& region) {
 
       CGAL_precondition(region.size() > 0);
       if (region.size() == 1) { // create new reference line and direction
-        m_line_of_best_fit = Line();
+        m_direction_of_best_fit = CGAL::NULL_VECTOR;
-        m_direction_of_best_fit = Vector();
       } else { // update reference line and direction
+        if (m_direction_of_best_fit == CGAL::NULL_VECTOR)
+          return false; // all points are equal
         CGAL_precondition(region.size() >= 2);
         std::tie(m_line_of_best_fit, m_direction_of_best_fit) =
           get_line_and_direction(region);
       }
+      return true;
     }
 
     /// @}

Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_polyline/One_ring_neighbor_query.h

@@ -60,7 +60,7 @@ namespace Polyline {
       \brief initializes all internal data structures.
 
       \param input_range
-      an instance of `InputRange`
+      an instance of `InputRange` with polyline vertices
 
       \pre `input_range.size() > 0`
     */
@@ -100,8 +100,8 @@ namespace Polyline {
       const std::size_t n = m_input_range.size();
       const std::size_t im = (query_index + n - 1) % n;
       const std::size_t ip = (query_index + 1) % n;
-      neighbors.push_back(im);
       neighbors.push_back(ip);
+      neighbors.push_back(im);
     }
 
     /// @}

Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_segment_set/Least_squares_line_fit_region.h

@@ -112,9 +112,7 @@ namespace Segment_set {
       const Segment& query_segment = get(m_segment_map, key);
       const Point& query_source = query_segment.source();
       const Point& query_target = query_segment.target();
-      CGAL_precondition(query_source != query_target);
       const Vector query_direction(query_source, query_target);
-      CGAL_precondition(query_direction != Vector());
 
       const FT squared_distance_to_fitted_line =
         get_max_squared_distance(query_segment);
@@ -139,7 +137,7 @@ namespace Segment_set {
       return (region.size() >= m_min_region_size);
     }
 
-    void update(const std::vector<std::size_t>& region) {
+    bool update(const std::vector<std::size_t>& region) {
 
       CGAL_precondition(region.size() > 0);
       if (region.size() == 1) { // create new reference line and direction
@@ -152,8 +150,9 @@ namespace Segment_set {
         const Segment& segment = get(m_segment_map, key);
         const Point& source = segment.source();
         const Point& target = segment.target();
-        CGAL_precondition(source != target);
+        if (source == target) return false;
 
+        CGAL_precondition(source != target);
         m_line_of_best_fit = Line(source, target);
         m_direction_of_best_fit = m_line_of_best_fit.to_vector();
 
@@ -162,6 +161,7 @@ namespace Segment_set {
         std::tie(m_line_of_best_fit, m_direction_of_best_fit) =
           get_line_and_direction(region);
       }
+      return true;
     }
 
     /// @}

Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_segment_set/Least_squares_line_fit_sorting.h

@@ -120,7 +120,12 @@ namespace Segment_set {
         neighbors.clear();
         m_neighbor_query(i, neighbors);
         neighbors.push_back(i);
-        m_scores[i] = m_segment_set_traits.create_line(
+        const auto& key = *(m_input_range.begin() + i);
+        const auto& segment = get(m_segment_map, key);
+        const auto& source = segment.source();
+        const auto& target = segment.target();
+        if (source == target) m_scores[i] = FT(0); // put it at the very back
+        else m_scores[i] = m_segment_set_traits.create_line(
           m_input_range, m_segment_map, neighbors).second;
       }
     }

Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h

@@ -147,7 +147,7 @@ namespace internal {
       const auto& item = get(item_map, key);
       items.push_back(iconverter(item));
     }
-    CGAL_postcondition(items.size() == region.size());
+    CGAL_precondition(items.size() == region.size());
 
     ILine_2 fitted_line;
     IPoint_2 fitted_centroid;
@@ -201,7 +201,7 @@ namespace internal {
       const auto& item = get(item_map, key);
       items.push_back(iconverter(item));
     }
-    CGAL_postcondition(items.size() == region.size());
+    CGAL_precondition(items.size() == region.size());
 
     ILine_3 fitted_line;
     IPoint_3 fitted_centroid;
@@ -260,7 +260,7 @@ namespace internal {
       const auto& item = get(item_map, key);
       items.push_back(iconverter(item));
     }
-    CGAL_postcondition(items.size() == region.size());
+    CGAL_precondition(items.size() == region.size());
 
     IPlane_3 fitted_plane;
     IPoint_3 fitted_centroid;
@@ -311,14 +311,14 @@ namespace internal {
 
       const auto hedge = halfedge(face, face_graph);
       const auto vertices = vertices_around_face(hedge, face_graph);
-      CGAL_postcondition(vertices.size() > 0);
+      CGAL_precondition(vertices.size() > 0);
 
       for (const auto vertex : vertices) {
         const auto& point = get(vertex_to_point_map, vertex);
         points.push_back(iconverter(point));
       }
     }
-    CGAL_postcondition(points.size() >= region.size());
+    CGAL_precondition(points.size() >= region.size());
 
     IPlane_3 fitted_plane;
     IPoint_3 fitted_centroid;

Shape_detection/todo.md

@@ -1,46 +1,4 @@
-- add region growing on segments
-- move all free functions to the same file
-- overload the line functions
-- do we need free functions here
 - update the docs
 - ---- submission ----
 - update the polyhedron demo
-- add new tests: polyline (2D, 3D), polyline with sorting (2D, 3D), free functions, randomness, strict tests on the data from PMP, tests similar to the basic_example
+- add new tests: polyline (2D, 3D), polyline with sorting (2D, 3D), free functions, randomness, strict tests on the data from PMP, tests similar to the basic_example, check polylines with equal points
-- ---- algorithm ----
-- create the polyline graph
-- create face to index map
-- create face to region map
-- iterate over all edges
-  - skip those, which have the same region neighbor
-  - preserve the boundary edges
-  - for each valid edge set its index as a neighbor for the source and target vertex
-  - use std map to map the index from the global range of vertices to the polyline vertices
-  - for each vertex save neighbor regions as well
-- create the neighbor query based on the saved neighbors
-- create the region type based on the neighbor regions
-  - skip corners (have more than 2 neighbors)
-  - split all other vertices into regions connecting corners
-  - corners will be returned as unassigned items
-- simplify the graph using the polyline regions
-  - iterate over all vertices in each region
-  - find the corner vertices for each region (must be two)
-  - if both corners are boundary, use region growing to detect linear segments
-  - insert boundary points of each linear segment as new corners
-  - map each corner to the region and take the opposite corner index (each region has max two corners)
-  - this creates new vertex range with neighbors
-  - while creating vertices also insert edges with neighbor regions
-  - based on the vertices and edges create oriented faces by possibly using sorting on Direction_2
-  - these faces are not necessary convex
-  - move each face vertex to the position of the intersection point between all neighbor region planes
-  - triangulate each face and tag exterior triangles beyond its boundaries, which are constraints in CDT
-- that gives the new simplified surface mesh
-
-- add region type for segments
-- the input graph can be split into polylines based on corners
-- regularize_face_selection_borders() or better apply graph cut on the computed regions to close holes
-- error in the graph cut can be controlled by the size to the ideal position of the corners
-
-- detect planes
-- regularize regions
-- decimate mesh
-- update corner positions