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 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..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 @@ -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; } @@ -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, - 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) @@ -1019,7 +1027,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 +1048,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; diff --git a/Shape_detection/test/Shape_detection/CMakeLists.txt b/Shape_detection/test/Shape_detection/CMakeLists.txt index 9e6719e15ce..cfe7db967d8 100644 --- a/Shape_detection/test/Shape_detection/CMakeLists.txt +++ b/Shape_detection/test/Shape_detection/CMakeLists.txt @@ -42,6 +42,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() + cgal_add_test(test_validity_sampled_data) endif() create_single_source_cgal_program( 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..d603e50f1ee --- /dev/null +++ b/Shape_detection/test/Shape_detection/test_validity_sampled_data.cpp @@ -0,0 +1,244 @@ +#include + +#define CGAL_RANSAC_EXPERIMENTAL_FIXES +#define USE_WEIGHTED_LEVELS +#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); +#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; +} + +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::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()); + 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; + + assert (nb_detected == ground_truth); + +#ifdef CGAL_TEST_SUITE + double timeout = 60; // 1 minute timeout + std::size_t nb_runs = 20; // +#else + double timeout = 120; // 2 minutes timeout + 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) + { + 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; + + // RANSAC should at least detect 75% of shapes + assert (detected_ransac[detected_ransac.size() / 2] > std::size_t(0.75 * ground_truth)); + +#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; +}