From 738763ae3272228cbeff0582a802c9dc63f29bbb Mon Sep 17 00:00:00 2001 From: Sven Oesau Date: Mon, 21 Oct 2024 14:26:30 +0200 Subject: [PATCH 1/3] cleaning up dependencies removing unused public KSP::Index replacing volume calculation using Delaunay_triangulation_3 --- .../include/CGAL/KSP/utils.h | 155 ------------------ .../include/CGAL/Kinetic_space_partition_3.h | 5 +- .../ksr_building.cpp | 2 +- .../ksr_parameters.cpp | 11 +- .../include/CGAL/KSR_3/Graphcut.h | 2 - .../CGAL/Kinetic_surface_reconstruction_3.h | 100 +++++++---- 6 files changed, 76 insertions(+), 199 deletions(-) diff --git a/Kinetic_space_partition/include/CGAL/KSP/utils.h b/Kinetic_space_partition/include/CGAL/KSP/utils.h index 64a2d9972965..5b43a1a28b58 100644 --- a/Kinetic_space_partition/include/CGAL/KSP/utils.h +++ b/Kinetic_space_partition/include/CGAL/KSP/utils.h @@ -29,17 +29,11 @@ #include // CGAL includes. -#include -#include -#include #include -#include #include #include #include -#include -#include #include // Boost includes. @@ -70,22 +64,6 @@ decltype(auto) distance(const Point_d& p, const Point_d& q) { return static_cast(CGAL::approximate_sqrt(sq_dist)); } -// Project 3D point onto 2D plane. -template -typename Kernel_traits::Kernel::Point_2 -point_2_from_point_3(const Point_3& point_3) { - return typename Kernel_traits::Kernel::Point_2( - point_3.x(), point_3.y()); -} - -// Get 3D point from a 2D point. -template -typename Kernel_traits::Kernel::Point_3 -point_3_from_point_2(const Point_2& point_2) { - return typename Kernel_traits::Kernel::Point_3( - point_2.x(), point_2.y(), typename Kernel_traits::Kernel::FT(0)); -} - // Normalize vector. template inline const Vector_d normalize(const Vector_d& v) { @@ -96,7 +74,6 @@ inline const Vector_d normalize(const Vector_d& v) { return v / static_cast(CGAL::approximate_sqrt(dot_product)); } - // Intersections. Used only in the 2D version. // For the 3D version, see conversions.h! template @@ -121,40 +98,6 @@ inline const ResultType intersection(const Type1& t1, const Type2& t2) { return out; } -// Get boundary points from a set of points. -template -void boundary_points_on_line_2( - const std::vector& input_range, - const std::vector& indices, - const Line_2& line, Point_2& p, Point_2& q) { - - using Traits = typename Kernel_traits::Kernel; - using FT = typename Traits::FT; - using Vector_2 = typename Traits::Vector_2; - - FT min_proj_value = (std::numeric_limits::max)(); - FT max_proj_value = -min_proj_value; - - const auto ref_vector = line.to_vector(); - const auto& ref_point = input_range[indices.front()]; - - for (const std::size_t index : indices) { - const auto& query = input_range[index]; - const auto point = line.projection(query); - const Vector_2 curr_vector(ref_point, point); - const FT value = CGAL::scalar_product(curr_vector, ref_vector); - - if (value < min_proj_value) { - min_proj_value = value; - p = point; - } - if (value > max_proj_value) { - max_proj_value = value; - q = point; - } - } -} - // Angles. // Converts radians to degrees. @@ -202,104 +145,6 @@ angle_2(const Direction_2& dir1, const Direction_2& dir2) { return CGAL::abs(convert_angle_2(angle_2)); } -// Classes. -template -class Indexer { -public: - std::size_t operator()(const IVertex& ivertex) { - const auto pair = m_indices.insert( - std::make_pair(ivertex, m_indices.size())); - const auto& item = pair.first; - const std::size_t idx = item->second; - return idx; - } - void clear() { m_indices.clear(); } - -private: - std::map m_indices; -}; - -template< - typename GeomTraits, - typename InputRange, - typename NeighborQuery> -class Estimate_normals_2 { - -public: - using Traits = GeomTraits; - using Input_range = InputRange; - using Neighbor_query = NeighborQuery; - - using Kernel = Traits; - using FT = typename Kernel::FT; - using Vector_2 = typename Kernel::Vector_2; - using Line_2 = typename Kernel::Line_2; - - using Indices = std::vector; - - using IK = CGAL::Exact_predicates_inexact_constructions_kernel; - using IPoint_2 = typename IK::Point_2; - using ILine_2 = typename IK::Line_2; - using Converter = CGAL::Cartesian_converter; - - Estimate_normals_2( - const Input_range& input_range, - const Neighbor_query& neighbor_query) : - m_input_range(input_range), - m_neighbor_query(neighbor_query) { - - CGAL_precondition(input_range.size() > 0); - } - - void get_normals(std::vector& normals) const { - - normals.clear(); - normals.reserve(m_input_range.size()); - - Indices neighbors; - for (std::size_t i = 0; i < m_input_range.size(); ++i) { - neighbors.clear(); - m_neighbor_query(i, neighbors); - const auto line = fit_line(neighbors); - auto normal = line.to_vector(); - normal = normal.perpendicular(CGAL::COUNTERCLOCKWISE); - normal = normalize(normal); - normals.push_back(normal); - } - CGAL_assertion(normals.size() == m_input_range.size()); - } - -private: - const Input_range& m_input_range; - const Neighbor_query& m_neighbor_query; - const Converter m_converter; - - const Line_2 fit_line(const Indices& indices) const { - CGAL_assertion(indices.size() > 0); - - std::vector points; - points.reserve(indices.size()); - for (const std::size_t index : indices) { - const auto& point = get(m_neighbor_query.point_map(), index); - points.push_back(m_converter(point)); - } - CGAL_assertion(points.size() == indices.size()); - - ILine_2 fitted_line; - IPoint_2 fitted_centroid; - CGAL::linear_least_squares_fitting_2( - points.begin(), points.end(), - fitted_line, fitted_centroid, - CGAL::Dimension_tag<0>()); - - const Line_2 line( - static_cast(fitted_line.a()), - static_cast(fitted_line.b()), - static_cast(fitted_line.c())); - return line; - } -}; - #endif } // namespace internal diff --git a/Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h b/Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h index 46c345a4c875..e149cbfdc729 100644 --- a/Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h +++ b/Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h @@ -71,8 +71,6 @@ class Kinetic_space_partition_3 { using Point_3 = typename Kernel::Point_3; - using Index = std::pair; - /*! identifies the support of a face in the exported linear cell complex, which is either an input polygon, identified by a non-negative number, a side of the bounding box in the rotated coordinate system or a face of the octree used to partition the scene. */ @@ -126,6 +124,8 @@ class Kinetic_space_partition_3 { using Triangle_2 = typename Kernel::Triangle_2; using Transform_3 = CGAL::Aff_transformation_3; + using Index = std::pair; + using Data_structure = KSP_3::internal::Data_structure; using IVertex = typename Data_structure::IVertex; @@ -730,6 +730,7 @@ class Kinetic_space_partition_3 { m_partition_nodes[i].m_data->face_to_volumes().clear(); } + if (m_parameters.verbose) std::cout << "ksp v: " << m_partition_nodes[0].m_data->vertices().size() << " f: " << m_partition_nodes[0].face2vertices.size() << " vol: " << m_volumes.size() << std::endl; return; diff --git a/Kinetic_surface_reconstruction/examples/Kinetic_surface_reconstruction/ksr_building.cpp b/Kinetic_surface_reconstruction/examples/Kinetic_surface_reconstruction/ksr_building.cpp index 50665db5ad9e..61494d81ced7 100644 --- a/Kinetic_surface_reconstruction/examples/Kinetic_surface_reconstruction/ksr_building.cpp +++ b/Kinetic_surface_reconstruction/examples/Kinetic_surface_reconstruction/ksr_building.cpp @@ -39,7 +39,7 @@ int main(const int, const char**) { std::vector vtx; std::vector > polylist; - std::vector lambdas{0.3, 0.5, 0.7, 0.8, 0.9, 0.95, 0.99}; + std::vector lambdas{0.5, 0.7, 0.8, 0.9}; bool non_empty = false; diff --git a/Kinetic_surface_reconstruction/examples/Kinetic_surface_reconstruction/ksr_parameters.cpp b/Kinetic_surface_reconstruction/examples/Kinetic_surface_reconstruction/ksr_parameters.cpp index 9ac52e26b9cc..9e624e3cc66f 100644 --- a/Kinetic_surface_reconstruction/examples/Kinetic_surface_reconstruction/ksr_parameters.cpp +++ b/Kinetic_surface_reconstruction/examples/Kinetic_surface_reconstruction/ksr_parameters.cpp @@ -167,13 +167,6 @@ int main(const int argc, const char** argv) { // Algorithm. KSR ksr(point_set, param); - FT max_d, max_dev; - std::size_t num; - ksr.estimate_detection_parameters(max_d, max_dev, num); - std::cout << "d: " << max_d << std::endl; - std::cout << "dev: " << max_dev << std::endl; - std::cout << "num: " << num << std::endl; - Timer timer; timer.start(); std::size_t num_shapes = ksr.detect_planar_shapes(param); @@ -206,7 +199,7 @@ int main(const int argc, const char** argv) { FT after_reconstruction = timer.time(); if (polylist.size() > 0) - CGAL::IO::write_polygon_soup("building_c_" + std::to_string(parameters.graphcut_lambda) + (parameters.use_ground ? "_g" : "_") + ".off", vtx, polylist); + CGAL::IO::write_polygon_soup("polylist_" + std::to_string(parameters.graphcut_lambda) + (parameters.use_ground ? "_g" : "_") + ".off", vtx, polylist); timer.stop(); const FT time = static_cast(timer.time()); @@ -230,7 +223,7 @@ int main(const int argc, const char** argv) { if (polylist.size() > 0) { non_empty = true; - CGAL::IO::write_polygon_soup("building_c_" + std::to_string(l) + (parameters.use_ground ? "_g" : "_") + ".off", vtx, polylist); + CGAL::IO::write_polygon_soup("polylist_" + std::to_string(l) + (parameters.use_ground ? "_g" : "_") + ".off", vtx, polylist); } } diff --git a/Kinetic_surface_reconstruction/include/CGAL/KSR_3/Graphcut.h b/Kinetic_surface_reconstruction/include/CGAL/KSR_3/Graphcut.h index 4c9453d4cb2d..6a3e129adcf2 100644 --- a/Kinetic_surface_reconstruction/include/CGAL/KSR_3/Graphcut.h +++ b/Kinetic_surface_reconstruction/include/CGAL/KSR_3/Graphcut.h @@ -21,8 +21,6 @@ #include #include #include -#include -#include #include // Internal includes. diff --git a/Kinetic_surface_reconstruction/include/CGAL/Kinetic_surface_reconstruction_3.h b/Kinetic_surface_reconstruction/include/CGAL/Kinetic_surface_reconstruction_3.h index 6cc639e6749e..05ce37d82adc 100644 --- a/Kinetic_surface_reconstruction/include/CGAL/Kinetic_surface_reconstruction_3.h +++ b/Kinetic_surface_reconstruction/include/CGAL/Kinetic_surface_reconstruction_3.h @@ -586,8 +586,6 @@ class Kinetic_surface_reconstruction_3 { typedef typename CDT::Finite_edges_iterator Finite_edges_iterator; typedef typename CDT::Finite_faces_iterator Finite_faces_iterator; - //using Visibility = KSR_3::Visibility; - using Index = typename KSP::Index; using Face_attribute = typename LCC::Base::template Attribute_descriptor<2>::type; using Volume_attribute = typename LCC::Base::template Attribute_descriptor<3>::type; @@ -1711,6 +1709,74 @@ class Kinetic_surface_reconstruction_3 { m_face_area_lcc[i] = m_face_area_lcc[i] * 2.0 * m_total_inliers / total_area; } + FT area(typename LCC::Dart_descriptor face_dart, Plane_3 &pl, std::vector *tris = nullptr) { + std::vector face; + From_exact from_exact; + + Dart_descriptor n = face_dart; + do { + face.push_back(from_exact(m_lcc.point(n))); + n = m_lcc.beta(n, 0); + } while (n != face_dart); + + pl = from_exact(m_lcc.template info<2>(face_dart).plane); + + Delaunay_2 tri; + for (const Point_3& p : face) + tri.insert(pl.to_2d(p)); + + FT face_area = 0; + + // Get area + for (auto fit = tri.finite_faces_begin(); fit != tri.finite_faces_end(); ++fit) { + const typename Kernel::Triangle_2 triangle( + fit->vertex(0)->point(), + fit->vertex(1)->point(), + fit->vertex(2)->point()); + face_area += triangle.area(); + if (tris) + tris->push_back(typename Kernel::Triangle_3(pl.to_3d(fit->vertex(0)->point()), pl.to_3d(fit->vertex(1)->point()), pl.to_3d(fit->vertex(2)->point()))); + } + + return face_area; + } + + FT volume(typename LCC::Dart_descriptor volume_dart) { + FT x = 0, y = 0, z = 0; + std::size_t count = 0; + From_exact from_exact; + + // Collect vertices to obtain point on the inside. + for (auto& fd : m_lcc.template one_dart_per_incident_cell<2, 3>(volume_dart)) { + typename LCC::Dart_descriptor fdh = m_lcc.dart_descriptor(fd); + + for (const auto& vd : m_lcc.template one_dart_per_incident_cell<0, 2>(fdh)) { + const auto &p = from_exact(m_lcc.point(m_lcc.dart_descriptor(vd))); + x += p.x(); + y += p.y(); + z += p.z(); + count++; + } + } + + Point_3 center(x / count, y / count, z / count); + + FT vol = 0; + // Second iteration for computing the area of each face and the volume spanned with the center point. + for (auto& fd : m_lcc.template one_dart_per_incident_cell<2, 3>(volume_dart)) { + typename LCC::Dart_descriptor fdh = m_lcc.dart_descriptor(fd); + + Plane_3 plane; + FT a = area(fdh, plane); + Vector_3 n = plane.orthogonal_vector(); + + FT distance = CGAL::abs((plane.point() - center) * n); + vol += distance * a / 3.0; + } + + return vol; + } + void count_volume_votes_lcc() { // const int debug_volume = -1; FT total_volume = 0; @@ -1744,8 +1810,6 @@ class Kinetic_surface_reconstruction_3 { const auto& point = get(m_point_map, p); const auto& normal = get(m_normal_map, p); -// count_points++; - for (std::size_t j = 0; j < idx; j++) { const Vector_3 vec(point, c[j]); const FT dot_product = vec * normal; @@ -1764,30 +1828,11 @@ class Kinetic_surface_reconstruction_3 { } } - for (auto &d : m_lcc.template one_dart_per_cell<3>()) { + for (auto& d : m_lcc.template one_dart_per_cell<3>()) { typename LCC::Dart_descriptor dh = m_lcc.dart_descriptor(d); - std::vector volume_vertices; - std::size_t volume_index = m_lcc.template info<3>(dh).volume_id; - - // Collect all vertices of volume to calculate volume - for (auto &fd : m_lcc.template one_dart_per_incident_cell<2, 3>(dh)) { - typename LCC::Dart_descriptor fdh = m_lcc.dart_descriptor(fd); - - for (const auto &vd : m_lcc.template one_dart_per_incident_cell<0, 2>(fdh)) - volume_vertices.push_back(from_exact(m_lcc.point(m_lcc.dart_descriptor(vd)))); - } - - Delaunay_3 tri; - for (const Point_3& p : volume_vertices) - tri.insert(p); - - m_volumes[volume_index] = FT(0); - for (auto cit = tri.finite_cells_begin(); cit != tri.finite_cells_end(); ++cit) { - const auto& tet = tri.tetrahedron(cit); - m_volumes[volume_index] += tet.volume(); - } + m_volumes[volume_index] = volume(dh); total_volume += m_volumes[volume_index]; } @@ -1795,9 +1840,6 @@ class Kinetic_surface_reconstruction_3 { // Normalize volumes for (FT& v : m_volumes) v /= total_volume; - -// for (std::size_t i = 0; i < m_volumes.size(); i++) -// std::cout << i << ": " << m_cost_matrix[0][i] << " o: " << m_cost_matrix[1][i] << " v: " << m_volumes[i] << std::endl; } template @@ -2002,8 +2044,6 @@ class Kinetic_surface_reconstruction_3 { } void map_points_to_faces(const std::size_t polygon_index, const std::vector& pts, std::vector > >& face_to_points) { - std::vector faces; - if (polygon_index >= m_kinetic_partition.input_planes().size()) assert(false); From 00c73407ce6cc50882ef93fd26939ebe96b5d05e Mon Sep 17 00:00:00 2001 From: Sven Oesau Date: Mon, 21 Oct 2024 15:24:07 +0200 Subject: [PATCH 2/3] Removing Triangulation_3 --- .../include/CGAL/Kinetic_surface_reconstruction_3.h | 3 --- .../package_info/Kinetic_surface_reconstruction/dependencies | 1 - 2 files changed, 4 deletions(-) diff --git a/Kinetic_surface_reconstruction/include/CGAL/Kinetic_surface_reconstruction_3.h b/Kinetic_surface_reconstruction/include/CGAL/Kinetic_surface_reconstruction_3.h index 05ce37d82adc..ab14c15f5b69 100644 --- a/Kinetic_surface_reconstruction/include/CGAL/Kinetic_surface_reconstruction_3.h +++ b/Kinetic_surface_reconstruction/include/CGAL/Kinetic_surface_reconstruction_3.h @@ -24,7 +24,6 @@ #include #include #include -#include #include #include #include @@ -560,8 +559,6 @@ class Kinetic_surface_reconstruction_3 { using Tds = CGAL::Triangulation_data_structure_2; using Delaunay_2 = CGAL::Delaunay_triangulation_2; - using Delaunay_3 = CGAL::Delaunay_triangulation_3; - typedef CGAL::Linear_cell_complex_traits<3, CGAL::Exact_predicates_exact_constructions_kernel> Traits; using LCC = CGAL::Linear_cell_complex_for_combinatorial_map<3, 3, Traits, typename KSP::Linear_cell_complex_min_items>; using Dart_descriptor = typename LCC::Dart_descriptor; diff --git a/Kinetic_surface_reconstruction/package_info/Kinetic_surface_reconstruction/dependencies b/Kinetic_surface_reconstruction/package_info/Kinetic_surface_reconstruction/dependencies index 485a6ca615ef..db4176cbbec0 100644 --- a/Kinetic_surface_reconstruction/package_info/Kinetic_surface_reconstruction/dependencies +++ b/Kinetic_surface_reconstruction/package_info/Kinetic_surface_reconstruction/dependencies @@ -46,4 +46,3 @@ Surface_mesh_segmentation TDS_2 TDS_3 Triangulation_2 -Triangulation_3 From 5ab26a97bfa1acddeb1527e8660c0d60d8d70ddf Mon Sep 17 00:00:00 2001 From: Sven Oesau Date: Mon, 21 Oct 2024 15:53:19 +0200 Subject: [PATCH 3/3] removing TDS_3 from dependencies --- .../package_info/Kinetic_surface_reconstruction/dependencies | 1 - 1 file changed, 1 deletion(-) diff --git a/Kinetic_surface_reconstruction/package_info/Kinetic_surface_reconstruction/dependencies b/Kinetic_surface_reconstruction/package_info/Kinetic_surface_reconstruction/dependencies index db4176cbbec0..a746754e7dcb 100644 --- a/Kinetic_surface_reconstruction/package_info/Kinetic_surface_reconstruction/dependencies +++ b/Kinetic_surface_reconstruction/package_info/Kinetic_surface_reconstruction/dependencies @@ -44,5 +44,4 @@ Stream_support Surface_mesh Surface_mesh_segmentation TDS_2 -TDS_3 Triangulation_2