From 60f85259cdab211688c48c03b7c64383ef6e9751 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 30 Nov 2020 11:16:39 +0100 Subject: [PATCH 01/10] Add validity test to shape detection --- .../test/Shape_detection/CMakeLists.txt | 2 + .../test/Shape_detection/test_validity.cpp | 196 ++++++++++++++++++ 2 files changed, 198 insertions(+) create mode 100644 Shape_detection/test/Shape_detection/test_validity.cpp diff --git a/Shape_detection/test/Shape_detection/CMakeLists.txt b/Shape_detection/test/Shape_detection/CMakeLists.txt index 9e6719e15ce..212d2ad86fe 100644 --- a/Shape_detection/test/Shape_detection/CMakeLists.txt +++ b/Shape_detection/test/Shape_detection/CMakeLists.txt @@ -15,6 +15,7 @@ include(CGAL_CreateSingleSourceCGALProgram) find_package(Eigen3 3.1.0 QUIET) # (3.1.0 or greater) include(CGAL_Eigen_support) if(EIGEN3_FOUND) + create_single_source_cgal_program("test_validity.cpp") create_single_source_cgal_program("test_region_growing_basic.cpp") create_single_source_cgal_program("test_region_growing_on_cube.cpp") create_single_source_cgal_program("test_region_growing_on_point_set_2.cpp") @@ -30,6 +31,7 @@ if(EIGEN3_FOUND) "test_region_growing_on_degenerated_mesh.cpp") foreach( target + test_validity test_region_growing_basic test_region_growing_on_cube test_region_growing_on_point_set_2 diff --git a/Shape_detection/test/Shape_detection/test_validity.cpp b/Shape_detection/test/Shape_detection/test_validity.cpp new file mode 100644 index 00000000000..b801fb53b91 --- /dev/null +++ b/Shape_detection/test/Shape_detection/test_validity.cpp @@ -0,0 +1,196 @@ +#include + +#include +#include +#include + +#include + +#include + +#include + +namespace SD = CGAL::Shape_detection; + +using Kernel = CGAL::Simple_cartesian; +using Point_3 = Kernel::Point_3; +using Vector_3 = Kernel::Vector_3; + +using Pwn = std::pair; +using Point_set = std::vector; +using Point_map = CGAL::First_of_pair_property_map; +using Normal_map = CGAL::Second_of_pair_property_map; + +using RG_query = SD::Point_set::Sphere_neighbor_query; +using RG_region = SD::Point_set::Least_squares_plane_fit_region; +using Region_growing = SD::Region_growing; + +using RANSAC_traits = SD::Efficient_RANSAC_traits; +using RANSAC = SD::Efficient_RANSAC; +using RANSAC_plane = SD::Plane; + +void test_random_planes(std::size_t nb_planes); + +int main() +{ + test_random_planes(1); + test_random_planes(10); + test_random_planes(100); + test_random_planes(1000); + test_random_planes(10000); + + return EXIT_SUCCESS; +} + +void test_random_planes(std::size_t nb_planes) +{ + CGAL::Random random; + std::cerr << "[TEST ON " << nb_planes << " RANDOM PLANES, SEED = " + << random.get_seed() << "] "; + + std::vector time_rg; + std::vector nb_detected_rg; + std::vector nb_unassigned_rg; + std::vector time_ransac; + std::vector nb_detected_ransac; + std::vector nb_unassigned_ransac; + + CGAL::Real_timer timer; + double timeout = 120.; // 2 minutes timeout + + timer.start(); + std::size_t nb_runs = 500; + for (std::size_t run = 0; run < nb_runs; ++ run) + { + std::size_t min_points = random.get_int(10, 200); + std::size_t max_points = random.get_int(200, 10000); + double cluster_epsilon = random.get_double(0.01, 10.); + double epsilon = random.get_double (0., cluster_epsilon / 10.); + double normal_threshold = random.get_double (0.75, 0.99); + + double domain_size = cluster_epsilon * std::sqrt(nb_planes) * min_points; + double spacing = 0.8 * cluster_epsilon / std::sqrt(2); // smaller than diagonal + double noise = 0.5 * epsilon; + + Point_set points; + for (std::size_t i = 0; i < nb_planes; ++ i) + { + // Generate random plane + Point_3 origin (random.get_double(-domain_size, domain_size), + random.get_double(-domain_size, domain_size), + random.get_double(-domain_size, domain_size)); + Vector_3 base1 (random.get_double(-1,1), random.get_double(-1,1), random.get_double(-1,1)); + base1 = base1 / std::sqrt(base1 * base1); + Vector_3 base2 (random.get_double(-1,1), random.get_double(-1,1), random.get_double(-1,1)); + base2 = base2 / std::sqrt(base2 * base2); + + Vector_3 normal = CGAL::cross_product (base1, base2); + normal = normal / std::sqrt (normal * normal); + + std::size_t nb_points = random.get_int (min_points, max_points); + + std::size_t nb_x = std::size_t(std::sqrt (double(nb_points))); + std::size_t nb_y = std::size_t(nb_points / double(nb_x)) + 1; + + for (std::size_t j = 0; j < nb_x; ++ j) + for (std::size_t k = 0; k < nb_y; ++ k) + points.emplace_back (origin + j * spacing * base1 + k * spacing * base2 + + normal * random.get_double(-noise/2, noise/2), + normal); + } + + // std::ofstream test("dump.pwn"); + // for (const auto& p : points) + // test << p.first << " " << p.second << std::endl; + CGAL::Real_timer t; + t.start(); + RG_query rg_query (points, cluster_epsilon); + RG_region rg_region (points, epsilon, normal_threshold, min_points); + Region_growing region_growing (points, rg_query, rg_region); + std::size_t nb_detected = 0; + std::size_t nb_unassigned = 0; + region_growing.detect (boost::make_function_output_iterator ([&](const auto&) { ++ nb_detected; })); + region_growing.unassigned_items (boost::make_function_output_iterator ([&](const auto&) { ++ nb_unassigned; })); + t.stop(); + + time_rg.push_back (t.time()); + nb_detected_rg.push_back (nb_detected); + nb_unassigned_rg.push_back (nb_unassigned / double(points.size())); + t.reset(); + + t.start(); + RANSAC ransac; + ransac.template add_shape_factory(); + ransac.set_input(points); + typename RANSAC::Parameters parameters; + parameters.probability = 0.05f; + parameters.min_points = min_points; + parameters.epsilon = epsilon; + parameters.cluster_epsilon = cluster_epsilon; + parameters.normal_threshold = normal_threshold; + ransac.detect(parameters); + t.stop(); + + time_ransac.push_back (t.time()); + nb_detected_ransac.push_back (ransac.shapes().size()); + nb_unassigned_ransac.push_back (ransac.number_of_unassigned_points() / double(points.size())); + +#if 0 + if (ransac.shapes().size() == 0) + { + std::cerr << "Detected 0 shapes with " << std::endl + << " * min points = " << min_points << std::endl + << " * epsilon = " << epsilon << std::endl + << " * cluster_epsilon = " << cluster_epsilon << std::endl + << " * normal_threshold = " << normal_threshold << std::endl; + + std::ofstream ofile("0shapes.ply", std::ios::binary); +// CGAL::set_binary_mode (ofile); + CGAL::write_ply_points (ofile, points, + CGAL::parameters::point_map(Point_map()). + normal_map(Normal_map())); + + exit(0); + } +#endif + + if (timer.time() > timeout) + { + nb_runs = run + 1; + break; + } + } + + std::cerr << "on " << nb_runs << " runs" << std::endl; + + + std::sort (time_ransac.begin(), time_ransac.end()); + std::sort (nb_detected_ransac.begin(), nb_detected_ransac.end()); + std::sort (nb_unassigned_ransac.begin(), nb_unassigned_ransac.end()); + std::sort (time_rg.begin(), time_rg.end()); + std::sort (nb_detected_rg.begin(), nb_detected_rg.end()); + std::sort (nb_unassigned_rg.begin(), nb_unassigned_rg.end()); + + std::cerr << " * Region Growing" << std::endl + << " - took between " << time_rg.front() << "s and " << time_rg.back() + << "s, median = " << time_rg[nb_runs / 2] << "s" << std::endl + << " - detected between " << nb_detected_rg.front() << " and " + << nb_detected_rg.back() << " planes (" << 100. * nb_detected_rg.front() / double(nb_planes) + << "% to " << 100. * nb_detected_rg.back() / double(nb_planes) << "%), median = " + << nb_detected_rg[nb_runs / 2] << " planes (" << 100. * nb_detected_rg[nb_runs / 2] / double(nb_planes) + << "%)" << std::endl + << " - left between " << 100. * nb_unassigned_rg.front() + << "% and " << 100. * nb_unassigned_rg.back() << "% of unassigned points, median = " + << 100. * nb_unassigned_rg[nb_runs / 2] << "%" << std::endl; + std::cerr << " * Efficient RANSAC" << std::endl + << " - took between " << time_ransac.front() << "s and " << time_ransac.back() + << "s, median = " << time_ransac[nb_runs / 2] << "s" << std::endl + << " - detected between " << nb_detected_ransac.front() << " and " + << nb_detected_ransac.back() << " planes (" << 100. * nb_detected_ransac.front() / double(nb_planes) + << "% to " << 100. * nb_detected_ransac.back() / double(nb_planes) << "%), median = " + << nb_detected_ransac[nb_runs / 2] << " planes (" << 100. * nb_detected_ransac[nb_runs / 2] / double(nb_planes) + << "%)" << std::endl + << " - left between " << 100. * nb_unassigned_ransac.front() + << "% and " << 100. * nb_unassigned_ransac.back() << "% of unassigned points, median = " + << 100. * nb_unassigned_ransac[nb_runs / 2] << "%" << std::endl; +} From 37d38da19014b72f7b76ddec67539ebb65b86571 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 1 Dec 2020 10:52:18 +0100 Subject: [PATCH 02/10] Add validity tests with sampled data --- .../test/Shape_detection/CMakeLists.txt | 15 ++ .../test/Shape_detection/test_validity.cpp | 21 +- .../test_validity_sampled_data.cpp | 228 ++++++++++++++++++ 3 files changed, 254 insertions(+), 10 deletions(-) create mode 100644 Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp diff --git a/Shape_detection/test/Shape_detection/CMakeLists.txt b/Shape_detection/test/Shape_detection/CMakeLists.txt index 212d2ad86fe..d5a4eca990a 100644 --- a/Shape_detection/test/Shape_detection/CMakeLists.txt +++ b/Shape_detection/test/Shape_detection/CMakeLists.txt @@ -44,6 +44,21 @@ if(EIGEN3_FOUND) target_link_libraries(${target} PUBLIC CGAL::Eigen_support) endforeach() + set(RANSAC_PROTO_DIR CACHE PATH "") + if (NOT RANSAC_PROTO_DIR STREQUAL "") + add_definitions(-DPOINTSWITHINDEX -DCGAL_TEST_RANSAC_PROTOTYPE) + include_directories(${RANSAC_PROTO_DIR}) + include_directories(${RANSAC_PROTO_DIR}/MiscLib/) + file(GLOB proto_src "${RANSAC_PROTO_DIR}/*.cpp") + file(GLOB proto_misc_src "${RANSAC_PROTO_DIR}/MiscLib/*.cpp") + add_library(libproto STATIC ${proto_src} ${proto_misc_src}) + add_executable(test_validity_sampled_data "test_validity_sampled_data.cpp") + target_link_libraries(test_validity_sampled_data libproto CGAL::CGAL CGAL::Eigen_support) + else() + add_executable(test_validity_sampled_data "test_validity_sampled_data.cpp") + target_link_libraries(test_validity_sampled_data CGAL::CGAL CGAL::Eigen_support) + endif() + endif() create_single_source_cgal_program( diff --git a/Shape_detection/test/Shape_detection/test_validity.cpp b/Shape_detection/test/Shape_detection/test_validity.cpp index b801fb53b91..3038082d850 100644 --- a/Shape_detection/test/Shape_detection/test_validity.cpp +++ b/Shape_detection/test/Shape_detection/test_validity.cpp @@ -1,5 +1,6 @@ #include +//#define CGAL_RANSAC_EXPERIMENTAL_FIXES #include #include #include @@ -33,9 +34,9 @@ void test_random_planes(std::size_t nb_planes); int main() { - test_random_planes(1); - test_random_planes(10); - test_random_planes(100); + // test_random_planes(1); + // test_random_planes(10); + // test_random_planes(100); test_random_planes(1000); test_random_planes(10000); @@ -56,20 +57,20 @@ void test_random_planes(std::size_t nb_planes) std::vector nb_unassigned_ransac; CGAL::Real_timer timer; - double timeout = 120.; // 2 minutes timeout + double timeout = 60.; // 1 minute timeout timer.start(); - std::size_t nb_runs = 500; + std::size_t nb_runs = 100; for (std::size_t run = 0; run < nb_runs; ++ run) { - std::size_t min_points = random.get_int(10, 200); - std::size_t max_points = random.get_int(200, 10000); + std::size_t min_points = random.get_int(10, 100); + std::size_t max_points = random.get_int(100, 1000); double cluster_epsilon = random.get_double(0.01, 10.); double epsilon = random.get_double (0., cluster_epsilon / 10.); double normal_threshold = random.get_double (0.75, 0.99); - double domain_size = cluster_epsilon * std::sqrt(nb_planes) * min_points; double spacing = 0.8 * cluster_epsilon / std::sqrt(2); // smaller than diagonal + double domain_size = spacing * std::sqrt(nb_planes) * std::sqrt(0.5 * (min_points + max_points)) * 0.5; double noise = 0.5 * epsilon; Point_set points; @@ -145,11 +146,11 @@ void test_random_planes(std::size_t nb_planes) << " * normal_threshold = " << normal_threshold << std::endl; std::ofstream ofile("0shapes.ply", std::ios::binary); -// CGAL::set_binary_mode (ofile); + CGAL::set_binary_mode (ofile); CGAL::write_ply_points (ofile, points, CGAL::parameters::point_map(Point_map()). normal_map(Normal_map())); - + ofile.close(); exit(0); } #endif diff --git a/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp b/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp new file mode 100644 index 00000000000..a194e9f43e9 --- /dev/null +++ b/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp @@ -0,0 +1,228 @@ +#include + +#define CGAL_RANSAC_EXPERIMENTAL_FIXES +#include +#include +#include + +#ifdef CGAL_TEST_RANSAC_PROTOTYPE +#include +#include +#include +#include +#endif + +#include +#include + +#include + +#include + +namespace SD = CGAL::Shape_detection; + +using Kernel = CGAL::Simple_cartesian; +using Point_3 = Kernel::Point_3; +using Vector_3 = Kernel::Vector_3; + +using Pwn = std::pair; +using Point_set = std::vector; +using Point_map = CGAL::First_of_pair_property_map; +using Normal_map = CGAL::Second_of_pair_property_map; + +using RG_query = SD::Point_set::Sphere_neighbor_query; +using RG_region = SD::Point_set::Least_squares_plane_fit_region; +using Region_growing = SD::Region_growing; + +using RANSAC_traits = SD::Efficient_RANSAC_traits; +using RANSAC = SD::Efficient_RANSAC; +using RANSAC_plane = SD::Plane; + +void test_copied_point_cloud (const Point_set& points, std::size_t nb); + +int main (int argc, char** argv) +{ + Point_set points; + const char* ifilename = (argc > 1) ? argv[1] : "data/point_set_3.xyz"; + std::ifstream ifile(ifilename); + + if (!ifile || + !CGAL::read_xyz_points( + ifile, + std::back_inserter(points), + CGAL::parameters::point_map(Point_map()). + normal_map(Normal_map()))) + { + std::cerr << "Reading error" << std::endl; + return EXIT_FAILURE; + } + + CGAL::cpp98::random_shuffle (points.begin(), points.end()); + + test_copied_point_cloud (points, 1); + test_copied_point_cloud (points, 2); + test_copied_point_cloud (points, 5); + test_copied_point_cloud (points, 10); + test_copied_point_cloud (points, 20); + test_copied_point_cloud (points, 50); + + return EXIT_SUCCESS; +} + +void test_copied_point_cloud (const Point_set& original_points, std::size_t nb) +{ + CGAL::Bbox_3 bbox = CGAL::bbox_3 + (CGAL::make_transform_iterator_from_property_map (original_points.begin(), Point_map()), + CGAL::make_transform_iterator_from_property_map (original_points.end(), Point_map())); + + std::cerr << "Ground truth = " << 6*nb*nb + 1 << " planes" << std::endl; + + Point_set points; + points.reserve (nb * nb * original_points.size()); + for (std::size_t x = 0; x < nb; ++ x) + for (std::size_t y = 0; y < nb; ++ y) + { + Vector_3 shift ((bbox.xmax() - bbox.xmin()) * x, (bbox.ymax() - bbox.ymin()) * y, 0); + for (std::size_t i = 0; i < original_points.size(); ++ i) + { + Vector_3 noise (CGAL::get_default_random().get_double(-0.0001, 0.0001), + CGAL::get_default_random().get_double(-0.0001, 0.0001), + CGAL::get_default_random().get_double(-0.0001, 0.0001)); + points.emplace_back (original_points[i].first + shift + noise, + original_points[i].second); + } + } + + bbox = CGAL::bbox_3 + (CGAL::make_transform_iterator_from_property_map (points.begin(), Point_map()), + CGAL::make_transform_iterator_from_property_map (points.end(), Point_map())); + + typename RANSAC::Parameters parameters; + parameters.probability = 0.01; + parameters.min_points = 100; + parameters.epsilon = 0.01; + parameters.cluster_epsilon = 0.01; + parameters.normal_threshold = 0.95; + + CGAL::Real_timer t; + t.start(); + RG_query rg_query (points, parameters.cluster_epsilon); + RG_region rg_region (points, parameters.epsilon, parameters.normal_threshold, parameters.min_points); + Region_growing region_growing (points, rg_query, rg_region); + std::size_t nb_detected = 0; + std::size_t nb_unassigned = 0; + region_growing.detect (boost::make_function_output_iterator ([&](const auto&) { ++ nb_detected; })); + region_growing.unassigned_items (boost::make_function_output_iterator ([&](const auto&) { ++ nb_unassigned; })); + t.stop(); + std::cerr << "Region Growing = " << nb_detected << " planes (" << 1000 * t.time() << "ms)" << std::endl; + + CGAL::Real_timer timer; + double timeout = 120; // 2 minutes timeout + timer.start(); + std::size_t nb_runs = 500; + std::vector detected_ransac; + std::vector times_ransac; + for (std::size_t run = 0; run < nb_runs; ++ run) + { + Point_set ordered_points = points; + + CGAL::Real_timer t; + t.start(); + RANSAC ransac; + ransac.template add_shape_factory(); + ransac.set_input(ordered_points); + ransac.detect(parameters); + t.stop(); + + detected_ransac.emplace_back (ransac.shapes().size()); + times_ransac.emplace_back (t.time() * 1000); + if (timer.time() > timeout) + { + nb_runs = run + 1; + break; + } + } + + std::sort (detected_ransac.begin(), detected_ransac.end()); + std::sort (times_ransac.begin(), times_ransac.end()); + std::cerr << "RANSAC = " << detected_ransac[detected_ransac.size() / 2] + << " planes (" << times_ransac[times_ransac.size() / 2] << "ms) (on " + << nb_runs << " runs, planes[" + << detected_ransac.front() << ";" << detected_ransac.back() << "], time[" + << times_ransac.front() << ";" << times_ransac.back() << "])" << std::endl; + +#ifdef CGAL_TEST_RANSAC_PROTOTYPE + { + CGAL::Real_timer timer; + double timeout = 120.; // 2 minute timeout + timer.start(); + std::size_t nb_runs = 500; + std::vector detected_ransac; + std::vector times_ransac; + for (std::size_t run = 0; run < nb_runs; ++ run) + { + PointCloud proto_points; + proto_points.reserve (points.size()); + + Point Pt; + for (std::size_t i = 0; i < points.size(); ++i) + { + Pt.pos[0] = static_cast(points[i].first.x()); + Pt.pos[1] = static_cast(points[i].first.y()); + Pt.pos[2] = static_cast(points[i].first.z()); + Pt.normal[0] = static_cast(points[i].second.x()); + Pt.normal[1] = static_cast(points[i].second.y()); + Pt.normal[2] = static_cast(points[i].second.z()); +#ifdef POINTSWITHINDEX + Pt.index = i; +#endif + proto_points.push_back(Pt); + } + + //manually set bounding box! + Vec3f cbbMin, cbbMax; + cbbMin[0] = static_cast(bbox.xmin()); + cbbMin[1] = static_cast(bbox.ymin()); + cbbMin[2] = static_cast(bbox.zmin()); + cbbMax[0] = static_cast(bbox.xmax()); + cbbMax[1] = static_cast(bbox.ymax()); + cbbMax[2] = static_cast(bbox.zmax()); + proto_points.setBBox(cbbMin, cbbMax); + + // Sets parameters for shape detection. + RansacShapeDetector::Options options; + options.m_epsilon = parameters.epsilon; + options.m_bitmapEpsilon = parameters.cluster_epsilon; + options.m_normalThresh = parameters.normal_threshold; + options.m_probability = parameters.probability; + options.m_minSupport = parameters.min_points; + + CGAL::Real_timer t; + t.start(); + RansacShapeDetector ransac (options); // the detector object + ransac.Add (new PlanePrimitiveShapeConstructor()); + MiscLib::Vector, std::size_t> > shapes; // stores the detected shapes + ransac.Detect (proto_points, 0, proto_points.size(), &shapes); + t.stop(); + + detected_ransac.emplace_back (shapes.size()); + times_ransac.emplace_back (t.time() * 1000); + if (timer.time() > timeout) + { + nb_runs = run + 1; + break; + } + } + + std::sort (detected_ransac.begin(), detected_ransac.end()); + std::sort (times_ransac.begin(), times_ransac.end()); + std::cerr << "RANSAC (proto) = " << detected_ransac[detected_ransac.size() / 2] + << " planes (" << times_ransac[times_ransac.size() / 2] << "ms) (on " + << nb_runs << " runs, planes[" + << detected_ransac.front() << ";" << detected_ransac.back() << "], time[" + << times_ransac.front() << ";" << times_ransac.back() << "])" << std::endl; + } +#endif + + std::cerr << std::endl; +} From d2d766a38e8173a726c23f47bbd97b39887e9e93 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 1 Dec 2020 11:50:14 +0100 Subject: [PATCH 03/10] Fix stop probability to take into account different shapes --- .../Efficient_RANSAC/Efficient_RANSAC.h | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h index 067dae11699..e4de8d1b03b 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h +++ b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h @@ -524,10 +524,10 @@ namespace CGAL { std::vector candidates; // Identifying minimum number of samples - std::size_t required_samples = 0; + m_required_samples = 0; for (std::size_t i = 0;i)(required_samples, tmp->minimum_sample_size()); + m_required_samples = (std::max)(m_required_samples, tmp->minimum_sample_size()); delete tmp; } @@ -569,7 +569,7 @@ namespace CGAL { select_random_octree_level(), indices, m_shape_index, - required_samples); + m_required_samples); if (callback && !callback(num_invalid / double(m_num_total_points))) return false; @@ -1019,7 +1019,8 @@ namespace CGAL { } inline FT stop_probability(std::size_t largest_candidate, std::size_t num_pts, std::size_t num_candidates, std::size_t octree_depth) const { - return (std::min)(std::pow((FT) 1.f - (FT) largest_candidate / FT(num_pts * octree_depth * 4), (int) num_candidates), (FT) 1); + return (std::min)(std::pow((FT) 1.f - (FT) largest_candidate + / (FT(num_pts) * (octree_depth+1) * (1 << (m_required_samples - 1))), (int) num_candidates), (FT) 1); } private: @@ -1039,6 +1040,7 @@ namespace CGAL { std::vector m_shape_index; std::size_t m_num_available_points; std::size_t m_num_total_points; + std::size_t m_required_samples; //give the index of the subset of point i std::vector m_index_subsets; From 4a3ad83f2181204a2465c781297258ca91746a8e Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 1 Dec 2020 11:51:33 +0100 Subject: [PATCH 04/10] Small optimization: pick several candidates at each loop --- .../Efficient_RANSAC/Efficient_RANSAC.h | 66 +++++++++++-------- 1 file changed, 37 insertions(+), 29 deletions(-) diff --git a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h index e4de8d1b03b..1da816c26a9 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h +++ b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h @@ -553,45 +553,52 @@ namespace CGAL { if (keep_searching) do { - // Generate candidates - //1. pick a point p1 randomly among available points - std::set indices; - bool done = false; - do { - do - first_sample = get_default_random()( - static_cast(m_num_available_points)); - while (m_shape_index[first_sample] != -1); + // Search (remaining_points / min_points) shapes (max 200 per iteration, min 1) + std::size_t search_number + = std::min(std::size_t(200), + std::max(std::size_t((m_num_available_points - num_invalid) / double(m_options.min_points)), + std::size_t(1))); + for (std::size_t nb = 0; nb < search_number; ++ nb) + { + // Generate candidates + //1. pick a point p1 randomly among available points + std::set indices; + bool done = false; + do { + do + first_sample = get_default_random()( + static_cast(m_num_available_points)); + while (m_shape_index[first_sample] != -1); - done = m_global_octree->drawSamplesFromCellContainingPoint( - get(m_point_pmap, - *(m_input_iterator_first + first_sample)), - select_random_octree_level(), - indices, - m_shape_index, - m_required_samples); + done = m_global_octree->drawSamplesFromCellContainingPoint( + get(m_point_pmap, + *(m_input_iterator_first + first_sample)), + select_random_octree_level(), + indices, + m_shape_index, + m_required_samples); - if (callback && !callback(num_invalid / double(m_num_total_points))) - return false; + if (callback && !callback(num_invalid / double(m_num_total_points))) + return false; - } while (m_shape_index[first_sample] != -1 || !done); + } while (m_shape_index[first_sample] != -1 || !done); - generated_candidates++; + generated_candidates++; - //add candidate for each type of primitives - for(typename std::vector::iterator it = - m_shape_factories.begin(); it != m_shape_factories.end(); it++) { + //add candidate for each type of primitives + for(typename std::vector::iterator it = + m_shape_factories.begin(); it != m_shape_factories.end(); it++) { if (callback && !callback(num_invalid / double(m_num_total_points))) return false; Shape *p = (Shape *) (*it)(); //compute the primitive and says if the candidate is valid p->compute(indices, - m_input_iterator_first, - m_traits, - m_point_pmap, - m_normal_pmap, - m_options.epsilon, - m_options.normal_threshold); + m_input_iterator_first, + m_traits, + m_point_pmap, + m_normal_pmap, + m_options.epsilon, + m_options.normal_threshold); if (p->is_valid()) { improve_bound(p, m_num_available_points - num_invalid, 1, 500); @@ -612,6 +619,7 @@ namespace CGAL { failed_candidates++; delete p; } + } } if (failed_candidates >= limit_failed_candidates) From 57f7466efed30a9c39927f8a07b73a4be115e142 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 1 Dec 2020 13:44:13 +0100 Subject: [PATCH 05/10] Fix RANSAC plugin --- .../Plugins/Point_set/Point_set_shape_detection_plugin.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_shape_detection_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_shape_detection_plugin.cpp index d71f08a2274..f75d9e95aeb 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_shape_detection_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_shape_detection_plugin.cpp @@ -553,7 +553,7 @@ private: op.min_points = dialog.min_points(); // Only extract shapes with a minimum number of points. op.epsilon = dialog.epsilon(); // maximum euclidean distance between point and shape. op.cluster_epsilon = dialog.cluster_epsilon(); // maximum euclidean distance between points to be clustered. - op.normal_threshold = dialog.normal_tolerance(); // normal_threshold < dot(surface_normal, point_normal); + op.normal_threshold = std::cos(CGAL_PI * dialog.normal_tolerance() / 180.); // normal_threshold < dot(surface_normal, point_normal); CGAL::Random rand(static_cast(time(0))); // Gets point set From e21d019c05dbbadaa60b12e0e88f88fd143637c8 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 1 Dec 2020 14:01:00 +0100 Subject: [PATCH 06/10] Finalize validity test --- .../test/Shape_detection/CMakeLists.txt | 2 - .../test/Shape_detection/test_validity.cpp | 197 ------------------ .../test_validity_sampled_data.cpp | 9 +- 3 files changed, 8 insertions(+), 200 deletions(-) delete mode 100644 Shape_detection/test/Shape_detection/test_validity.cpp diff --git a/Shape_detection/test/Shape_detection/CMakeLists.txt b/Shape_detection/test/Shape_detection/CMakeLists.txt index d5a4eca990a..49f314edeb8 100644 --- a/Shape_detection/test/Shape_detection/CMakeLists.txt +++ b/Shape_detection/test/Shape_detection/CMakeLists.txt @@ -15,7 +15,6 @@ include(CGAL_CreateSingleSourceCGALProgram) find_package(Eigen3 3.1.0 QUIET) # (3.1.0 or greater) include(CGAL_Eigen_support) if(EIGEN3_FOUND) - create_single_source_cgal_program("test_validity.cpp") create_single_source_cgal_program("test_region_growing_basic.cpp") create_single_source_cgal_program("test_region_growing_on_cube.cpp") create_single_source_cgal_program("test_region_growing_on_point_set_2.cpp") @@ -31,7 +30,6 @@ if(EIGEN3_FOUND) "test_region_growing_on_degenerated_mesh.cpp") foreach( target - test_validity test_region_growing_basic test_region_growing_on_cube test_region_growing_on_point_set_2 diff --git a/Shape_detection/test/Shape_detection/test_validity.cpp b/Shape_detection/test/Shape_detection/test_validity.cpp deleted file mode 100644 index 3038082d850..00000000000 --- a/Shape_detection/test/Shape_detection/test_validity.cpp +++ /dev/null @@ -1,197 +0,0 @@ -#include - -//#define CGAL_RANSAC_EXPERIMENTAL_FIXES -#include -#include -#include - -#include - -#include - -#include - -namespace SD = CGAL::Shape_detection; - -using Kernel = CGAL::Simple_cartesian; -using Point_3 = Kernel::Point_3; -using Vector_3 = Kernel::Vector_3; - -using Pwn = std::pair; -using Point_set = std::vector; -using Point_map = CGAL::First_of_pair_property_map; -using Normal_map = CGAL::Second_of_pair_property_map; - -using RG_query = SD::Point_set::Sphere_neighbor_query; -using RG_region = SD::Point_set::Least_squares_plane_fit_region; -using Region_growing = SD::Region_growing; - -using RANSAC_traits = SD::Efficient_RANSAC_traits; -using RANSAC = SD::Efficient_RANSAC; -using RANSAC_plane = SD::Plane; - -void test_random_planes(std::size_t nb_planes); - -int main() -{ - // test_random_planes(1); - // test_random_planes(10); - // test_random_planes(100); - test_random_planes(1000); - test_random_planes(10000); - - return EXIT_SUCCESS; -} - -void test_random_planes(std::size_t nb_planes) -{ - CGAL::Random random; - std::cerr << "[TEST ON " << nb_planes << " RANDOM PLANES, SEED = " - << random.get_seed() << "] "; - - std::vector time_rg; - std::vector nb_detected_rg; - std::vector nb_unassigned_rg; - std::vector time_ransac; - std::vector nb_detected_ransac; - std::vector nb_unassigned_ransac; - - CGAL::Real_timer timer; - double timeout = 60.; // 1 minute timeout - - timer.start(); - std::size_t nb_runs = 100; - for (std::size_t run = 0; run < nb_runs; ++ run) - { - std::size_t min_points = random.get_int(10, 100); - std::size_t max_points = random.get_int(100, 1000); - double cluster_epsilon = random.get_double(0.01, 10.); - double epsilon = random.get_double (0., cluster_epsilon / 10.); - double normal_threshold = random.get_double (0.75, 0.99); - - double spacing = 0.8 * cluster_epsilon / std::sqrt(2); // smaller than diagonal - double domain_size = spacing * std::sqrt(nb_planes) * std::sqrt(0.5 * (min_points + max_points)) * 0.5; - double noise = 0.5 * epsilon; - - Point_set points; - for (std::size_t i = 0; i < nb_planes; ++ i) - { - // Generate random plane - Point_3 origin (random.get_double(-domain_size, domain_size), - random.get_double(-domain_size, domain_size), - random.get_double(-domain_size, domain_size)); - Vector_3 base1 (random.get_double(-1,1), random.get_double(-1,1), random.get_double(-1,1)); - base1 = base1 / std::sqrt(base1 * base1); - Vector_3 base2 (random.get_double(-1,1), random.get_double(-1,1), random.get_double(-1,1)); - base2 = base2 / std::sqrt(base2 * base2); - - Vector_3 normal = CGAL::cross_product (base1, base2); - normal = normal / std::sqrt (normal * normal); - - std::size_t nb_points = random.get_int (min_points, max_points); - - std::size_t nb_x = std::size_t(std::sqrt (double(nb_points))); - std::size_t nb_y = std::size_t(nb_points / double(nb_x)) + 1; - - for (std::size_t j = 0; j < nb_x; ++ j) - for (std::size_t k = 0; k < nb_y; ++ k) - points.emplace_back (origin + j * spacing * base1 + k * spacing * base2 - + normal * random.get_double(-noise/2, noise/2), - normal); - } - - // std::ofstream test("dump.pwn"); - // for (const auto& p : points) - // test << p.first << " " << p.second << std::endl; - CGAL::Real_timer t; - t.start(); - RG_query rg_query (points, cluster_epsilon); - RG_region rg_region (points, epsilon, normal_threshold, min_points); - Region_growing region_growing (points, rg_query, rg_region); - std::size_t nb_detected = 0; - std::size_t nb_unassigned = 0; - region_growing.detect (boost::make_function_output_iterator ([&](const auto&) { ++ nb_detected; })); - region_growing.unassigned_items (boost::make_function_output_iterator ([&](const auto&) { ++ nb_unassigned; })); - t.stop(); - - time_rg.push_back (t.time()); - nb_detected_rg.push_back (nb_detected); - nb_unassigned_rg.push_back (nb_unassigned / double(points.size())); - t.reset(); - - t.start(); - RANSAC ransac; - ransac.template add_shape_factory(); - ransac.set_input(points); - typename RANSAC::Parameters parameters; - parameters.probability = 0.05f; - parameters.min_points = min_points; - parameters.epsilon = epsilon; - parameters.cluster_epsilon = cluster_epsilon; - parameters.normal_threshold = normal_threshold; - ransac.detect(parameters); - t.stop(); - - time_ransac.push_back (t.time()); - nb_detected_ransac.push_back (ransac.shapes().size()); - nb_unassigned_ransac.push_back (ransac.number_of_unassigned_points() / double(points.size())); - -#if 0 - if (ransac.shapes().size() == 0) - { - std::cerr << "Detected 0 shapes with " << std::endl - << " * min points = " << min_points << std::endl - << " * epsilon = " << epsilon << std::endl - << " * cluster_epsilon = " << cluster_epsilon << std::endl - << " * normal_threshold = " << normal_threshold << std::endl; - - std::ofstream ofile("0shapes.ply", std::ios::binary); - CGAL::set_binary_mode (ofile); - CGAL::write_ply_points (ofile, points, - CGAL::parameters::point_map(Point_map()). - normal_map(Normal_map())); - ofile.close(); - exit(0); - } -#endif - - if (timer.time() > timeout) - { - nb_runs = run + 1; - break; - } - } - - std::cerr << "on " << nb_runs << " runs" << std::endl; - - - std::sort (time_ransac.begin(), time_ransac.end()); - std::sort (nb_detected_ransac.begin(), nb_detected_ransac.end()); - std::sort (nb_unassigned_ransac.begin(), nb_unassigned_ransac.end()); - std::sort (time_rg.begin(), time_rg.end()); - std::sort (nb_detected_rg.begin(), nb_detected_rg.end()); - std::sort (nb_unassigned_rg.begin(), nb_unassigned_rg.end()); - - std::cerr << " * Region Growing" << std::endl - << " - took between " << time_rg.front() << "s and " << time_rg.back() - << "s, median = " << time_rg[nb_runs / 2] << "s" << std::endl - << " - detected between " << nb_detected_rg.front() << " and " - << nb_detected_rg.back() << " planes (" << 100. * nb_detected_rg.front() / double(nb_planes) - << "% to " << 100. * nb_detected_rg.back() / double(nb_planes) << "%), median = " - << nb_detected_rg[nb_runs / 2] << " planes (" << 100. * nb_detected_rg[nb_runs / 2] / double(nb_planes) - << "%)" << std::endl - << " - left between " << 100. * nb_unassigned_rg.front() - << "% and " << 100. * nb_unassigned_rg.back() << "% of unassigned points, median = " - << 100. * nb_unassigned_rg[nb_runs / 2] << "%" << std::endl; - std::cerr << " * Efficient RANSAC" << std::endl - << " - took between " << time_ransac.front() << "s and " << time_ransac.back() - << "s, median = " << time_ransac[nb_runs / 2] << "s" << std::endl - << " - detected between " << nb_detected_ransac.front() << " and " - << nb_detected_ransac.back() << " planes (" << 100. * nb_detected_ransac.front() / double(nb_planes) - << "% to " << 100. * nb_detected_ransac.back() / double(nb_planes) << "%), median = " - << nb_detected_ransac[nb_runs / 2] << " planes (" << 100. * nb_detected_ransac[nb_runs / 2] / double(nb_planes) - << "%)" << std::endl - << " - left between " << 100. * nb_unassigned_ransac.front() - << "% and " << 100. * nb_unassigned_ransac.back() << "% of unassigned points, median = " - << 100. * nb_unassigned_ransac[nb_runs / 2] << "%" << std::endl; -} diff --git a/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp b/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp index a194e9f43e9..85ee4de534f 100644 --- a/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp +++ b/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp @@ -1,6 +1,7 @@ #include #define CGAL_RANSAC_EXPERIMENTAL_FIXES +#define USE_WEIGHTED_LEVELS #include #include #include @@ -75,7 +76,8 @@ void test_copied_point_cloud (const Point_set& original_points, std::size_t nb) (CGAL::make_transform_iterator_from_property_map (original_points.begin(), Point_map()), CGAL::make_transform_iterator_from_property_map (original_points.end(), Point_map())); - std::cerr << "Ground truth = " << 6*nb*nb + 1 << " planes" << std::endl; + std::size_t ground_truth = 6*nb*nb+1; + std::cerr << "Ground truth = " << ground_truth << " planes" << std::endl; Point_set points; points.reserve (nb * nb * original_points.size()); @@ -116,6 +118,8 @@ void test_copied_point_cloud (const Point_set& original_points, std::size_t nb) t.stop(); std::cerr << "Region Growing = " << nb_detected << " planes (" << 1000 * t.time() << "ms)" << std::endl; + assert (nb_detected == ground_truth); + CGAL::Real_timer timer; double timeout = 120; // 2 minutes timeout timer.start(); @@ -151,6 +155,9 @@ void test_copied_point_cloud (const Point_set& original_points, std::size_t nb) << detected_ransac.front() << ";" << detected_ransac.back() << "], time[" << times_ransac.front() << ";" << times_ransac.back() << "])" << std::endl; + // RANSAC should at least detect 75% of shapes + assert (detected_ransac[detected_ransac.size() / ] > std::size_t(0.75 * ground_truth)); + #ifdef CGAL_TEST_RANSAC_PROTOTYPE { CGAL::Real_timer timer; From 9294461c64c88142e4c8391af128b968a50043e6 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 1 Dec 2020 14:56:17 +0100 Subject: [PATCH 07/10] Call ctest on new validity test Co-authored-by: Laurent Rineau --- Shape_detection/test/Shape_detection/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Shape_detection/test/Shape_detection/CMakeLists.txt b/Shape_detection/test/Shape_detection/CMakeLists.txt index 49f314edeb8..cfe7db967d8 100644 --- a/Shape_detection/test/Shape_detection/CMakeLists.txt +++ b/Shape_detection/test/Shape_detection/CMakeLists.txt @@ -56,7 +56,7 @@ if(EIGEN3_FOUND) add_executable(test_validity_sampled_data "test_validity_sampled_data.cpp") target_link_libraries(test_validity_sampled_data CGAL::CGAL CGAL::Eigen_support) endif() - + cgal_add_test(test_validity_sampled_data) endif() create_single_source_cgal_program( From 09b61c13d9be3680d2f62a7d5498bc9585715d76 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 2 Dec 2020 09:22:48 +0100 Subject: [PATCH 08/10] Fix assertion --- .../test/Shape_detection/test_validity_sampled_data.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp b/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp index 85ee4de534f..c146fb6d626 100644 --- a/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp +++ b/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp @@ -156,7 +156,7 @@ void test_copied_point_cloud (const Point_set& original_points, std::size_t nb) << times_ransac.front() << ";" << times_ransac.back() << "])" << std::endl; // RANSAC should at least detect 75% of shapes - assert (detected_ransac[detected_ransac.size() / ] > std::size_t(0.75 * ground_truth)); + assert (detected_ransac[detected_ransac.size() / 2] > std::size_t(0.75 * ground_truth)); #ifdef CGAL_TEST_RANSAC_PROTOTYPE { From 873aa26a8f33834513131ddaf8e8bdfd6b48e4b1 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Thu, 3 Dec 2020 09:27:01 +0100 Subject: [PATCH 09/10] Avoid testsuite timeouts --- .../Shape_detection/test_validity_sampled_data.cpp | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp b/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp index c146fb6d626..d603e50f1ee 100644 --- a/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp +++ b/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp @@ -64,8 +64,10 @@ int main (int argc, char** argv) test_copied_point_cloud (points, 2); test_copied_point_cloud (points, 5); test_copied_point_cloud (points, 10); +#ifndef CGAL_TEST_SUITE // Disable tests too large for testsuite test_copied_point_cloud (points, 20); test_copied_point_cloud (points, 50); +#endif return EXIT_SUCCESS; } @@ -120,10 +122,17 @@ void test_copied_point_cloud (const Point_set& original_points, std::size_t nb) assert (nb_detected == ground_truth); - CGAL::Real_timer timer; +#ifdef CGAL_TEST_SUITE + double timeout = 60; // 1 minute timeout + std::size_t nb_runs = 20; // +#else double timeout = 120; // 2 minutes timeout - timer.start(); std::size_t nb_runs = 500; +#endif + + CGAL::Real_timer timer; + timer.start(); + std::vector detected_ransac; std::vector times_ransac; for (std::size_t run = 0; run < nb_runs; ++ run) From b2573a7b38be6711a8853cdd6bb2df6cb8b25fd4 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Thu, 3 Dec 2020 09:34:31 +0100 Subject: [PATCH 10/10] Fix Windows min/max garbage --- .../Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h index 1da816c26a9..beff95c50e9 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h +++ b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h @@ -555,9 +555,9 @@ namespace CGAL { do { // Search (remaining_points / min_points) shapes (max 200 per iteration, min 1) std::size_t search_number - = std::min(std::size_t(200), - std::max(std::size_t((m_num_available_points - num_invalid) / double(m_options.min_points)), - std::size_t(1))); + = (std::min)(std::size_t(200), + (std::max)(std::size_t((m_num_available_points - num_invalid) / double(m_options.min_points)), + std::size_t(1))); for (std::size_t nb = 0; nb < search_number; ++ nb) { // Generate candidates