From 60511740bfabf16c44d50daf17b17384f85a7aeb Mon Sep 17 00:00:00 2001 From: Dmitry Anisimov Date: Tue, 2 Feb 2021 16:43:57 +0100 Subject: [PATCH] better future points --- .../include/CGAL/KSR_3/Data_structure.h | 206 +++++++++++++----- .../include/CGAL/KSR_3/Initializer.h | 6 +- .../include/CGAL/KSR_3/Polygon_splitter.h | 4 +- .../include/CGAL/KSR_3/Support_plane.h | 25 ++- .../CGAL/Kinetic_shape_reconstruction_3.h | 98 ++++----- 5 files changed, 222 insertions(+), 117 deletions(-) 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 2c3cf9e8f36..3994ad37cf1 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h @@ -284,6 +284,34 @@ public: return m_map_volumes; } + void precompute_iedge_data() { + + for (std::size_t i = 0; i < number_of_support_planes(); ++i) { + auto& unique_iedges = support_plane(i).unique_iedges(); + CGAL_assertion(unique_iedges.size() > 0); + + auto& iedges = this->iedges(i); + auto& ibboxes = this->ibboxes(i); + auto& isegments = this->isegments(i); + + iedges.clear(); + iedges.reserve(unique_iedges.size()); + std::copy(unique_iedges.begin(), unique_iedges.end(), std::back_inserter(iedges)); + unique_iedges.clear(); + + ibboxes.clear(); + isegments.clear(); + + ibboxes.reserve(iedges.size()); + isegments.reserve(iedges.size()); + + for (const auto& iedge : iedges) { + isegments.push_back(segment_2(i, iedge)); + ibboxes.push_back(isegments.back().bbox()); + } + } + } + void set_limit_lines() { m_limit_lines.clear(); @@ -611,19 +639,19 @@ public: for (std::size_t i = 0; i < n; ++i) { const auto& iplanes = m_intersection_graph.intersected_planes(intersections[i].first); for (const std::size_t sp_idx : iplanes) { - support_plane(sp_idx).iedges().erase(intersections[i].first); + support_plane(sp_idx).unique_iedges().erase(intersections[i].first); } const auto edges = m_intersection_graph.split_edge( intersections[i].first, vertices[i]); const auto& iplanes_1 = m_intersection_graph.intersected_planes(edges.first); for (const std::size_t sp_idx : iplanes_1) { - support_plane(sp_idx).iedges().insert(edges.first); + support_plane(sp_idx).unique_iedges().insert(edges.first); } const auto& iplanes_2 = m_intersection_graph.intersected_planes(edges.second); for (const std::size_t sp_idx : iplanes_2) { - support_plane(sp_idx).iedges().insert(edges.second); + support_plane(sp_idx).unique_iedges().insert(edges.second); } const auto new_edge = m_intersection_graph.add_edge( @@ -631,8 +659,8 @@ public: m_intersection_graph.intersected_planes(new_edge).insert(common_planes_idx[i]); m_intersection_graph.set_line(new_edge, map_lines_idx[common_planes_idx[i]]); - support_plane(support_plane_idx).iedges().insert(new_edge); - support_plane(common_planes_idx[i]).iedges().insert(new_edge); + support_plane(support_plane_idx).unique_iedges().insert(new_edge); + support_plane(common_planes_idx[i]).unique_iedges().insert(new_edge); } } @@ -660,7 +688,7 @@ public: } support_plane(support_plane_idx).set_iedge(vertices[i], vertices[(i + 1) % 4], iedge); - support_plane(support_plane_idx).iedges().insert(iedge); + support_plane(support_plane_idx).unique_iedges().insert(iedge); } } @@ -1113,7 +1141,7 @@ public: m_intersection_graph.set_line(iedge, line_idx); for (const auto support_plane_idx : support_planes_idx) { - support_plane(support_plane_idx).iedges().insert(iedge); + support_plane(support_plane_idx).unique_iedges().insert(iedge); } } } @@ -1134,9 +1162,26 @@ public: return m_intersection_graph.incident_edges(ivertex); } - const std::set& iedges(const std::size_t support_plane_idx) const { + const std::vector& iedges(const std::size_t support_plane_idx) const { return support_plane(support_plane_idx).iedges(); } + std::vector& iedges(const std::size_t support_plane_idx) { + return support_plane(support_plane_idx).iedges(); + } + + const std::vector& isegments(const std::size_t support_plane_idx) const { + return support_plane(support_plane_idx).isegments(); + } + std::vector& isegments(const std::size_t support_plane_idx) { + return support_plane(support_plane_idx).isegments(); + } + + const std::vector& ibboxes(const std::size_t support_plane_idx) const { + return support_plane(support_plane_idx).ibboxes(); + } + std::vector& ibboxes(const std::size_t support_plane_idx) { + return support_plane(support_plane_idx).ibboxes(); + } const std::set& intersected_planes(const IEdge& iedge) const { return m_intersection_graph.intersected_planes(iedge); @@ -1761,9 +1806,9 @@ public: Point_2 future_point_a, future_point_b; Vector_2 future_direction_a, future_direction_b; compute_future_point_and_direction( - pvertex_p, source_p, pvertex, prev, future_point_a, future_direction_a); + target_p, pvertex_p, source_p, pvertex, prev, future_point_a, future_direction_a); compute_future_point_and_direction( - pvertex_p, target_p, pvertex, next, future_point_b, future_direction_b); + source_p, pvertex_p, target_p, pvertex, next, future_point_b, future_direction_b); CGAL_assertion(future_direction_a * future_direction_b < FT(0)); CGAL_assertion(future_direction_a != Vector_2()); CGAL_assertion(future_direction_b != Vector_2()); @@ -1893,7 +1938,7 @@ public: } compute_future_point_and_direction( - pvertex_p, target_p, pvertex, pthird, future_point, future_direction); + source_p, pvertex_p, target_p, pvertex, pthird, future_point, future_direction); CGAL_assertion(future_direction != Vector_2()); direction(pvertex) = future_direction; @@ -1927,7 +1972,7 @@ public: if (m_verbose) std::cout << "- swap source and target" << std::endl; compute_future_point_and_direction( - pother_p, target_p, pother, pthird, future_point, future_direction); + source_p, pother_p, target_p, pother, pthird, future_point, future_direction); CGAL_assertion(future_direction != Vector_2()); direction(pother) = future_direction; @@ -2059,7 +2104,7 @@ public: Point_2 future_point; Vector_2 future_direction; compute_future_point_and_direction( - pother_p, target_p, pother, pthird, future_point, future_direction); + source_p, pother_p, target_p, pother, pthird, future_point, future_direction); CGAL_assertion(future_direction != Vector_2()); if (target_pface == null_pface()) { // in case we have 1 pface @@ -2104,7 +2149,7 @@ public: support_plane(pother).set_point(pother.second, future_point); connect(pother, iedge); - CGAL_assertion_msg(false, "TODO: TRANSFER 2, ADD NEW FUTURE POINTS AND DIRECTIONS!"); + // CGAL_assertion_msg(false, "TODO: TRANSFER 2, ADD NEW FUTURE POINTS AND DIRECTIONS!"); } if (m_verbose) { @@ -2400,7 +2445,7 @@ public: // std::cout << "- opoint: " << to_3d(pvertex.first, opoint) << std::endl; CGAL_assertion_msg(ipoint != opoint, "TODO: BACK, HANDLE ZERO LENGTH IEDGE!"); compute_future_point_and_direction( - ipoint, opoint, back, prev, future_point, future_direction); + ipoint, ipoint, opoint, back, prev, future_point, future_direction); CGAL_assertion(future_direction != Vector_2()); // Crop the pvertex. @@ -2515,7 +2560,7 @@ public: // std::cout << "- opoint: " << to_3d(pvertex.first, opoint) << std::endl; CGAL_assertion_msg(ipoint != opoint, "TODO: FRONT, HANDLE ZERO LENGTH IEDGE!"); compute_future_point_and_direction( - ipoint, opoint, front, next, future_point, future_direction); + ipoint, ipoint, opoint, front, next, future_point, future_direction); CGAL_assertion(future_direction != Vector_2()); // Crop the pvertex. @@ -2642,14 +2687,14 @@ public: // std::cout << "- opoint1: " << to_3d(pvertex.first, opoint) << std::endl; CGAL_assertion_msg(ipoint != opoint, "TODO: OPEN, HANDLE ZERO LENGTH IEDGE!"); compute_future_point_and_direction( - ipoint, opoint, front, next, future_point_a, future_direction_a); + ipoint, ipoint, opoint, front, next, future_point_a, future_direction_a); CGAL_assertion(future_direction_a != Vector_2()); opoint = point_2(pvertex.first, opposite(crossed_iedges.back().first, ivertex)); // std::cout << "- opoint2: " << to_3d(pvertex.first, opoint) << std::endl; CGAL_assertion_msg(ipoint != opoint, "TODO: OPEN, HANDLE ZERO LENGTH IEDGE!"); compute_future_point_and_direction( - ipoint, opoint, back, prev, future_point_b, future_direction_b); + ipoint, ipoint, opoint, back, prev, future_point_b, future_direction_b); CGAL_assertion(future_direction_b != Vector_2()); // Crop the pvertex. @@ -2816,6 +2861,13 @@ public: std::cout << "- adding new pface: " << std::endl; } + std::cout << "idx: " << idx << std::endl; + for (const auto& pvertex : pvertices) { + if (pvertex != null_pvertex()) { + std::cout << "pv: " << point_3(pvertex) << std::endl; + } + } + const auto& pv1 = pvertices[idx]; CGAL_assertion(pv1 != null_pvertex()); if (m_verbose) { @@ -2858,7 +2910,7 @@ public: std::cout << "- opoint: " << to_3d(pvertex.first, opoint) << std::endl; CGAL_assertion_msg(ipoint != opoint, "TODO: CREATE PVERTEX, HANDLE ZERO LENGTH IEDGE!"); compute_future_point_and_direction( - ipoint, opoint, pother, pthird, future_point, future_direction); + ipoint, ipoint, opoint, pother, pthird, future_point, future_direction); CGAL_assertion(future_direction != Vector_2()); const auto propagated = add_pvertex(pvertex.first, future_point); @@ -3810,6 +3862,7 @@ public: } const bool check_integrity( + const bool is_initialized = true, const bool check_simplicity = false, const bool check_convexity = false) const { @@ -3820,15 +3873,33 @@ public: return false; } - for (const auto& iedge : this->iedges(i)) { - const auto& iplanes = this->intersected_planes(iedge); - if (iplanes.find(i) == iplanes.end()) { + if (is_initialized) { + const auto& iedges = this->iedges(i); + CGAL_assertion(iedges.size() > 0); + for (const auto& iedge : iedges) { + const auto& iplanes = this->intersected_planes(iedge); + if (iplanes.find(i) == iplanes.end()) { - const std::string msg = "ERROR: SUPPORT PLANE " + std::to_string(i) + - " IS INTERSECTED BY " + str(iedge) + - " BUT IT CLAIMS IT DOES NOT INTERSECT IT!"; - CGAL_assertion_msg(false, msg.c_str()); - return false; + const std::string msg = "ERROR: SUPPORT PLANE " + std::to_string(i) + + " IS INTERSECTED BY " + str(iedge) + + " BUT IT CLAIMS IT DOES NOT INTERSECT IT!"; + CGAL_assertion_msg(false, msg.c_str()); + return false; + } + } + } else { + const auto& iedges = support_plane(i).unique_iedges(); + CGAL_assertion(iedges.size() > 0); + for (const auto& iedge : iedges) { + const auto& iplanes = this->intersected_planes(iedge); + if (iplanes.find(i) == iplanes.end()) { + + const std::string msg = "ERROR: SUPPORT PLANE " + std::to_string(i) + + " IS INTERSECTED BY " + str(iedge) + + " BUT IT CLAIMS IT DOES NOT INTERSECT IT!"; + CGAL_assertion_msg(false, msg.c_str()); + return false; + } } } } @@ -3837,14 +3908,28 @@ public: const auto& iplanes = this->intersected_planes(iedge); for (const auto support_plane_idx : iplanes) { - const auto& sp_iedges = this->iedges(support_plane_idx); - if (sp_iedges.find(iedge) == sp_iedges.end()) { + if (is_initialized) { + const auto& sp_iedges = this->iedges(support_plane_idx); + CGAL_assertion(sp_iedges.size() > 0); + if (std::find(sp_iedges.begin(), sp_iedges.end(), iedge) == sp_iedges.end()) { - const std::string msg = "ERROR: IEDGE " + str(iedge) + - " INTERSECTS SUPPORT PLANE " + std::to_string(support_plane_idx) + - " BUT IT CLAIMS IT IS NOT INTERSECTED BY IT!"; - CGAL_assertion_msg(false, msg.c_str()); - return false; + const std::string msg = "ERROR: IEDGE " + str(iedge) + + " INTERSECTS SUPPORT PLANE " + std::to_string(support_plane_idx) + + " BUT IT CLAIMS IT IS NOT INTERSECTED BY IT!"; + CGAL_assertion_msg(false, msg.c_str()); + return false; + } + } else { + const auto& sp_iedges = support_plane(support_plane_idx).unique_iedges(); + CGAL_assertion(sp_iedges.size() > 0); + if (sp_iedges.find(iedge) == sp_iedges.end()) { + + const std::string msg = "ERROR: IEDGE " + str(iedge) + + " INTERSECTS SUPPORT PLANE " + std::to_string(support_plane_idx) + + " BUT IT CLAIMS IT IS NOT INTERSECTED BY IT!"; + CGAL_assertion_msg(false, msg.c_str()); + return false; + } } } } @@ -3948,8 +4033,6 @@ public: } // First, traverse only boundary volumes. - // TODO: SORT HERE BY PFACE AREA! - // Actually, we should sort by both number of edges and area! bool is_found_new_volume = false; std::size_t volume_size = 0; int num_volumes = 0; @@ -3985,8 +4068,6 @@ public: } } - // TODO: SORT HERE BY PFACE AREA! - // Actually, we should sort by both number of edges and area! std::sort(other_pfaces.begin(), other_pfaces.end(), [&](const PFace& pface1, const PFace& pface2) -> bool { const auto pedges1 = pedges_of_pface(pface1); @@ -4610,10 +4691,22 @@ private: *************************************/ const std::pair > compute_future_point( - const Point_2& q0, const Point_2& q1, + const Point_2& source, const Point_2& query, const Point_2& target, const PVertex& pv0, const PVertex& pv1) const { + auto q0 = query; + const auto& q1 = target; + const FT tol = KSR::tolerance(); + const FT sq_dist = CGAL::squared_distance(query, target); + if (sq_dist < tol) { + if (m_verbose) { + std::cout << "- warning: query is almost equal to target" << std::endl; + std::cout << "- replacing query with source" << std::endl; + } + q0 = source; + } + const auto q2 = point_2(pv0, m_current_time + FT(1)); const auto q3 = point_2(pv1, m_current_time + FT(1)); @@ -4636,7 +4729,7 @@ private: const FT sq_d1 = (x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0); CGAL_assertion(sq_d1 >= FT(0)); - if (m_verbose) std::cout << "sq d1: " << sq_d1 << std::endl; + // if (m_verbose) std::cout << "sq d1: " << sq_d1 << std::endl; CGAL_assertion_msg(sq_d1 >= tol, "TODO: FUTURE POINT, HANDLE ZERO CURRENT EDGE!"); const FT a2 = a1 / sq_d1, b2 = b1 / sq_d1, c2 = c1 / sq_d1; @@ -4683,23 +4776,22 @@ private: } const Vector_2 compute_future_direction( - const Point_2& q0, const Point_2& q1, + const Point_2& source, const Point_2& target, const PVertex& pv0, const PVertex& pv1) const { auto pinit = point_2(pv0); - const Line_2 iedge_line(q0, q1); + const Line_2 iedge_line(source, target); pinit = iedge_line.projection(pinit); - const auto res = compute_future_point(q0, q1, pv0, pv1); + const auto res = compute_future_point(source, source, target, pv0, pv1); const auto& future_point = res.first; - CGAL_assertion_msg(future_point != pinit, - "TODO: ZERO LENGTH FUTURE DIRECTION!"); + CGAL_assertion_msg(future_point != pinit, "ERROR: ZERO LENGTH FUTURE DIRECTION!"); const Vector_2 future_direction(pinit, future_point); if (m_verbose) { std::cout.precision(20); - // std::cout << "- future point: " << to_3d(pv0.first, future_point) << std::endl; - // std::cout << "- future direction: " << future_direction << std::endl; + // std::cout << "- future point: " << to_3d(pv0.first, future_point) << std::endl; + // std::cout << "- future direction: " << future_direction << std::endl; } // CGAL_assertion_msg(false, "TODO: COMPUTE FUTURE DIRECTION!"); @@ -4707,24 +4799,24 @@ private: } void compute_future_point_and_direction( - const Point_2& q0, const Point_2& q1, + const Point_2& source, const Point_2& query, const Point_2& target, const PVertex& pv0, const PVertex& pv1, Point_2& future_point, Vector_2& future_direction) const { - const auto& pinit = q0; - const auto res = compute_future_point(q0, q1, pv0, pv1); + const auto& pinit = query; + const auto res = compute_future_point(source, query, target, pv0, pv1); future_point = res.first; if (m_verbose) { - std::cout << - "- w1: " << res.second.first << ";" << - " w2: " << res.second.second << std::endl; + std::cout << "-" << + " w1: " << res.second.first << ";" << + " w2: " << res.second.second << std::endl; std::cout << "- future point: " << to_3d(pv0.first, future_point) << std::endl; } - CGAL_assertion_msg(future_point != pinit, - "TODO: ZERO LENGTH FUTURE DIRECTION!"); - CGAL_assertion_msg(res.second.first <= FT(1), "TODO: W1, WRONG ORIENTATION!"); - CGAL_assertion_msg(res.second.second >= FT(0), "TODO: W2, WRONG ORIENTATION!"); + + CGAL_assertion_msg(future_point != pinit, "ERROR: ZERO LENGTH FUTURE DIRECTION!"); + CGAL_assertion_msg(res.second.first <= FT(1), "ERROR: W1, WRONG ORIENTATION!"); + CGAL_assertion_msg(res.second.second >= FT(0), "ERROR: W2, WRONG ORIENTATION!"); future_direction = Vector_2(pinit, future_point); future_point = pinit - m_current_time * future_direction; diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h index 75e8b61e8b9..1b2b8db70a6 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h @@ -103,9 +103,9 @@ public: // KSR_3::dump_segmented_edges(m_data, "init"); } - CGAL_assertion(m_data.check_integrity()); + CGAL_assertion(m_data.check_integrity(false, true, true)); make_polygons_intersection_free(); - CGAL_assertion(m_data.check_integrity()); + CGAL_assertion(m_data.check_integrity(false, true, true)); set_k_intersections(k); if (m_verbose) std::cout << "done" << std::endl; @@ -135,7 +135,7 @@ public: m_data.convert(ds); m_data.clear(); - CGAL_assertion(ds.check_integrity()); + CGAL_assertion(ds.check_integrity(false, true, true)); CGAL_assertion(ds.check_bbox()); } diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Polygon_splitter.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Polygon_splitter.h index 36f9f1e0a89..c115e147ba5 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Polygon_splitter.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Polygon_splitter.h @@ -312,7 +312,7 @@ private: } // Then, add intersection vertices + constraints. - const auto& iedges = m_data.iedges(support_plane_idx); + const auto& iedges = m_data.support_plane(support_plane_idx).unique_iedges(); for (const auto& iedge : iedges) { const auto source = m_data.source(iedge); const auto target = m_data.target(iedge); @@ -666,7 +666,7 @@ private: std::vector& bbox) const { CGAL_assertion(support_plane_idx >= 6); - const auto& iedges = m_data.iedges(support_plane_idx); + const auto& iedges = m_data.support_plane(support_plane_idx).unique_iedges(); std::vector points; points.reserve(iedges.size() * 2); 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 1ebf7a5be0e..15d2cdbd5d7 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Support_plane.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Support_plane.h @@ -52,6 +52,7 @@ public: using Mesh = CGAL::Surface_mesh; using Intersection_graph = KSR_3::Intersection_graph; + using Bbox_2 = CGAL::Bbox_2; using IVertex = typename Intersection_graph::Vertex_descriptor; using IEdge = typename Intersection_graph::Edge_descriptor; @@ -84,7 +85,10 @@ private: F_uint_map k_map; V_original_map v_original_map; V_time_map v_time_map; - std::set iedges; + std::set unique_iedges; + std::vector iedges; + std::vector isegments; + std::vector ibboxes; unsigned int k; }; @@ -294,10 +298,10 @@ public: sp.data().k_map[fi] = m_data->k_map[face]; } - sp.data().iedges.clear(); - for (const auto& iedge : m_data->iedges) { + sp.data().unique_iedges.clear(); + for (const auto& iedge : m_data->unique_iedges) { CGAL_assertion(iedge != IG::null_iedge()); - sp.data().iedges.insert(emap.at(iedge)); + sp.data().unique_iedges.insert(emap.at(iedge)); } } @@ -525,8 +529,17 @@ public: return (m_data->direction[vi] == CGAL::NULL_VECTOR); } - const std::set& iedges() const { return m_data->iedges; } - std::set& iedges() { return m_data->iedges; } + const std::set& unique_iedges() const { return m_data->unique_iedges; } + std::set& unique_iedges() { return m_data->unique_iedges; } + + const std::vector& iedges() const { return m_data->iedges; } + std::vector& iedges() { return m_data->iedges; } + + const std::vector& isegments() const { return m_data->isegments; } + std::vector& isegments() { return m_data->isegments; } + + const std::vector& ibboxes() const { return m_data->ibboxes; } + std::vector& ibboxes() { return m_data->ibboxes; } const Point_2 to_2d(const Point_3& point) const { return m_data->plane.to_2d(point); 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 d98cddb10a7..777811d0364 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h +++ b/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h @@ -32,6 +32,7 @@ #include #include #include +#include // Internal includes. #include @@ -74,6 +75,7 @@ private: using Polygon_mesh = CGAL::Surface_mesh; using Vertex_index = typename Polygon_mesh::Vertex_index; + using Timer = CGAL::Real_timer; private: const bool m_debug; @@ -106,6 +108,7 @@ public: const PolygonMap polygon_map, const NamedParameters& np) { + Timer timer; const unsigned int k = parameters::choose_parameter( parameters::get_parameter(np, internal_np::k_intersections), 1); unsigned int n = parameters::choose_parameter( @@ -149,11 +152,16 @@ public: std::cout << "* reorient: " << is_reorient << std::endl; } + timer.reset(); + timer.start(); const FT time_step = static_cast(m_initializer.initialize( input_range, polygon_map, k, CGAL::to_double(enlarge_bbox_ratio), reorient)); m_initializer.convert(m_data); m_data.set_limit_lines(); - CGAL_assertion(m_data.check_integrity()); + m_data.precompute_iedge_data(); + CGAL_assertion(m_data.check_integrity(true, true, true)); + timer.stop(); + const double time_to_initialize = timer.time(); if (k == 0) { CGAL_warning_msg(k > 0, @@ -177,6 +185,8 @@ public: std::cout << "* propagation started" << std::endl; } + timer.reset(); + timer.start(); std::size_t num_iterations = 0; m_min_time = FT(0); m_max_time = time_step; @@ -202,6 +212,8 @@ public: return false; } } + timer.stop(); + const double time_to_partition = timer.time(); if (m_verbose) { if (m_verbose && !m_debug) std::cout << std::endl; @@ -209,22 +221,34 @@ public: std::cout << "* number of events: " << global_iteration << std::endl; } + timer.reset(); + timer.start(); if (m_verbose) std::cout << std::endl << "--- FINALIZING PARTITION:" << std::endl; if (m_debug) dump(m_data, "jiter-final-a-result"); m_data.finalize(); if (m_verbose) std::cout << "* checking final mesh integrity ..."; - CGAL_assertion(m_data.check_integrity()); + CGAL_assertion(m_data.check_integrity(true, true, true)); if (m_verbose) std::cout << " done" << std::endl; if (m_debug) dump(m_data, "jiter-final-b-result"); - // std::cout << std::endl << "CLEANING SUCCESS!" << std::endl << std::endl; // exit(EXIT_SUCCESS); - if (m_verbose) std::cout << "* getting volumes ..." << std::endl; m_data.create_polyhedra(); + timer.stop(); + const double time_to_finalize = timer.time(); if (m_verbose) { std::cout << "* found " << m_data.number_of_volumes(-1) << " volumes" << std::endl; } + + if (m_verbose) std::cout << std::endl << "--- TIMING (sec.):" << std::endl; + const double total_time = + time_to_initialize + time_to_partition + time_to_finalize; + if (m_verbose) { + std::cout << "* initialization: " << time_to_initialize << std::endl; + std::cout << "* partition: " << time_to_partition << std::endl; + std::cout << "* finalization: " << time_to_finalize << std::endl; + std::cout << "* total time: " << total_time << std::endl; + } return true; } @@ -638,12 +662,10 @@ private: m_data.update_positions(m_max_time); bool still_running = false; - - std::vector iedges; - std::vector segments; - std::vector bboxes; for (std::size_t i = 0; i < m_data.number_of_support_planes(); ++i) { - initialize_search_structures(i, iedges, segments, bboxes); + const auto& iedges = m_data.iedges(i); + const auto& segments = m_data.isegments(i); + const auto& bboxes = m_data.ibboxes(i); for (const auto pvertex : m_data.pvertices(i)) { if (compute_events_of_pvertex(pvertex, iedges, segments, bboxes)) { still_running = true; @@ -654,37 +676,16 @@ private: return still_running; } - void initialize_search_structures( - const std::size_t i, - std::vector& iedges, - std::vector& segments, - std::vector& bboxes) { - - iedges.clear(); - segments.clear(); - bboxes.clear(); - - // To get random access, copy in vector (suboptimal to do this - // all the time, maybe this should be done once and for all and - // replace the set). - iedges.reserve(m_data.iedges(i).size()); - std::copy(m_data.iedges(i).begin(), m_data.iedges(i).end(), std::back_inserter(iedges)); - - // Precompute segments and bboxes. - segments.reserve(iedges.size()); - bboxes.reserve(iedges.size()); - for (const auto& iedge : iedges) { - segments.push_back(m_data.segment_2(i, iedge)); - bboxes.push_back(segments.back().bbox()); - } - } - const bool compute_events_of_pvertex( const PVertex& pvertex, const std::vector& iedges, const std::vector& segments, const std::vector& bboxes) { + CGAL_assertion(iedges.size() > 0); + CGAL_assertion(iedges.size() == segments.size()); + CGAL_assertion(iedges.size() == bboxes.size()); + std::cout.precision(20); if (m_data.is_frozen(pvertex)) { return false; @@ -1099,6 +1100,15 @@ private: "ERROR: CONSTRAINED PVERTEX MEETS IEDGE! WHAT IS WRONG?"); } + void apply_event_pvertices_meet_ivertex( + const PVertex& pvertex, const PVertex& pother, const IVertex& ivertex, const Event& event) { + + CGAL_assertion( m_data.has_iedge(pvertex)); + CGAL_assertion(!m_data.has_iedge(pother)); + CGAL_assertion_msg(false, + "ERROR: PVERTICES MEET IVERTEX! IT SHOULD NOT EVER HAPPEN!"); + } + // VALID EVENTS! void apply_event_pvertex_meets_ivertex( const PVertex& pvertex, const IVertex& ivertex, const Event& event) { @@ -1137,14 +1147,6 @@ private: // CGAL_assertion_msg(false, "TODO: PVERTEX MEETS IVERTEX!"); } - void apply_event_pvertices_meet_ivertex( - const PVertex& pvertex, const PVertex& pother, const IVertex& ivertex, const Event& event) { - - CGAL_assertion( m_data.has_iedge(pvertex)); - CGAL_assertion(!m_data.has_iedge(pother)); - CGAL_assertion_msg(false, "TODO: PVERTICES MEET IVERTEX!"); - } - void apply_event_unconstrained_pvertex_meets_ivertex( const PVertex& pvertex, const IVertex& ivertex, const Event& event) { @@ -1295,7 +1297,7 @@ private: const bool check_pvertex_meets_iedge_global_k( const PVertex& pvertex, const IEdge& iedge) { - if (m_verbose) { + if (m_debug) { std::cout << "- k intersections before: " << m_data.k(pvertex.first) << std::endl; } @@ -1336,7 +1338,7 @@ private: const bool check_pedge_meets_iedge_global_k( const PVertex& pvertex, const PVertex& pother, const IEdge& iedge) { - if (m_verbose) { + if (m_debug) { std::cout << "- k intersections before: " << m_data.k(pvertex.first) << std::endl; } @@ -1388,13 +1390,11 @@ private: m_min_time = m_data.current_time(); m_data.update_positions(m_max_time); - std::vector iedges; - std::vector segments; - std::vector bboxes; - const auto& pfront = pvertices.front(); CGAL_assertion(pfront != Data_structure::null_pvertex()); - initialize_search_structures(pfront.first, iedges, segments, bboxes); + const auto& iedges = m_data.iedges(pfront.first); + const auto& segments = m_data.isegments(pfront.first); + const auto& bboxes = m_data.ibboxes(pfront.first); for (const auto& pvertex : pvertices) { if (pvertex == Data_structure::null_pvertex()) continue;