diff --git a/Shape_detection/benchmark/Shape_detection/benchmark_region_growing_on_point_set_2.cpp b/Shape_detection/benchmark/Shape_detection/benchmark_region_growing_on_point_set_2.cpp index d2f20137bae..d58cb1fb0db 100644 --- a/Shape_detection/benchmark/Shape_detection/benchmark_region_growing_on_point_set_2.cpp +++ b/Shape_detection/benchmark/Shape_detection/benchmark_region_growing_on_point_set_2.cpp @@ -8,7 +8,7 @@ #include // CGAL includes. -#include +#include #include #include @@ -32,7 +32,7 @@ using Neighbor_query = SD::Point_set::Sphere_neighbor_query; using Region_growing = SD::Region_growing; -using Timer = CGAL::Timer; +using Timer = CGAL::Real_timer; using Region = std::vector; void benchmark_region_growing_on_point_set_2( diff --git a/Shape_detection/benchmark/Shape_detection/benchmark_region_growing_on_point_set_3.cpp b/Shape_detection/benchmark/Shape_detection/benchmark_region_growing_on_point_set_3.cpp index 60b03393c8a..9f3121e8924 100644 --- a/Shape_detection/benchmark/Shape_detection/benchmark_region_growing_on_point_set_3.cpp +++ b/Shape_detection/benchmark/Shape_detection/benchmark_region_growing_on_point_set_3.cpp @@ -8,7 +8,7 @@ #include // CGAL includes. -#include +#include #include #include @@ -34,7 +34,7 @@ using Neighbor_query = SD::Point_set::K_neighbor_query; using Region_growing = SD::Region_growing; -using Timer = CGAL::Timer; +using Timer = CGAL::Real_timer; using Region = std::vector; void create_input_range( diff --git a/Shape_detection/doc/Shape_detection/Concepts/RegionType.h b/Shape_detection/doc/Shape_detection/Concepts/RegionType.h index 06242426799..4b89b0c7084 100644 --- a/Shape_detection/doc/Shape_detection/Concepts/RegionType.h +++ b/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` diff --git a/Shape_detection/doc/Shape_detection/PackageDescription.txt b/Shape_detection/doc/Shape_detection/PackageDescription.txt index c75a154171f..ffe08f8e0b5 100644 --- a/Shape_detection/doc/Shape_detection/PackageDescription.txt +++ b/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` - `CGAL::Shape_detection::Point_set::Sphere_neighbor_query` + - `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` diff --git a/Shape_detection/doc/Shape_detection/Shape_detection.txt b/Shape_detection/doc/Shape_detection/Shape_detection.txt index 9b90999d8a0..31fcd504a12 100644 --- a/Shape_detection/doc/Shape_detection/Shape_detection.txt +++ b/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 diff --git a/Shape_detection/examples/Shape_detection/data/circles.ply b/Shape_detection/examples/Shape_detection/data/circles.ply index e8767879066..a9a9c15aa3c 100644 Binary files a/Shape_detection/examples/Shape_detection/data/circles.ply and b/Shape_detection/examples/Shape_detection/data/circles.ply differ diff --git a/Shape_detection/examples/Shape_detection/data/spheres.ply b/Shape_detection/examples/Shape_detection/data/spheres.ply index 13901fd744f..6eb225f1b02 100644 Binary files a/Shape_detection/examples/Shape_detection/data/spheres.ply and b/Shape_detection/examples/Shape_detection/data/spheres.ply differ diff --git a/Shape_detection/examples/Shape_detection/include/utils.h b/Shape_detection/examples/Shape_detection/include/utils.h index af523ef86a2..93e726cd698 100644 --- a/Shape_detection/examples/Shape_detection/include/utils.h +++ b/Shape_detection/examples/Shape_detection/include/utils.h @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include diff --git a/Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp b/Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp index 149dc9d9da3..ae87852c3e4 100644 --- a/Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp +++ b/Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp @@ -1,123 +1,120 @@ -#include -#include -#include - +#include "include/utils.h" #include #include - +#include #include #include - #include -#include - -using Kernel = CGAL::Simple_cartesian; -using Point_3 = Kernel::Point_3; -using Vector_3 = Kernel::Vector_3; -using Point_2 = Kernel::Point_2; +// Typedefs. +using Kernel = CGAL::Simple_cartesian; +using FT = Kernel::FT; +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; using Point_set_2 = CGAL::Point_set_3; -using Point_map = Point_set_2::Point_map; +using Point_set_3 = CGAL::Point_set_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; +using Region_type = CGAL::Shape_detection::Point_set::Least_squares_circle_fit_region; +using Sorting = CGAL::Shape_detection::Point_set::Least_squares_circle_fit_sorting; +using Region_growing = CGAL::Shape_detection::Region_growing; -using Neighbor_query = Shape_detection::K_neighbor_query - ; -using Region_type = Shape_detection::Least_squares_circle_fit_region - ; -using Sorting = Shape_detection::Least_squares_circle_fit_sorting - ; -using Region_growing = CGAL::Shape_detection::Region_growing - ; +int main(int argc, char** argv) { -int main (int argc, char** argv) -{ - std::ifstream ifile (argc > 1 ? argv[1] : "data/circles.ply"); - Point_set_3 points3; - ifile >> points3; + // Load ply data either from a local folder or a user-provided file. + const bool is_default_input = argc > 1 ? false : true; + std::ifstream in(is_default_input ? "data/circles.ply" : argv[1]); - std::cerr << points3.size() << " points read" << std::endl; - - // Input should have normals - assert (points3.has_normal_map()); - - 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())); + CGAL::IO::set_ascii_mode(in); + if (!in) { + std::cerr << "ERROR: cannot read the input file!" << std::endl; + return EXIT_FAILURE; } - // Default parameters for data/circles.ply - const std::size_t k = 12; - const double max_distance = 0.01; - const double max_angle = 10.; + 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 FT max_distance = FT(1) / FT(100); + 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:: - k_neighbors(k). - point_map(points.point_map())); + point_set_2, CGAL::parameters::k_neighbors(k).point_map(point_set_2.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()). - normal_map(points.normal_map())); + point_map(point_set_2.point_map()). + normal_map(point_set_2.normal_map())); - // Sort indices + // Sort indices. Sorting sorting( - points, neighbor_query, CGAL::parameters:: - point_map(points.point_map())); + point_set_2, neighbor_query, CGAL::parameters::point_map(point_set_2.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 - red = points3.add_property_map("red", 0).first, - green = points3.add_property_map("green", 0).first, - blue = points3.add_property_map("blue", 0).first; + red = point_set_3.add_property_map("red" , 0).first, + green = point_set_3.add_property_map("green", 0).first, + blue = point_set_3.add_property_map("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& region) { - std::size_t nb_circles = 0; - CGAL::Real_timer timer; - timer.start(); - region_growing.detect - (boost::make_function_output_iterator - ([&](const std::vector& region) - { - // Assign a random color to each region - unsigned char r = static_cast(random.get_int(64, 192)); - unsigned char g = static_cast(random.get_int(64, 192)); - unsigned char b = static_cast(random.get_int(64, 192)); - for (const std::size_t& idx : region) - { - red[idx] = r; + // Assign a random color to each region. + const unsigned char r = static_cast(random.get_int(64, 192)); + const unsigned char g = static_cast(random.get_int(64, 192)); + const unsigned char b = static_cast(random.get_int(64, 192)); + for (const std::size_t idx : region) { + red[idx] = r; green[idx] = g; - blue[idx] = b; + blue[idx] = b; } - ++ nb_circles; - })); - timer.stop(); - - std::cerr << nb_circles << " circles detected in " - << timer.time() << " seconds" << std::endl; - - // Save in colored_circles.ply - std::ofstream out ("colored_circles.ply"); - CGAL::IO::set_binary_mode (out); - out << points3; + ++num_circles; + } + ) + ); + std::cout << "* number of found circles: " << num_circles << std::endl; + assert(is_default_input && num_circles == 0); + // 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; } diff --git a/Shape_detection/examples/Shape_detection/region_growing_cylinders_on_point_set_3.cpp b/Shape_detection/examples/Shape_detection/region_growing_cylinders_on_point_set_3.cpp index f73d30efe08..ba1b82a43f2 100644 --- a/Shape_detection/examples/Shape_detection/region_growing_cylinders_on_point_set_3.cpp +++ b/Shape_detection/examples/Shape_detection/region_growing_cylinders_on_point_set_3.cpp @@ -1,103 +1,99 @@ -#include -#include -#include -#include - +#include "include/utils.h" #include #include - +#include #include #include - #include - -using Kernel = CGAL::Simple_cartesian; -using Point_3 = Kernel::Point_3; +// Typedefs. +using Kernel = CGAL::Simple_cartesian; +using FT = Kernel::FT; +using Point_3 = Kernel::Point_3; using Vector_3 = Kernel::Vector_3; -using Point_set = CGAL::Point_set_3; -using Point_map = typename Point_set::Point_map; +using Point_set = CGAL::Point_set_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; +using Region_type = CGAL::Shape_detection::Point_set::Least_squares_cylinder_fit_region; +using Region_growing = CGAL::Shape_detection::Region_growing; -using Neighbor_query = Shape_detection::K_neighbor_query - ; -using Region_type = Shape_detection::Least_squares_cylinder_fit_region - ; -using Region_growing = CGAL::Shape_detection::Region_growing - ; +int main(int argc, char** argv) { -int main (int argc, char** argv) -{ - std::ifstream ifile (argc > 1 ? argv[1] : "data/cube.pwn"); - Point_set points; - ifile >> points; + // Load ply data either from a local folder or a user-provided file. + const bool is_default_input = argc > 1 ? false : true; + std::ifstream in(is_default_input ? "data/cube.pwn" : argv[1]); - 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 - assert (points.has_normal_map()); + Point_set point_set; + 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 - const std::size_t k = 24; - const double max_distance = 0.05; - const double max_angle = 5.; + // Default parameter values for the data file cuble.pwn. + const std::size_t k = 24; + const FT max_distance = FT(1) / FT(20); + 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:: - k_neighbors(k). - point_map(points.point_map())); + point_set, CGAL::parameters::k_neighbors(k).point_map(point_set.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()). - normal_map(points.normal_map())); + point_map(point_set.point_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 - red = points.add_property_map("red", 0).first, - green = points.add_property_map("green", 0).first, - blue = points.add_property_map("blue", 0).first; + red = point_set.add_property_map("red" , 0).first, + green = point_set.add_property_map("green", 0).first, + blue = point_set.add_property_map("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& region) { - std::size_t nb_cylinders = 0; - CGAL::Real_timer timer; - timer.start(); - region_growing.detect - (boost::make_function_output_iterator - ([&](const std::vector& region) - { - // Assign a random color to each region - unsigned char r = static_cast(random.get_int(64, 192)); - unsigned char g = static_cast(random.get_int(64, 192)); - unsigned char b = static_cast(random.get_int(64, 192)); - for (const std::size_t& idx : region) - { - red[idx] = r; + // Assign a random color to each region. + const unsigned char r = static_cast(random.get_int(64, 192)); + const unsigned char g = static_cast(random.get_int(64, 192)); + const unsigned char b = static_cast(random.get_int(64, 192)); + for (const std::size_t idx : region) { + red[idx] = r; green[idx] = g; - blue[idx] = b; + blue[idx] = b; } - ++ nb_cylinders; - })); - timer.stop(); - - std::cerr << nb_cylinders << " cylinders detected in " - << timer.time() << " seconds" << std::endl; - - // Save in colored_cylinders.ply - std::ofstream out ("colored_cylinders.ply"); - CGAL::IO::set_binary_mode (out); - out << points; + ++num_cylinders; + } + ) + ); + std::cout << "* number of found cylinders: " << num_cylinders << std::endl; + assert(is_default_input && num_cylinders == 0); + // 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; } diff --git a/Shape_detection/examples/Shape_detection/region_growing_lines_on_point_set_2.cpp b/Shape_detection/examples/Shape_detection/region_growing_lines_on_point_set_2.cpp index a1862c82e60..e3fbc9bf791 100644 --- a/Shape_detection/examples/Shape_detection/region_growing_lines_on_point_set_2.cpp +++ b/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 > 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( point_set_2, regions, fullpath); return EXIT_SUCCESS; diff --git a/Shape_detection/examples/Shape_detection/region_growing_planes_on_point_set_3.cpp b/Shape_detection/examples/Shape_detection/region_growing_planes_on_point_set_3.cpp index 8d7f719aec9..40c4a20fb6f 100644 --- a/Shape_detection/examples/Shape_detection/region_growing_planes_on_point_set_3.cpp +++ b/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(); diff --git a/Shape_detection/examples/Shape_detection/region_growing_planes_on_polygon_mesh.cpp b/Shape_detection/examples/Shape_detection/region_growing_planes_on_polygon_mesh.cpp index ef2de51a703..414eb2d4a7a 100644 --- a/Shape_detection/examples/Shape_detection/region_growing_planes_on_polygon_mesh.cpp +++ b/Shape_detection/examples/Shape_detection/region_growing_planes_on_polygon_mesh.cpp @@ -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 > 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; diff --git a/Shape_detection/examples/Shape_detection/region_growing_spheres_on_point_set_3.cpp b/Shape_detection/examples/Shape_detection/region_growing_spheres_on_point_set_3.cpp index 3badba367e1..f5293a6c7a0 100644 --- a/Shape_detection/examples/Shape_detection/region_growing_spheres_on_point_set_3.cpp +++ b/Shape_detection/examples/Shape_detection/region_growing_spheres_on_point_set_3.cpp @@ -1,103 +1,99 @@ -#include -#include -#include - +#include "include/utils.h" #include #include - +#include #include #include - #include -#include - -using Kernel = CGAL::Simple_cartesian; -using Point_3 = Kernel::Point_3; +// Typedefs. +using Kernel = CGAL::Simple_cartesian; +using FT = Kernel::FT; +using Point_3 = Kernel::Point_3; using Vector_3 = Kernel::Vector_3; -using Point_set = CGAL::Point_set_3; -using Point_map = typename Point_set::Point_map; +using Point_set = CGAL::Point_set_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; +using Region_type = CGAL::Shape_detection::Point_set::Least_squares_sphere_fit_region; +using Region_growing = CGAL::Shape_detection::Region_growing; -using Neighbor_query = Shape_detection::K_neighbor_query - ; -using Region_type = Shape_detection::Least_squares_sphere_fit_region - ; -using Region_growing = CGAL::Shape_detection::Region_growing - ; +int main(int argc, char** argv) { -int main (int argc, char** argv) -{ - std::ifstream ifile (argc > 1 ? argv[1] : "data/spheres.ply"); - Point_set points; - ifile >> points; + // Load ply data either from a local folder or a user-provided file. + const bool is_default_input = argc > 1 ? false : true; + std::ifstream in(is_default_input ? "data/spheres.ply" : argv[1]); - 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 - assert (points.has_normal_map()); + Point_set point_set; + 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 - const std::size_t k = 12; - const double max_distance = 0.01; - const double max_angle = 10.; + // Default parameter values for the data file spheres.ply. + const std::size_t k = 12; + const FT max_distance = FT(1) / FT(100); + 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:: - k_neighbors(k). - point_map(points.point_map())); + point_set, CGAL::parameters::k_neighbors(k).point_map(point_set.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()). - normal_map(points.normal_map())); + point_map(point_set.point_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 - red = points.add_property_map("red", 0).first, - green = points.add_property_map("green", 0).first, - blue = points.add_property_map("blue", 0).first; + red = point_set.add_property_map("red" , 0).first, + green = point_set.add_property_map("green", 0).first, + blue = point_set.add_property_map("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& region) { - std::size_t nb_spheres = 0; - CGAL::Real_timer timer; - timer.start(); - region_growing.detect - (boost::make_function_output_iterator - ([&](const std::vector& region) - { - // Assign a random color to each region - unsigned char r = static_cast(random.get_int(64, 192)); - unsigned char g = static_cast(random.get_int(64, 192)); - unsigned char b = static_cast(random.get_int(64, 192)); - for (const std::size_t& idx : region) - { - red[idx] = r; + // Assign a random color to each region. + const unsigned char r = static_cast(random.get_int(64, 192)); + const unsigned char g = static_cast(random.get_int(64, 192)); + const unsigned char b = static_cast(random.get_int(64, 192)); + for (const std::size_t idx : region) { + red[idx] = r; green[idx] = g; - blue[idx] = b; + blue[idx] = b; } - ++ nb_spheres; - })); - timer.stop(); - - std::cerr << nb_spheres << " spheres detected in " - << timer.time() << " seconds" << std::endl; - - // Save in colored_spheres.ply - std::ofstream out ("colored_spheres.ply"); - CGAL::IO::set_binary_mode (out); - out << points; + ++num_spheres; + } + ) + ); + std::cout << "* number of found spheres: " << num_spheres << std::endl; + assert(is_default_input && num_spheres == 0); + // 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; } diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_region.h index 76b19110fef..ba196b6eae9 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_region.h +++ b/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 - if (indices.size() < 6) + // computed circle means something. + 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 (m_radius < m_min_radius || m_radius > m_max_radius) + // If radius is out of bound, nothing fits, early ending. + 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,24 +324,25 @@ 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); - FT distance_to_circle = CGAL::abs (distance_to_center - m_radius); + const FT distance_to_center = m_sqrt (sq_dist); + 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; normal = normal / m_sqrt (sq_norm); - Vector_2 ray (m_center, query_point); + Vector_2 ray(m_center, 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; } @@ -375,17 +379,14 @@ namespace Point_set { bool update(const std::vector& region) { CGAL_precondition(region.size() > 0); - auto unary_function - = [&](const std::size_t& idx) -> const Point_2& - { - return get (m_point_map, *(m_input_range.begin() + idx)); - }; + auto unary_function = [&](const std::size_t& idx) -> const Point_2& { + return get(m_point_map, *(m_input_range.begin() + idx)); + }; - internal::circle_fit - (make_range(boost::make_transform_iterator - (region.begin(), unary_function), - boost::make_transform_iterator - (region.end(), unary_function)), + internal::create_circle_2( + make_range( + boost::make_transform_iterator(region.begin(), unary_function), + boost::make_transform_iterator(region.end(), unary_function)), m_sqrt, m_squared_distance_2, m_center, m_radius); return true; } diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_sorting.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_sorting.h index d4facdfb28b..567e0f67ffc 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_sorting.h +++ b/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)) - { - // Score is min squared distance to sphere + if (internal::create_circle_2( + points, sqrt, squared_distance_2, fitted_center, fitted_radius)) { + + // Score is min squared distance to sphere. m_scores[i] = Local_FT(0); - for (const Local_point_2& p : points) - m_scores[i] += abs (sqrt(squared_distance_2(p, fitted_center)) - fitted_radius); + for (const Local_point_2& p : points) { + m_scores[i] += abs(sqrt(squared_distance_2(p, fitted_center)) - fitted_radius); + } + } else { + m_scores[i] = Local_FT(std::numeric_limits::max()); } - else - m_scores[i] = Local_FT(std::numeric_limits::infinity()); } } }; diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_region.h index acd5d452464..f99917467fc 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_region.h +++ b/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 - if (indices.size() < 6) + // computed cylinder means something. + 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 (m_radius < m_min_radius || m_radius > m_max_radius) + // If radius is out of bound, nothing fits, early ending. + 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); - FT distance_to_cylinder = CGAL::abs (distance_to_center - m_radius); + const FT distance_to_center = m_sqrt(sq_dist); + 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); - Vector_3 ray (proj, 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); + 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 - std::vector& aregion - = const_cast&>(region); - cpp98::random_shuffle (aregion.begin(), aregion.end()); + // Shuffle to avoid always picking 2 close points. + std::vector& aregion = + const_cast&>(region); + cpp98::random_shuffle(aregion.begin(), aregion.end()); using VT = typename std::iterator_traits::value_type; - auto unary_function - = [&](const std::size_t& idx) -> VT - { - return *(m_input_range.begin() + idx); - }; + auto unary_function = [&](const std::size_t& idx) -> VT { + return *(m_input_range.begin() + idx); + }; - internal::cylinder_fit - (make_range (boost::make_transform_iterator - (region.begin(), unary_function), - boost::make_transform_iterator - (region.end(), unary_function)), - m_point_map, m_normal_map, m_sqrt, m_squared_distance_3, - m_axis, m_radius); + internal::create_cylinder_3( + make_range( + boost::make_transform_iterator(region.begin(), unary_function), + boost::make_transform_iterator(region.end(), unary_function)), + m_point_map, m_normal_map, m_sqrt, m_squared_distance_3, m_axis, m_radius); return true; } diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_sorting.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_sorting.h index 302074274a1..d3e11973dda 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_sorting.h +++ b/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::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(), - sqrt, squared_distance_3, fitted_line, fitted_radius)) - { - // Score is min squared distance to cylinder + if (internal::create_cylinder_3( + points, Local_point_map(), Local_normal_map(), + sqrt, squared_distance_3, fitted_line, fitted_radius)) { + + // Score is min squared distance to cylinder. m_scores[i] = Local_FT(0); - for (const Local_pwn& pwn : points) - m_scores[i] += abs (sqrt(squared_distance_3(pwn.first, fitted_line)) - fitted_radius); + for (const Local_pwn& pwn : points) { + m_scores[i] += abs(sqrt(squared_distance_3(pwn.first, fitted_line)) - fitted_radius); + } + } else { + m_scores[i] = Local_FT(std::numeric_limits::max()); } - else - m_scores[i] = Local_FT(std::numeric_limits::infinity()); } } }; diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_region.h index 70da23f88b2..a2dd3652833 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_region.h +++ b/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 - if (indices.size() < 6) + // computed sphere means something. + 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 (m_radius < m_min_radius || m_radius > m_max_radius) + // If radius is out of bound, nothing fits, early ending. + 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,24 +329,25 @@ 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); - FT distance_to_sphere = CGAL::abs (distance_to_center - m_radius); + const FT distance_to_center = m_sqrt(sq_dist); + 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; normal = normal / m_sqrt (sq_norm); - Vector_3 ray (m_center, query_point); + Vector_3 ray(m_center, 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; } @@ -380,17 +384,14 @@ namespace Point_set { bool update(const std::vector& region) { CGAL_precondition(region.size() > 0); - auto unary_function - = [&](const std::size_t& idx) -> const Point_3& - { - return get (m_point_map, *(m_input_range.begin() + idx)); - }; + auto unary_function = [&](const std::size_t& idx) -> const Point_3& { + return get (m_point_map, *(m_input_range.begin() + idx)); + }; - internal::sphere_fit - (make_range (boost::make_transform_iterator - (region.begin(), unary_function), - boost::make_transform_iterator - (region.end(), unary_function)), + internal::create_sphere_3( + make_range( + boost::make_transform_iterator(region.begin(), unary_function), + boost::make_transform_iterator(region.end(), unary_function)), m_sqrt, m_squared_distance_3, m_center, m_radius); return true; } diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_sorting.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_sorting.h index 317ca38d617..bad9a8f605c 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_sorting.h +++ b/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::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)) - { - // Score is min squared distance to sphere + if (internal::create_sphere_3( + points, sqrt, squared_distance_3, fitted_center, fitted_radius)) { + + // Score is min squared distance to sphere. m_scores[i] = Local_FT(0); - for (const Local_point_3& p : points) - m_scores[i] += abs (sqrt(squared_distance_3(p, fitted_center)) - fitted_radius); + for (const Local_point_3& p : points) { + m_scores[i] += abs(sqrt(squared_distance_3(p, fitted_center)) - fitted_radius); + } + } else { + m_scores[i] = Local_FT(std::numeric_limits::max()); } - else - m_scores[i] = Local_FT(std::numeric_limits::infinity()); } } }; diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h index ef83fb846ba..8815c068844 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h +++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h @@ -346,240 +346,272 @@ namespace internal { return std::make_pair(plane, static_cast(score)); } -template -bool circle_fit (const PointRange& points, - const Sqrt& sqrt, - const Squared_distance_2& squared_distance_2, - Point_2& center, FT& radius) -{ - using Diagonalize_traits = Eigen_diagonalize_traits; - using Covariance_matrix = typename Diagonalize_traits::Covariance_matrix; - using Vector = typename Diagonalize_traits::Vector; - using Matrix = typename Diagonalize_traits::Matrix; + template< + typename PointRange, + typename Sqrt, + 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) { - // Use bbox to compute diagonalization with smaller coordinates to - // avoid loss of precision when inverting large coordinates - Bbox_2 bbox = bbox_2 (points.begin(), points.end()); + using Diagonalize_traits = Eigen_diagonalize_traits; + using Covariance_matrix = typename Diagonalize_traits::Covariance_matrix; - // Circle least squares fitting, - // 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(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 ri = x^2 + y^2 - // -> Minimize Sum(D + B*xi + C*yi + ri)^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 = -0.5 * eigenvector(1) / eigenvector(3) ; -0.5 * eigenvector(2) / eigenvector(3) - Covariance_matrix A - = { FT(0), FT(0), FT(0), FT(0), FT(0), - FT(0), FT(0), FT(0), FT(0), FT(0) }; + using Vector = typename Diagonalize_traits::Vector; + using Matrix = typename Diagonalize_traits::Matrix; - A[0] = static_cast(points.size()); - for (const Point_2& p : points) - { - FT x = p.x() - bbox.xmin(); - FT y = p.y() - bbox.ymin(); - FT r = x*x + y*y; - A[1] += x; - A[2] += y; - A[3] += r; - A[4] += x * x; - A[5] += x * y; - A[6] += x * r; - A[7] += y * y; - A[8] += y * r; - A[9] += r * r; + // Use bbox to compute diagonalization with smaller coordinates to + // avoid loss of precision when inverting large coordinates. + const Bbox_2 bbox = bbox_2(points.begin(), points.end()); + + // Circle least squares fitting: + // 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(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; and ri = x^2 + y^2. + // -> Minimize sum(D + B * xi + C * yi + ri)^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 = [ + // -0.5 * eigenvector(1) / eigenvector(3); + // -0.5 * eigenvector(2) / eigenvector(3); ] + Covariance_matrix A = { + FT(0), FT(0), FT(0), FT(0), FT(0), + FT(0), FT(0), FT(0), FT(0), FT(0) + }; + + A[0] = static_cast(points.size()); + for (const Point_2& p : points) { + + const FT x = p.x() - bbox.xmin(); + const FT y = p.y() - bbox.ymin(); + const FT r = x * x + y * y; + + A[1] += x; + A[2] += y; + A[3] += r; + A[4] += x * x; + A[5] += x * y; + A[6] += x * r; + A[7] += y * y; + A[8] += y * r; + A[9] += r * r; + } + + Vector eigenvalues = { + 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) + }; + Diagonalize_traits:: + diagonalize_selfadjoint_covariance_matrix(A, eigenvalues, eigenvectors); + + // Perfect line case, no circle can be fitted. + if (eigenvectors[3] == 0) { + return false; + } + + // Other cases. + 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) { + radius += sqrt(squared_distance_2(p, center)); + } + radius /= static_cast(points.size()); + return true; } - Vector eigenvalues = { 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) }; + template< + typename PointRange, + typename Sqrt, + 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) { - Diagonalize_traits::diagonalize_selfadjoint_covariance_matrix - (A, eigenvalues, eigenvectors); + using Diagonalize_traits = Eigen_diagonalize_traits; + using Covariance_matrix = typename Diagonalize_traits::Covariance_matrix; - // Perfect line case, no circle can be fitted - if (eigenvectors[3] == 0) - return false; + using Vector = typename Diagonalize_traits::Vector; + using Matrix = typename Diagonalize_traits::Matrix; - center = Point_2 (bbox.xmin() - FT(0.5) * (eigenvectors[1] / eigenvectors[3]), - bbox.ymin() - FT(0.5) * (eigenvectors[2] / eigenvectors[3])); + // Use bbox to compute diagonalization with smaller coordinates to + // avoid loss of precision when inverting large coordinates. + const Bbox_3 bbox = bbox_3(points.begin(), points.end()); - radius = FT(0); - for (const Point_2& p : points) - radius += sqrt (squared_distance_2 (p, center)); - radius /= points.size(); + // Sphere least squares fitting. + // See create_circle_2() above for more details on computation. + 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) + }; - return true; -} + A[0] = static_cast(points.size()); + for (const Point_3& p : points) { -template -bool sphere_fit (const PointRange& points, - const Sqrt& sqrt, - const Squared_distance_3& squared_distance_3, - Point_3& center, FT& radius) -{ - using Diagonalize_traits = Eigen_diagonalize_traits; - using Covariance_matrix = typename Diagonalize_traits::Covariance_matrix; - using Vector = typename Diagonalize_traits::Vector; - using Matrix = typename Diagonalize_traits::Matrix; + const FT x = p.x() - bbox.xmin(); + const FT y = p.y() - bbox.ymin(); + const FT z = p.z() - bbox.zmin(); + const FT r = x * x + y * y + z * z; - // Use bbox to compute diagonalization with smaller coordinates to - // avoid loss of precision when inverting large coordinates - Bbox_3 bbox = bbox_3 (points.begin(), points.end()); + A[1] += x; + A[2] += y; + A[3] += z; + A[4] += r; + A[5] += x * x; + A[6] += x * y; + A[7] += x * z; + A[8] += x * r; + A[9] += y * y; + A[10] += y * z; + A[11] += y * r; + A[12] += z * z; + A[13] += z * r; + A[14] += r * r; + } - // Sphere least squares fitting - // (see Least_square_circle_fit_region for details about computation) - 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) }; + Vector eigenvalues = { + 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) + }; + Diagonalize_traits:: + diagonalize_selfadjoint_covariance_matrix(A, eigenvalues, eigenvectors); - A[0] = static_cast(points.size()); - for (const Point_3& p : points) - { - FT x = p.x() - bbox.xmin(); - FT y = p.y() - bbox.ymin(); - FT z = p.z() - bbox.zmin(); - FT r = x*x + y*y + z*z; - A[1] += x; - A[2] += y; - A[3] += z; - A[4] += r; - A[5] += x * x; - A[6] += x * y; - A[7] += x * z; - A[8] += x * r; - A[9] += y * y; - A[10] += y * z; - A[11] += y * r; - A[12] += z * z; - A[13] += z * r; - A[14] += r * r; + // Perfect plane case, no sphere can be fitted. + if (eigenvectors[4] == 0) { + return false; + } + + // Other cases. + const FT half = FT(1) / FT(2); + 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) { + radius += sqrt(squared_distance_3(p, center)); + } + radius /= static_cast(points.size()); + return true; } - Vector eigenvalues = { 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) }; + template< + typename PointRange, + typename PointMap, + typename NormalMap, + 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) { - Diagonalize_traits::diagonalize_selfadjoint_covariance_matrix - (A, eigenvalues, eigenvectors); + using Point_3 = typename boost::property_traits::value_type; + using Vector_3 = typename boost::property_traits::value_type; - // Perfect plane case, no sphere can be fitted - if (eigenvectors[4] == 0) - return false; + const Point_3& ref = get(point_map, *(points.begin())); - center = Point_3 (bbox.xmin() - FT(0.5) * (eigenvectors[1] / eigenvectors[4]), - bbox.ymin() - FT(0.5) * (eigenvectors[2] / eigenvectors[4]), - bbox.zmin() - FT(0.5) * (eigenvectors[3] / eigenvectors[4])); + radius = FT(0); + std::size_t nb = 0; + Vector_3 mean_axis = CGAL::NULL_VECTOR; + Point_3 point_on_axis = ORIGIN; - radius = FT(0); - for (const Point_3& p : points) - radius += sqrt (squared_distance_3 (p, center)); - radius /= points.size(); + for (std::size_t i = 0; i < points.size() - 1; ++i) { - return true; -} + 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(1) / FT(100)) { + continue; + } + axis = axis / sqrt(axis * axis); -template -bool cylinder_fit (const PointRange& points, - PointMap point_map, NormalMap normal_map, - const Sqrt& sqrt, - const Squared_distance_3& squared_distance_3, - Line_3& axis, FT& radius) -{ - using Point_3 = typename boost::property_traits::value_type; - using Vector_3 = typename boost::property_traits::value_type; + const Point_3& p0 = get(point_map, *std::next(points.begin(), i)); + const Point_3& p1 = get(point_map, *std::next(points.begin(), i + 1)); - const Point_3& ref = get(point_map, *(points.begin())); + Vector_3 xdir = v0 - axis * (v0 * axis); + xdir = xdir / sqrt(xdir * xdir); - Vector_3 mean_axis = CGAL::NULL_VECTOR; - std::size_t nb = 0; - radius = FT(0); - Point_3 point_on_axis = ORIGIN; - 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)) - continue; - axis = axis / sqrt(axis * axis); + Vector_3 ydir = cross_product(axis, xdir); + ydir = ydir / sqrt (ydir * ydir); - const Point_3& p0 = get(point_map, *std::next(points.begin(), i)); - const Point_3& p1 = get(point_map, *std::next(points.begin(), i + 1)); + const FT v1x = v1 * ydir; + const FT v1y = -v1 * xdir; - Vector_3 xdir = v0 - axis * (v0 * axis); - xdir = xdir / sqrt (xdir * xdir); + Vector_3 d(p0, p1); + const FT ox = xdir * d; + const FT oy = ydir * d; + const FT ldist = v1x * ox + v1y * oy; - Vector_3 ydir = cross_product (axis, xdir); - ydir = ydir / sqrt (ydir * ydir); + FT r = ldist / v1x; + Point_3 point = p0 + xdir * r; + const Line_3 line(point, axis); + point = line.projection(ref); + point_on_axis = barycenter(point_on_axis, static_cast(nb), point, FT(1)); + r += abs(r); - FT v1x = v1 * ydir; - FT v1y = -v1 * xdir; + if (nb != 0 && axis * mean_axis < 0) { + axis = -axis; + } + mean_axis = mean_axis + axis; + ++nb; + } - Vector_3 d (p0, p1); + if (nb == 0) { + return false; + } - FT ox = xdir * d; - FT oy = ydir * d; - FT ldist = v1x * ox + v1y * oy; + mean_axis = mean_axis / sqrt(mean_axis * mean_axis); + axis = Line_3(point_on_axis, mean_axis); + radius /= static_cast(nb); - 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(nb), point, FT(1)); - - radius += abs(radius); - - if (nb != 0 && axis * mean_axis < 0) - axis = -axis; - - mean_axis = mean_axis + axis; - ++ nb; + radius = FT(0); + for (const auto& p : points) { + const Point_3& p0 = get(point_map, p); + radius += sqrt(squared_distance_3(p0, axis)); + } + radius /= static_cast(points.size()); + return true; } - 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 = FT(0); - for (const auto& p : points) - { - const Point_3& p0 = get(point_map, p); - radius += sqrt(squared_distance_3(p0, axis)); - } - radius /= points.size(); - - return true; -} - } // namespace internal } // namespace Shape_detection } // namespace CGAL diff --git a/Shape_detection/package_info/Shape_detection/long_description.txt b/Shape_detection/package_info/Shape_detection/long_description.txt index e49e558f3d3..b9cf7d3c5bc 100644 --- a/Shape_detection/package_info/Shape_detection/long_description.txt +++ b/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, -detecting planes in a 3D point set, detecting planes on a polygon mesh, detecting lines +CGAL provides several instances of the algorithm: detecting lines and circles in a 2D point set, +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. diff --git a/Shape_detection/test/Shape_detection/test_region_growing_on_point_set_2_with_sorting.cpp b/Shape_detection/test/Shape_detection/test_region_growing_on_point_set_2_with_sorting.cpp index 7d7674650c6..fd140f316a1 100644 --- a/Shape_detection/test/Shape_detection/test_region_growing_on_point_set_2_with_sorting.cpp +++ b/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(argc, argv, 62, 66) // && - // test(argc, argv, 196, 200) + test(argc, argv, 62, 66) && + test(argc, argv, 196, 200) ) ? EXIT_SUCCESS : EXIT_FAILURE ); } diff --git a/Shape_detection/test/Shape_detection/test_region_growing_on_point_set_3_with_sorting.cpp b/Shape_detection/test/Shape_detection/test_region_growing_on_point_set_3_with_sorting.cpp index bac9f101e26..9baede56975 100644 --- a/Shape_detection/test/Shape_detection/test_region_growing_on_point_set_3_with_sorting.cpp +++ b/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 (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::infinity(); + const double max_radius = std::numeric_limits::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::infinity(); + const double max_radius = std::numeric_limits::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; }