From f741e3516d287d9681523ca595b36459046bb992 Mon Sep 17 00:00:00 2001 From: Sven Oesau Date: Fri, 27 Jan 2023 19:21:15 +0100 Subject: [PATCH] addition of Kinetic_shape_partitioning_Traits_3 separation of Kinetic_shape_reconstruction_3 into Kinetic_shape_partitioning_3 and Kinetic_shape_reconstruction_3 (WIP) adding Linear_cell_complex to dependencies renaming a few parameters in KSR::Parameters_3 removal and renaming of named parameters one pass on the documentation --- .../Kinetic_shape_reconstruction/dependencies | 2 + .../kinetic_precomputed_shapes.cpp | 26 +- .../kinetic_random_shapes.cpp | 21 +- .../kinetic_reconstruction.cpp | 76 +- .../include/CGAL/KSR/conversions.h | 97 -- .../include/CGAL/KSR/parameters.h | 12 +- .../include/CGAL/KSR_3/Data_structure.h | 73 +- .../include/CGAL/KSR_3/FacePropagation.h | 25 +- .../include/CGAL/KSR_3/Finalizer.h | 11 +- .../include/CGAL/KSR_3/Graphcut.h | 22 +- .../include/CGAL/KSR_3/Initializer.h | 66 +- .../include/CGAL/KSR_3/Intersection_graph.h | 21 +- .../include/CGAL/KSR_3/Reconstruction.h | 51 +- .../include/CGAL/KSR_3/Support_plane.h | 74 +- .../include/CGAL/KSR_3/Visibility.h | 36 +- .../CGAL/Kinetic_shape_partitioning_3.h | 1000 +++++++++++++ .../CGAL/Kinetic_shape_reconstruction_3.h | 1299 ++++++----------- .../kinetic_2d_stress_test.cpp | 12 +- .../kinetic_3d_test_all.cpp | 100 +- .../internal/parameters_interface.h | 8 +- 20 files changed, 1785 insertions(+), 1247 deletions(-) delete mode 100644 Kinetic_shape_reconstruction/include/CGAL/KSR/conversions.h create mode 100644 Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_partitioning_3.h diff --git a/Kinetic_shape_reconstruction/doc/Kinetic_shape_reconstruction/dependencies b/Kinetic_shape_reconstruction/doc/Kinetic_shape_reconstruction/dependencies index 8c8126899c5..0f8c15ab6cb 100644 --- a/Kinetic_shape_reconstruction/doc/Kinetic_shape_reconstruction/dependencies +++ b/Kinetic_shape_reconstruction/doc/Kinetic_shape_reconstruction/dependencies @@ -1,3 +1,5 @@ Manual Kernel_23 BGL +Linear_cell_complex +Surface_mesh diff --git a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes.cpp b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes.cpp index e287f1668a4..9a8faf02baa 100644 --- a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes.cpp +++ b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes.cpp @@ -1,7 +1,7 @@ #include #include #include -#include +#include #include #include #include @@ -12,6 +12,8 @@ using SCD = CGAL::Simple_cartesian; using EPICK = CGAL::Exact_predicates_inexact_constructions_kernel; using EPECK = CGAL::Exact_predicates_exact_constructions_kernel; +using Traits = typename CGAL::Kinetic_shape_partitioning_Traits_3, CGAL::Identity_property_map >; + using Kernel = EPICK; using FT = typename Kernel::FT; using Point_2 = typename Kernel::Point_2; @@ -20,7 +22,7 @@ using Segment_3 = typename Kernel::Segment_3; using Triangle_2 = typename Kernel::Triangle_2; using Surface_mesh = CGAL::Surface_mesh; -using KSR = CGAL::Kinetic_shape_reconstruction_3; +using KSR = CGAL::Kinetic_shape_partitioning_3; using Timer = CGAL::Real_timer; struct Polygon_map { @@ -81,10 +83,8 @@ int main(const int argc, const char** argv) { const bool verbose = true; const bool debug = false; const unsigned int k = (argc > 2 ? std::atoi(argv[2]) : 1); // intersections - const unsigned int subdiv = 0; const double eratio = 1.1; const bool orient = false; - const bool use_hm = false; // Algorithm. KSR ksr(verbose, debug); @@ -92,24 +92,24 @@ int main(const int argc, const char** argv) { Timer timer; timer.start(); - const bool is_ksr_success = ksr.partition( + bool is_ksr_success = ksr.initialize( input_faces, polygon_map, CGAL::parameters:: - k_intersections(k). - n_subdivisions(subdiv). - enlarge_bbox_ratio(eratio). - reorient(orient). - use_hybrid_mode(use_hm)); + bbox_dilation_ratio(eratio). + reorient_bbox(orient)); + + if (is_ksr_success) + is_ksr_success = ksr.partition(k); + assert(is_ksr_success); const std::string success = is_ksr_success ? "SUCCESS" : "FAILED"; timer.stop(); const FT time = static_cast(timer.time()); - const std::size_t num_events = ksr.number_of_events(); // Output. const int support_plane_idx = -1; const int num_support_planes = ksr.number_of_support_planes(); assert(num_support_planes > 6); - assert(ksr.support_plane_index(0) == 6); + assert(ksr.support_plane_index(0) == 6);/* // Vertices. const std::size_t num_vertices = ksr.number_of_vertices(support_plane_idx); @@ -218,7 +218,7 @@ int main(const int argc, const char** argv) { // CGAL::IO::write_PLY(output_file_support_plane, support_planes[i]); // output_file_support_plane.close(); // } - // std::cout << "* partition support planes exported successfully" << std::endl; + // std::cout << "* partition support planes exported successfully" << std::endl;*/ std::cout << std::endl << "3D KINETIC " << success << " in " << time << " seconds!" << std::endl << std::endl; diff --git a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_random_shapes.cpp b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_random_shapes.cpp index 0fbf256b5b1..2b8a9f768f6 100644 --- a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_random_shapes.cpp +++ b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_random_shapes.cpp @@ -6,7 +6,7 @@ #include #include #include -#include +#include #include #include #include @@ -36,7 +36,9 @@ using IPoint_3 = typename EPICK::Point_3; using IPolygon_3 = std::vector; using IPolygon_3_map = CGAL::Identity_property_map; -using KSR = CGAL::Kinetic_shape_reconstruction_3; +using Traits = typename CGAL::Kinetic_shape_partitioning_Traits_3, CGAL::Identity_property_map >; + +using KSP = CGAL::Kinetic_shape_partitioning_3; const std::vector box_vertices_to_faces(const int i) { const int _vertices_to_faces[8][3] = { @@ -403,14 +405,19 @@ int main(const int argc, const char** argv) { assert(input_polygons.size() == rnd_polygons.size()); // Algorithm. - KSR ksr(true, false); + KSP ksp(true, false); const IPolygon_3_map polygon_map; const unsigned int k = (argc > 3 ? std::atoi(argv[3]) : 1); std::cout << "* input k: " << k << std::endl; - const bool is_ksr_success = ksr.partition( - input_polygons, polygon_map, CGAL::parameters::k_intersections(k)); - assert(is_ksr_success); - const std::string success = is_ksr_success ? "SUCCESS" : "FAILED"; + + bool is_ksp_success = ksp.initialize( + input_polygons, polygon_map); + + if (is_ksp_success) + ksp.partition(k); + + assert(is_ksp_success); + const std::string success = is_ksp_success ? "SUCCESS" : "FAILED"; std::cout << std::endl << "3D KINETIC " << success << "!" << std::endl << std::endl; return EXIT_SUCCESS; diff --git a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_reconstruction.cpp b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_reconstruction.cpp index adddb2929a1..660f7d58088 100644 --- a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_reconstruction.cpp +++ b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_reconstruction.cpp @@ -10,6 +10,7 @@ #include "include/Terminal_parser.h" using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using EPECK = CGAL::Exact_predicates_exact_constructions_kernel; using FT = typename Kernel::FT; using Point_3 = typename Kernel::Point_3; using Vector_3 = typename Kernel::Vector_3; @@ -17,12 +18,14 @@ using Segment_3 = typename Kernel::Segment_3; using Point_set = CGAL::Point_set_3; using Point_map = typename Point_set::Point_map; -using Vector_map = typename Point_set::Vector_map; +using Normal_map = typename Point_set::Vector_map; using Label_map = typename Point_set:: template Property_map; using Semantic_map = CGAL::KSR::Semantic_from_label_map; using Region_map = typename Point_set:: template Property_map; -using KSR = CGAL::Kinetic_shape_reconstruction_3; +using Traits = typename CGAL::Kinetic_shape_partitioning_Traits_3; + +using KSR = CGAL::Kinetic_shape_reconstruction_3; using Parameters = CGAL::KSR::All_parameters; using Terminal_parser = CGAL::KSR::Terminal_parser; @@ -111,7 +114,6 @@ int main(const int argc, const char** argv) { } else { std::ifstream input_file(parameters.data, std::ios_base::binary); - bool f = input_file.is_open(); input_file >> point_set; input_file.close(); } @@ -129,18 +131,6 @@ int main(const int argc, const char** argv) { std::cout << "verbose " << parameters.verbose << std::endl; std::cout << "debug " << parameters.debug << std::endl; - // Define a map from a user-defined label to the semantic label. - const Label_map label_map = point_set. template property_map("label").first; - const bool is_defined = point_set. template property_map("label").second; - const Semantic_map semantic_map( - label_map, - is_defined, - parameters.gi, - parameters.bi, - parameters.ii, - parameters.vi, - parameters.verbose); - // Algorithm. KSR ksr(parameters.verbose, parameters.debug); @@ -149,13 +139,38 @@ int main(const int argc, const char** argv) { Timer timer; timer.start(); - bool is_ksr_success; + std::size_t num_shapes = ksr.detect_planar_shapes(point_set, + CGAL::parameters::distance_threshold(parameters.distance_threshold) + .angle_threshold(parameters.angle_threshold) + .k_neighbors(parameters.k_neighbors) + .min_region_size(parameters.min_region_size)); + + std::cout << num_shapes << " detected planar shapes" << std::endl; + + num_shapes = ksr.regularize_shapes(CGAL::parameters::regularize_parallelism(true).regularize_coplanarity(true).regularize_axis_symmetry(false).regularize_orthogonality(false)); + + std::cout << num_shapes << " detected planar shapes after regularization" << std::endl; + + bool is_ksr_success = ksr.initialize_partitioning(); + + if (!is_ksr_success) { + std::cout << "Initializing kinetic partitioning failed!" << std::endl; + return 1; + } + + is_ksr_success = ksr.partition(parameters.k_intersections); + + if (!is_ksr_success) { + std::cout << "Initializing kinetic partitioning failed!" << std::endl; + return 2; + } + + ksr.setup_energyterms(); + + ksr.reconstruct(parameters.graphcut_beta); +/* if (is_segmented) is_ksr_success = ksr.reconstruct( - point_set, - point_set.point_map(), - point_set.normal_map(), - semantic_map, region_map, CGAL::parameters:: k_neighbors(parameters.k_neighbors). @@ -167,10 +182,6 @@ int main(const int argc, const char** argv) { graphcut_beta(parameters.graphcut_beta)); else is_ksr_success = ksr.reconstruct( - point_set, - point_set.point_map(), - point_set.normal_map(), - semantic_map, base, CGAL::parameters:: k_neighbors(parameters.k_neighbors). @@ -179,35 +190,38 @@ int main(const int argc, const char** argv) { min_region_size(parameters.min_region_size). regularize(parameters.regularize). k_intersections(parameters.k_intersections). - graphcut_beta(parameters.graphcut_beta)); + graphcut_beta(parameters.graphcut_beta));*/ assert(is_ksr_success); const std::string success = is_ksr_success ? "SUCCESS" : "FAILED"; timer.stop(); const FT time = static_cast(timer.time()); + const KSR::KSP& ksp = ksr.partitioning(); + // Output. CGAL::Linear_cell_complex_for_combinatorial_map<3, 3> lcc; - ksr.get_linear_cell_complex(lcc); + ksp.get_linear_cell_complex(lcc); +/* // Vertices. std::vector all_vertices; - ksr.output_partition_vertices( + ksp.output_partition_vertices( std::back_inserter(all_vertices), -1); // Edges. std::vector all_edges; - ksr.output_partition_edges( + ksp.output_partition_edges( std::back_inserter(all_edges), -1); // Faces. std::vector< std::vector > all_faces; - ksr.output_partition_faces( + ksp.output_partition_faces( std::back_inserter(all_faces), -1, 6); // Model. std::vector output_vertices; std::vector< std::vector > output_faces; - ksr.output_reconstructed_model( + ksp.output_reconstructed_model( std::back_inserter(output_vertices), std::back_inserter(output_faces)); const std::size_t num_vertices = output_vertices.size(); @@ -251,7 +265,7 @@ int main(const int argc, const char** argv) { return EXIT_FAILURE; } output_file_model.close(); - std::cout << "* the reconstructed model exported successfully" << std::endl; + std::cout << "* the reconstructed model exported successfully" << std::endl;*/ std::cout << std::endl << "3D KINETIC RECONSTRUCTION " << success << " in " << time << " seconds!" << std::endl << std::endl; diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR/conversions.h b/Kinetic_shape_reconstruction/include/CGAL/KSR/conversions.h deleted file mode 100644 index 1ca0862d0d0..00000000000 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR/conversions.h +++ /dev/null @@ -1,97 +0,0 @@ -// Copyright (c) 2021 GeometryFactory SARL (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// -// $URL$ -// $Id$ -// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial -// -// -// Author(s) : Simon Giraudot, Dmitry Anisimov - -#ifndef CGAL_KSR_CONVERSIONS_H -#define CGAL_KSR_CONVERSIONS_H - -// #include - -// STL includes. -#include -#include - -// CGAL includes. -#include -#include -#include -#include - -namespace CGAL { -namespace KSR { - -template -class Kinetic_traits_3 { - - using IK = GeomTraits; // assume GT is inexact here - - using IT = typename GeomTraits::FT; - using Point_2 = typename GeomTraits::Point_2; - using Line_2 = typename GeomTraits::Line_2; - - // TODO: This is very naive way of doing this. We should compare IT and ET - // and, in case they are the same or IT is already exact, abort conversion! - using EK = CGAL::Exact_predicates_exact_constructions_kernel; - using ET = typename EK::FT; - - using IK_to_EK = CGAL::Cartesian_converter; - using EK_to_IK = CGAL::Cartesian_converter; - -public: - Kinetic_traits_3() : - m_use_hybrid_mode(false) { } - - inline const Point_2 intersection(const Line_2& t1, const Line_2& t2) const { - - Point_2 out; - const bool is_intersection_found = intersection(t1, t2, out); - CGAL_assertion(is_intersection_found); - return out; - } - - template - inline bool intersection( - const Type1& t1, const Type2& t2, ResultType& result) const { - - const auto inter = intersection_impl(t1, t2); - if (!inter) return false; - if (const ResultType* typed_inter = boost::get(&*inter)) { - result = *typed_inter; - return true; - } - return false; - } - -private: - const bool m_use_hybrid_mode; - IK_to_EK m_inexact_to_exact; - EK_to_IK m_exact_to_inexact; - - template - decltype(auto) intersection_impl(const Type1& t1, const Type2& t2) const { - - // if (!m_use_hybrid_mode) { - return CGAL::intersection(t1, t2); - // } else { // convert to exact - - // // TODO: It does not compile with EPECK as input kernel. - // const auto exact_t1 = m_inexact_to_exact(t1); - // const auto exact_t2 = m_inexact_to_exact(t2); - // const auto exact_result = CGAL::intersection(exact_t1, exact_t2); - // return m_exact_to_inexact(exact_result); - // } - } -}; - -} // namespace KSR -} // namespace CGAL - -#endif // CGAL_KSR_CONVERSIONS_H diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR/parameters.h b/Kinetic_shape_reconstruction/include/CGAL/KSR/parameters.h index 6c9cad59b4f..7a4315d1b5b 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR/parameters.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR/parameters.h @@ -24,19 +24,19 @@ struct Parameters_3 { unsigned int k = 1; // k intersections unsigned int n = 0; // n subdivisions, not implemented yet - FT enlarge_bbox_ratio = FT(11) / FT(10); // ratio to enlarge bbox - FT distance_tolerance = FT(0.005) / FT(10); // distance tolerance between planes + FT bbox_dilation_ratio = FT(11) / FT(10); // ratio to enlarge bbox + FT angle_tolerance = FT(5); + FT distance_tolerance = FT(5) / FT(10); // distance tolerance between planes - bool reorient = false; // true - optimal bounding box, false - axis aligned + bool reorient_bbox = false; // true - optimal bounding box, false - axis aligned // All files are saved in the current build directory. bool verbose = true; // print basic verbose information bool debug = false; // print all steps and substeps + export initial and final configurations - bool export_all = false; // export all intermediate configurations and events // See also global tolerance inside utils.h! - Parameters_3(const bool v = true, const bool d = false, const bool e = false) : - verbose(v), debug(d), export_all(e) { } + Parameters_3(const bool v = true, const bool d = false) : + verbose(v), debug(d) { } }; } // namespace KSR diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h index 0bee15d2ed2..3c1e57d22d3 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h @@ -24,7 +24,6 @@ #include #include #include -#include #include #include @@ -35,40 +34,40 @@ namespace KSR_3 { #ifdef DOXYGEN_RUNNING #else -template +template class Data_structure { public: - using Kernel = GeomTraits; - using Intersection_Kernel = Intersection_Traits; + using Kernel = typename Traits::Kernel; + using Intersection_Kernel = typename Traits::Intersection_Kernel; - using Support_plane = KSR_3::Support_plane; - using Intersection_graph = KSR_3::Intersection_graph; + using Support_plane = KSR_3::Support_plane; + using Intersection_graph = KSR_3::Intersection_graph; using FaceEvent = typename Support_plane::FaceEvent; private: - using FT = typename Kernel::FT; - using Point_2 = typename Kernel::Point_2; + using FT = typename Traits::FT; + using Point_2 = typename Traits::Point_2; using IkPoint_2 = typename Intersection_Kernel::Point_2; - using Point_3 = typename Kernel::Point_3; + using Point_3 = typename Traits::Point_3; using IkPoint_3 = typename Intersection_Kernel::Point_3; - using Segment_2 = typename Kernel::Segment_2; + using Segment_2 = typename Traits::Segment_2; using IkSegment_2 = typename Intersection_Kernel::Segment_2; - using Segment_3 = typename Kernel::Segment_3; + using Segment_3 = typename Traits::Segment_3; using IkSegment_3 = typename Intersection_Kernel::Segment_3; - using Vector_2 = typename Kernel::Vector_2; - using Direction_2 = typename Kernel::Direction_2; + using Vector_2 = typename Traits::Vector_2; + using Direction_2 = typename Traits::Direction_2; using IkDirection_2 = typename Intersection_Kernel::Direction_2; using Triangle_2 = typename Kernel::Triangle_2; - using Line_2 = typename Kernel::Line_2; + using Line_2 = typename Traits::Line_2; using IkLine_2 = typename Intersection_Kernel::Line_2; - using Plane_3 = typename Kernel::Plane_3; + using Plane_3 = typename Traits::Plane_3; using Polygon_2 = CGAL::Polygon_2; using Parameters = KSR::Parameters_3; - using To_exact = CGAL::Cartesian_converter; - using From_exact = CGAL::Cartesian_converter; + using To_exact = typename Traits::To_exact; + using From_exact = typename Traits::From_exact; public: using Mesh = typename Support_plane::Mesh; @@ -185,6 +184,7 @@ public: struct Volume_cell { std::vector pfaces; + std::vector faces;// Indices into m_face2vertices in m_data. std::vector pface_oriented_outwards; std::vector neighbors; std::set pvertices; @@ -218,8 +218,6 @@ public: } }; - using Kinetic_traits = KSR::Kinetic_traits_3; - private: std::map< std::pair, Point_2> m_points; std::map< std::pair, Vector_2> m_directions; @@ -230,18 +228,29 @@ private: From_exact from_exact; const Parameters& m_parameters; - Kinetic_traits m_kinetic_traits; std::vector m_volumes; - std::map > m_map_volumes; + std::vector > m_face2vertices; + std::map > m_face2volumes; + std::map > m_face2input_polygon; std::map m_input_polygon_map; Reconstructed_model m_reconstructed_model; + template + inline bool intersection( + const Type1& t1, const Type2& t2, ResultType& result) const { + + const auto inter = CGAL::intersection(t1, t2); + if (!inter) return false; + if (const ResultType* typed_inter = boost::get(&*inter)) { + result = *typed_inter; + return true; + } + return false; + } + public: - Data_structure(const Parameters& parameters) : - m_parameters(parameters), - m_kinetic_traits() - { } + Data_structure(const Parameters& parameters) : m_parameters(parameters), to_exact(), from_exact() { } /******************************* ** INITIALIZATION ** @@ -254,7 +263,7 @@ public: m_intersection_graph.clear(); m_volumes.clear(); - m_map_volumes.clear(); + m_face2volumes.clear(); m_input_polygon_map.clear(); m_reconstructed_model.clear(); } @@ -299,8 +308,8 @@ public: ** GENERAL ** ********************************/ - std::map >& pface_neighbors() { return m_map_volumes; } - const std::map >& pface_neighbors() const { return m_map_volumes; } + std::map >& pface_neighbors() { return m_face2volumes; } + const std::map >& pface_neighbors() const { return m_face2volumes; } const std::vector& support_planes() const { return m_support_planes; } std::vector& support_planes() { return m_support_planes; } @@ -568,6 +577,9 @@ public: std::vector& volumes() { return m_volumes; } const std::vector& volumes() const { return m_volumes; } + const std::vector face(std::size_t face_index) const { return &m_face2vertices[face_index]; } + const Point_3& vertex(std::size_t vertex_index) const { return &m_face2vertices[face_index]; } + Reconstructed_model& reconstructed_model() { return m_reconstructed_model; } const Reconstructed_model& reconstructed_model() const { return m_reconstructed_model; } @@ -604,7 +616,7 @@ public: const PointRange& polygon, const bool is_bbox) { const Support_plane new_support_plane( - polygon, is_bbox, m_parameters.distance_tolerance, number_of_support_planes()); + polygon, is_bbox, m_parameters.distance_tolerance, m_parameters.angle_tolerance, number_of_support_planes()); std::size_t support_plane_idx = KSR::no_element(); for (std::size_t i = 0; i < number_of_support_planes(); ++i) { @@ -644,9 +656,8 @@ public: for (const auto iedge : all_iedges) { const auto segment = segment_3(iedge); typename Intersection_graph::Edge_property* p = (typename Intersection_graph::Edge_property*)iedge.get_property(); - if (!m_kinetic_traits.intersection(plane, segment, point)) { + if (!intersection(plane, segment, point)) continue; - } const auto isource = source(iedge); const auto itarget = target(iedge); diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/FacePropagation.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/FacePropagation.h index 6743b9871ab..df1cfd2752d 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/FacePropagation.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/FacePropagation.h @@ -29,21 +29,22 @@ namespace KSR_3 { #ifdef DOXYGEN_RUNNING #else -template +template class FacePropagation { public: - using Kernel = GeomTraits; + using Kernel = typename Traits::Kernel; + using Intersection_Kernel = typename Traits::Intersection_Kernel; private: - using FT = typename Kernel::FT; - using Point_2 = typename Kernel::Point_2; - using Vector_2 = typename Kernel::Vector_2; - using Segment_2 = typename Kernel::Segment_2; - using Direction_2 = typename Kernel::Direction_2; - using Line_2 = typename Kernel::Line_2; + using FT = typename Traits::FT; + using Point_2 = typename Traits::Point_2; + using Vector_2 = typename Traits::Vector_2; + using Segment_2 = typename Traits::Segment_2; + using Direction_2 = typename Traits::Direction_2; + using Line_2 = typename Traits::Line_2; - using Data_structure = KSR_3::Data_structure; + using Data_structure = KSR_3::Data_structure; using IVertex = typename Data_structure::IVertex; using IEdge = typename Data_structure::IEdge; @@ -57,7 +58,6 @@ private: using Face_index = typename Data_structure::Face_index; using Parameters = KSR::Parameters_3; - using Kinetic_traits = KSR::Kinetic_traits_3; using FaceEvent = typename Data_structure::Support_plane::FaceEvent; @@ -69,11 +69,11 @@ private: public: FacePropagation(Data_structure& data, const Parameters& parameters) : - m_data(data), m_parameters(parameters), m_kinetic_traits(), + m_data(data), m_parameters(parameters), m_min_time(-FT(1)), m_max_time(-FT(1)) { } - const std::pair propagate() { + const std::pair propagate(std::size_t k) { std::size_t num_queue_calls = 0; std::size_t num_events = 0; @@ -97,7 +97,6 @@ public: private: Data_structure& m_data; const Parameters& m_parameters; - Kinetic_traits m_kinetic_traits; FT m_min_time; FT m_max_time; diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Finalizer.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Finalizer.h index 4e9a2832344..7c4ff1ddad4 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Finalizer.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Finalizer.h @@ -20,7 +20,6 @@ #include #include #include -#include #include @@ -30,11 +29,11 @@ namespace KSR_3 { #ifdef DOXYGEN_RUNNING #else -template +template class Finalizer { public: - using Kernel = GeomTraits; + using Kernel = typename Traits; private: using FT = typename Kernel::FT; @@ -48,7 +47,7 @@ private: using Direction_2 = typename Kernel::Direction_2; using Tetrahedron_3 = typename Kernel::Tetrahedron_3; - using Data_structure = KSR_3::Data_structure; + using Data_structure = KSR_3::Data_structure; using IVertex = typename Data_structure::IVertex; using IEdge = typename Data_structure::IEdge; @@ -91,11 +90,10 @@ private: }; using Parameters = KSR::Parameters_3; - using Kinetic_traits = KSR::Kinetic_traits_3; public: Finalizer(Data_structure& data, const Parameters& parameters) : - m_data(data), m_parameters(parameters), m_kinetic_traits() + m_data(data), m_parameters(parameters) { } void create_polyhedra() { @@ -124,7 +122,6 @@ public: private: Data_structure& m_data; const Parameters& m_parameters; - Kinetic_traits m_kinetic_traits; /******************************* ** EXTRACTING VOLUMES ** diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Graphcut.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Graphcut.h index 0f5cb253ee6..d17529038ba 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Graphcut.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Graphcut.h @@ -38,20 +38,20 @@ namespace KSR_3 { #ifdef DOXYGEN_RUNNING #else - template + template class Graphcut { public: - using Kernel = GeomTraits; - using Intersection_Kernel = Intersection_Traits; + using Kernel = typename Traits::Kernel; + using Intersection_Kernel = typename Traits::Intersection_Kernel; - using FT = typename Kernel::FT; - using Point_3 = typename Kernel::Point_3; - using Triangle_2 = typename Kernel::Triangle_2; - using Triangle_3 = typename Kernel::Triangle_3; + using FT = typename Traits::FT; + using Point_3 = typename Traits::Point_3; + using Triangle_2 = typename Traits::Triangle_2; + using Triangle_3 = typename Traits::Triangle_3; using Indices = std::vector; - using Data_structure = KSR_3::Data_structure; + using Data_structure = KSR_3::Data_structure; using Mesh = typename Data_structure::Mesh; using Volume_cell = typename Data_structure::Volume_cell; using PFace = typename Data_structure::PFace; @@ -60,8 +60,8 @@ namespace KSR_3 { using Delaunay_2 = CGAL::Delaunay_triangulation_2; using Delaunay_3 = CGAL::Delaunay_triangulation_3; - //using Converter = CGAL::Cartesian_converter; - using From_exact = CGAL::Cartesian_converter; + + using From_exact = typename Traits::From_exact; struct Wrapper { PFace pface; @@ -333,7 +333,7 @@ namespace KSR_3 { for (std::size_t i = 0; i < 6; i++) { cost_matrix[0][i] = 100000000; } - cost_matrix[0][0] = -100000000; + //cost_matrix[0][0] = -100000000; for (std::size_t i = 0; i < 6; i++) { cost_matrix[1][i] = -cost_matrix[0][i]; diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h index 6af2113f9c7..798061b12b1 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h @@ -28,7 +28,6 @@ #include #include #include -#include #include @@ -40,25 +39,25 @@ namespace KSR_3 { #ifdef DOXYGEN_RUNNING #else -template +template class Initializer { public: - using Kernel = GeomTraits; - using Intersection_Kernel = Intersection_Traits; + using Kernel = typename Traits::Kernel; + using Intersection_Kernel = typename Traits::Intersection_Kernel; private: - using FT = typename Kernel::FT; - using Point_2 = typename Kernel::Point_2; - using Point_3 = typename Kernel::Point_3; - using Vector_2 = typename Kernel::Vector_2; - using Segment_2 = typename Kernel::Segment_2; - using Segment_3 = typename Kernel::Segment_3; - using Line_2 = typename Kernel::Line_2; - using Transform_3 = typename Kernel::Aff_transformation_3; - using Direction_2 = typename Kernel::Direction_2; + using FT = typename Traits::FT; + using Point_2 = typename Traits::Point_2; + using Point_3 = typename Traits::Point_3; + using Vector_2 = typename Traits::Vector_2; + using Segment_2 = typename Traits::Segment_2; + using Segment_3 = typename Traits::Segment_3; + using Line_2 = typename Traits::Line_2; + using Transform_3 = typename Traits::Transform_3; + using Direction_2 = typename Traits::Direction_2; - using Data_structure = KSR_3::Data_structure; + using Data_structure = KSR_3::Data_structure; using Support_plane = typename Data_structure::Support_plane; using IEdge = typename Data_structure::IEdge; using IFace = typename Data_structure::IFace; @@ -68,23 +67,21 @@ private: using IVertex = typename Data_structure::IVertex; - using To_exact = CGAL::Cartesian_converter; - using From_exact = CGAL::Cartesian_converter; + using To_exact = typename Traits::To_exact; + using From_exact = typename Traits::From_exact; using Bbox_3 = CGAL::Bbox_3; using OBB_traits = CGAL::Oriented_bounding_box_traits_3; using Planar_shape_type = KSR::Planar_shape_type; using Parameters = KSR::Parameters_3; - using Kinetic_traits = KSR::Kinetic_traits_3; using Timer = CGAL::Real_timer; public: Initializer(Data_structure& data, const Parameters& parameters) : m_data(data), m_parameters(parameters), - m_merge_type(Planar_shape_type::CONVEX_HULL), - m_kinetic_traits() + m_merge_type(Planar_shape_type::CONVEX_HULL) { } template< @@ -98,8 +95,8 @@ public: std::array bbox; create_bounding_box( input_range, polygon_map, - m_parameters.enlarge_bbox_ratio, - m_parameters.reorient, bbox); + m_parameters.bbox_dilation_ratio, + m_parameters.reorient_bbox, bbox); const double time_to_bbox = timer.time(); @@ -179,7 +176,6 @@ private: Data_structure& m_data; const Parameters& m_parameters; const Planar_shape_type m_merge_type; - Kinetic_traits m_kinetic_traits; template< typename InputRange, @@ -1000,9 +996,9 @@ private: const auto& qj = merged[jp]; const Segment_2 segment(pj, qj); Point_2 inter; - const bool is_intersected = - m_kinetic_traits.intersection(segment, edge, inter); - if (is_intersected) return false; + const bool is_intersected = intersection(segment, edge, inter); + if (is_intersected) + return false; } } return true; @@ -1118,13 +1114,8 @@ private: typename Intersection_Kernel::Point_2 point; typename Intersection_Kernel::Segment_3 seg_a(m_data.point_3(it_a->second.first), m_data.point_3(it_a->second.second)); typename Intersection_Kernel::Segment_3 seg_b(m_data.point_3(it_b->second.first), m_data.point_3(it_b->second.second)); - if (!m_kinetic_traits.intersection( - m_data.to_2d(common_plane_idx, seg_a), - m_data.to_2d(common_plane_idx, seg_b), - point)) { - + if (!intersection(m_data.to_2d(common_plane_idx, seg_a), m_data.to_2d(common_plane_idx, seg_b), point)) continue; - } crossed_vertices.push_back( m_data.add_ivertex(m_data.to_3d(common_plane_idx, point), union_set)); @@ -1201,6 +1192,19 @@ private: } } } + + template + inline bool intersection( + const Type1& t1, const Type2& t2, ResultType& result) const { + + const auto inter = CGAL::intersection(t1, t2); + if (!inter) return false; + if (const ResultType* typed_inter = boost::get(&*inter)) { + result = *typed_inter; + return true; + } + return false; + } }; #endif //DOXYGEN_RUNNING diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Intersection_graph.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Intersection_graph.h index c2e9692b78e..44b320cf64e 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Intersection_graph.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Intersection_graph.h @@ -31,19 +31,18 @@ namespace KSR_3 { #ifdef DOXYGEN_RUNNING #else -template +template class Intersection_graph { public: - using Kernel = GeomTraits; - using Interxsection_Kernel = Intersection_Traits; + using Kernel = typename Traits::Kernel; + using Intersection_Kernel = typename Traits::Intersection_Kernel; - using FT = typename Intersection_Traits::FT; - using Point_2 = typename Intersection_Traits::Point_2; - using Point_3 = typename Intersection_Traits::Point_3; - using Segment_3 = typename Intersection_Traits::Segment_3; - using Line_3 = typename Intersection_Traits::Line_3; - using Polygon_2 = typename CGAL::Polygon_2; + using Point_2 = typename Traits::IK_Point_2; + using Point_3 = typename Traits::IK_Point_3; + using Segment_3 = typename Traits::IK_Segment_3; + using Line_3 = typename Traits::IK_Line_3; + using Polygon_2 = typename CGAL::Polygon_2; using Inexact_FT = typename Kernel::FT; @@ -106,7 +105,7 @@ public: Face_property(std::size_t support_plane_idx) : support_plane(support_plane_idx), part_of_partition(false) {} std::size_t support_plane; bool part_of_partition; - CGAL::Polygon_2 poly; + CGAL::Polygon_2 poly; std::vector pts; std::vector edges; std::vector vertices; @@ -371,7 +370,7 @@ public: void set_crossed(const Edge_descriptor& edge, std::size_t sp_idx) { m_graph[edge].crossed.insert(sp_idx); } }; -template std::size_t Intersection_graph::Edge_property::edge_counter = 0; +template std::size_t Intersection_graph::Edge_property::edge_counter = 0; #endif //DOXYGEN_RUNNING diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Reconstruction.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Reconstruction.h index 99a90bf535f..76bd9963dad 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Reconstruction.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Reconstruction.h @@ -46,8 +46,11 @@ namespace CGAL { \tparam IntersectionKernel must be a model of `Kernel`. Is used for the creation of the intersection graph. An exact kernel is suggested. + + \tparam InputRange + must be a model of `ConstRange` whose iterator type is `RandomAccessIterator`. */ -template +template class Kinetic_reconstruction_3 { public: /*! @@ -267,34 +270,30 @@ public: #else namespace KSR_3 { -template< -typename InputRange, -typename PointMap, -typename VectorMap, -typename SemanticMap, -typename GeomTraits, -typename Intersection_Kernel> +template class Reconstruction { public: - using Input_range = InputRange; - using Point_map = PointMap; - using Vector_map = VectorMap; + using Kernel = typename Traits::Kernel; + using Intersection_Kernel = typename Traits::Intersection_Kernel; + using Input_range = typename Traits::Input_range; + using Point_map = typename Traits::Point_map; + using Vector_map = typename Traits::Normal_map; using Semantic_map = SemanticMap; - using Kernel = GeomTraits; private: - using FT = typename Kernel::FT; - using Point_2 = typename Kernel::Point_2; - using Line_2 = typename Kernel::Line_2; - using Point_3 = typename Kernel::Point_3; - using Plane_3 = typename Kernel::Plane_3; - using Vector_2 = typename Kernel::Vector_2; - using Vector_3 = typename Kernel::Vector_3; - using Segment_2 = typename Kernel::Segment_2; - using Segment_3 = typename Kernel::Segment_3; + using FT = typename Traits::FT; + using Point_2 = typename Traits::Point_2; + using Line_2 = typename Traits::Line_2; + using Point_3 = typename Traits::Point_3; + using Plane_3 = typename Traits::Plane_3; + using Vector_2 = typename Traits::Vector_2; + using Vector_3 = typename Traits::Vector_3; + using Segment_2 = typename Traits::Segment_2; + using Segment_3 = typename Traits::Segment_3; - using Data_structure = KSR_3::Data_structure; + using Data_structure = KSR_3::Data_structure; using PFace = typename Data_structure::PFace; using Point_map_3 = KSR::Item_property_map; @@ -309,7 +308,7 @@ private: using Polygon_3 = std::vector; using Polygon_map = CGAL::Identity_property_map; - using From_EK = CGAL::Cartesian_converter; + using From_EK = typename Traits::From_exact; struct Vertex_info { FT z = FT(0); }; struct Face_info { }; @@ -360,8 +359,8 @@ private: Region_growing; using Visibility_label = KSR::Visibility_label; - using Visibility = KSR_3::Visibility; - using Graphcut = KSR_3::Graphcut; + using Visibility = KSR_3::Visibility; + using Graphcut = KSR_3::Graphcut; public: @@ -954,7 +953,7 @@ private: for (const auto& indices : result) { region.clear(); for (const std::size_t index : indices) { - region.push_back(input_range[index]); + region.push_back(index); } regions.push_back(region); } diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Support_plane.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Support_plane.h index 616b3a0d7cb..3bf3ebe97e4 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Support_plane.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Support_plane.h @@ -29,30 +29,30 @@ namespace KSR_3 { #ifdef DOXYGEN_RUNNING #else -template +template class Support_plane { public: - using Kernel = GeomTraits; - using Intersection_Kernel = Intersection_Traits; + using Kernel = typename Traits::Kernel; + using Intersection_Kernel = typename Traits::Intersection_Kernel; using To_exact = CGAL::Cartesian_converter; using From_exact = CGAL::Cartesian_converter; - using FT = typename Kernel::FT; - using Point_2 = typename Kernel::Point_2; - using Point_3 = typename Kernel::Point_3; - using Vector_2 = typename Kernel::Vector_2; - using Vector_3 = typename Kernel::Vector_3; - using Direction_2 = typename Kernel::Direction_2; - using Segment_2 = typename Kernel::Segment_2; - using Segment_3 = typename Kernel::Segment_3; - using Line_2 = typename Kernel::Line_2; - using Line_3 = typename Kernel::Line_3; - using Plane_3 = typename Kernel::Plane_3; - using Triangle_2 = typename Kernel::Triangle_2; + using FT = typename Traits::FT; + using Point_2 = typename Traits::Point_2; + using Point_3 = typename Traits::Point_3; + using Vector_2 = typename Traits::Vector_2; + using Vector_3 = typename Traits::Vector_3; + using Direction_2 = typename Traits::Direction_2; + using Segment_2 = typename Traits::Segment_2; + using Segment_3 = typename Traits::Segment_3; + using Line_2 = typename Traits::Line_2; + using Line_3 = typename Traits::Line_3; + using Plane_3 = typename Traits::Plane_3; + using Triangle_2 = typename Traits::Triangle_2; using Mesh = CGAL::Surface_mesh; - using Intersection_graph = KSR_3::Intersection_graph; + using Intersection_graph = KSR_3::Intersection_graph; using Bbox_2 = CGAL::Bbox_2; using IVertex = typename Intersection_graph::Vertex_descriptor; @@ -113,8 +113,11 @@ private: std::vector original_vectors; std::vector original_directions; std::vector original_rays; - int k; + FT distance_tolerance; + FT angle_tolerance; + + int k; }; std::shared_ptr m_data; @@ -126,7 +129,7 @@ public: } template - Support_plane(const PointRange& polygon, const bool is_bbox, const FT distance_tolerance, std::size_t idx) : + Support_plane(const PointRange& polygon, const bool is_bbox, const FT distance_tolerance, const FT angle_tolerance, std::size_t idx) : m_data(std::make_shared()) { To_exact to_EK; @@ -159,6 +162,7 @@ public: m_data->exact_plane = to_EK(m_data->plane); m_data->is_bbox = is_bbox; m_data->distance_tolerance = distance_tolerance; + m_data->angle_tolerance = angle_tolerance; std::vector tris(points.size() - 2); for (std::size_t i = 2; i < points.size(); i++) { @@ -279,6 +283,10 @@ public: return m_data->distance_tolerance; } + FT angle_tolerance() const { + return m_data->angle_tolerance; + } + void clear_pfaces() { m_data->mesh.clear(); add_property_maps(); @@ -756,14 +764,14 @@ public: } }; -template -bool operator==(const Support_plane& a, const Support_plane& b) { +template +bool operator==(const Support_plane& a, const Support_plane& b) { if (a.is_bbox() || b.is_bbox()) { return false; } - using FT = typename Kernel::FT; + using FT = typename Traits::FT; const auto& planea = a.plane(); const auto& planeb = b.plane(); @@ -779,24 +787,15 @@ bool operator==(const Support_plane& a, const Suppo // return false; // } - const FT vtol = FT(5); // degrees // TODO: We should put it as a parameter. FT aval = approximate_angle(va, vb); CGAL_assertion(aval >= FT(0) && aval <= FT(180)); - if (aval >= FT(90)) aval = FT(180) - aval; + if (aval >= FT(90)) + aval = FT(180) - aval; - // std::cout << "aval: " << aval << " : " << vtol << std::endl; - if (aval >= vtol) { + if (aval >= a.angle_tolerance()) { return false; } - // Are the planes coplanar? - // const FT ptol = KSR::point_tolerance(); - // const auto pa = planea.point(); - // const auto pb = planeb.projection(pa); - // const FT bval = KSR::distance(pa, pb); - - // TODO: We should put it as a parameter. - const FT ptol = a.distance_tolerance(); const auto pa1 = a.to_3d(a.centroid()); const auto pb1 = planeb.projection(pa1); const auto pb2 = b.to_3d(b.centroid()); @@ -807,14 +806,9 @@ bool operator==(const Support_plane& a, const Suppo const FT bval = (CGAL::max)(bval1, bval2); CGAL_assertion(bval >= FT(0)); - // if (bval < ptol) { - // std::cout << "2 " << pa << " " << pb << std::endl; - // std::cout << "bval: " << bval << " : " << ptol << std::endl; - // } + if (bval >= a.distance_tolerance()) + return false; - // std::cout << "bval: " << bval << " : " << ptol << std::endl; - if (bval >= ptol) return false; - // std::cout << "- found coplanar planes" << std::endl; return true; } diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Visibility.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Visibility.h index 1612ff4a0c6..c2087dd7743 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Visibility.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Visibility.h @@ -36,47 +36,42 @@ namespace KSR_3 { #ifdef DOXYGEN_RUNNING #else - template< - typename GeomTraits, - typename Intersection_Traits, - typename PointMap_3, - typename VectorMap_3> + template class Visibility { public: - using Kernel = GeomTraits; - using Intersection_Kernel = Intersection_Traits; - using Point_map_3 = PointMap_3; - using Vector_map_3 = VectorMap_3; + using Kernel = typename Traits::Kernel; + using Intersection_Kernel = typename Traits::Intersection_Kernel; + using Point_map_3 = typename Traits::Point_map; + using Vector_map_3 = typename Traits::Normal_map; - using FT = typename Kernel::FT; - using Point_3 = typename Kernel::Point_3; - using Vector_3 = typename Kernel::Vector_3; + using FT = typename Traits::FT; + using Point_3 = typename Traits::Point_3; + using Vector_3 = typename Traits::Vector_3; using Indices = std::vector; - using Data_structure = KSR_3::Data_structure; + using Data_structure = KSR_3::Data_structure; using PFace = typename Data_structure::PFace; using Volume_cell = typename Data_structure::Volume_cell; using Delaunay_3 = CGAL::Delaunay_triangulation_3; using Generator = CGAL::Random_points_in_tetrahedron_3; - using From_EK = CGAL::Cartesian_converter; + using From_EK = typename Traits::From_exact; using Visibility_label = KSR::Visibility_label; Visibility( const Data_structure& data, - const std::map& pface_points, + const std::map& face2points, const Point_map_3& point_map_3, const Vector_map_3& normal_map_3) : m_data(data), - m_pface_points(pface_points), + m_face2points(face2points), m_point_map_3(point_map_3), m_normal_map_3(normal_map_3), - m_num_samples(0), - m_random(0) { - CGAL_assertion(m_pface_points.size() > 0); + m_num_samples(0) { + CGAL_assertion(m_face2points.size() > 0); m_inliers = 0; for (const auto pair : m_pface_points) { m_inliers += pair.second.size(); @@ -116,12 +111,11 @@ namespace KSR_3 { private: const Data_structure& m_data; - const std::map& m_pface_points; + const std::map& m_face2points; const Point_map_3& m_point_map_3; const Vector_map_3& m_normal_map_3; mutable std::size_t m_num_samples; mutable std::size_t m_inliers; - Random m_random; void estimate_volume_label(Volume_cell& volume) const { diff --git a/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_partitioning_3.h b/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_partitioning_3.h new file mode 100644 index 00000000000..f1ce1f01f7e --- /dev/null +++ b/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_partitioning_3.h @@ -0,0 +1,1000 @@ +// Copyright (c) 2019 GeometryFactory SARL (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// +// Author(s) : Simon Giraudot, Dmitry Anisimov, Sven Oesau + +#ifndef CGAL_KINETIC_SHAPE_PARTITIONING_3_H +#define CGAL_KINETIC_SHAPE_PARTITIONING_3_H + +// #include + +#include"Kinetic_shape_partitioning_Traits.h" + +// Boost includes. +#include +#include + +// CGAL includes. +#include +#include +#include +#include +#include + +#include +#include +#include + +// Internal includes. +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace CGAL { + +#ifdef DOXYGEN_NS +/*! +* \ingroup PkgKineticPartition + \brief Creates the kinetic partitioning of the bounding box. + + \tparam Kernel + must be a model of `Kernel`. Is used for non-critical calculations. + + \tparam IntersectionKernel + must be a model of `Kernel`. Is used for the creation of the intersection graph. An exact kernel is suggested. +*/ +template +class Kinetic_partitioning_3 { +public: + + + /*! + \brief Propagates the kinetic polygons in the initialized partition. + + \param k + maximum number of allowed intersections for each input polygon before its expansion stops. + + @return + success of kinetic partitioning. + + \pre `successful initialization` + */ + /// @} + + bool partition(std::size_t k); + +} +#else + +/*! +* \ingroup PkgKineticPartition + \brief Creates the kinetic partitioning of the bounding box. + + \tparam Kernel + must be a model of `Kernel`. Is used for non-critical calculations. + + \tparam IntersectionKernel + must be a model of `Kernel`. Is used for the creation of the intersection graph. An exact kernel is suggested. +*/ +template +class Kinetic_shape_partitioning_3 { + +public: + using Kernel = typename Traits::Kernel; + using Intersection_Kernel = typename Traits::Intersection_Kernel; + + using Point_3 = typename Traits::Point_3; + +private: + using FT = typename Traits::FT; + + using Data_structure = KSR_3::Data_structure; + + using IVertex = typename Data_structure::IVertex; + using IEdge = typename Data_structure::IEdge; + + using From_exact = typename Traits::From_exact; + + using Initializer = KSR_3::Initializer; + using Propagation = KSR_3::FacePropagation; + using Finalizer = KSR_3::Finalizer; + + using Polygon_mesh = CGAL::Surface_mesh; + using Timer = CGAL::Real_timer; + using Parameters = KSR::Parameters_3; + +private: + Parameters m_parameters; + Data_structure m_data; + std::size_t m_num_events; + +public: + /// \name Initialization + /// @{ + /*! + \brief Constructs an empty kinetic shape partitioning object. + + \param verbose + prints information about the creation into std::cout. + + \param debug + writes intermediate results into files. + + */ + Kinetic_shape_partitioning_3( + const bool verbose = true, + const bool debug = false) : + m_parameters(verbose, debug), // use true here to export all steps + m_data(m_parameters), + m_num_events(0) + { } + + /*! + \brief Initializes the kinetic partitioning of the bounding box. + + \tparam InputRange + must be a model of `ConstRange` whose iterator type is `RandomAccessIterator`. + + \tparam PolygonMap + contains index ranges to form polygons from InputRange + + \tparam NamedParameters + a sequence of \ref bgl_namedparameters "Named Parameters" + + \param input_range + an instance of `InputRange` with 3D points and corresponding 3D normal vectors + + \param polygon_map + a range of polygons defined by a range of indices into input_range + + \param np + a sequence of \ref bgl_namedparameters "Named Parameters" + among the ones listed below + + @return + success. The initialization fails if the input data is empty. + + \pre `input_range.size() > 0 and polygon_map.size() > 0` + + \cgalNamedParamsBegin + \cgalParamNBegin{reorient_bbox} + \cgalParamDescription{Use the oriented bounding box instead of the axis-aligned bounding box.} + \cgalParamType{bool} + \cgalParamDefault{false} + \cgalParamNEnd + \cgalParamNBegin{bbox_dilation_ratio} + \cgalParamDescription{Factor for extension of the bounding box of the input data to be used for the partitioning.} + \cgalParamType{FT} + \cgalParamDefault{1.1} + \cgalParamNEnd + \cgalParamNBegin{angle_tolerance} + \cgalParamDescription{The tolerance angle to snap the planes of two input polygons into one plane.} + \cgalParamType{FT} + \cgalParamDefault{5} + \cgalParamNEnd + \cgalParamNBegin{distance_tolerance} + \cgalParamDescription{The tolerance distance to snap the planes of two input polygons into one plane.} + \cgalParamType{FT} + \cgalParamDefault{0.5} + \cgalParamNEnd + \cgalNamedParamsEnd + */ + + template< + typename InputRange, + typename PolygonMap, + typename NamedParameters = parameters::Default_named_parameters> + bool initialize( + const InputRange& input_range, + const PolygonMap polygon_map, + const NamedParameters& np = CGAL::parameters::default_values()) { + + Timer timer; + m_parameters.bbox_dilation_ratio = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::bbox_dilation_ratio), FT(11) / FT(10)); + m_parameters.angle_tolerance = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::angle_tolerance), FT(5) / FT(10)); + m_parameters.distance_tolerance = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::distance_tolerance), FT(5) / FT(10)); + m_parameters.reorient_bbox = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::reorient_bbox), false); + + std::cout.precision(20); + if (input_range.size() == 0) { + CGAL_warning_msg(input_range.size() > 0, + "Warning: Your input is empty!"); + return false; + } + + if (m_parameters.bbox_dilation_ratio < FT(1)) { + CGAL_warning_msg(m_parameters.bbox_dilation_ratio >= FT(1), + "Warning: You set enlarge_bbox_ratio < 1.0! The valid range is [1.0, +inf). Setting to 1.0!"); + m_parameters.bbox_dilation_ratio = FT(1); + } + + if (m_parameters.verbose) { + const unsigned int num_blocks = static_cast(std::pow(m_parameters.n + 1, 3)); + //const std::string is_reorient = (m_parameters.reorient ? "true" : "false"); + + std::cout << std::endl << "--- PARTITION OPTIONS: " << std::endl; + std::cout << "* enlarge bbox ratio: " << m_parameters.bbox_dilation_ratio << std::endl; + } + + if (m_parameters.verbose) { + std::cout << std::endl << "--- INITIALIZING PARTITION:" << std::endl; + + // Initialization. + timer.reset(); + timer.start(); + } + + m_data.clear(); + + Initializer initializer(m_data, m_parameters); + initializer.initialize(input_range, polygon_map); + + // Timing. + if (m_parameters.verbose) { + timer.stop(); + const double time_to_initialize = timer.time(); + std::cout << "* initialization time: " << time_to_initialize << std::endl; + } + /* + + if (m_parameters.k == 0) { // for k = 0, we skip propagation + CGAL_warning_msg(m_parameters.k > 0, + "WARNING: YOU SET K TO 0! THAT MEANS NO PROPAGATION! THE VALID VALUES ARE {1,2,...}. INTERSECT AND RETURN!"); + return false; + } + + if (m_parameters.verbose) { + std::cout << std::endl << "--- RUNNING THE QUEUE:" << std::endl; + std::cout << "* propagation started" << std::endl; + } + + // Propagation. + timer.reset(); + timer.start(); + std::size_t num_queue_calls = 0; + + Propagation propagation(m_data, m_parameters); + std::tie(num_queue_calls, m_num_events) = propagation.propagate(); + + timer.stop(); + const double time_to_propagate = timer.time(); + + if (m_parameters.verbose) { + std::cout << "* propagation finished" << std::endl; + std::cout << "* number of queue calls: " << num_queue_calls << std::endl; + std::cout << "* number of events handled: " << m_num_events << std::endl; + } + + if (m_parameters.verbose) { + std::cout << std::endl << "--- FINALIZING PARTITION:" << std::endl; + } + + // Finalization. + timer.reset(); + timer.start(); + if (m_parameters.debug) + dump(m_data, "final-" + m_parameters.k); + + Finalizer finalizer(m_data, m_parameters); + + if (m_parameters.verbose) + std::cout << "* checking final mesh integrity ..."; + + CGAL_assertion(m_data.check_integrity(true, true, true)); + + if (m_parameters.verbose) + std::cout << " done" << std::endl; + + if (m_parameters.verbose) + std::cout << "* getting volumes ..." << std::endl; + + finalizer.create_polyhedra(); + timer.stop(); + const double time_to_finalize = timer.time(); + + if (m_parameters.verbose) { + std::cout << "* found all together " << m_data.number_of_volumes() << " volumes" << std::endl; + + for (std::size_t i = 0; i < m_data.number_of_support_planes(); i++) { + dump_2d_surface_mesh(m_data, i, "final-surface-mesh-" + std::to_string(i)); + } + }*/ + + //std::cout << "* propagation: " << time_to_propagate << std::endl; + //std::cout << "* finalization: " << time_to_finalize << std::endl; + //std::cout << "* total time: " << total_time << std::endl; + + return true; + } + /*! + \brief Propagates the kinetic polygons in the initialized partition. + + \param k + maximum number of allowed intersections for each input polygon before its expansion stops. + + @return + success of kinetic partitioning. It can fail if the input data is empty, the partitioning has not been initialized or k is 0. + + \pre `successful initialization` + */ + bool partition(std::size_t k) { + Timer timer; + std::cout.precision(20); + + // Already initialized? + if (m_data.number_of_support_planes() < 6) { + std::cout << "Kinetic partitioning not initialized or empty. Number of support planes: " << m_data.number_of_support_planes() << std::endl; + + return false; + } + + if (k == 0) { // for k = 0, we skip propagation + std::cout << "k needs to be a positive number" << std::endl; + + return false; + } + + if (m_parameters.verbose) { + std::cout << std::endl << "--- RUNNING THE QUEUE:" << std::endl; + std::cout << "* propagation started" << std::endl; + } + + // Propagation. + timer.reset(); + timer.start(); + std::size_t num_queue_calls = 0; + + Propagation propagation(m_data, m_parameters); + std::tie(num_queue_calls, m_num_events) = propagation.propagate(k); + + timer.stop(); + const double time_to_propagate = timer.time(); + + if (m_parameters.verbose) { + std::cout << "* propagation finished" << std::endl; + std::cout << "* number of queue calls: " << num_queue_calls << std::endl; + std::cout << "* number of events handled: " << m_num_events << std::endl; + } + + if (m_parameters.verbose) { + std::cout << std::endl << "--- FINALIZING PARTITION:" << std::endl; + } + + // Finalization. + timer.reset(); + timer.start(); + if (m_parameters.debug) + dump(m_data, "final-" + m_parameters.k); + + Finalizer finalizer(m_data, m_parameters); + + if (m_parameters.verbose) + std::cout << "* checking final mesh integrity ..."; + + CGAL_assertion(m_data.check_integrity(true, true, true)); + + if (m_parameters.verbose) + std::cout << " done" << std::endl; + + if (m_parameters.verbose) + std::cout << "* getting volumes ..." << std::endl; + + finalizer.create_polyhedra(); + timer.stop(); + const double time_to_finalize = timer.time(); + + if (m_parameters.verbose) + std::cout << "* found all together " << m_data.number_of_volumes() << " volumes" << std::endl; + + if (m_parameters.debug) + for (std::size_t i = 0; i < m_data.number_of_support_planes(); i++) + dump_2d_surface_mesh(m_data, i, "final-surface-mesh-" + std::to_string(i)); + + // Timing. + if (m_parameters.verbose) { + std::cout << std::endl << "--- TIMING (sec.):" << std::endl; + + std::cout << "* propagation: " << time_to_propagate << std::endl; + std::cout << "* finalization: " << time_to_finalize << std::endl; + } + + return true; + } + + /// @} + /* + template< + typename InputRange, + typename PointMap, + typename VectorMap, + typename SemanticMap, + typename NamedParameters> + bool reconstruct( + const InputRange& input_range, + const PointMap point_map, + const VectorMap normal_map, + const SemanticMap semantic_map, + const std::string file_name, + const NamedParameters& np) { + + using Reconstruction = KSR_3::Reconstruction< + InputRange, PointMap, VectorMap, SemanticMap, Kernel, Intersection_Kernel>; + + Reconstruction reconstruction( + input_range, point_map, normal_map, semantic_map, m_data, m_parameters.verbose, m_parameters.debug); + + bool success = reconstruction.detect_planar_shapes(file_name, np); + if (!success) { + CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, DETECTING PLANAR SHAPES FAILED!"); + return false; + } + // exit(EXIT_SUCCESS); + + success = reconstruction.regularize_planar_shapes(np); + if (!success) { + CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, REGULARIZATION FAILED!"); + return false; + } + // exit(EXIT_SUCCESS); + + success = partition( + reconstruction.planar_shapes(), reconstruction.polygon_map(), np); + if (!success) { + CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, PARTITION FAILED!"); + return false; + } + // exit(EXIT_SUCCESS); + + success = reconstruction.compute_model(np); + if (!success) { + CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, COMPUTING MODEL FAILED!"); + return false; + } + return success; + } + + template< + typename InputRange, + typename PointMap, + typename VectorMap, + typename SemanticMap, + typename RegionMap, + typename NamedParameters> + bool reconstruct( + const InputRange& input_range, + const PointMap point_map, + const VectorMap normal_map, + const SemanticMap semantic_map, + const RegionMap region_map, + const NamedParameters& np) { + + using Reconstruction = KSR_3::Reconstruction< + InputRange, PointMap, VectorMap, SemanticMap, Kernel, Intersection_Kernel>; + + Reconstruction reconstruction( + input_range, point_map, normal_map, semantic_map, m_data, m_parameters.verbose, m_parameters.debug); + + bool success = reconstruction.planar_shapes_from_map(region_map, np); + if (!success) { + CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, DETECTING PLANAR SHAPES FAILED!"); + return false; + } + // exit(EXIT_SUCCESS); + + //success = reconstruction.regularize_planar_shapes(np); + if (!success) { + CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, REGULARIZATION FAILED!"); + return false; + } + // exit(EXIT_SUCCESS); + + success = partition( + reconstruction.planar_shapes(), reconstruction.polygon_map(), np); + if (!success) { + CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, PARTITION FAILED!"); + return false; + } + // exit(EXIT_SUCCESS); + + success = reconstruction.compute_model(np); + if (!success) { + CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, COMPUTING MODEL FAILED!"); + return false; + } + return success; + } +*/ + + /******************************* + ** Access ** + ********************************/ + + /// \name Access + /// @{ + /*! + \brief Number of support planes in the kinetic partitioning. They originate from the planes of the input polygons and the bounding box. + + @return + number of vertices. + + \pre `successful partitioning` + */ + + int number_of_support_planes() const { + return static_cast(m_data.number_of_support_planes()); + } + + /*! + \brief Number of vertices in the kinetic partitioning. + + @return + number of vertices. + + \pre `successful partitioning` + */ + std::size_t number_of_vertices() const { + return 0; + } + + /*! + \brief Number of convex faces in the kinetic partitioning. + + @return + number of convex faces. + + \pre `successful partitioning` + */ + std::size_t number_of_faces() const { + return 0; + } + + /*! + \brief Number of convex volumes created by the kinetic partitioning. + + @return + number of convex volumes. + + \pre `successful partitioning` + */ + std::size_t number_of_volumes() const { + return m_data.volumes().size(); + } + + /*! + \brief Point vector for mapping vertex indices to positions. + + @return + vector of points. + + \pre `successful partitioning` + */ + const std::vector& vertices() const; + + /*! + \brief Vertex indices of convex face. + + \param face_index + index of the query face. + + @return + vector of vertex indices. + + \pre `successful partitioning` + */ + const std::vector& vertices(std::size_t face_index) const; + + /*! + \brief Face indices of the convex volume. + + \param volume_index + index of the query volume. + + @return + vector of face indices. + + \pre `successful partitioning` + */ + const std::vector& face(std::size_t volume_index) const { + CGAL_assertion(m_data.number_of_volumes() > volume_index); + return m_data.volumes()[volume_index].faces; + } + + /*! + \brief Indices of adjacent volumes. Negative indices correspond to the empty spaces behind the sides of the bounding box. + + \param face_index + index of the query face. + + @return + pair of adjacent volumes. + + -1 zmin + -2 ymin + -3 xmax + -4 ymax + -5 xmin + -6 zmax + + \pre `successful partitioning` + */ + const std::pair& neighbors(std::size_t face_index) const { + CGAL_assertion(m_data.number_of_volumes() > volume_index); + return m_data. + } + + /*! + \brief Retrieves the support plane generated from the input polygon. + + \param input_polygon_index + index of the input polygon. + + @return + index into polygon_map provided on initialization. + + \pre `successful partitioning` + */ + std::size_t support_plane_index(const std::size_t input_polygon_index) const { + const int support_plane_idx = m_data.support_plane_index(input_polygon_index); + CGAL_assertion(support_plane_idx >= 6); + return support_plane_idx; + } + + /*! + \brief Creates a linear cell complex from the kinetic partitioning. + + \tparam LCC + Linear_cell_complex_for_combinatorial_map<3, 3,...> + The dimension of the combinatorial map and the dimension of the ambient space have to be 3. + + \param lcc + an instance of LCC + + \pre `successful partitioning` + */ + + template + void get_linear_cell_complex(LCC& lcc) const { + using LCC_Kernel = typename LCC::Traits; + CGAL::Cartesian_converter conv; + lcc.clear(); + + std::vector used_vertices(m_data.igraph().number_of_vertices(), false); + std::vector remap(m_data.igraph().number_of_vertices(), -1); + std::vector mapped_vertices; + mapped_vertices.reserve(m_data.igraph().number_of_vertices()); + + for (const auto& volume : m_data.volumes()) { + for (const auto& vertex : volume.pvertices) { + CGAL_assertion(m_data.has_ivertex(vertex)); + IVertex ivertex = m_data.ivertex(vertex); + if (remap[ivertex] == -1) { + remap[ivertex] = static_cast(mapped_vertices.size()); + mapped_vertices.push_back(conv(m_data.point_3(ivertex))); + } + } + } + + CGAL::Linear_cell_complex_incremental_builder_3 ib(lcc); + for (const auto& p : mapped_vertices) + ib.add_vertex(p); + + for (const auto& vol : m_data.volumes()) { + ib.begin_surface(); + for (std::size_t i = 0; i < vol.pfaces.size(); i++) { + auto vertex_range = m_data.pvertices_of_pface(vol.pfaces[i]); + ib.begin_facet(); + if (vol.pface_oriented_outwards[i]) { + typename Data_structure::PVertex_of_pface_iterator it = vertex_range.begin(); + while (it != vertex_range.end()) { + CGAL_assertion(m_data.has_ivertex(*it)); + IVertex ivertex = m_data.ivertex(*it); + ib.add_vertex_to_facet(static_cast(remap[ivertex])); + it++; + } + } + else { + typename Data_structure::PVertex_of_pface_iterator it = vertex_range.end()--; + do { + CGAL_assertion(m_data.has_ivertex(*it)); + IVertex ivertex = m_data.ivertex(*it); + ib.add_vertex_to_facet(static_cast(remap[ivertex])); + it--; + if (it == vertex_range.begin()) + break; + } while (true); + } + ib.end_facet(); + } + ib.end_surface(); + } + + lcc.display_characteristics(std::cout) << std::endl; + } + + /// @} + + /* + template + VertexOutputIterator output_partition_vertices( + VertexOutputIterator vertices, const int support_plane_idx = -1) const { + From_exact from_EK; + + CGAL_assertion(support_plane_idx < number_of_support_planes()); + if (support_plane_idx >= number_of_support_planes()) return vertices; + if (support_plane_idx < 0) { + const auto all_ivertices = m_data.ivertices(); + for (const auto ivertex : all_ivertices) { + *(vertices++) = from_EK(m_data.point_3(ivertex)); + } + return vertices; + } + + CGAL_assertion(support_plane_idx >= 0); + const std::size_t sp_idx = static_cast(support_plane_idx); + const auto all_pvertices = m_data.pvertices(sp_idx); + for (const auto pvertex : all_pvertices) { + CGAL_assertion(m_data.has_ivertex(pvertex)); + const auto ivertex = m_data.ivertex(pvertex); + *(vertices++) = from_EK(m_data.point_3(ivertex)); + } + return vertices; + } + + template + EdgeOutputIterator output_partition_edges( + EdgeOutputIterator edges, const int support_plane_idx = -1) const { + From_exact from_EK; + + CGAL_assertion(support_plane_idx < number_of_support_planes()); + if (support_plane_idx >= number_of_support_planes()) return edges; + if (support_plane_idx < 0) { + const auto all_iedges = m_data.iedges(); + for (const auto iedge : all_iedges) { + *(edges++) = from_EK(m_data.segment_3(iedge)); + } + return edges; + } + + CGAL_assertion(support_plane_idx >= 0); + const std::size_t sp_idx = static_cast(support_plane_idx); + const auto all_pedges = m_data.pedges(sp_idx); + for (const auto pedge : all_pedges) { + CGAL_assertion(m_data.has_iedge(pedge)); + const auto iedge = m_data.iedge(pedge); + *(edges++) = from_EK(m_data.segment_3(iedge)); + } + return edges; + } + + template + FaceOutputIterator output_partition_faces( + FaceOutputIterator faces, + const int support_plane_idx = -1, + const int begin = 0) const { + + KSR::Indexer indexer; + CGAL_assertion(support_plane_idx < number_of_support_planes()); + if (support_plane_idx >= number_of_support_planes()) return faces; + if (support_plane_idx < 0) { + const auto all_ivertices = m_data.ivertices(); + for (const auto ivertex : all_ivertices) indexer(ivertex); + for (int i = begin; i < number_of_support_planes(); ++i) { + const std::size_t sp_idx = static_cast(i); + output_partition_faces(faces, indexer, sp_idx); + } + return faces; + } + + CGAL_assertion(support_plane_idx >= 0); + const std::size_t sp_idx = static_cast(support_plane_idx); + const auto all_pvertices = m_data.pvertices(sp_idx); + for (const auto pvertex : all_pvertices) { + CGAL_assertion(m_data.has_ivertex(pvertex)); + const auto ivertex = m_data.ivertex(pvertex); + indexer(ivertex); + } + return output_partition_faces(faces, indexer, sp_idx); + } + + void output_support_plane( + Polygon_mesh& polygon_mesh, const int support_plane_idx) const { + From_exact from_EK; + + polygon_mesh.clear(); + CGAL_assertion(support_plane_idx >= 0); + if (support_plane_idx < 0) return; + CGAL_assertion(support_plane_idx < number_of_support_planes()); + if (support_plane_idx >= number_of_support_planes()) return; + const std::size_t sp_idx = static_cast(support_plane_idx); + + std::vector vertices; + std::vector map_vertices; + + map_vertices.clear(); + const auto all_pvertices = m_data.pvertices(sp_idx); + for (const auto pvertex : all_pvertices) { + CGAL_assertion(m_data.has_ivertex(pvertex)); + const auto ivertex = m_data.ivertex(pvertex); + + if (map_vertices.size() <= pvertex.second) + map_vertices.resize(pvertex.second + 1); + map_vertices[pvertex.second] = + polygon_mesh.add_vertex(from_EK(m_data.point_3(ivertex))); + } + + const auto all_pfaces = m_data.pfaces(sp_idx); + for (const auto pface : all_pfaces) { + vertices.clear(); + const auto pvertices = m_data.pvertices_of_pface(pface); + for (const auto pvertex : pvertices) { + vertices.push_back(map_vertices[pvertex.second]); + } + polygon_mesh.add_face(vertices); + } + } + + template + VolumeOutputIterator output_partition_volumes( + VolumeOutputIterator volumes) const { + for (std::size_t i = 0; i < m_data.number_of_volumes(); ++i) { + output_partition_volume(volumes, i); + } + return volumes; + } + + template + VolumeOutputIterator output_partition_volume( + VolumeOutputIterator volumes, const std::size_t volume_index) const { + + CGAL_assertion(volume_index < number_of_volumes()); + if (volume_index >= number_of_volumes()) return volumes; + + std::vector vertices; + std::vector< std::vector > faces; + output_partition_volume( + std::back_inserter(vertices), std::back_inserter(faces), volume_index); + CGAL::Polygon_mesh_processing::orient_polygon_soup(vertices, faces); + + Polygon_mesh polygon_mesh; + CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh( + vertices, faces, polygon_mesh); + *(volumes++) = polygon_mesh; + return volumes; + } + + template + void output_partition_volume( + VertexOutputIterator vertices, FaceOutputIterator faces, + const std::size_t volume_index) const { + From_exact from_EK; + + CGAL_assertion(volume_index < number_of_volumes()); + if (volume_index >= number_of_volumes()) return; + + const auto& volume = m_data.volumes()[volume_index]; + + std::size_t num_vertices = 0; + KSR::Indexer indexer; + + std::vector face; + const auto& pfaces = volume.pfaces; + for (const auto& pface : pfaces) { + face.clear(); + const auto pvertices = m_data.pvertices_of_pface(pface); + for (const auto pvertex : pvertices) { + + CGAL_assertion(m_data.has_ivertex(pvertex)); + const auto ivertex = m_data.ivertex(pvertex); + const std::size_t idx = indexer(ivertex); + + if (idx == num_vertices) { + *(vertices++) = from_EK(m_data.point_3(ivertex)); + ++num_vertices; + } + face.push_back(idx); + } + *(faces++) = face; + } + } + + template + void output_reconstructed_model( + VertexOutputIterator vertices, FaceOutputIterator faces) const { + From_exact from_EK; + + const auto& model = m_data.reconstructed_model(); + CGAL_assertion(model.pfaces.size() > 0); + + std::size_t num_vertices = 0; + KSR::Indexer indexer; + + std::vector face; + const auto& pfaces = model.pfaces; + for (const auto& pface : pfaces) { + face.clear(); + const auto pvertices = m_data.pvertices_of_pface(pface); + for (const auto pvertex : pvertices) { + + CGAL_assertion(m_data.has_ivertex(pvertex)); + const auto ivertex = m_data.ivertex(pvertex); + const std::size_t idx = indexer(ivertex); + + if (idx == num_vertices) { + *(vertices++) = from_EK(m_data.point_3(ivertex)); + ++num_vertices; + } + face.push_back(idx); + } + *(faces++) = face; + } + } + + void output_reconstructed_model(Polygon_mesh& polygon_mesh) const { + + std::vector vertices; + std::vector< std::vector > faces; + output_reconstructed_model( + std::back_inserter(vertices), std::back_inserter(faces)); + CGAL::Polygon_mesh_processing::orient_polygon_soup(vertices, faces); + polygon_mesh.clear(); + CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh( + vertices, faces, polygon_mesh); + } +*/ + + /******************************* + ** MEMORY ** + ********************************/ + + void clear() { + m_data.clear(); + m_num_events = 0; + } + +private: + /* + + template + FaceOutputIterator output_partition_faces( + FaceOutputIterator faces, KSR::Indexer& indexer, + const std::size_t sp_idx) const { + + std::vector face; + const auto all_pfaces = m_data.pfaces(sp_idx); + for (const auto pface : all_pfaces) { + face.clear(); + const auto pvertices = m_data.pvertices_of_pface(pface); + for (const auto pvertex : pvertices) { + CGAL_assertion(m_data.has_ivertex(pvertex)); + const auto ivertex = m_data.ivertex(pvertex); + const std::size_t idx = indexer(ivertex); + face.push_back(idx); + } + *(faces++) = face; + } + return faces; + }*/ +}; + +#endif //DOXYGEN_NS + +} // namespace CGAL + +#endif // CGAL_KINETIC_SHAPE_PARTITIONING_3_H diff --git a/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h b/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h index 91456c789c2..623e4a7498a 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h +++ b/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h @@ -1,4 +1,4 @@ -// Copyright (c) 2019 GeometryFactory SARL (France). +// Copyright (c) 2019 GeometryFactory Sarl (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -8,84 +8,285 @@ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // -// Author(s) : Simon Giraudot, Dmitry Anisimov +// Author(s) : Simon Giraudot, Dmitry Anisimov, Sven Oesau #ifndef CGAL_KINETIC_SHAPE_RECONSTRUCTION_3_H #define CGAL_KINETIC_SHAPE_RECONSTRUCTION_3_H -// #include - -// Boost includes. -#include -#include - -// CGAL includes. -#include -#include -#include -#include -#include - -#include -#include -#include - -// Internal includes. -#include +//#include +#include #include -#include -#include -#include -#include -#include -#include - -namespace CGAL { - -#ifdef DOXYGEN_NS +namespace CGAL +{ /*! * \ingroup PkgKineticPartition - \brief Creates the kinetic partitioning of the bounding box. + \brief Piece-wise linear reconstruction via inside/outside labeling of a kinetic partition using graph cut. - \tparam Kernel + \tparam GeomTraits must be a model of `Kernel`. Is used for non-critical calculations. - \tparam IntersectionKernel + \tparam IntersectionTraits must be a model of `Kernel`. Is used for the creation of the intersection graph. An exact kernel is suggested. -*/ -template -class Kinetic_partitioning_3 { -public: - /// \name Initialization - /// @{ - /*! - \brief Initializes the kinetic partitioning of the bounding box. \tparam InputRange - must be a model of `ConstRange` whose iterator type is `RandomAccessIterator`. + must be a model of `ConstRange` whose iterator type is `RandomAccessIterator`. +*/ +template +class Kinetic_shape_reconstruction_3 { +public: + using Input_range = typename Traits::Input_range; + using Kernel = typename Traits::Kernel; + using Intersection_Kernel = typename Traits::Intersection_Kernel; - \tparam PolygonMap - contains index ranges to form polygons from InputRange + using FT = typename Kernel::FT; - \tparam NamedParameters - a sequence of \ref bgl_namedparameters "Named Parameters" + using Point_2 = typename Traits::Point_2; + using Point_3 = typename Traits::Point_3; + using Plane_3 = typename Traits::Plane_3; - \param input_range - an instance of `InputRange` with 3D points and corresponding 3D normal vectors + using Point_map = typename Traits::Point_map; + using Normal_map = NormalMap; - \param polygon_map - a range of polygons defined by a range of indices into input_range + using Indices = std::vector; + using Polygon_3 = std::vector; + + using KSP = Kinetic_shape_partitioning_3; + + using Mesh = Surface_mesh; + + using Neighbor_query_3 = CGAL::Shape_detection::Point_set::K_neighbor_query; + using Planar_region = CGAL::Shape_detection::Point_set::Least_squares_plane_fit_region; + using Planar_sorting = CGAL::Shape_detection::Point_set::Least_squares_plane_fit_sorting; + using Region_growing_3 = CGAL::Shape_detection::Region_growing; + + Kinetic_shape_reconstruction_3(bool verbose, bool debug) : m_kinetic_partition(verbose, debug) { + } + /*! + \brief Detects shapes in the provided point cloud + + \tparam InputRange + must be a model of `ConstRange` whose iterator type is `RandomAccessIterator`. + + \tparam NamedParameters + a sequence of \ref bgl_namedparameters "Named Parameters" + + \param input_range + an instance of `InputRange` with 3D points and corresponding 3D normal vectors + + \param np + an instance of `NamedParameters`. + + \cgalNamedParamsBegin + \cgalParamNBegin{point_map} + \cgalParamDescription{a property map associating points to the elements of `input_range`} + \cgalParamDefault{`PointMap()`} + \cgalParamNEnd + \cgalParamNBegin{normal_map} + \cgalParamDescription{a property map associating normals to the elements of `input_range`} + \cgalParamDefault{`NormalMap()`} + \cgalParamNEnd + \cgalParamNBegin{k_neighbors} + \cgalParamDescription{the number of returned neighbors per each query point} + \cgalParamType{`std::size_t`} + \cgalParamDefault{12} + \cgalParamNEnd + \cgalParamNBegin{maximum_distance} + \cgalParamDescription{maximum distance from a point to a planar shape} + \cgalParamType{`GeomTraits::FT`} + \cgalParamDefault{1} + \cgalParamNEnd + \cgalParamNBegin{maximum_angle} + \cgalParamDescription{maximum angle in degrees between the normal of a point and the plane normal} + \cgalParamType{`GeomTraits::FT`} + \cgalParamDefault{25 degrees} + \cgalParamNEnd + \cgalParamNBegin{minimum_region_size} + \cgalParamDescription{minimum number of 3D points a region must have} + \cgalParamType{`std::size_t`} + \cgalParamDefault{5} + \cgalParamNEnd + \cgalNamedParamsEnd + + */ + template< + typename InputRange, + typename CGAL_NP_TEMPLATE_PARAMETERS> + std::size_t detect_planar_shapes( + InputRange input_range, + const CGAL_NP_CLASS& np = parameters::default_values()) { + + if (m_verbose) + std::cout << std::endl << "--- DETECTING PLANAR SHAPES: " << std::endl; + + m_planes.clear(); + m_polygons.clear(); + m_region_map.clear(); + + m_point_map = Point_set_processing_3_np_helper::get_point_map(input_range, np); + m_normal_map = Point_set_processing_3_np_helper::get_normal_map(input_range, np); + + create_planar_shapes(np); + + CGAL_assertion(m_planes.size() == m_polygons.size()); + CGAL_assertion(m_polygons.size() == m_region_map.size()); + +// if (m_debug) +// KSR_3::dump_polygons("detected-planar-shapes"); + + if (m_polygons.size() == 0) { + if (m_verbose) + std::cout << "* no planar shapes found" << std::endl; + return false; + } + return true; + } + + /*! + \brief Regularizes detected planar shapes by using `CGAL::Shape_regularization::Planes::regularize_planes` and merging coplanar planes afterwards. + + \tparam NamedParameters + a sequence of \ref bgl_namedparameters "Named Parameters" + + \param np + an instance of `NamedParameters`. + + \cgalNamedParamsBegin + \cgalParamNBegin{point_map} + \cgalParamDescription{a property map associating points to the elements of `input_range` that has been passed to detect_planar_shapes} + \cgalParamType{bool} + \cgalParamDefault{false} + \cgalParamNEnd + \cgalParamNBegin{normal_map} + \cgalParamDescription{a property map associating normals to the elements of `input_range` that has been passed to detect_planar_shapes} + \cgalParamType{bool} + \cgalParamDefault{false} + \cgalParamNEnd + \cgalParamNBegin{maximum_angle} + \cgalParamDescription{maximum allowed angle in degrees between plane normals used for parallelism, orthogonality, and axis symmetry} + \cgalParamType{FT} + \cgalParamDefault{25 degrees} + \cgalParamNEnd + \cgalParamNBegin{maximum_offset} + \cgalParamDescription{maximum allowed orthogonal distance between two parallel planes such that they are considered to be coplanar} + \cgalParamType{FT} + \cgalParamDefault{0.01} + \cgalParamNEnd + \cgalParamNBegin{regularize_parallelity} + \cgalParamDescription{indicates whether parallelism should be regularized or not} + \cgalParamType{bool} + \cgalParamDefault{true} + \cgalParamNEnd + \cgalParamNBegin{regularize_orthogonality} + \cgalParamDescription{indicates whether orthogonality should be regularized or not} + \cgalParamType{bool} + \cgalParamDefault{true} + \cgalParamNEnd + \cgalParamNBegin{regularize_coplanarity} + \cgalParamDescription{indicates whether coplanarity should be regularized or not} + \cgalParamType{bool} + \cgalParamDefault{true} + \cgalParamNEnd + \cgalNamedParamsEnd + + */ + template + std::size_t regularize_shapes( + const NamedParameters& np) { + + /*if (m_verbose) + std::cout << std::endl << "--- REGULARIZING PLANAR SHAPES: " << std::endl; + + const bool regularize = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::regularize), false); + if (!regularize) { + if (m_verbose) std::cout << "* user-defined, skipping" << std::endl; + return true; + } + + if (m_polygons.size() == 0) { + if (m_verbose) std::cout << "* no planes found, skipping" << std::endl; + return false; + } + + // Regularize. + + std::vector planes; + std::vector regions; + create_planes_and_regions(planes, regions); + + CGAL_assertion(planes.size() > 0); + CGAL_assertion(planes.size() == regions.size()); + + Plane_map plane_map; + Point_to_plane_map point_to_plane_map(m_input_range, regions); + CGAL::Shape_regularization::Planes::regularize_planes( + m_input_range, + planes, + plane_map, + point_to_plane_map, + true, true, true, false, + max_accepted_angle, + max_distance_to_plane, + symmetry_axis); + + const std::size_t num_polygons = m_polygons.size(); + + m_planes.clear(); + m_polygons.clear(); + m_region_map.clear(); + for (std::size_t i = 0; i < regions.size(); ++i) { + const auto& plane = planes[i]; + const auto& region = regions[i]; + + const std::size_t shape_idx = add_planar_shape(region, plane); + CGAL_assertion(shape_idx != std::size_t(-1)); + m_region_map[shape_idx] = region; + } + CGAL_assertion(m_polygons.size() == num_polygons); + CGAL_assertion(m_polygons.size() == m_planes.size()); + CGAL_assertion(m_polygons.size() == m_region_map.size()); + + if (m_verbose) + std::cout << "* num regularized planes: " << m_planes.size() << std::endl; + + + if (m_debug) + dump_polygons("regularized-planar-shapes");*/ + return true; + } + + /*! + \brief Retrieves the detected shapes. + + @return + vector with a plane equation for each detected planar shape. + + \pre `successful shape detection` + */ + const std::vector& detected_planar_shapes(std::vector >& indices) { + return m_planes; + } + + /*! + \brief Retrieves the indices of detected shapes. + + @return + vector with a plane equation for each detected planar shape. + + \pre `successful shape detection` + */ + const std::vector >& detected_planar_shape_indices() { + return m_planar_regions; + } + + /*! + \brief initializes the kinetic partitioning. \param np a sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below - @return - success - - \pre `input_range.size() > 0 and polygon_map.size() > 0` - \cgalNamedParamsBegin \cgalParamNBegin{reorient_bbox} \cgalParamDescription{Use the oriented bounding box instead of the axis-aligned bounding box.} @@ -108,846 +309,262 @@ public: \cgalParamDefault{5} \cgalParamNEnd \cgalNamedParamsEnd + + \pre `successful shape detection` */ - template< - typename InputRange, - typename PolygonMap, - typename NamedParameters> - bool initialization( - const InputRange input_range, - const PolygonMap polygon_map, - const NamedParameters& np); + template + bool initialize_partitioning(const CGAL_NP_CLASS& np = parameters::default_values()) { + if (m_polygons.size() == 0) { + std::cout << "No planar shapes available to create kinetic partitioning." << std::endl; + return false; + } + + using Polygon_3 = std::vector; + using Polygon_map = CGAL::Identity_property_map; + + return m_kinetic_partition.initialize(m_polygons, Polygon_map()); + } /*! \brief Propagates the kinetic polygons in the initialized partition. - \param k - maximum number of allowed intersections for each input polygon before its expansion stops. + \param k + maximum number of allowed intersections for each input polygon before its expansion stops. - @return - success of kinetic partitioning. + @return + success of kinetic partitioning. - \pre `successful initialization` + \pre `successful initialization` */ - /// @} - - bool partition(std::size_t k); - - /// \name Access - /// @{ + bool partition(std::size_t k) { + return m_kinetic_partition.partition(k); + } /*! - \brief Number of vertices in the kinetic partitioning. + \brief Access to the kinetic partitioning. - @return - number of vertices. + @return + created kinetic partitioning data structure - \pre `successful partitioning` + \pre `successful partitioning` */ - std::size_t number_of_vertices() const; + const Kinetic_shape_partitioning_3& partitioning() const { + return m_kinetic_partition; + } /*! - \brief Number of convex faces in the kinetic partitioning. + \brief Creates the visibility (data-) and regularity energy terms from the input point cloud and the kinetic partitioning. - @return - number of convex faces. + @return + success. - \pre `successful partitioning` + \pre `successful initialization` */ - std::size_t number_of_faces() const; + bool setup_energyterms() { + if (m_kinetic_partition.number_of_volumes() == 0) { + if (m_verbose) std::cout << "Kinetic partition is not constructed or does not have volumes" << std::endl; + return false; + } + +/* + if (m_verbose) std::cout << "* computing visibility ... "; + std::map face2points; + assign_points_to_pfaces(face2points); + const Visibility visibility( + m_data, face2points, m_point_map_3, m_normal_map_3); + + CGAL_assertion(m_data.volumes().size() > 0); + visibility.compute(m_data.volumes()); + //dump_visibility("visibility/visibility", pface_points); + + if (m_verbose) { + std::cout << "done" << std::endl; + std::cout << "* applying graphcut ... "; + } + + const FT beta = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::graphcut_beta), FT(1) / FT(2));*/ + + return false; + } /*! - \brief Number of convex volumes created by the kinetic partitioning. + \brief Provides the data and regularity energy terms for reconstruction via graph-cut. - @return - number of convex volumes. + \param edges + contains a vector of pairs of volume indices. Indicates which volumes should be connected in the graph cut formulation. - \pre `successful partitioning` + \param edge_costs + contains the cost for each edge specified in `edges` for two labels with different labels. For equal labels, the cost is 0. Needs to be index compatible to the `edges` parameter. + + \param cost_matrix + provides the cost of a label for each volume cell. The first index corresponds to the label and the second index corresponds to the volume index. + + @return + fails if the dimensions of parameters does not match the kinetic partitioning. + + \pre `successful initialization` */ - std::size_t number_of_volumes() const; + template + bool setup_energyterms( + const std::vector< std::pair >& edges, + const std::vector& edge_costs, + const std::vector< std::vector >& cost_matrix); /*! - \brief Point vector for mapping vertex indices to positions. + \brief Propagates the kinetic polygons in the initialized partition. - @return - vector of points. + \param lambda + trades the impact of the data term for impact of the regularization term. Should be in the range [0, 1). - \pre `successful partitioning` + @return + success of reconstruction. + + \pre `successful initialization` */ - const std::vector& vertices() const; + bool reconstruct(FT lambda) { + return false; + } /*! - \brief Vertex indices of convex face. + \brief Provides the reconstructed surface mesh - \param face_index - index of the query face. + \param mesh + a mesh object to store the reconstructed surface. - @return - vector of vertex indices. - - \pre `successful partitioning` + \pre `successful reconstruction` */ - const std::vector& vertices(std::size_t face_index) const; - - /*! - \brief Face indices of the convex volume. - - \param volume_index - index of the query volume. - - @return - vector of face indices. - - \pre `successful partitioning` - */ - const std::vector& face(std::size_t volume_index) const; - - /*! - \brief Indices of adjacent volumes. Negative indices correspond to the empty spaces behind the sides of the bounding box. - - \param face_index - index of the query face. - - @return - pair of adjacent volumes. - - -1 zmin - -2 ymin - -3 xmax - -4 ymax - -5 xmin - -6 zmax - - \pre `successful partitioning` - */ - const std::pair& neighbors(std::size_t face_index) const; - - /*! - \brief Retrieves the input polygon this face originates from. - - \param face_index - index of the query face. - - @return - index into polygon_map provided on initialization. - - \pre `successful partitioning` - */ - const std::size_t input_polygon(std::size_t face_index) const; - - /*! - \brief Creates a linear cell complex from the kinetic partitioning. - - \tparam LCC - Linear_cell_complex_for_combinatorial_map<3, 3,...> - The dimension of the combinatorial map and the dimension of the ambient space have to be 3. - - \param lcc - an instance of LCC - - \pre `successful partitioning` - */ - template - void get_linear_cell_complex(LCC& lcc) const; - - /// @} -} -#else - -template -class Kinetic_shape_reconstruction_3 { - -public: - using Kernel = GeomTraits; - using LCC = Linear_cell_complex_for_combinatorial_map<3, 3>; + void output_reconstructed_model(Mesh& mesh); private: - using FT = typename Kernel::FT; - using Point_3 = typename Kernel::Point_3; + bool m_verbose; + bool m_debug; - using Data_structure = KSR_3::Data_structure; + Input_range m_points; + Point_map m_point_map; + Normal_map m_normal_map; - using IVertex = typename Data_structure::IVertex; - using IEdge = typename Data_structure::IEdge; + std::vector > m_planar_regions; + std::map m_region_map; - using From_exact = CGAL::Cartesian_converter; + std::vector m_planes; + std::vector m_polygons; + KSP m_kinetic_partition; - using Initializer = KSR_3::Initializer; - using Propagation = KSR_3::FacePropagation; - using Finalizer = KSR_3::Finalizer; + std::size_t add_convex_hull_shape( + const std::vector& region, const Plane_3& plane) { - using Polygon_mesh = CGAL::Surface_mesh; - using Vertex_index = typename Polygon_mesh::Vertex_index; - using Timer = CGAL::Real_timer; - using Parameters = KSR::Parameters_3; + std::vector points; + points.reserve(region.size()); + for (const std::size_t idx : region) { + CGAL_assertion(idx < m_points.size()); + const auto& p = get(m_point_map, idx); + const auto q = plane.projection(p); + const auto point = plane.to_2d(q); + points.push_back(point); + } + CGAL_assertion(points.size() == region.size()); -private: - Parameters m_parameters; - Data_structure m_data; - std::size_t m_num_events; + std::vector ch; + CGAL::convex_hull_2(points.begin(), points.end(), std::back_inserter(ch)); -public: - Kinetic_shape_reconstruction_3( - const bool verbose = true, - const bool debug = false) : - m_parameters(verbose, debug, false), // use true here to export all steps - m_data(m_parameters), - m_num_events(0) - { } - - template< - typename InputRange, - typename PolygonMap, - typename NamedParameters> - bool partition( - const InputRange & input_range, - const PolygonMap polygon_map, - const NamedParameters & np) { - - Timer timer; - m_parameters.k = parameters::choose_parameter( - parameters::get_parameter(np, internal_np::k_intersections), 1); - m_parameters.n = parameters::choose_parameter( - parameters::get_parameter(np, internal_np::n_subdivisions), 0); - m_parameters.enlarge_bbox_ratio = parameters::choose_parameter( - parameters::get_parameter(np, internal_np::enlarge_bbox_ratio), FT(11) / FT(10)); - m_parameters.distance_tolerance = parameters::choose_parameter( - parameters::get_parameter(np, internal_np::distance_tolerance), FT(5) / FT(10)); - m_parameters.reorient = parameters::choose_parameter( - parameters::get_parameter(np, internal_np::reorient), false); - - std::cout.precision(20); - if (input_range.size() == 0) { - CGAL_warning_msg(input_range.size() > 0, - "WARNING: YOUR INPUT IS EMPTY! RETURN WITH NO CHANGE!"); - return false; + std::vector polygon; + for (const auto& p : ch) { + const auto point = plane.to_3d(p); + polygon.push_back(point); } - if (m_parameters.n != 0) { - CGAL_assertion_msg(false, "TODO: IMPLEMENT KINETIC SUBDIVISION!"); - if (m_parameters.n > 3) { - CGAL_warning_msg(m_parameters.n <= 3, - "WARNING: DOES IT MAKE SENSE TO HAVE MORE THAN 64 INPUT BLOCKS? SETTING N TO 3!"); - m_parameters.n = 3; + const std::size_t shape_idx = m_polygons.size(); + m_polygons.push_back(polygon); + m_planes.push_back(plane); + return shape_idx; + } + + template + void create_planar_shapes(const NamedParameters& np) { + + if (m_points.size() < 3) { + if (m_verbose) std::cout << "* no points found, skipping" << std::endl; + return; + } + if (m_verbose) std::cout << "* getting planar shapes using region growing" << std::endl; + + // Parameters. + const std::size_t k = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::k_neighbors), 12); + const FT max_distance_to_plane = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::distance_threshold), FT(1)); + const FT max_accepted_angle = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::angle_threshold), FT(15)); + const std::size_t min_region_size = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::min_region_size), 50); + + // Region growing. + Neighbor_query_3 neighbor_query(m_points, k, m_point_map); + + Planar_region planar_region(m_points, + max_distance_to_plane, max_accepted_angle, min_region_size, + m_point_map, m_normal_map); + + Planar_sorting sorting( + m_points, neighbor_query, m_point_map); + sorting.sort(); + + std::vector result; + Region_growing_3 region_growing( + m_points, neighbor_query, planar_region, sorting.seed_map()); + region_growing.detect(std::back_inserter(result)); + + // Convert indices. + m_planar_regions.clear(); + m_planar_regions.reserve(result.size()); + + Indices region; + for (const auto& indices : result) { + region.clear(); + for (const std::size_t index : indices) { + region.push_back(index); } + m_planar_regions.push_back(region); + const auto plane = fit_plane(region); + const std::size_t shape_idx = add_convex_hull_shape(region, plane); + CGAL_assertion(shape_idx != std::size_t(-1)); + m_region_map[shape_idx] = region; } + CGAL_assertion(m_planar_regions.size() == result.size()); - if (m_parameters.enlarge_bbox_ratio < FT(1)) { - CGAL_warning_msg(m_parameters.enlarge_bbox_ratio >= FT(1), - "WARNING: YOU SET ENLARGE_BBOX_RATIO < 1.0! THE VALID RANGE IS [1.0, +INF). SETTING TO 1.0!"); - m_parameters.enlarge_bbox_ratio = FT(1); - } - - if (m_parameters.verbose) { - const unsigned int num_blocks = static_cast(std::pow(m_parameters.n + 1, 3)); - const std::string is_reorient = (m_parameters.reorient ? "true" : "false"); - - std::cout << std::endl << "--- PARTITION OPTIONS: " << std::endl; - std::cout << "* number of intersections k: " << m_parameters.k << std::endl; - std::cout << "* number of subdivisions per bbox side: " << m_parameters.n << std::endl; - std::cout << "* number of subdivision blocks: " << num_blocks << std::endl; - std::cout << "* enlarge bbox ratio: " << m_parameters.enlarge_bbox_ratio << std::endl; - std::cout << "* reorient: " << is_reorient << std::endl; - } - - if (m_parameters.verbose) { - std::cout << std::endl << "--- INITIALIZING PARTITION:" << std::endl; - } - - // Initialization. - timer.reset(); - timer.start(); - m_data.clear(); - - Initializer initializer(m_data, m_parameters); - initializer.initialize(input_range, polygon_map); - - timer.stop(); - const double time_to_initialize = timer.time(); - - // if (m_parameters.verbose) { - // std::cout << std::endl << "* initialization (sec.): " << time_to_initialize << std::endl; - // std::cout << "INITIALIZATION SUCCESS!" << std::endl << std::endl; - // } - // exit(EXIT_SUCCESS); - - // Output planes. - // for (std::size_t i = 6; i < m_data.number_of_support_planes(); ++i) { - // std::cout << m_data.support_plane(i).plane() << std::endl; - // } - - if (m_parameters.k == 0) { // for k = 0, we skip propagation - CGAL_warning_msg(m_parameters.k > 0, - "WARNING: YOU SET K TO 0! THAT MEANS NO PROPAGATION! THE VALID VALUES ARE {1,2,...}. INTERSECT AND RETURN!"); - return false; - } - - if (m_parameters.verbose) { - std::cout << std::endl << "--- RUNNING THE QUEUE:" << std::endl; - std::cout << "* propagation started" << std::endl; - } - - // Propagation. - timer.reset(); - timer.start(); - std::size_t num_queue_calls = 0; - - Propagation propagation(m_data, m_parameters); - std::tie(num_queue_calls, m_num_events) = propagation.propagate(); - - timer.stop(); - const double time_to_propagate = timer.time(); - - if (m_parameters.verbose) { - std::cout << "* propagation finished" << std::endl; - std::cout << "* number of queue calls: " << num_queue_calls << std::endl; - std::cout << "* number of events handled: " << m_num_events << std::endl; - } - - if (m_parameters.verbose) { - std::cout << std::endl << "--- FINALIZING PARTITION:" << std::endl; - } - - // Finalization. - timer.reset(); - timer.start(); - if (m_parameters.debug) - dump(m_data, "final-" + m_parameters.k); - - Finalizer finalizer(m_data, m_parameters); - - if (m_parameters.verbose) - std::cout << "* checking final mesh integrity ..."; - - CGAL_assertion(m_data.check_integrity(true, true, true)); - - if (m_parameters.verbose) - std::cout << " done" << std::endl; - - if (m_parameters.verbose) - std::cout << "* getting volumes ..." << std::endl; - - finalizer.create_polyhedra(); - timer.stop(); - const double time_to_finalize = timer.time(); - - if (m_parameters.verbose) { - std::cout << "* found all together " << m_data.number_of_volumes() << " volumes" << std::endl; - - for (std::size_t i = 0; i < m_data.number_of_support_planes(); i++) { - dump_2d_surface_mesh(m_data, i, "final-surface-mesh-" + std::to_string(i)); - } - } - - // Timing. - if (m_parameters.verbose) { - std::cout << std::endl << "--- TIMING (sec.):" << std::endl; - } - const double total_time = - time_to_initialize + time_to_propagate + time_to_finalize; - - std::cout << "* initialization: " << time_to_initialize << std::endl; - std::cout << "* propagation: " << time_to_propagate << std::endl; - std::cout << "* finalization: " << time_to_finalize << std::endl; - std::cout << "* total time: " << total_time << std::endl; - - return true; + if (m_verbose) + std::cout << "* found " << m_polygons.size() << " planar shapes" << std::endl; } - template< - typename InputRange, - typename PointMap, - typename VectorMap, - typename SemanticMap, - typename NamedParameters> - bool reconstruct( - const InputRange& input_range, - const PointMap point_map, - const VectorMap normal_map, - const SemanticMap semantic_map, - const std::string file_name, - const NamedParameters& np) { + const Plane_3 fit_plane(const std::vector& region) const { - using Reconstruction = KSR_3::Reconstruction< - InputRange, PointMap, VectorMap, SemanticMap, Kernel, Intersection_Kernel>; - - Reconstruction reconstruction( - input_range, point_map, normal_map, semantic_map, m_data, m_parameters.verbose, m_parameters.debug); - - bool success = reconstruction.detect_planar_shapes(file_name, np); - if (!success) { - CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, DETECTING PLANAR SHAPES FAILED!"); - return false; + std::vector points; + points.reserve(region.size()); + for (const std::size_t idx : region) { + CGAL_assertion(idx < m_points.size()); + points.push_back(get(m_point_map, idx)); } - // exit(EXIT_SUCCESS); - - success = reconstruction.regularize_planar_shapes(np); - if (!success) { - CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, REGULARIZATION FAILED!"); - return false; - } - // exit(EXIT_SUCCESS); - - success = partition( - reconstruction.planar_shapes(), reconstruction.polygon_map(), np); - if (!success) { - CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, PARTITION FAILED!"); - return false; - } - // exit(EXIT_SUCCESS); - - success = reconstruction.compute_model(np); - if (!success) { - CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, COMPUTING MODEL FAILED!"); - return false; - } - return success; - } - - template< - typename InputRange, - typename PointMap, - typename VectorMap, - typename SemanticMap, - typename RegionMap, - typename NamedParameters> - bool reconstruct( - const InputRange& input_range, - const PointMap point_map, - const VectorMap normal_map, - const SemanticMap semantic_map, - const RegionMap region_map, - const NamedParameters& np) { - - using Reconstruction = KSR_3::Reconstruction< - InputRange, PointMap, VectorMap, SemanticMap, Kernel, Intersection_Kernel>; - - Reconstruction reconstruction( - input_range, point_map, normal_map, semantic_map, m_data, m_parameters.verbose, m_parameters.debug); - - bool success = reconstruction.planar_shapes_from_map(region_map, np); - if (!success) { - CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, DETECTING PLANAR SHAPES FAILED!"); - return false; - } - // exit(EXIT_SUCCESS); - - //success = reconstruction.regularize_planar_shapes(np); - if (!success) { - CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, REGULARIZATION FAILED!"); - return false; - } - // exit(EXIT_SUCCESS); - - success = partition( - reconstruction.planar_shapes(), reconstruction.polygon_map(), np); - if (!success) { - CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, PARTITION FAILED!"); - return false; - } - // exit(EXIT_SUCCESS); - - success = reconstruction.compute_model(np); - if (!success) { - CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, COMPUTING MODEL FAILED!"); - return false; - } - return success; - } - - /******************************* - ** STATISTICS ** - ********************************/ - - std::size_t number_of_events() const { - return m_num_events; - } - - int number_of_support_planes() const { - return static_cast(m_data.number_of_support_planes()); - } - - std::size_t number_of_vertices(const int support_plane_idx = -1) const { - - CGAL_assertion(support_plane_idx < number_of_support_planes()); - if (support_plane_idx >= number_of_support_planes()) return std::size_t(-1); - if (support_plane_idx < 0) { - return m_data.igraph().number_of_vertices(); - } - - CGAL_assertion(support_plane_idx >= 0); - const std::size_t sp_idx = static_cast(support_plane_idx); - return static_cast(m_data.mesh(sp_idx).number_of_vertices()); - } - - std::size_t number_of_edges(const int support_plane_idx = -1) const { - - CGAL_assertion(support_plane_idx < number_of_support_planes()); - if (support_plane_idx >= number_of_support_planes()) return std::size_t(-1); - if (support_plane_idx < 0) { - return m_data.igraph().number_of_edges(); - } - - CGAL_assertion(support_plane_idx >= 0); - const std::size_t sp_idx = static_cast(support_plane_idx); - return static_cast(m_data.mesh(sp_idx).number_of_edges()); - } - - std::size_t number_of_faces(const int support_plane_idx = -1) const { - - CGAL_assertion(support_plane_idx < number_of_support_planes()); - if (support_plane_idx >= number_of_support_planes()) return std::size_t(-1); - if (support_plane_idx < 0) { - std::size_t num_all_faces = 0; - for (int i = 0; i < number_of_support_planes(); ++i) { - const std::size_t num_faces = static_cast( - m_data.mesh(static_cast(i)).number_of_faces()); - num_all_faces += num_faces; - } - return num_all_faces; - } - - CGAL_assertion(support_plane_idx >= 0); - const std::size_t sp_idx = static_cast(support_plane_idx); - const std::size_t num_faces = static_cast( - m_data.mesh(sp_idx).number_of_faces()); - return num_faces; - } - - std::size_t number_of_volumes() const { - return m_data.volumes().size(); - } - - int support_plane_index(const std::size_t polygon_index) const { - const int support_plane_idx = m_data.support_plane_index(polygon_index); - CGAL_assertion(support_plane_idx >= 6); - return support_plane_idx; - } - - template - void get_linear_cell_complex(LCC &lcc) const { - using LCC_Kernel = typename LCC::Traits; - CGAL::Cartesian_converter conv; - lcc.clear(); - - std::vector used_vertices(m_data.igraph().number_of_vertices(), false); - std::vector remap(m_data.igraph().number_of_vertices(), -1); - std::vector mapped_vertices; - mapped_vertices.reserve(m_data.igraph().number_of_vertices()); - - for (const auto& volume : m_data.volumes()) { - for (const auto& vertex : volume.pvertices) { - CGAL_assertion(m_data.has_ivertex(vertex)); - IVertex ivertex = m_data.ivertex(vertex); - if (remap[ivertex] == -1) { - remap[ivertex] = static_cast(mapped_vertices.size()); - mapped_vertices.push_back(conv(m_data.point_3(ivertex))); - } - } - } - - CGAL::Linear_cell_complex_incremental_builder_3 ib(lcc); - for (const auto& p : mapped_vertices) - ib.add_vertex(p); - - for (const auto& vol : m_data.volumes()) { - ib.begin_surface(); - for (std::size_t i = 0; i < vol.pfaces.size(); i++) { - auto vertex_range = m_data.pvertices_of_pface(vol.pfaces[i]); - ib.begin_facet(); - if (vol.pface_oriented_outwards[i]) { - typename Data_structure::PVertex_of_pface_iterator it = vertex_range.begin(); - while (it != vertex_range.end()) { - CGAL_assertion(m_data.has_ivertex(*it)); - IVertex ivertex = m_data.ivertex(*it); - ib.add_vertex_to_facet(static_cast(remap[ivertex])); - it++; - } - } - else { - typename Data_structure::PVertex_of_pface_iterator it = vertex_range.end()--; - do { - CGAL_assertion(m_data.has_ivertex(*it)); - IVertex ivertex = m_data.ivertex(*it); - ib.add_vertex_to_facet(static_cast(remap[ivertex])); - it--; - if (it == vertex_range.begin()) - break; - } while (true); - } - ib.end_facet(); - } - ib.end_surface(); - } - - lcc.display_characteristics(std::cout) << std::endl; - } - - /******************************* - ** OUTPUT ** - ********************************/ - - template - VertexOutputIterator output_partition_vertices( - VertexOutputIterator vertices, const int support_plane_idx = -1) const { - From_exact from_EK; - - CGAL_assertion(support_plane_idx < number_of_support_planes()); - if (support_plane_idx >= number_of_support_planes()) return vertices; - if (support_plane_idx < 0) { - const auto all_ivertices = m_data.ivertices(); - for (const auto ivertex : all_ivertices) { - *(vertices++) = from_EK(m_data.point_3(ivertex)); - } - return vertices; - } - - CGAL_assertion(support_plane_idx >= 0); - const std::size_t sp_idx = static_cast(support_plane_idx); - const auto all_pvertices = m_data.pvertices(sp_idx); - for (const auto pvertex : all_pvertices) { - CGAL_assertion(m_data.has_ivertex(pvertex)); - const auto ivertex = m_data.ivertex(pvertex); - *(vertices++) = from_EK(m_data.point_3(ivertex)); - } - return vertices; - } - - template - EdgeOutputIterator output_partition_edges( - EdgeOutputIterator edges, const int support_plane_idx = -1) const { - From_exact from_EK; - - CGAL_assertion(support_plane_idx < number_of_support_planes()); - if (support_plane_idx >= number_of_support_planes()) return edges; - if (support_plane_idx < 0) { - const auto all_iedges = m_data.iedges(); - for (const auto iedge : all_iedges) { - *(edges++) = from_EK(m_data.segment_3(iedge)); - } - return edges; - } - - CGAL_assertion(support_plane_idx >= 0); - const std::size_t sp_idx = static_cast(support_plane_idx); - const auto all_pedges = m_data.pedges(sp_idx); - for (const auto pedge : all_pedges) { - CGAL_assertion(m_data.has_iedge(pedge)); - const auto iedge = m_data.iedge(pedge); - *(edges++) = from_EK(m_data.segment_3(iedge)); - } - return edges; - } - - template - FaceOutputIterator output_partition_faces( - FaceOutputIterator faces, - const int support_plane_idx = -1, - const int begin = 0) const { - - KSR::Indexer indexer; - CGAL_assertion(support_plane_idx < number_of_support_planes()); - if (support_plane_idx >= number_of_support_planes()) return faces; - if (support_plane_idx < 0) { - const auto all_ivertices = m_data.ivertices(); - for (const auto ivertex : all_ivertices) indexer(ivertex); - for (int i = begin; i < number_of_support_planes(); ++i) { - const std::size_t sp_idx = static_cast(i); - output_partition_faces(faces, indexer, sp_idx); - } - return faces; - } - - CGAL_assertion(support_plane_idx >= 0); - const std::size_t sp_idx = static_cast(support_plane_idx); - const auto all_pvertices = m_data.pvertices(sp_idx); - for (const auto pvertex : all_pvertices) { - CGAL_assertion(m_data.has_ivertex(pvertex)); - const auto ivertex = m_data.ivertex(pvertex); - indexer(ivertex); - } - return output_partition_faces(faces, indexer, sp_idx); - } - - void output_support_plane( - Polygon_mesh& polygon_mesh, const int support_plane_idx) const { - From_exact from_EK; - - polygon_mesh.clear(); - CGAL_assertion(support_plane_idx >= 0); - if (support_plane_idx < 0) return; - CGAL_assertion(support_plane_idx < number_of_support_planes()); - if (support_plane_idx >= number_of_support_planes()) return; - const std::size_t sp_idx = static_cast(support_plane_idx); - - std::vector vertices; - std::vector map_vertices; - - map_vertices.clear(); - const auto all_pvertices = m_data.pvertices(sp_idx); - for (const auto pvertex : all_pvertices) { - CGAL_assertion(m_data.has_ivertex(pvertex)); - const auto ivertex = m_data.ivertex(pvertex); - - if (map_vertices.size() <= pvertex.second) - map_vertices.resize(pvertex.second + 1); - map_vertices[pvertex.second] = - polygon_mesh.add_vertex(from_EK(m_data.point_3(ivertex))); - } - - const auto all_pfaces = m_data.pfaces(sp_idx); - for (const auto pface : all_pfaces) { - vertices.clear(); - const auto pvertices = m_data.pvertices_of_pface(pface); - for (const auto pvertex : pvertices) { - vertices.push_back(map_vertices[pvertex.second]); - } - polygon_mesh.add_face(vertices); - } - } - - template - VolumeOutputIterator output_partition_volumes( - VolumeOutputIterator volumes) const { - for (std::size_t i = 0; i < m_data.number_of_volumes(); ++i) { - output_partition_volume(volumes, i); - } - return volumes; - } - - template - VolumeOutputIterator output_partition_volume( - VolumeOutputIterator volumes, const std::size_t volume_index) const { - - CGAL_assertion(volume_index < number_of_volumes()); - if (volume_index >= number_of_volumes()) return volumes; - - std::vector vertices; - std::vector< std::vector > faces; - output_partition_volume( - std::back_inserter(vertices), std::back_inserter(faces), volume_index); - CGAL::Polygon_mesh_processing::orient_polygon_soup(vertices, faces); - - Polygon_mesh polygon_mesh; - CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh( - vertices, faces, polygon_mesh); - *(volumes++) = polygon_mesh; - return volumes; - } - - template - void output_partition_volume( - VertexOutputIterator vertices, FaceOutputIterator faces, - const std::size_t volume_index) const { - - CGAL_assertion(volume_index < number_of_volumes()); - if (volume_index >= number_of_volumes()) return; - - const auto& volume = m_data.volumes()[volume_index]; - - std::size_t num_vertices = 0; - KSR::Indexer indexer; - - std::vector face; - const auto& pfaces = volume.pfaces; - for (const auto& pface : pfaces) { - face.clear(); - const auto pvertices = m_data.pvertices_of_pface(pface); - for (const auto pvertex : pvertices) { - - CGAL_assertion(m_data.has_ivertex(pvertex)); - const auto ivertex = m_data.ivertex(pvertex); - const std::size_t idx = indexer(ivertex); - - if (idx == num_vertices) { - *(vertices++) = m_data.point_3(ivertex); - ++num_vertices; - } - face.push_back(idx); - } - *(faces++) = face; - } - } - - template - void output_reconstructed_model( - VertexOutputIterator vertices, FaceOutputIterator faces) const { - From_exact from_EK; - - const auto& model = m_data.reconstructed_model(); - CGAL_assertion(model.pfaces.size() > 0); - - std::size_t num_vertices = 0; - KSR::Indexer indexer; - - std::vector face; - const auto& pfaces = model.pfaces; - for (const auto& pface : pfaces) { - face.clear(); - const auto pvertices = m_data.pvertices_of_pface(pface); - for (const auto pvertex : pvertices) { - - CGAL_assertion(m_data.has_ivertex(pvertex)); - const auto ivertex = m_data.ivertex(pvertex); - const std::size_t idx = indexer(ivertex); - - if (idx == num_vertices) { - *(vertices++) = from_EK(m_data.point_3(ivertex)); - ++num_vertices; - } - face.push_back(idx); - } - *(faces++) = face; - } - } - - void output_reconstructed_model(Polygon_mesh& polygon_mesh) const { - - std::vector vertices; - std::vector< std::vector > faces; - output_reconstructed_model( - std::back_inserter(vertices), std::back_inserter(faces)); - CGAL::Polygon_mesh_processing::orient_polygon_soup(vertices, faces); - polygon_mesh.clear(); - CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh( - vertices, faces, polygon_mesh); - } - - /******************************* - ** MEMORY ** - ********************************/ - - void clear() { - m_data.clear(); - m_num_events = 0; - } - -private: - - template - FaceOutputIterator output_partition_faces( - FaceOutputIterator faces, KSR::Indexer& indexer, - const std::size_t sp_idx) const { - - std::vector face; - const auto all_pfaces = m_data.pfaces(sp_idx); - for (const auto pface : all_pfaces) { - face.clear(); - const auto pvertices = m_data.pvertices_of_pface(pface); - for (const auto pvertex : pvertices) { - CGAL_assertion(m_data.has_ivertex(pvertex)); - const auto ivertex = m_data.ivertex(pvertex); - const std::size_t idx = indexer(ivertex); - face.push_back(idx); - } - *(faces++) = face; - } - return faces; + CGAL_assertion(points.size() == region.size()); + + Plane_3 fitted_plane; + Point_3 fitted_centroid; + CGAL::linear_least_squares_fitting_3( + points.begin(), points.end(), + fitted_plane, fitted_centroid, + CGAL::Dimension_tag<0>()); + + const Plane_3 plane( + static_cast(fitted_plane.a()), + static_cast(fitted_plane.b()), + static_cast(fitted_plane.c()), + static_cast(fitted_plane.d())); + return plane; } }; -#endif //DOXYGEN_RUNNING } // namespace CGAL + #endif // CGAL_KINETIC_SHAPE_RECONSTRUCTION_3_H diff --git a/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_2d_stress_test.cpp b/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_2d_stress_test.cpp index acaafd05bd6..49342cbe410 100644 --- a/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_2d_stress_test.cpp +++ b/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_2d_stress_test.cpp @@ -1,14 +1,7 @@ #include -using SC = CGAL::Simple_cartesian; - #include -using EPECK = CGAL::Exact_predicates_exact_constructions_kernel; - #include -using EPICK = CGAL::Exact_predicates_inexact_constructions_kernel; - #include -using EPECK_to_EPICK = CGAL::Cartesian_converter; #include #include @@ -16,6 +9,11 @@ using EPECK_to_EPICK = CGAL::Cartesian_converter; #include #include +using SC = CGAL::Simple_cartesian; +using EPECK = CGAL::Exact_predicates_exact_constructions_kernel; +using EPICK = CGAL::Exact_predicates_inexact_constructions_kernel; +using EPECK_to_EPICK = CGAL::Cartesian_converter; + using FT = typename EPECK::FT; using Point_2 = typename EPECK::Point_2; using Vector_2 = typename EPECK::Vector_2; diff --git a/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_3d_test_all.cpp b/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_3d_test_all.cpp index 1ecbaa233e2..9aa0cb06e2d 100644 --- a/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_3d_test_all.cpp +++ b/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_3d_test_all.cpp @@ -1,7 +1,7 @@ #include #include #include -#include +#include #include #include #include @@ -12,6 +12,8 @@ using EPICK = CGAL::Exact_predicates_inexact_constructions_kernel; using EPECK = CGAL::Exact_predicates_exact_constructions_kernel; using Timer = CGAL::Real_timer; +using Traits = typename CGAL::Kinetic_shape_partitioning_Traits_3, CGAL::Identity_property_map >; + template struct Polygon_map { @@ -48,11 +50,11 @@ bool run_test( std::vector< std::vector >& all_times, std::size_t& num_tests) { - using Point_3 = typename Traits::Point_3; - using Segment_3 = typename Traits::Segment_3; + using Point_3 = typename Traits::Kernel::Point_3; + using Segment_3 = typename Traits::Kernel::Segment_3; using Surface_mesh = CGAL::Surface_mesh; - using KSR = CGAL::Kinetic_shape_reconstruction_3; + using KSP = CGAL::Kinetic_shape_partitioning_3; ++num_tests; std::string baseDir = "C:/dev/kinetic/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/"; @@ -83,34 +85,35 @@ bool run_test( double time = 0.0; for (std::size_t iter = 0; iter < num_iters; ++iter) { std::cout << std::endl << "--ITERATION #" << iter + 1 << " BEGIN!" << std::endl; - KSR ksr(true, true); // first verbose, second debug + KSP ksp(true, false); // first verbose, second debug // Running KSR. Timer timer; timer.start(); - const bool is_ksr_success = ksr.partition( - input_faces, polygon_map, CGAL::parameters::k_intersections(k)); - assert(is_ksr_success); - if (!is_ksr_success) return false; + + bool is_ksp_success = ksp.initialize( + input_faces, polygon_map); + + if (is_ksp_success) + ksp.partition(k); + + assert(is_ksp_success); + if (!is_ksp_success) return false; timer.stop(); time += timer.time(); // Testing results. - const std::size_t num_events = ksr.number_of_events(); - assert(num_events > 0); - const int num_support_planes = ksr.number_of_support_planes(); + const int num_support_planes = ksp.number_of_support_planes(); - const int num_vertices = static_cast(ksr.number_of_vertices()); - const int num_edges = static_cast(ksr.number_of_edges()); - const int num_faces = static_cast(ksr.number_of_faces()); - const int num_volumes = static_cast(ksr.number_of_volumes()); + const int num_vertices = static_cast(ksp.number_of_vertices()); + const int num_faces = static_cast(ksp.number_of_faces()); + const int num_volumes = static_cast(ksp.number_of_volumes()); std::cout << std::endl << "--RESULTS: "; std::cout << num_support_planes << ","; std::cout << num_vertices << ","; - std::cout << num_edges << ","; std::cout << num_faces << ","; std::cout << num_volumes << std::endl; @@ -136,14 +139,14 @@ bool run_test( if (num_volumes < results[5]) return false;*/ CGAL::Linear_cell_complex_for_combinatorial_map<3, 3> lcc; - ksr.get_linear_cell_complex(lcc); - + ksp.get_linear_cell_complex(lcc); +/* std::vector output_vertices; - ksr.output_partition_vertices( + ksp.output_partition_vertices( std::back_inserter(output_vertices)); assert(static_cast(num_vertices) == output_vertices.size()); if (static_cast(num_vertices) != output_vertices.size()) return false; - /* + std::vector output_edges; ksr.output_partition_edges( @@ -163,18 +166,16 @@ bool run_test( assert(static_cast(num_volumes) == output_volumes.size()); if (static_cast(num_volumes) != output_volumes.size()) return false;*/ - ksr.clear(); - assert(ksr.number_of_support_planes() == 0); - assert(ksr.number_of_vertices() == 0); - assert(ksr.number_of_edges() == 0); - assert(ksr.number_of_faces() == 0); - assert(ksr.number_of_volumes() == 0); + ksp.clear(); + assert(ksp.number_of_support_planes() == 0); + assert(ksp.number_of_vertices() == 0); + assert(ksp.number_of_faces() == 0); + assert(ksp.number_of_volumes() == 0); - if (ksr.number_of_support_planes() != 0) return false; - if (ksr.number_of_vertices() != 0) return false; - if (ksr.number_of_edges() != 0) return false; - if (ksr.number_of_faces() != 0) return false; - if (ksr.number_of_volumes() != 0) return false; + if (ksp.number_of_support_planes() != 0) return false; + if (ksp.number_of_vertices() != 0) return false; + if (ksp.number_of_faces() != 0) return false; + if (ksp.number_of_volumes() != 0) return false; std::cout << std::endl << "--ITERATION #" << iter + 1 << " END!" << std::endl; } @@ -202,19 +203,18 @@ void run_all_tests() { // Number of allowed intersections k. std::vector ks; for (unsigned int k = 1; k <= 6; ++k) { - //ks.push_back(k); + ks.push_back(k); } - ks.push_back(3); - //results = { 9,1,28,56,35,6 }; - //run_test("data/stress-test-1/test-8-rnd-polygons-3-4.off", ks, num_iters, results, all_times, num_tests); - //results = { 16,1,133,315,212,34 }; - //run_test("data/real-data-test/test-10-polygons.off", ks, num_iters, results, all_times, num_tests); - //results = { 10,1,37,77,46,6 }; - //run_test("data/stress-test-4/test-4-rnd-polygons-4-6.off", ks, num_iters, results, all_times, num_tests); - //results = { 10,1,37,77,46,6 }; - //assert(run_test("data/edge-case-test/test-box.off", ks, num_iters, results, all_times, num_tests)); - //results = {7,1,12,20,11,2}; - //assert(run_test("data/edge-case-test/test-flat-bbox-xy-split.off", ks, num_iters, results, all_times, num_tests)); + results = { 9,1,28,56,35,6 }; + run_test("data/stress-test-1/test-8-rnd-polygons-3-4.off", ks, num_iters, results, all_times, num_tests); + results = { 16,1,133,315,212,34 }; + run_test("data/real-data-test/test-10-polygons.off", ks, num_iters, results, all_times, num_tests); + results = { 10,1,37,77,46,6 }; + run_test("data/stress-test-4/test-4-rnd-polygons-4-6.off", ks, num_iters, results, all_times, num_tests); + results = { 10,1,37,77,46,6 }; + run_test("data/edge-case-test/test-box.off", ks, num_iters, results, all_times, num_tests); + results = {7,1,12,20,11,2}; + run_test("data/edge-case-test/test-flat-bbox-xy-split.off", ks, num_iters, results, all_times, num_tests); // Edge case tests. @@ -419,12 +419,14 @@ void run_all_tests() { } } - const auto kernel_name = boost::typeindex::type_id().pretty_name(); + const auto kernel_name = boost::typeindex::type_id().pretty_name(); + const auto intersection_kernel_name = boost::typeindex::type_id().pretty_name(); if (num_tests != 0) { - std::cout << std::endl << kernel_name << + std::cout << std::endl << kernel_name << " with " << intersection_kernel_name << " intersections" << ": ALL " << num_tests << " TESTS SUCCESS!" << std::endl << std::endl; - } else { - std::cout << std::endl << kernel_name << + } + else { + std::cout << std::endl << kernel_name << " with " << intersection_kernel_name << " intersections" << ": ALL " << num_tests << " TESTS FAILED!" << std::endl << std::endl; } } @@ -438,6 +440,6 @@ int main(const int /* argc */, const char** /* argv */) { // Passes all tests except for those when // intersections lead to accumulated errors. - run_all_tests(); + run_all_tests(); return EXIT_SUCCESS; } diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h index d00bddcf6f1..bd5c92a3671 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -251,12 +251,10 @@ CGAL_add_named_parameter(overlap_t, overlap, overlap) CGAL_add_named_parameter(maximum_normal_deviation_t, maximum_normal_deviation, maximum_normal_deviation) // kinetic parameters -CGAL_add_named_parameter(k_intersections_t, k_intersections, k_intersections) -CGAL_add_named_parameter(n_subdivisions_t, n_subdivisions, n_subdivisions) -CGAL_add_named_parameter(enlarge_bbox_ratio_t, enlarge_bbox_ratio, enlarge_bbox_ratio) -CGAL_add_named_parameter(reorient_t, reorient, reorient) -CGAL_add_named_parameter(use_hybrid_mode_t, use_hybrid_mode, use_hybrid_mode) +CGAL_add_named_parameter(bbox_dilation_ratio_t, bbox_dilation_ratio, bbox_dilation_ratio) +CGAL_add_named_parameter(reorient_bbox_t, reorient_bbox, reorient_bbox) CGAL_add_named_parameter(distance_tolerance_t, distance_tolerance, distance_tolerance) +CGAL_add_named_parameter(angle_tolerance_t, angle_tolerance, angle_tolerance) // region growing CGAL_add_named_parameter(k_neighbors_t, k_neighbors, k_neighbors)