fixed examples + some other small stuff

Shape_detection/benchmark/Shape_detection/benchmark_region_growing_on_point_set_2.cpp

@@ -8,7 +8,7 @@
 #include <iterator>
 
 // CGAL includes.
-#include <CGAL/Timer.h>
+#include <CGAL/Real_timer.h>
 #include <CGAL/property_map.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
@@ -32,7 +32,7 @@ using Neighbor_query = SD::Point_set::Sphere_neighbor_query<Kernel, Input_range,
 using Region_type    = SD::Point_set::Least_squares_line_fit_region<Kernel, Input_range, Point_map, Normal_map>;
 using Region_growing = SD::Region_growing<Input_range, Neighbor_query, Region_type>;
 
-using Timer  = CGAL::Timer;
+using Timer  = CGAL::Real_timer;
 using Region = std::vector<std::size_t>;
 
 void benchmark_region_growing_on_point_set_2(

Shape_detection/benchmark/Shape_detection/benchmark_region_growing_on_point_set_3.cpp

@@ -8,7 +8,7 @@
 #include <iterator>
 
 // CGAL includes.
-#include <CGAL/Timer.h>
+#include <CGAL/Real_timer.h>
 #include <CGAL/property_map.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
@@ -34,7 +34,7 @@ using Neighbor_query = SD::Point_set::K_neighbor_query<Kernel, Input_range, Poin
 using Region_type    = SD::Point_set::Least_squares_plane_fit_region<Kernel, Input_range, Point_map, Normal_map>;
 using Region_growing = SD::Region_growing<Input_range, Neighbor_query, Region_type>;
 
-using Timer  = CGAL::Timer;
+using Timer  = CGAL::Real_timer;
 using Region = std::vector<std::size_t>;
 
 void create_input_range(

Shape_detection/doc/Shape_detection/Concepts/RegionType.h

@@ -10,7 +10,10 @@ this region.
 
 \cgalHasModel
 - `CGAL::Shape_detection::Point_set::Least_squares_line_fit_region`
+- `CGAL::Shape_detection::Point_set::Least_squares_circle_fit_region`
 - `CGAL::Shape_detection::Point_set::Least_squares_plane_fit_region`
+- `CGAL::Shape_detection::Point_set::Least_squares_sphere_fit_region`
+- `CGAL::Shape_detection::Point_set::Least_squares_cylinder_fit_region`
 - `CGAL::Shape_detection::Segment_set::Least_squares_line_fit_region`
 - `CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_region`
 - `CGAL::Shape_detection::Polyline::Least_squares_line_fit_region`

Shape_detection/doc/Shape_detection/PackageDescription.txt

@@ -118,6 +118,7 @@ and the Region Growing approach for detecting shapes in a set of arbitrary items
 ### Point Set ###
 - `CGAL::Shape_detection::Point_set::K_neighbor_query<GeomTraits, InputRange, PointMap>`
 - `CGAL::Shape_detection::Point_set::Sphere_neighbor_query<GeomTraits, InputRange, PointMap>`
+
 - `CGAL::Shape_detection::Point_set::Least_squares_line_fit_region<GeomTraits, InputRange, PointMap, NormalMap>`
 - `CGAL::Shape_detection::Point_set::Least_squares_circle_fit_region<GeomTraits, InputRange, PointMap, NormalMap>`
 - `CGAL::Shape_detection::Point_set::Least_squares_plane_fit_region<GeomTraits, InputRange, PointMap, NormalMap>`

Shape_detection/doc/Shape_detection/Shape_detection.txt

@@ -105,6 +105,7 @@ A lower probability provides a higher reliability and determinism at the cost of
 
 
 \subsection Shape_detection_RANSACExamples Examples
+
 The main class `Shape_detection::Efficient_RANSAC` takes a template parameter `Shape_detection::Efficient_RANSAC_traits` that defines the geometric types and input format.
 Property maps provide a means to interface with the user-specific data structures.
 The first parameter of the `Shape_detection::Efficient_RANSAC_traits` class is the common `Kernel`.
@@ -373,6 +374,7 @@ to the quality of the least squares sphere fit applied to the neighbors of each
 - `CGAL::Shape_detection::Point_set::Least_squares_cylinder_fit_sorting` - indices of 3D input points are sorted with respect
 to the quality of the least squares cylinder fit applied to the neighbors of each point.
 
+
 \subsubsection Shape_detection_RegionGrowingPoints_parameters Parameters
 
 The classes in Section \ref Shape_detection_RegionGrowingPoints "Region Growing On Point Set" depend on a few parameters
@@ -469,6 +471,7 @@ The following example shows another similar example, this time detecting (infini
 
 \cgalExample{Shape_detection/region_growing_cylinders_on_point_set_3.cpp}
 
+
 \subsubsection Shape_detection_RegionGrowingPoints_performance Performance
 
 The main parameter that affects the region growing algorithm on a point set is the neighborhood size

Shape_detection/examples/Shape_detection/include/utils.h

@@ -16,6 +16,7 @@
 #include <CGAL/memory.h>
 #include <CGAL/IO/PLY.h>
 #include <CGAL/IO/OFF.h>
+#include <CGAL/Real_timer.h>
 #include <CGAL/property_map.h>
 #include <CGAL/IO/write_ply_points.h>
 #include <CGAL/IO/Color.h>

Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp

@@ -1,123 +1,120 @@
-#include <CGAL/Real_timer.h>
+#include "include/utils.h"
-#include <CGAL/Random.h>
-#include <CGAL/Simple_cartesian.h>
-
 #include <CGAL/Point_set_3.h>
 #include <CGAL/Point_set_3/IO.h>
-
+#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Shape_detection/Region_growing/Region_growing.h>
 #include <CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h>
-
 #include <boost/iterator/function_output_iterator.hpp>
 
-#include <fstream>
+// Typedefs.
-
 using Kernel   = CGAL::Simple_cartesian<double>;
-using Point_3 = Kernel::Point_3;
+using FT       = Kernel::FT;
-using Vector_3 = Kernel::Vector_3;
 using Point_2  = Kernel::Point_2;
+using Point_3  = Kernel::Point_3;
 using Vector_2 = Kernel::Vector_2;
+using Vector_3 = Kernel::Vector_3;
 
-using Point_set_3 = CGAL::Point_set_3<Point_3, Vector_3>;
 using Point_set_2 = CGAL::Point_set_3<Point_2, Vector_2>;
+using Point_set_3 = CGAL::Point_set_3<Point_3, Vector_3>;
+
 using Point_map  = Point_set_2::Point_map;
 using Normal_map = Point_set_2::Vector_map;
 
-namespace Shape_detection = CGAL::Shape_detection::Point_set;
+using Neighbor_query = CGAL::Shape_detection::Point_set::K_neighbor_query<Kernel, Point_set_2, Point_map>;
+using Region_type    = CGAL::Shape_detection::Point_set::Least_squares_circle_fit_region<Kernel, Point_set_2, Point_map, Normal_map>;
+using Sorting        = CGAL::Shape_detection::Point_set::Least_squares_circle_fit_sorting<Kernel, Point_set_2, Neighbor_query, Point_map>;
+using Region_growing = CGAL::Shape_detection::Region_growing<Point_set_2, Neighbor_query, Region_type, typename Sorting::Seed_map>;
 
-using Neighbor_query = Shape_detection::K_neighbor_query
+int main(int argc, char** argv) {
-  <Kernel, Point_set_2, Point_map>;
-using Region_type = Shape_detection::Least_squares_circle_fit_region
-  <Kernel, Point_set_2, Point_map, Normal_map>;
-using Sorting = Shape_detection::Least_squares_circle_fit_sorting
-  <Kernel, Point_set_2, Neighbor_query, Point_map>;
-using Region_growing = CGAL::Shape_detection::Region_growing
-  <Point_set_2, Neighbor_query, Region_type, typename Sorting::Seed_map>;
 
-int main (int argc, char** argv)
+  // Load ply data either from a local folder or a user-provided file.
-{
+  const bool is_default_input = argc > 1 ? false : true;
-  std::ifstream ifile (argc > 1 ? argv[1] : "data/circles.ply");
+  std::ifstream in(is_default_input ? "data/circles.ply" : argv[1]);
-  Point_set_3 points3;
-  ifile >> points3;
 
-  std::cerr << points3.size() << " points read" << std::endl;
+  CGAL::IO::set_ascii_mode(in);
-
+  if (!in) {
-  // Input should have normals
+    std::cerr << "ERROR: cannot read the input file!" << std::endl;
-  assert (points3.has_normal_map());
+    return EXIT_FAILURE;
-
-  Point_set_2 points;
-  points.add_normal_map();
-  for (Point_set_3::Index idx : points3)
-  {
-    const Point_3& p = points3.point(idx);
-    const Vector_3& n = points3.normal(idx);
-    points.insert (Point_2 (p.x(), p.y()), Vector_2 (n.x(), n.y()));
   }
 
-  // Default parameters for data/circles.ply
+  Point_set_3 point_set_3;
+  in >> point_set_3;
+  in.close();
+  std::cout << "* number of input points: " << point_set_3.size() << std::endl;
+  assert(is_default_input && point_set_3.size() == 1101);
+  assert(point_set_3.has_normal_map()); // input should have normals
+
+  // Create a 2D point set.
+  Point_set_2 point_set_2;
+  point_set_2.add_normal_map();
+  for (const auto& idx : point_set_3) {
+    const Point_3& point = point_set_3.point(idx);
+    const Vector_3& normal = point_set_3.normal(idx);
+    point_set_2.insert(
+      Point_2(point.x(), point.y()),
+      Vector_2(normal.x(), normal.y()));
+  }
+
+  // Default parameter values for the data file circles.ply.
   const std::size_t k               = 12;
-  const double max_distance = 0.01;
+  const FT          max_distance    = FT(1) / FT(100);
-  const double max_angle = 10.;
+  const FT          max_angle       = FT(10);
   const std::size_t min_region_size = 20;
 
+  // Create instances of the classes Neighbor_query and Region_type.
   Neighbor_query neighbor_query(
-    points, CGAL::parameters::
+    point_set_2, CGAL::parameters::k_neighbors(k).point_map(point_set_2.point_map()));
-    k_neighbors(k).
-    point_map(points.point_map()));
 
   Region_type region_type(
-    points, CGAL::parameters::
+    point_set_2,
+    CGAL::parameters::
     maximum_distance(max_distance).
     maximum_angle(max_angle).
     minimum_region_size(min_region_size).
-    point_map(points.point_map()).
+    point_map(point_set_2.point_map()).
-    normal_map(points.normal_map()));
+    normal_map(point_set_2.normal_map()));
 
-  // Sort indices
+  // Sort indices.
   Sorting sorting(
-    points, neighbor_query, CGAL::parameters::
+    point_set_2, neighbor_query, CGAL::parameters::point_map(point_set_2.point_map()));
-    point_map(points.point_map()));
   sorting.sort();
 
+  // Create an instance of the region growing class.
   Region_growing region_growing(
-    points, neighbor_query, region_type, sorting.seed_map());
+    point_set_2, neighbor_query, region_type, sorting.seed_map());
 
-  // Add maps to get colored output
+  // Add maps to get a colored output.
   Point_set_3::Property_map<unsigned char>
-    red = points3.add_property_map<unsigned char>("red", 0).first,
+    red   = point_set_3.add_property_map<unsigned char>("red"  , 0).first,
-    green = points3.add_property_map<unsigned char>("green", 0).first,
+    green = point_set_3.add_property_map<unsigned char>("green", 0).first,
-    blue = points3.add_property_map<unsigned char>("blue", 0).first;
+    blue  = point_set_3.add_property_map<unsigned char>("blue" , 0).first;
 
+  // Run the algorithm.
   CGAL::Random random;
+  std::size_t num_circles = 0;
+  region_growing.detect(
+    boost::make_function_output_iterator(
+      [&](const std::vector<std::size_t>& region) {
 
-  std::size_t nb_circles = 0;
+        // Assign a random color to each region.
-  CGAL::Real_timer timer;
+        const unsigned char r = static_cast<unsigned char>(random.get_int(64, 192));
-  timer.start();
+        const unsigned char g = static_cast<unsigned char>(random.get_int(64, 192));
-  region_growing.detect
+        const unsigned char b = static_cast<unsigned char>(random.get_int(64, 192));
-    (boost::make_function_output_iterator
+        for (const std::size_t idx : region) {
-     ([&](const std::vector<std::size_t>& region)
-      {
-        // Assign a random color to each region
-        unsigned char r = static_cast<unsigned char>(random.get_int(64, 192));
-        unsigned char g = static_cast<unsigned char>(random.get_int(64, 192));
-        unsigned char b = static_cast<unsigned char>(random.get_int(64, 192));
-        for (const std::size_t& idx : region)
-        {
           red[idx]   = r;
           green[idx] = g;
           blue[idx]  = b;
         }
-        ++ nb_circles;
+        ++num_circles;
-      }));
+      }
-  timer.stop();
+    )
-
+  );
-  std::cerr << nb_circles << " circles detected in "
+  std::cout << "* number of found circles: " << num_circles << std::endl;
-            << timer.time() << " seconds" << std::endl;
+  assert(is_default_input && num_circles == 0);
-
-  // Save in colored_circles.ply
-  std::ofstream out ("colored_circles.ply");
-  CGAL::IO::set_binary_mode (out);
-  out << points3;
 
+  // Save regions to a file.
+  std::ofstream out("circles_point_set_2.ply");
+  CGAL::IO::set_ascii_mode(out);
+  out << point_set_3;
   return EXIT_SUCCESS;
 }

Shape_detection/examples/Shape_detection/region_growing_cylinders_on_point_set_3.cpp

@@ -1,18 +1,14 @@
-#include <fstream>
+#include "include/utils.h"
-#include <CGAL/Real_timer.h>
-#include <CGAL/Random.h>
-#include <CGAL/Simple_cartesian.h>
-
 #include <CGAL/Point_set_3.h>
 #include <CGAL/Point_set_3/IO.h>
-
+#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Shape_detection/Region_growing/Region_growing.h>
 #include <CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h>
-
 #include <boost/iterator/function_output_iterator.hpp>
 
-
+// Typedefs.
 using Kernel   = CGAL::Simple_cartesian<double>;
+using FT       = Kernel::FT;
 using Point_3  = Kernel::Point_3;
 using Vector_3 = Kernel::Vector_3;
 
@@ -20,84 +16,84 @@ using Point_set = CGAL::Point_set_3<Point_3>;
 using Point_map  = typename Point_set::Point_map;
 using Normal_map = typename Point_set::Vector_map;
 
-namespace Shape_detection = CGAL::Shape_detection::Point_set;
+using Neighbor_query = CGAL::Shape_detection::Point_set::K_neighbor_query<Kernel, Point_set, Point_map>;
+using Region_type    = CGAL::Shape_detection::Point_set::Least_squares_cylinder_fit_region<Kernel, Point_set, Point_map, Normal_map>;
+using Region_growing = CGAL::Shape_detection::Region_growing<Point_set, Neighbor_query, Region_type>;
 
-using Neighbor_query = Shape_detection::K_neighbor_query
+int main(int argc, char** argv) {
-  <Kernel, Point_set, Point_map>;
-using Region_type = Shape_detection::Least_squares_cylinder_fit_region
-  <Kernel, Point_set, Point_map, Normal_map>;
-using Region_growing = CGAL::Shape_detection::Region_growing
-  <Point_set, Neighbor_query, Region_type>;
 
-int main (int argc, char** argv)
+  // Load ply data either from a local folder or a user-provided file.
-{
+  const bool is_default_input = argc > 1 ? false : true;
-  std::ifstream ifile (argc > 1 ? argv[1] : "data/cube.pwn");
+  std::ifstream in(is_default_input ? "data/cube.pwn" : argv[1]);
-  Point_set points;
-  ifile >> points;
 
-  std::cerr << points.size() << " points read" << std::endl;
+  CGAL::IO::set_ascii_mode(in);
+  if (!in) {
+    std::cerr << "ERROR: cannot read the input file!" << std::endl;
+    return EXIT_FAILURE;
+  }
 
-  // Input should have normals
+  Point_set point_set;
-  assert (points.has_normal_map());
+  in >> point_set;
+  in.close();
+  std::cout << "* number of input points: " << point_set.size() << std::endl;
+  assert(is_default_input && point_set.size() == 50000);
+  assert(point_set.has_normal_map()); // input should have normals
 
-  // Default parameters for data/cube.pwn
+  // Default parameter values for the data file cuble.pwn.
   const std::size_t k               = 24;
-  const double max_distance = 0.05;
+  const FT          max_distance    = FT(1) / FT(20);
-  const double max_angle = 5.;
+  const FT          max_angle       = FT(5);
   const std::size_t min_region_size = 200;
 
+  // Create instances of the classes Neighbor_query and Region_type.
   Neighbor_query neighbor_query(
-    points, CGAL::parameters::
+    point_set, CGAL::parameters::k_neighbors(k).point_map(point_set.point_map()));
-    k_neighbors(k).
-    point_map(points.point_map()));
 
   Region_type region_type(
-    points, CGAL::parameters::
+    point_set,
+    CGAL::parameters::
     maximum_distance(max_distance).
     maximum_angle(max_angle).
     minimum_region_size(min_region_size).
-    point_map(points.point_map()).
+    point_map(point_set.point_map()).
-    normal_map(points.normal_map()));
+    normal_map(point_set.normal_map()));
 
+  // Create an instance of the region growing class.
   Region_growing region_growing(
-    points, neighbor_query, region_type);
+    point_set, neighbor_query, region_type);
 
-  // Add maps to get colored output
+  // Add maps to get a colored output.
   Point_set::Property_map<unsigned char>
-    red = points.add_property_map<unsigned char>("red", 0).first,
+    red   = point_set.add_property_map<unsigned char>("red"  , 0).first,
-    green = points.add_property_map<unsigned char>("green", 0).first,
+    green = point_set.add_property_map<unsigned char>("green", 0).first,
-    blue = points.add_property_map<unsigned char>("blue", 0).first;
+    blue  = point_set.add_property_map<unsigned char>("blue" , 0).first;
 
+  // Run the algorithm.
   CGAL::Random random;
+  std::size_t num_cylinders = 0;
+  region_growing.detect(
+    boost::make_function_output_iterator(
+      [&](const std::vector<std::size_t>& region) {
 
-  std::size_t nb_cylinders = 0;
+        // Assign a random color to each region.
-  CGAL::Real_timer timer;
+        const unsigned char r = static_cast<unsigned char>(random.get_int(64, 192));
-  timer.start();
+        const unsigned char g = static_cast<unsigned char>(random.get_int(64, 192));
-  region_growing.detect
+        const unsigned char b = static_cast<unsigned char>(random.get_int(64, 192));
-    (boost::make_function_output_iterator
+        for (const std::size_t idx : region) {
-     ([&](const std::vector<std::size_t>& region)
-      {
-        // Assign a random color to each region
-        unsigned char r = static_cast<unsigned char>(random.get_int(64, 192));
-        unsigned char g = static_cast<unsigned char>(random.get_int(64, 192));
-        unsigned char b = static_cast<unsigned char>(random.get_int(64, 192));
-        for (const std::size_t& idx : region)
-        {
           red[idx]   = r;
           green[idx] = g;
           blue[idx]  = b;
         }
-        ++ nb_cylinders;
+        ++num_cylinders;
-      }));
+      }
-  timer.stop();
+    )
-
+  );
-  std::cerr << nb_cylinders << " cylinders detected in "
+  std::cout << "* number of found cylinders: " << num_cylinders << std::endl;
-            << timer.time() << " seconds" << std::endl;
+  assert(is_default_input && num_cylinders == 0);
-
-  // Save in colored_cylinders.ply
-  std::ofstream out ("colored_cylinders.ply");
-  CGAL::IO::set_binary_mode (out);
-  out << points;
 
+  // Save regions to a file.
+  std::ofstream out("cylinders_point_set_3.ply");
+  CGAL::IO::set_ascii_mode(out);
+  out << point_set;
   return EXIT_SUCCESS;
 }

Shape_detection/examples/Shape_detection/region_growing_lines_on_point_set_2.cpp

@@ -40,7 +40,7 @@ int main(int argc, char *argv[]) {
   std::cout << "* number of input points: " << point_set_2.size() << std::endl;
   assert(is_default_input && point_set_2.size() == 3634);
 
-  // Default parameter values for the data file point_set_2.xyz.
+  // Default parameter values for the data file buildings_outline.xyz.
   const FT          sphere_radius   = FT(5);
   const FT          max_distance    = FT(45) / FT(10);
   const FT          max_angle       = FT(45);
@@ -64,11 +64,11 @@ int main(int argc, char *argv[]) {
   // Run the algorithm.
   std::vector< std::vector<std::size_t> > regions;
   region_growing.detect(std::back_inserter(regions));
-  std::cout << "* number of found regions: " << regions.size() << std::endl;
+  std::cout << "* number of found lines: " << regions.size() << std::endl;
   assert(is_default_input && regions.size() == 65);
 
   // Save regions to a file.
-  const std::string fullpath = (argc > 2 ? argv[2] : "regions_point_set_2.ply");
+  const std::string fullpath = (argc > 2 ? argv[2] : "lines_point_set_2.ply");
   utils::save_point_regions_2<Kernel, Point_set_2, Point_map>(
     point_set_2, regions, fullpath);
   return EXIT_SUCCESS;

Shape_detection/examples/Shape_detection/region_growing_planes_on_point_set_3.cpp

@@ -40,7 +40,7 @@ int main(int argc, char *argv[]) {
   std::cout << "* number of input points: " << input_range.size() << std::endl;
   assert(is_default_input && input_range.size() == 8075);
 
-  // Default parameter values for the data file point_set_3.xyz.
+  // Default parameter values for the data file building.xyz.
   const std::size_t k               = 12;
   const FT          max_distance    = FT(2);
   const FT          max_angle       = FT(20);
@@ -73,11 +73,11 @@ int main(int argc, char *argv[]) {
     output_range, number_of_regions);
   region_growing.detect(
     boost::make_function_output_iterator(inserter));
-  std::cout << "* number of found regions: " << number_of_regions << std::endl;
+  std::cout << "* number of found planes: " << number_of_regions << std::endl;
   assert(is_default_input && number_of_regions == 7);
 
   // Save regions to a file.
-  const std::string fullpath = (argc > 2 ? argv[2] : "regions_point_set_3.ply");
+  const std::string fullpath = (argc > 2 ? argv[2] : "planes_point_set_3.ply");
   std::ofstream out(fullpath);
   out << output_range;
   out.close();

Shape_detection/examples/Shape_detection/region_growing_planes_on_polygon_mesh.cpp

@@ -48,7 +48,7 @@ int main(int argc, char *argv[]) {
   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 polygon_mesh.off.
+  // 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;
@@ -75,12 +75,12 @@ int main(int argc, char *argv[]) {
   // Run the algorithm.
   std::vector< std::vector<std::size_t> > regions;
   region_growing.detect(std::back_inserter(regions));
-  std::cout << "* number of found regions: " << regions.size() << std::endl;
+  std::cout << "* number of found planes: " << regions.size() << std::endl;
   assert(is_default_input && regions.size() == 326);
 
   // Save regions to a file only if it is stored in CGAL::Surface_mesh.
   #if defined(USE_SURFACE_MESH)
-    const std::string fullpath = (argc > 2 ? argv[2] : "regions_polygon_mesh.ply");
+    const std::string fullpath = (argc > 2 ? argv[2] : "planes_polygon_mesh.ply");
     utils::save_polygon_mesh_regions(polygon_mesh, regions, fullpath);
   #endif
   return EXIT_SUCCESS;

Shape_detection/examples/Shape_detection/region_growing_spheres_on_point_set_3.cpp

@@ -1,18 +1,14 @@
-#include <CGAL/Real_timer.h>
+#include "include/utils.h"
-#include <CGAL/Random.h>
-#include <CGAL/Simple_cartesian.h>
-
 #include <CGAL/Point_set_3.h>
 #include <CGAL/Point_set_3/IO.h>
-
+#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Shape_detection/Region_growing/Region_growing.h>
 #include <CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h>
-
 #include <boost/iterator/function_output_iterator.hpp>
 
-#include <fstream>
+// Typedefs.
-
 using Kernel   = CGAL::Simple_cartesian<double>;
+using FT       = Kernel::FT;
 using Point_3  = Kernel::Point_3;
 using Vector_3 = Kernel::Vector_3;
 
@@ -20,84 +16,84 @@ using Point_set = CGAL::Point_set_3<Point_3>;
 using Point_map  = typename Point_set::Point_map;
 using Normal_map = typename Point_set::Vector_map;
 
-namespace Shape_detection = CGAL::Shape_detection::Point_set;
+using Neighbor_query = CGAL::Shape_detection::Point_set::K_neighbor_query<Kernel, Point_set, Point_map>;
+using Region_type    = CGAL::Shape_detection::Point_set::Least_squares_sphere_fit_region<Kernel, Point_set, Point_map, Normal_map>;
+using Region_growing = CGAL::Shape_detection::Region_growing<Point_set, Neighbor_query, Region_type>;
 
-using Neighbor_query = Shape_detection::K_neighbor_query
+int main(int argc, char** argv) {
-  <Kernel, Point_set, Point_map>;
-using Region_type = Shape_detection::Least_squares_sphere_fit_region
-  <Kernel, Point_set, Point_map, Normal_map>;
-using Region_growing = CGAL::Shape_detection::Region_growing
-  <Point_set, Neighbor_query, Region_type>;
 
-int main (int argc, char** argv)
+  // Load ply data either from a local folder or a user-provided file.
-{
+  const bool is_default_input = argc > 1 ? false : true;
-  std::ifstream ifile (argc > 1 ? argv[1] : "data/spheres.ply");
+  std::ifstream in(is_default_input ? "data/spheres.ply" : argv[1]);
-  Point_set points;
-  ifile >> points;
 
-  std::cerr << points.size() << " points read" << std::endl;
+  CGAL::IO::set_ascii_mode(in);
+  if (!in) {
+    std::cerr << "ERROR: cannot read the input file!" << std::endl;
+    return EXIT_FAILURE;
+  }
 
-  // Input should have normals
+  Point_set point_set;
-  assert (points.has_normal_map());
+  in >> point_set;
+  in.close();
+  std::cout << "* number of input points: " << point_set.size() << std::endl;
+  assert(is_default_input && point_set.size() == 5969);
+  assert(point_set.has_normal_map()); // input should have normals
 
-  // Default parameters for data/spheres.ply
+  // Default parameter values for the data file spheres.ply.
   const std::size_t k               = 12;
-  const double max_distance = 0.01;
+  const FT          max_distance    = FT(1) / FT(100);
-  const double max_angle = 10.;
+  const FT          max_angle       = FT(10);
   const std::size_t min_region_size = 50;
 
+  // Create instances of the classes Neighbor_query and Region_type.
   Neighbor_query neighbor_query(
-    points, CGAL::parameters::
+    point_set, CGAL::parameters::k_neighbors(k).point_map(point_set.point_map()));
-    k_neighbors(k).
-    point_map(points.point_map()));
 
   Region_type region_type(
-    points, CGAL::parameters::
+    point_set,
+    CGAL::parameters::
     maximum_distance(max_distance).
     maximum_angle(max_angle).
     minimum_region_size(min_region_size).
-    point_map(points.point_map()).
+    point_map(point_set.point_map()).
-    normal_map(points.normal_map()));
+    normal_map(point_set.normal_map()));
 
+  // Create an instance of the region growing class.
   Region_growing region_growing(
-    points, neighbor_query, region_type);
+    point_set, neighbor_query, region_type);
 
-  // Add maps to get colored output
+  // Add maps to get a colored output.
   Point_set::Property_map<unsigned char>
-    red = points.add_property_map<unsigned char>("red", 0).first,
+    red   = point_set.add_property_map<unsigned char>("red"  , 0).first,
-    green = points.add_property_map<unsigned char>("green", 0).first,
+    green = point_set.add_property_map<unsigned char>("green", 0).first,
-    blue = points.add_property_map<unsigned char>("blue", 0).first;
+    blue  = point_set.add_property_map<unsigned char>("blue" , 0).first;
 
+  // Run the algorithm.
   CGAL::Random random;
+  std::size_t num_spheres = 0;
+  region_growing.detect(
+    boost::make_function_output_iterator(
+      [&](const std::vector<std::size_t>& region) {
 
-  std::size_t nb_spheres = 0;
+        // Assign a random color to each region.
-  CGAL::Real_timer timer;
+        const unsigned char r = static_cast<unsigned char>(random.get_int(64, 192));
-  timer.start();
+        const unsigned char g = static_cast<unsigned char>(random.get_int(64, 192));
-  region_growing.detect
+        const unsigned char b = static_cast<unsigned char>(random.get_int(64, 192));
-    (boost::make_function_output_iterator
+        for (const std::size_t idx : region) {
-     ([&](const std::vector<std::size_t>& region)
-      {
-        // Assign a random color to each region
-        unsigned char r = static_cast<unsigned char>(random.get_int(64, 192));
-        unsigned char g = static_cast<unsigned char>(random.get_int(64, 192));
-        unsigned char b = static_cast<unsigned char>(random.get_int(64, 192));
-        for (const std::size_t& idx : region)
-        {
           red[idx]   = r;
           green[idx] = g;
           blue[idx]  = b;
         }
-        ++ nb_spheres;
+        ++num_spheres;
-      }));
+      }
-  timer.stop();
+    )
-
+  );
-  std::cerr << nb_spheres << " spheres detected in "
+  std::cout << "* number of found spheres: " << num_spheres << std::endl;
-            << timer.time() << " seconds" << std::endl;
+  assert(is_default_input && num_spheres == 0);
-
-  // Save in colored_spheres.ply
-  std::ofstream out ("colored_spheres.ply");
-  CGAL::IO::set_binary_mode (out);
-  out << points;
 
+  // Save regions to a file.
+  std::ofstream out("spheres_point_set_3.ply");
+  CGAL::IO::set_ascii_mode(out);
+  out << point_set;
   return EXIT_SUCCESS;
 }

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

@@ -303,17 +303,20 @@ namespace Point_set {
       CGAL_precondition(query_index < m_input_range.size());
 
       // First, we need to integrate at least 6 points so that the
-      // computed circle means something
+      // computed circle means something.
-      if (indices.size() < 6)
+      if (indices.size() < 6) {
         return true;
+      }
 
       // TODO: Why do we get so many nan in this class?
-      if (std::isnan(m_radius))
+      if (std::isnan(m_radius)) {
         return false;
+      }
 
-      // If radius is out of bound, nothing fits, early ending
+      // If radius is out of bound, nothing fits, early ending.
-      if (m_radius < m_min_radius || m_radius > m_max_radius)
+      if (m_radius < m_min_radius || m_radius > m_max_radius) {
         return false;
+      }
 
       const auto& key = *(m_input_range.begin() + query_index);
       const Point_2& query_point = get(m_point_map, key);
@@ -321,11 +324,12 @@ namespace Point_set {
 
       const FT sq_dist = m_squared_distance_2(query_point, m_center);
       if (std::isnan(sq_dist)) return false;
-      FT distance_to_center = m_sqrt (sq_dist);
+      const FT distance_to_center = m_sqrt (sq_dist);
-      FT distance_to_circle = CGAL::abs (distance_to_center - m_radius);
+      const FT distance_to_circle = CGAL::abs (distance_to_center - m_radius);
 
-      if (distance_to_circle > m_distance_threshold)
+      if (distance_to_circle > m_distance_threshold) {
         return false;
+      }
 
       const FT sq_norm = normal * normal;
       if (std::isnan(sq_norm)) return false;
@@ -336,9 +340,9 @@ namespace Point_set {
       if (std::isnan(sq_ray)) return false;
       ray = ray / m_sqrt (sq_ray);
 
-      if (CGAL::abs (normal * ray) < m_cos_value_threshold)
+      if (CGAL::abs(normal * ray) < m_cos_value_threshold) {
         return false;
-
+      }
       return true;
     }
 
@@ -375,17 +379,14 @@ namespace Point_set {
     bool update(const std::vector<std::size_t>& region) {
 
       CGAL_precondition(region.size() > 0);
-      auto unary_function
+      auto unary_function = [&](const std::size_t& idx) -> const Point_2& {
-        = [&](const std::size_t& idx) -> const Point_2&
-          {
         return get(m_point_map, *(m_input_range.begin() + idx));
       };
 
-      internal::circle_fit
+      internal::create_circle_2(
-        (make_range(boost::make_transform_iterator
+        make_range(
-                    (region.begin(), unary_function),
+          boost::make_transform_iterator(region.begin(), unary_function),
-                    boost::make_transform_iterator
+          boost::make_transform_iterator(region.end(), unary_function)),
-                    (region.end(), unary_function)),
         m_sqrt, m_squared_distance_2, m_center, m_radius);
       return true;
     }

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

@@ -231,15 +231,17 @@ namespace Point_set {
         Local_point_2 fitted_center;
         Local_FT fitted_radius;
 
-        if (internal::circle_fit (points, sqrt, squared_distance_2, fitted_center, fitted_radius))
+        if (internal::create_circle_2(
-        {
+          points, sqrt, squared_distance_2, fitted_center, fitted_radius)) {
-          // Score is min squared distance to sphere
+
+          // Score is min squared distance to sphere.
           m_scores[i] = Local_FT(0);
-          for (const Local_point_2& p : points)
+          for (const Local_point_2& p : points) {
             m_scores[i] += abs(sqrt(squared_distance_2(p, fitted_center)) - fitted_radius);
           }
-        else
+        } else {
-          m_scores[i] = Local_FT(std::numeric_limits<double>::infinity());
+          m_scores[i] = Local_FT(std::numeric_limits<double>::max());
+        }
       }
     }
   };

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

@@ -305,17 +305,20 @@ namespace Point_set {
       CGAL_precondition(query_index < m_input_range.size());
 
       // First, we need to integrate at least 6 points so that the
-      // computed cylinder means something
+      // computed cylinder means something.
-      if (indices.size() < 6)
+      if (indices.size() < 6) {
         return true;
+      }
 
       // TODO: Why do we get so many nan in this class?
-      if (std::isnan(m_radius))
+      if (std::isnan(m_radius)) {
         return false;
+      }
 
-      // If radius is out of bound, nothing fits, early ending
+      // If radius is out of bound, nothing fits, early ending.
-      if (m_radius < m_min_radius || m_radius > m_max_radius)
+      if (m_radius < m_min_radius || m_radius > m_max_radius) {
         return false;
+      }
 
       const auto& key = *(m_input_range.begin() + query_index);
       const Point_3& query_point = get(m_point_map, key);
@@ -326,25 +329,26 @@ namespace Point_set {
       if (m_axis.to_vector() == Vector_3(0, 0, 0)) return false;
       const FT sq_dist = m_squared_distance_3(query_point, m_axis);
       if (std::isnan(sq_dist)) return false;
-      FT distance_to_center = m_sqrt (sq_dist);
+      const FT distance_to_center = m_sqrt(sq_dist);
-      FT distance_to_cylinder = CGAL::abs (distance_to_center - m_radius);
+      const FT distance_to_cylinder = CGAL::abs(distance_to_center - m_radius);
 
-      if (distance_to_cylinder > m_distance_threshold)
+      if (distance_to_cylinder > m_distance_threshold) {
         return false;
+      }
 
       const FT sq_norm = normal * normal;
       if (std::isnan(sq_norm)) return false;
       normal = normal / m_sqrt (sq_norm);
 
-      Point_3 proj = m_axis.projection(query_point);
+      const Point_3 proj = m_axis.projection(query_point);
       Vector_3 ray(proj, query_point);
       const FT sq_ray = ray * ray;
       if (std::isnan(sq_ray)) return false;
       ray = ray / m_sqrt(sq_ray);
 
-      if (CGAL::abs (normal * ray) < m_cos_value_threshold)
+      if (CGAL::abs(normal * ray) < m_cos_value_threshold) {
         return false;
-
+      }
       return true;
     }
 
@@ -382,25 +386,21 @@ namespace Point_set {
 
       CGAL_precondition(region.size() > 0);
 
-      // Shuffle to avoid always picking 2 close points
+      // Shuffle to avoid always picking 2 close points.
-      std::vector<std::size_t>& aregion
+      std::vector<std::size_t>& aregion =
-        = const_cast<std::vector<std::size_t>&>(region);
+        const_cast<std::vector<std::size_t>&>(region);
       cpp98::random_shuffle(aregion.begin(), aregion.end());
 
       using VT = typename std::iterator_traits<typename InputRange::const_iterator>::value_type;
-      auto unary_function
+      auto unary_function = [&](const std::size_t& idx) -> VT {
-        = [&](const std::size_t& idx) -> VT
-          {
         return *(m_input_range.begin() + idx);
       };
 
-      internal::cylinder_fit
+      internal::create_cylinder_3(
-        (make_range (boost::make_transform_iterator
+        make_range(
-                    (region.begin(), unary_function),
+          boost::make_transform_iterator(region.begin(), unary_function),
-                    boost::make_transform_iterator
+          boost::make_transform_iterator(region.end(), unary_function)),
-                    (region.end(), unary_function)),
+        m_point_map, m_normal_map, m_sqrt, m_squared_distance_3, m_axis, m_radius);
-        m_point_map, m_normal_map, m_sqrt, m_squared_distance_3,
-        m_axis, m_radius);
       return true;
     }
 

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

@@ -236,16 +236,14 @@ namespace Point_set {
       typename internal::Get_sqrt<Local_traits>::Sqrt sqrt;
       typename Local_traits::Compute_squared_distance_3 squared_distance_3;
 
-      for (std::size_t i = 0; i < m_input_range.size(); ++i)
+      for (std::size_t i = 0; i < m_input_range.size(); ++i) {
-      {
 
         neighbors.clear();
         m_neighbor_query(i, neighbors);
         neighbors.push_back(i);
 
         points.clear();
-        for (std::size_t j = 0; j < neighbors.size(); ++j)
+        for (std::size_t j = 0; j < neighbors.size(); ++j) {
-        {
           CGAL_precondition(neighbors[j] < m_input_range.size());
 
           const auto& key = *(m_input_range.begin() + neighbors[j]);
@@ -257,16 +255,18 @@ namespace Point_set {
         Local_line_3 fitted_line;
         Local_FT fitted_radius;
 
-        if (internal::cylinder_fit (points, Local_point_map(), Local_normal_map(),
+        if (internal::create_cylinder_3(
-                                    sqrt, squared_distance_3, fitted_line, fitted_radius))
+          points, Local_point_map(), Local_normal_map(),
-        {
+          sqrt, squared_distance_3, fitted_line, fitted_radius)) {
-          // Score is min squared distance to cylinder
+
+          // Score is min squared distance to cylinder.
           m_scores[i] = Local_FT(0);
-          for (const Local_pwn& pwn : points)
+          for (const Local_pwn& pwn : points) {
             m_scores[i] += abs(sqrt(squared_distance_3(pwn.first, fitted_line)) - fitted_radius);
           }
-        else
+        } else {
-          m_scores[i] = Local_FT(std::numeric_limits<double>::infinity());
+          m_scores[i] = Local_FT(std::numeric_limits<double>::max());
+        }
       }
     }
   };

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

@@ -308,17 +308,20 @@ namespace Point_set {
       CGAL_precondition(query_index < m_input_range.size());
 
       // First, we need to integrate at least 6 points so that the
-      // computed sphere means something
+      // computed sphere means something.
-      if (indices.size() < 6)
+      if (indices.size() < 6) {
         return true;
+      }
 
       // TODO: Why do we get so many nan in this class?
-      if (std::isnan(m_radius))
+      if (std::isnan(m_radius)) {
         return false;
+      }
 
-      // If radius is out of bound, nothing fits, early ending
+      // If radius is out of bound, nothing fits, early ending.
-      if (m_radius < m_min_radius || m_radius > m_max_radius)
+      if (m_radius < m_min_radius || m_radius > m_max_radius) {
         return false;
+      }
 
       const auto& key = *(m_input_range.begin() + query_index);
       const Point_3& query_point = get(m_point_map, key);
@@ -326,11 +329,12 @@ namespace Point_set {
 
       const FT sq_dist = m_squared_distance_3(query_point, m_center);
       if (std::isnan(sq_dist)) return false;
-      FT distance_to_center = m_sqrt (sq_dist);
+      const FT distance_to_center = m_sqrt(sq_dist);
-      FT distance_to_sphere = CGAL::abs (distance_to_center - m_radius);
+      const FT distance_to_sphere = CGAL::abs(distance_to_center - m_radius);
 
-      if (distance_to_sphere > m_distance_threshold)
+      if (distance_to_sphere > m_distance_threshold) {
         return false;
+      }
 
       const FT sq_norm = normal * normal;
       if (std::isnan(sq_norm)) return false;
@@ -341,9 +345,9 @@ namespace Point_set {
       if (std::isnan(sq_ray)) return false;
       ray = ray / m_sqrt (sq_ray);
 
-      if (CGAL::abs (normal * ray) < m_cos_value_threshold)
+      if (CGAL::abs(normal * ray) < m_cos_value_threshold) {
         return false;
-
+      }
       return true;
     }
 
@@ -380,17 +384,14 @@ namespace Point_set {
     bool update(const std::vector<std::size_t>& region) {
 
       CGAL_precondition(region.size() > 0);
-      auto unary_function
+      auto unary_function = [&](const std::size_t& idx) -> const Point_3& {
-        = [&](const std::size_t& idx) -> const Point_3&
-          {
         return get (m_point_map, *(m_input_range.begin() + idx));
       };
 
-      internal::sphere_fit
+      internal::create_sphere_3(
-        (make_range (boost::make_transform_iterator
+        make_range(
-                    (region.begin(), unary_function),
+          boost::make_transform_iterator(region.begin(), unary_function),
-                    boost::make_transform_iterator
+          boost::make_transform_iterator(region.end(), unary_function)),
-                    (region.end(), unary_function)),
         m_sqrt, m_squared_distance_3, m_center, m_radius);
       return true;
     }

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

@@ -212,16 +212,14 @@ namespace Point_set {
       typename internal::Get_sqrt<Local_traits>::Sqrt sqrt;
       typename Local_traits::Compute_squared_distance_3 squared_distance_3;
 
-      for (std::size_t i = 0; i < m_input_range.size(); ++i)
+      for (std::size_t i = 0; i < m_input_range.size(); ++i) {
-      {
 
         neighbors.clear();
         m_neighbor_query(i, neighbors);
         neighbors.push_back(i);
 
         points.clear();
-        for (std::size_t j = 0; j < neighbors.size(); ++j)
+        for (std::size_t j = 0; j < neighbors.size(); ++j) {
-        {
           CGAL_precondition(neighbors[j] < m_input_range.size());
 
           const auto& key = *(m_input_range.begin() + neighbors[j]);
@@ -232,15 +230,17 @@ namespace Point_set {
         Local_point_3 fitted_center;
         Local_FT fitted_radius;
 
-        if (internal::sphere_fit (points, sqrt, squared_distance_3, fitted_center, fitted_radius))
+        if (internal::create_sphere_3(
-        {
+          points, sqrt, squared_distance_3, fitted_center, fitted_radius)) {
-          // Score is min squared distance to sphere
+
+          // Score is min squared distance to sphere.
           m_scores[i] = Local_FT(0);
-          for (const Local_point_3& p : points)
+          for (const Local_point_3& p : points) {
             m_scores[i] += abs(sqrt(squared_distance_3(p, fitted_center)) - fitted_radius);
           }
-        else
+        } else {
-          m_scores[i] = Local_FT(std::numeric_limits<double>::infinity());
+          m_scores[i] = Local_FT(std::numeric_limits<double>::max());
+        }
       }
     }
   };

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

@@ -346,50 +346,57 @@ namespace internal {
     return std::make_pair(plane, static_cast<FT>(score));
   }
 
-template <typename PointRange,
+  template<
-          typename Sqrt, typename Squared_distance_2,
+  typename PointRange,
-          typename Point_2, typename FT>
+  typename Sqrt,
-bool circle_fit (const PointRange& points,
+  typename Squared_distance_2,
+  typename Point_2,
+  typename FT>
+  bool create_circle_2(
+    const PointRange& points,
     const Sqrt& sqrt,
     const Squared_distance_2& squared_distance_2,
-                 Point_2& center, FT& radius)
+    Point_2& center, FT& radius) {
-{
+
     using Diagonalize_traits = Eigen_diagonalize_traits<FT, 4>;
     using Covariance_matrix = typename Diagonalize_traits::Covariance_matrix;
+
     using Vector = typename Diagonalize_traits::Vector;
     using Matrix = typename Diagonalize_traits::Matrix;
 
     // Use bbox to compute diagonalization with smaller coordinates to
-  // avoid loss of precision when inverting large coordinates
+    // avoid loss of precision when inverting large coordinates.
-  Bbox_2 bbox = bbox_2 (points.begin(), points.end());
+    const Bbox_2 bbox = bbox_2(points.begin(), points.end());
 
-  // Circle least squares fitting,
+    // Circle least squares fitting:
-  // Circle of center (a,b) and radius R
+    // Circle of center (a, b) and radius R:
     // Ri = sqrt((xi - a)^2 + (yi - b)^2)
-  // Minimize Sum(Ri^2 - R^2)^2
+    // Minimize sum(Ri^2 - R^2)^2
-  // -> Minimize Sum(xi^2 + yi^2 − 2 a*xi − 2 b*yi + a^2 + b^2 − R^2)^2
+    // -> Minimize sum(xi^2 + yi^2 − 2 * a * xi − 2 * b * yi + a^2 + b^2 − R^2)^2
-  // let B=-2a ; C=-2b; D= a^2 + b^2 - R^2
+    // Let B = -2a; C = -2b; D = a^2 + b^2 - R^2; and ri = x^2 + y^2.
-  // let ri = x^2 + y^2
+    // -> Minimize sum(D + B * xi + C * yi + ri)^2
-  // -> Minimize Sum(D + B*xi + C*yi + ri)^2
+    // -> Minimize sum(1 + B / D * xi + C / D * yi + ri / D)^2
-  // -> Minimize Sum(1 + B/D*xi + C/D*yi + ri/D)^2
     // -> system of linear equations
     // -> diagonalize matrix
     //    NB   x   y   r
     //        xx  xy  xr
     //            yy  yr
     //                rr
-  //
+    // -> center coordinates = [
-  // -> center coordinates = -0.5 * eigenvector(1) / eigenvector(3) ; -0.5 * eigenvector(2) / eigenvector(3)
+    // -0.5 * eigenvector(1) / eigenvector(3);
-  Covariance_matrix A
+    // -0.5 * eigenvector(2) / eigenvector(3); ]
-    = { FT(0), FT(0), FT(0), FT(0), FT(0),
+    Covariance_matrix A = {
-        FT(0), FT(0), FT(0), FT(0), FT(0) };
+      FT(0), FT(0), FT(0), FT(0), FT(0),
+      FT(0), FT(0), FT(0), FT(0), FT(0)
+    };
 
     A[0] = static_cast<FT>(points.size());
-  for (const Point_2& p : points)
+    for (const Point_2& p : points) {
-  {
+
-    FT x = p.x() - bbox.xmin();
+      const FT x = p.x() - bbox.xmin();
-    FT y = p.y() - bbox.ymin();
+      const FT y = p.y() - bbox.ymin();
-    FT r = x*x + y*y;
+      const FT r = x * x + y * y;
+
       A[1] += x;
       A[2] += y;
       A[3] += r;
@@ -401,61 +408,75 @@ bool circle_fit (const PointRange& points,
       A[9] += r * r;
     }
 
-  Vector eigenvalues = { FT(0), FT(0), FT(0), FT(0) };
+    Vector eigenvalues = {
-  Matrix eigenvectors = { FT(0), FT(0), FT(0), FT(0),
+      FT(0), FT(0), FT(0), FT(0)
+    };
+    Matrix eigenvectors = {
       FT(0), FT(0), FT(0), FT(0),
       FT(0), FT(0), FT(0), FT(0),
-                          FT(0), FT(0), FT(0), FT(0) };
+      FT(0), FT(0), FT(0), FT(0),
+      FT(0), FT(0), FT(0), FT(0)
+    };
+    Diagonalize_traits::
+      diagonalize_selfadjoint_covariance_matrix(A, eigenvalues, eigenvectors);
 
-  Diagonalize_traits::diagonalize_selfadjoint_covariance_matrix
+    // Perfect line case, no circle can be fitted.
-    (A, eigenvalues, eigenvectors);
+    if (eigenvectors[3] == 0) {
-
-  // Perfect line case, no circle can be fitted
-  if (eigenvectors[3] == 0)
       return false;
+    }
 
-  center = Point_2 (bbox.xmin() - FT(0.5) * (eigenvectors[1] / eigenvectors[3]),
+    // Other cases.
-                    bbox.ymin() - FT(0.5) * (eigenvectors[2] / eigenvectors[3]));
+    const FT half = FT(1) / FT(2);
+    center = Point_2(
+      bbox.xmin() - half * (eigenvectors[1] / eigenvectors[3]),
+      bbox.ymin() - half * (eigenvectors[2] / eigenvectors[3]));
 
     radius = FT(0);
-  for (const Point_2& p : points)
+    for (const Point_2& p : points) {
       radius += sqrt(squared_distance_2(p, center));
-  radius /= points.size();
+    }
-
+    radius /= static_cast<FT>(points.size());
     return true;
   }
 
-template <typename PointRange,
+  template<
-          typename Sqrt, typename Squared_distance_3,
+  typename PointRange,
-          typename Point_3, typename FT>
+  typename Sqrt,
-bool sphere_fit (const PointRange& points,
+  typename Squared_distance_3,
+  typename Point_3,
+  typename FT>
+  bool create_sphere_3(
+    const PointRange& points,
     const Sqrt& sqrt,
     const Squared_distance_3& squared_distance_3,
-                 Point_3& center, FT& radius)
+    Point_3& center, FT& radius) {
-{
+
     using Diagonalize_traits = Eigen_diagonalize_traits<FT, 5>;
     using Covariance_matrix = typename Diagonalize_traits::Covariance_matrix;
+
     using Vector = typename Diagonalize_traits::Vector;
     using Matrix = typename Diagonalize_traits::Matrix;
 
     // Use bbox to compute diagonalization with smaller coordinates to
-  // avoid loss of precision when inverting large coordinates
+    // avoid loss of precision when inverting large coordinates.
-  Bbox_3 bbox = bbox_3 (points.begin(), points.end());
+    const Bbox_3 bbox = bbox_3(points.begin(), points.end());
 
-  // Sphere least squares fitting
+    // Sphere least squares fitting.
-  // (see Least_square_circle_fit_region for details about computation)
+    // See create_circle_2() above for more details on computation.
-  Covariance_matrix A
+    Covariance_matrix A = {
-    = { FT(0), FT(0), FT(0), FT(0), FT(0),
       FT(0), FT(0), FT(0), FT(0), FT(0),
-        FT(0), FT(0), FT(0), FT(0), FT(0) };
+      FT(0), FT(0), FT(0), FT(0), FT(0),
+      FT(0), FT(0), FT(0), FT(0), FT(0)
+    };
 
     A[0] = static_cast<FT>(points.size());
-  for (const Point_3& p : points)
+    for (const Point_3& p : points) {
-  {
+
-    FT x = p.x() - bbox.xmin();
+      const FT x = p.x() - bbox.xmin();
-    FT y = p.y() - bbox.ymin();
+      const FT y = p.y() - bbox.ymin();
-    FT z = p.z() - bbox.zmin();
+      const FT z = p.z() - bbox.zmin();
-    FT r = x*x + y*y + z*z;
+      const FT r = x * x + y * y + z * z;
+
       A[1] += x;
       A[2] += y;
       A[3] += z;
@@ -472,59 +493,75 @@ bool sphere_fit (const PointRange& points,
       A[14] += r * r;
     }
 
-  Vector eigenvalues = { FT(0), FT(0), FT(0), FT(0), FT(0) };
+    Vector eigenvalues = {
-  Matrix eigenvectors = { FT(0), FT(0), FT(0), FT(0), FT(0),
+      FT(0), FT(0), FT(0), FT(0), FT(0)
+    };
+    Matrix eigenvectors = {
       FT(0), FT(0), FT(0), FT(0), FT(0),
       FT(0), FT(0), FT(0), FT(0), FT(0),
       FT(0), FT(0), FT(0), FT(0), FT(0),
-                          FT(0), FT(0), FT(0), FT(0), FT(0) };
+      FT(0), FT(0), FT(0), FT(0), FT(0),
+      FT(0), FT(0), FT(0), FT(0), FT(0)
+    };
+    Diagonalize_traits::
+      diagonalize_selfadjoint_covariance_matrix(A, eigenvalues, eigenvectors);
 
-  Diagonalize_traits::diagonalize_selfadjoint_covariance_matrix
+    // Perfect plane case, no sphere can be fitted.
-    (A, eigenvalues, eigenvectors);
+    if (eigenvectors[4] == 0) {
-
-  // Perfect plane case, no sphere can be fitted
-  if (eigenvectors[4] == 0)
       return false;
+    }
 
-  center = Point_3 (bbox.xmin() - FT(0.5) * (eigenvectors[1] / eigenvectors[4]),
+    // Other cases.
-                    bbox.ymin() - FT(0.5) * (eigenvectors[2] / eigenvectors[4]),
+    const FT half = FT(1) / FT(2);
-                    bbox.zmin() - FT(0.5) * (eigenvectors[3] / eigenvectors[4]));
+    center = Point_3(
+      bbox.xmin() - half * (eigenvectors[1] / eigenvectors[4]),
+      bbox.ymin() - half * (eigenvectors[2] / eigenvectors[4]),
+      bbox.zmin() - half * (eigenvectors[3] / eigenvectors[4]));
 
     radius = FT(0);
-  for (const Point_3& p : points)
+    for (const Point_3& p : points) {
       radius += sqrt(squared_distance_3(p, center));
-  radius /= points.size();
+    }
-
+    radius /= static_cast<FT>(points.size());
     return true;
   }
 
-template <typename PointRange, typename PointMap, typename NormalMap,
+  template<
-          typename Sqrt, typename Squared_distance_3,
+  typename PointRange,
-          typename Line_3, typename FT>
+  typename PointMap,
-bool cylinder_fit (const PointRange& points,
+  typename NormalMap,
-                   PointMap point_map, NormalMap normal_map,
+  typename Sqrt,
+  typename Squared_distance_3,
+  typename Line_3,
+  typename FT>
+  bool create_cylinder_3(
+    const PointRange& points,
+    PointMap point_map,
+    NormalMap normal_map,
     const Sqrt& sqrt,
     const Squared_distance_3& squared_distance_3,
-                   Line_3& axis, FT& radius)
+    Line_3& axis, FT& radius) {
-{
+
     using Point_3 = typename boost::property_traits<PointMap>::value_type;
     using Vector_3 = typename boost::property_traits<NormalMap>::value_type;
 
     const Point_3& ref = get(point_map, *(points.begin()));
 
-  Vector_3 mean_axis = CGAL::NULL_VECTOR;
-  std::size_t nb = 0;
     radius = FT(0);
+    std::size_t nb = 0;
+    Vector_3 mean_axis = CGAL::NULL_VECTOR;
     Point_3 point_on_axis = ORIGIN;
-  for (std::size_t i = 0; i < points.size() - 1; ++ i)
+
-  {
+    for (std::size_t i = 0; i < points.size() - 1; ++i) {
+
       Vector_3 v0 = get(normal_map, *std::next(points.begin(), i));
       v0 = v0 / sqrt(v0 * v0);
       Vector_3 v1 = get(normal_map, *std::next(points.begin(), i + 1));
       v1 = v1 / sqrt(v1 * v1);
       Vector_3 axis = cross_product(v0, v1);
-    if (sqrt(axis.squared_length()) < (FT)(0.01))
+      if (sqrt(axis.squared_length()) < FT(1) / FT(100)) {
         continue;
+      }
       axis = axis / sqrt(axis * axis);
 
       const Point_3& p0 = get(point_map, *std::next(points.begin(), i));
@@ -536,47 +573,42 @@ bool cylinder_fit (const PointRange& points,
       Vector_3 ydir = cross_product(axis, xdir);
       ydir = ydir / sqrt (ydir * ydir);
 
-    FT v1x = v1 * ydir;
+      const FT v1x =  v1 * ydir;
-    FT v1y = -v1 * xdir;
+      const FT v1y = -v1 * xdir;
 
       Vector_3 d(p0, p1);
+      const FT ox = xdir * d;
+      const FT oy = ydir * d;
+      const FT ldist = v1x * ox + v1y * oy;
 
-    FT ox = xdir * d;
+      FT r = ldist / v1x;
-    FT oy = ydir * d;
+      Point_3 point = p0 + xdir * r;
-    FT ldist = v1x * ox + v1y * oy;
+      const Line_3 line(point, axis);
-
-    FT radius = ldist / v1x;
-
-    Point_3 point = p0 + xdir * radius;
-    Line_3 line (point, axis);
       point = line.projection(ref);
-
       point_on_axis = barycenter(point_on_axis, static_cast<FT>(nb), point, FT(1));
+      r += abs(r);
 
-    radius += abs(radius);
+      if (nb != 0 && axis * mean_axis < 0) {
-
-    if (nb != 0 && axis * mean_axis < 0)
         axis = -axis;
-
+      }
       mean_axis = mean_axis + axis;
       ++nb;
     }
 
-  if (nb == 0)
+    if (nb == 0) {
       return false;
+    }
 
     mean_axis = mean_axis / sqrt(mean_axis * mean_axis);
     axis = Line_3(point_on_axis, mean_axis);
-  radius /= nb;
+    radius /= static_cast<FT>(nb);
 
     radius = FT(0);
-  for (const auto& p : points)
+    for (const auto& p : points) {
-  {
       const Point_3& p0 = get(point_map, p);
       radius += sqrt(squared_distance_3(p0, axis));
     }
-  radius /= points.size();
+    radius /= static_cast<FT>(points.size());
-
     return true;
   }
 

Shape_detection/package_info/Shape_detection/long_description.txt

@@ -9,7 +9,7 @@ Additional shapes can be detected, given a custom shape class by the user.
 
 In addition to the RANSAC algorithm, the generic Region Growing algorithm is implemented.
 This algorithm can detect shapes in a set of arbitrary items. In particular,
-CGAL provides five instances of the algorithm: detecting lines in a 2D point set,
+CGAL provides several instances of the algorithm: detecting lines and circles in a 2D point set,
-detecting planes in a 3D point set, detecting planes on a polygon mesh, detecting lines
+detecting planes, spheres, and cylinders in a 3D point set, detecting planes on a polygon mesh, detecting lines
 in a 2D/3D segment set, and detecting lines on a 2D/3D polyline.
 Other custom instances can be easily added by the user.

Shape_detection/test/Shape_detection/test_region_growing_on_point_set_2_with_sorting.cpp

@@ -107,8 +107,8 @@ bool test (int argc, char** argv, const std::size_t minr, const std::size_t maxr
 int main(int argc, char *argv[]) {
 
   return ((
-    test<Line_region, Line_sorting>(argc, argv, 62, 66) // &&
+    test<Line_region, Line_sorting>(argc, argv, 62, 66) &&
-    // test<Circle_region, Circle_sorting>(argc, argv, 196, 200)
+    test<Circle_region, Circle_sorting>(argc, argv, 196, 200)
   ) ? EXIT_SUCCESS : EXIT_FAILURE );
 
 }

Shape_detection/test/Shape_detection/test_region_growing_on_point_set_3_with_sorting.cpp

@@ -138,7 +138,6 @@ int main(int argc, char *argv[]) {
   if (!success)
     return EXIT_FAILURE;
 
-  /*
   success =
     test<Sphere_region, Sphere_sorting>
     (argc, argv,
@@ -154,7 +153,7 @@ int main(int argc, char *argv[]) {
        const std::size_t min_region_size = 50;
        // No constraint on radius
        const double min_radius = 0.;
-       const double max_radius = std::numeric_limits<double>::infinity();
+       const double max_radius = std::numeric_limits<double>::max();
        return Sphere_region
          (input_range, tolerance, max_angle, min_region_size,
           min_radius, max_radius,
@@ -182,7 +181,7 @@ int main(int argc, char *argv[]) {
        const std::size_t min_region_size = 200;
        // No constraint on radius
        const double min_radius = 0.;
-       const double max_radius = std::numeric_limits<double>::infinity();
+       const double max_radius = std::numeric_limits<double>::max();
        return Cylinder_region
          (input_range, tolerance, max_angle, min_region_size,
           min_radius, max_radius,
@@ -193,7 +192,7 @@ int main(int argc, char *argv[]) {
         return (r.size() > 2 && r.size() < 30);
       });
   if (!success)
-    return EXIT_FAILURE; */
+    return EXIT_FAILURE;
 
   return EXIT_SUCCESS;
 }