From 6db0179b2d0f21985df20ed65b8e033f36bc866f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 12 Nov 2025 15:44:32 +0100 Subject: [PATCH 01/16] add experimental functions to compute the kernel of a triangle mesh --- .../CGAL/Polygon_mesh_processing/clip.h | 50 ++++- .../Polygon_mesh_processing/internal/kernel.h | 203 ++++++++++++++++++ .../refine_with_plane.h | 10 + .../Polygon_mesh_processing/CMakeLists.txt | 1 + .../test_mesh_kernel.cpp | 90 ++++++++ .../internal/parameters_interface.h | 1 + 6 files changed, 345 insertions(+), 10 deletions(-) create mode 100644 Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h create mode 100644 Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h index 79e449530ee..eef48d12119 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h @@ -265,6 +265,24 @@ bool close(PolygonMesh& pm, VertexPointMap vpm, typename Traits::Vector_3 plane_ #endif +template +void naive_close(PolygonMesh& pm, Visitor& visitor) +{ + using halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + + std::vector< halfedge_descriptor > border_cycles; + extract_boundary_cycles(pm, std::back_inserter(border_cycles)); + std::vector< Bbox_3 > bboxes; + + for (halfedge_descriptor h : border_cycles) + { + visitor.before_face_copy(boost::graph_traits::null_face(), pm, pm); + Euler::fill_hole(h, pm); + visitor.after_face_copy(boost::graph_traits::null_face(), pm, face(h, pm), pm); + } +} + template ::halfedge_descriptor; @@ -1026,11 +1044,12 @@ bool clip(PolygonMesh& pm, const bool allow_self_intersections = choose_parameter(get_parameter(np, internal_np::allow_self_intersections), false); bool triangulate = !choose_parameter(get_parameter(np, internal_np::do_not_triangulate_faces), false); - + constexpr bool traits_supports_cdt2 = !internal::Has_member_Does_not_support_CDT2::value; + const bool used_for_kernel = choose_parameter(get_parameter(np, internal_np::used_for_kernel), false); auto vos = get(dynamic_vertex_property_t(), pm); auto ecm = get(dynamic_edge_property_t(), pm, false); - if (triangulate && !is_triangle_mesh(pm)) + if (traits_supports_cdt2 && triangulate && !is_triangle_mesh(pm)) triangulate = false; refine_with_plane(pm, plane, parameters::vertex_oriented_side_map(vos) @@ -1082,15 +1101,26 @@ bool clip(PolygonMesh& pm, remove_connected_components(pm, ccs_to_remove, fcc); - if (clip_volume) + if constexpr (traits_supports_cdt2) { - //TODO: add in the traits construct_orthogonal_vector - if (triangulate) - internal::close_and_triangulate(pm, vpm, plane.orthogonal_vector(), visitor); - else - if (!internal::close(pm, vpm, plane.orthogonal_vector(), visitor)) - internal::close_and_triangulate(pm, vpm, plane.orthogonal_vector(), visitor); + if (clip_volume) + { + if (!used_for_kernel) + { + //TODO: add in the traits construct_orthogonal_vector + if (triangulate) + internal::close_and_triangulate(pm, vpm, plane.orthogonal_vector(), visitor); + else + if (!internal::close(pm, vpm, plane.orthogonal_vector(), visitor)) + internal::close_and_triangulate(pm, vpm, plane.orthogonal_vector(), visitor); + } + else + internal::naive_close(pm, visitor); + } } + else + if (clip_volume) + internal::naive_close(pm, visitor); return true; } diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h new file mode 100644 index 00000000000..10c51b99ca9 --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h @@ -0,0 +1,203 @@ +// Copyright (c) 2016-2025 GeometryFactory (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) : Sebastien Loriot + +#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_KERNEL_H +#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_KERNEL_H + +#include +#include +#include +#include + + +namespace CGAL { +namespace Polygon_mesh_processing { +namespace experimental { + +template +struct Three_point_cut_plane_traits +{ + using FT = typename Kernel::FT; + using Plane_3 = std::array; + using Point_3 = typename Kernel::Point_3; + + struct Does_not_support_CDT2{}; + + struct Oriented_side_3 + { + Oriented_side operator()(const Plane_3& plane, const Point_3& p) const + { + return orientation(plane[0], plane[1], plane[2], p); + } + }; + + struct Construct_plane_line_intersection_point_3 + { + Point_3 operator()(const Plane_3& plane, const Point_3& p, const Point_3& q) + { + typename Kernel::Construct_plane_line_intersection_point_3 construction; + typename Kernel::Plane_3 k_plane(plane[0], plane[1], plane[2]); + + return construction(k_plane, p, q); + } + }; + + Oriented_side_3 oriented_side_3_object() const + { + return Oriented_side_3(); + } + + Construct_plane_line_intersection_point_3 construct_plane_line_intersection_point_3_object() const + { + return Construct_plane_line_intersection_point_3(); + } + +#ifndef CGAL_PLANE_CLIP_DO_NOT_USE_BOX_INTERSECTION_D +// for does self-intersect + using Segment_3 = typename Kernel::Segment_3; + using Triangle_3 = typename Kernel::Triangle_3; + using Construct_segment_3 = typename Kernel::Construct_segment_3; + using Construct_triangle_3 =typename Kernel::Construct_triangle_3; + using Do_intersect_3 = typename Kernel::Do_intersect_3; + Construct_segment_3 construct_segment_3_object() const { return Construct_segment_3(); } + Construct_triangle_3 construct_triangle_3_object() const { return Construct_triangle_3(); } + Do_intersect_3 do_intersect_3_object() const { return Do_intersect_3(); } +#endif +}; + + +template +TriangleMesh +kernel(const TriangleMesh& pm, + const NamedParameters& np = parameters::default_values()) +{ + // TODO: bench if with EPECK we can directly use Kernel::Plane_3 + // TODO: TriangleMesh as output is not correct since it is actually a PolygonMesh + using parameters::choose_parameter; + using parameters::get_parameter; + + using GT = typename GetGeomTraits::type; + auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_const_property_map(vertex_point, pm)); + // GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits)); + using Point_3 = typename GT::Point_3; + using parameters::choose_parameter; + using parameters::get_parameter; + + + + CGAL::Bbox_3 bb3 = bbox(pm, np); + TriangleMesh kernel; + CGAL::make_hexahedron(Point_3(bb3.xmax(),bb3.ymin(),bb3.zmin()), Point_3(bb3.xmax(),bb3.ymax(),bb3.zmin()), Point_3(bb3.xmin(),bb3.ymax(),bb3.zmin()), Point_3(bb3.xmin(),bb3.ymin(),bb3.zmin()), + Point_3(bb3.xmin(),bb3.ymin(),bb3.zmax()), Point_3(bb3.xmax(),bb3.ymin(),bb3.zmax()), Point_3(bb3.xmax(),bb3.ymax(),bb3.zmax()), Point_3(bb3.xmin(),bb3.ymax(),bb3.zmax()), + kernel); + + Three_point_cut_plane_traits kgt; + for (auto f : faces(pm)) + { + auto h = halfedge(f, pm); + auto plane = make_array( get(vpm,source(h, pm)), + get(vpm,target(h, pm)), + get(vpm,target(next(h, pm), pm)) ); + clip(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(kgt).do_not_triangulate_faces(true).used_for_kernel(true)); + if (is_empty(kernel)) break; + } + + return kernel; +}; + +template +TriangleMesh +kernel_using_chull(const TriangleMesh& pm, + const NamedParameters& np = parameters::default_values()) +{ + using parameters::choose_parameter; + using parameters::get_parameter; + + using GT = typename GetGeomTraits::type; + auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_const_property_map(vertex_point, pm)); + // GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits)); + using parameters::choose_parameter; + using parameters::get_parameter; + + std::vector planes; + planes.reserve(faces(pm).size()); + + + for (auto f : faces(pm)) + { + auto h = halfedge(f, pm); + planes.emplace_back( get(vpm,source(h, pm)), + get(vpm,target(h, pm)), + get(vpm,target(next(h, pm), pm)) ); + } + + TriangleMesh kernel; + + // if no point inside the intersection is provided, one + // will be automatically found using linear programming + CGAL::halfspace_intersection_3(planes.begin(), + planes.end(), + kernel ); + + return kernel; +} + +template +TriangleMesh +kernel_using_chull_and_constructions(const TriangleMesh& pm, + const NamedParameters& np = parameters::default_values()) +{ + using parameters::choose_parameter; + using parameters::get_parameter; + + using GT = typename GetGeomTraits::type; + auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_const_property_map(vertex_point, pm)); + // GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits)); + using parameters::choose_parameter; + using parameters::get_parameter; + + std::vector planes; + planes.reserve(faces(pm).size()); + + + for (auto f : faces(pm)) + { + auto h = halfedge(f, pm); + planes.emplace_back( get(vpm,source(h, pm)), + get(vpm,target(h, pm)), + get(vpm,target(next(h, pm), pm)) ); + } + + TriangleMesh kernel; + + // if no point inside the intersection is provided, one + // will be automatically found using linear programming + CGAL::halfspace_intersection_with_constructions_3(planes.begin(), + planes.end(), + kernel ); + + return kernel; +} + + + + + +} } } // end of CGAL::Polygon_mesh_processing::experimental + +#endif // CGAL_POLYGON_MESH_PROCESSING_INTERNAL_KERNEL_H \ No newline at end of file diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine_with_plane.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine_with_plane.h index 03101901e83..bfce4ff7714 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine_with_plane.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine_with_plane.h @@ -24,10 +24,18 @@ #ifndef CGAL_PLANE_CLIP_DO_NOT_USE_BOX_INTERSECTION_D #include #endif +#include namespace CGAL { namespace Polygon_mesh_processing { +namespace internal +{ +BOOST_MPL_HAS_XXX_TRAIT_NAMED_DEF(Has_member_Does_not_support_CDT2, + Does_not_support_CDT2, + false) +} // internal namespace + #ifndef DOXYGEN_RUNNING template struct Default_cut_visitor @@ -74,6 +82,8 @@ struct Orthogonal_cut_plane_traits using Plane_3 = std::pair; using Point_3 = typename Kernel::Point_3; + struct Does_not_support_CDT2{}; + struct Oriented_side_3 { Oriented_side operator()(const Plane_3& plane, const Point_3& p) const diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 6fcadef5628..2950d6d0509 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -74,6 +74,7 @@ create_single_source_cgal_program("test_corefinement_nm_bo.cpp") create_single_source_cgal_program("test_corefinement_cavities.cpp") create_single_source_cgal_program("issue_8730.cpp") create_single_source_cgal_program("issue_7164.cpp") +create_single_source_cgal_program("test_mesh_kernel.cpp") # create_single_source_cgal_program("test_pmp_repair_self_intersections.cpp") find_package(Eigen3 QUIET) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp new file mode 100644 index 00000000000..f5fef238e45 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp @@ -0,0 +1,90 @@ +#include +#include +#include +#include +#include + +#include + +using K=CGAL::Exact_predicates_inexact_constructions_kernel; +using EK=CGAL::Exact_predicates_exact_constructions_kernel; +using Mesh=CGAL::Surface_mesh; +using EMesh=CGAL::Surface_mesh; + +namespace PMP=CGAL::Polygon_mesh_processing; + +void test_traits() +{ + Mesh m; + std::ifstream(CGAL::data_file_path("meshes/elephant.off")) >> m; + + std::pair p(0,0.05); + using Traits = PMP::Orthogonal_cut_plane_traits; + + PMP::clip(m, p, CGAL::parameters::geom_traits(Traits())); + + std::ofstream("clipped.off") << m; +} + +void test_kernel() +{ + CGAL::Real_timer timer; + Mesh m; + std::ifstream(CGAL::data_file_path("meshes/blobby.off")) >> m; + + timer.start(); + Mesh kernel = PMP::experimental::kernel(m); + timer.stop(); + + std::ofstream("kernel.off") << kernel; + std::cout << "test_kernel done in " << timer.time() << "\n"; +} + +void test_exact_kernel() +{ + CGAL::Real_timer timer; + EMesh m; + std::ifstream(CGAL::data_file_path("meshes/blobby.off")) >> m; + timer.start(); + EMesh kernel = PMP::experimental::kernel(m); + timer.stop(); + + std::ofstream("ekernel.off") << kernel; + std::cout << "test_exact_kernel done in " << timer.time() << "\n"; +} + +void test_kernel_with_chull() +{ + CGAL::Real_timer timer; + Mesh m; + std::ifstream(CGAL::data_file_path("meshes/blobby.off")) >> m; + + timer.start(); + Mesh kernel = PMP::experimental::kernel_using_chull(m); + timer.stop(); + + std::ofstream("kernel_with_chull.off") << kernel; + std::cout << "test_kernel_with_chull done in " << timer.time() << "\n"; +} +void test_kernel_with_chull_and_constructions() +{ + CGAL::Real_timer timer; + Mesh m; + std::ifstream(CGAL::data_file_path("meshes/blobby.off")) >> m; + + timer.start(); + Mesh kernel = PMP::experimental::kernel_using_chull_and_constructions(m); + timer.stop(); + + std::ofstream("kernel_with_chull_and_constructions.off") << kernel; + std::cout << "test_kernel_with_chull_and_constructions done in " << timer.time() << "\n"; +} + +int main() +{ + test_traits(); + test_kernel(); + test_exact_kernel(); + test_kernel_with_chull(); + test_kernel_with_chull_and_constructions(); +} 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 73c54858bf2..6f4164d8d56 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -200,6 +200,7 @@ CGAL_add_named_parameter(force_filtering_t, force_filtering, force_filtering) //internal CGAL_add_named_parameter(weight_calculator_t, weight_calculator, weight_calculator) CGAL_add_named_parameter(use_bool_op_to_clip_surface_t, use_bool_op_to_clip_surface, use_bool_op_to_clip_surface) +CGAL_add_named_parameter(used_for_kernel_t, used_for_kernel, used_for_kernel) // List of named parameters used in the Point Set Processing package CGAL_add_named_parameter(query_point_t, query_point_map, query_point_map) From 8467e71f55ca9aa6e5985d3de05c9ce76094ded9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 12 Nov 2025 16:18:44 +0100 Subject: [PATCH 02/16] add an early exit option --- .../Polygon_mesh_processing/internal/kernel.h | 39 +++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h index 10c51b99ca9..04a0850cb0f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h @@ -109,6 +109,45 @@ kernel(const TriangleMesh& pm, auto plane = make_array( get(vpm,source(h, pm)), get(vpm,target(h, pm)), get(vpm,target(next(h, pm), pm)) ); + +#ifdef CGAL_USE_OPTI_WITH_BBOX + auto pred = kgt.oriented_side_3_object(); + auto gbox = bbox(kernel); + std::array corners = CGAL::make_array(Point_3(bb3.xmax(),bb3.ymin(),bb3.zmin()), + Point_3(bb3.xmax(),bb3.ymax(),bb3.zmin()), + Point_3(bb3.xmin(),bb3.ymax(),bb3.zmin()), + Point_3(bb3.xmin(),bb3.ymin(),bb3.zmin()), + Point_3(bb3.xmin(),bb3.ymin(),bb3.zmax()), + Point_3(bb3.xmax(),bb3.ymin(),bb3.zmax()), + Point_3(bb3.xmax(),bb3.ymax(),bb3.zmax()), + Point_3(bb3.xmin(),bb3.ymax(),bb3.zmax())); + int i=0; + auto first_ori=pred(plane, corners[i]); + while(++i!=8 && first_ori==ON_ORIENTED_BOUNDARY) + first_ori=pred(plane, corners[i]); + + if (i==8) continue; + bool all_the_same=true; + for (;i<8;++i) + { + auto other_ori=pred(plane, corners[i]); + if (other_ori!=ON_ORIENTED_BOUNDARY && other_ori!=first_ori) + { + all_the_same=false; + break; + } + } + + if (all_the_same) + { + if (first_ori==ON_NEGATIVE_SIDE) continue; + else + { + return TriangleMesh(); + } + } +#endif + clip(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(kgt).do_not_triangulate_faces(true).used_for_kernel(true)); if (is_empty(kernel)) break; } From 7df4982a194f6d49ffd9d430436b9d11217c0e43 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 13 Nov 2025 12:21:43 +0100 Subject: [PATCH 03/16] kernel is empty if not a topological sphere --- .../include/CGAL/Polygon_mesh_processing/internal/kernel.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h index 04a0850cb0f..1164e8e93b9 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h @@ -94,6 +94,12 @@ kernel(const TriangleMesh& pm, using parameters::choose_parameter; using parameters::get_parameter; + //TODO: what do we do with a mesh that is not closed? + //TODO: what do we do if the input is not a triangle mesh? + + if (vertices(pm).size() - edges(pm).size() + faces(pm).size() != 2) + return TriangleMesh(); + CGAL::Bbox_3 bb3 = bbox(pm, np); From 02ca6fe09bde7dcac8a5a31e53f1b0e38f330f44 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 14 Nov 2025 10:06:44 +0100 Subject: [PATCH 04/16] handle case of empty kernel --- .../Polygon_mesh_processing/internal/kernel.h | 24 +++++++++++-------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h index 1164e8e93b9..c3ccedc204a 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h @@ -189,13 +189,15 @@ kernel_using_chull(const TriangleMesh& pm, get(vpm,target(next(h, pm), pm)) ); } + auto origin_opt = CGAL::halfspace_intersection_interior_point_3(planes.begin(), planes.end()); + TriangleMesh kernel; - // if no point inside the intersection is provided, one - // will be automatically found using linear programming - CGAL::halfspace_intersection_3(planes.begin(), - planes.end(), - kernel ); + if (origin_opt.has_value()) + CGAL::halfspace_intersection_3(planes.begin(), + planes.end(), + kernel, + origin_opt); return kernel; } @@ -228,13 +230,15 @@ kernel_using_chull_and_constructions(const TriangleMesh& pm, get(vpm,target(next(h, pm), pm)) ); } + auto origin_opt = CGAL::halfspace_intersection_interior_point_3(planes.begin(), planes.end()); + TriangleMesh kernel; - // if no point inside the intersection is provided, one - // will be automatically found using linear programming - CGAL::halfspace_intersection_with_constructions_3(planes.begin(), - planes.end(), - kernel ); + if (origin_opt.has_value()) + CGAL::halfspace_intersection_with_constructions_3(planes.begin(), + planes.end(), + kernel, + origin_opt); return kernel; } From 5954d07a8b08b514a9df0cd5d57800a11d9a5a16 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 14 Nov 2025 10:09:29 +0100 Subject: [PATCH 05/16] update the test to accept a filename as input --- .../test_mesh_kernel.cpp | 27 ++++++++++--------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp index f5fef238e45..1a6e6cdebad 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp @@ -26,11 +26,11 @@ void test_traits() std::ofstream("clipped.off") << m; } -void test_kernel() +void test_kernel(std::string fname) { CGAL::Real_timer timer; Mesh m; - std::ifstream(CGAL::data_file_path("meshes/blobby.off")) >> m; + std::ifstream(fname) >> m; timer.start(); Mesh kernel = PMP::experimental::kernel(m); @@ -40,11 +40,11 @@ void test_kernel() std::cout << "test_kernel done in " << timer.time() << "\n"; } -void test_exact_kernel() +void test_exact_kernel(std::string fname) { CGAL::Real_timer timer; EMesh m; - std::ifstream(CGAL::data_file_path("meshes/blobby.off")) >> m; + std::ifstream(fname) >> m; timer.start(); EMesh kernel = PMP::experimental::kernel(m); timer.stop(); @@ -53,11 +53,11 @@ void test_exact_kernel() std::cout << "test_exact_kernel done in " << timer.time() << "\n"; } -void test_kernel_with_chull() +void test_kernel_with_chull(std::string fname) { CGAL::Real_timer timer; Mesh m; - std::ifstream(CGAL::data_file_path("meshes/blobby.off")) >> m; + std::ifstream(fname) >> m; timer.start(); Mesh kernel = PMP::experimental::kernel_using_chull(m); @@ -66,11 +66,11 @@ void test_kernel_with_chull() std::ofstream("kernel_with_chull.off") << kernel; std::cout << "test_kernel_with_chull done in " << timer.time() << "\n"; } -void test_kernel_with_chull_and_constructions() +void test_kernel_with_chull_and_constructions(std::string fname) { CGAL::Real_timer timer; Mesh m; - std::ifstream(CGAL::data_file_path("meshes/blobby.off")) >> m; + std::ifstream(fname) >> m; timer.start(); Mesh kernel = PMP::experimental::kernel_using_chull_and_constructions(m); @@ -80,11 +80,12 @@ void test_kernel_with_chull_and_constructions() std::cout << "test_kernel_with_chull_and_constructions done in " << timer.time() << "\n"; } -int main() +int main(int argc, char** argv) { test_traits(); - test_kernel(); - test_exact_kernel(); - test_kernel_with_chull(); - test_kernel_with_chull_and_constructions(); + const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"); + test_kernel(filename); + test_exact_kernel(filename); + test_kernel_with_chull(filename); + test_kernel_with_chull_and_constructions(filename); } From fd54a4107fc4abff317d4e3ead37fafff0a9ba09 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 14 Nov 2025 10:22:11 +0100 Subject: [PATCH 06/16] handle all file types --- .../test_mesh_kernel.cpp | 24 +++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp index 1a6e6cdebad..0aedae81a73 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp @@ -30,7 +30,11 @@ void test_kernel(std::string fname) { CGAL::Real_timer timer; Mesh m; - std::ifstream(fname) >> m; + if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m)) + { + std::cerr << "ERROR: cannot read " << fname << "\n"; + exit(1); + } timer.start(); Mesh kernel = PMP::experimental::kernel(m); @@ -44,7 +48,11 @@ void test_exact_kernel(std::string fname) { CGAL::Real_timer timer; EMesh m; - std::ifstream(fname) >> m; + if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m)) + { + std::cerr << "ERROR: cannot read " << fname << "\n"; + exit(1); + } timer.start(); EMesh kernel = PMP::experimental::kernel(m); timer.stop(); @@ -57,7 +65,11 @@ void test_kernel_with_chull(std::string fname) { CGAL::Real_timer timer; Mesh m; - std::ifstream(fname) >> m; + if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m)) + { + std::cerr << "ERROR: cannot read " << fname << "\n"; + exit(1); + } timer.start(); Mesh kernel = PMP::experimental::kernel_using_chull(m); @@ -70,7 +82,11 @@ void test_kernel_with_chull_and_constructions(std::string fname) { CGAL::Real_timer timer; Mesh m; - std::ifstream(fname) >> m; + if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m)) + { + std::cerr << "ERROR: cannot read " << fname << "\n"; + exit(1); + } timer.start(); Mesh kernel = PMP::experimental::kernel_using_chull_and_constructions(m); From 3b96bc9dc1baf93ee57d2ff84c427aef9afb412a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Mon, 17 Nov 2025 10:53:25 +0100 Subject: [PATCH 07/16] do not construct plane --- .../include/CGAL/Polygon_mesh_processing/internal/kernel.h | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h index c3ccedc204a..e060a6cbbc4 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h @@ -45,9 +45,7 @@ struct Three_point_cut_plane_traits Point_3 operator()(const Plane_3& plane, const Point_3& p, const Point_3& q) { typename Kernel::Construct_plane_line_intersection_point_3 construction; - typename Kernel::Plane_3 k_plane(plane[0], plane[1], plane[2]); - - return construction(k_plane, p, q); + return construction(plane[0], plane[1], plane[2], p, q); } }; From 15095827f3d893dc09a86f11b9de118bee596c49 Mon Sep 17 00:00:00 2001 From: lvalque Date: Wed, 19 Nov 2025 18:04:02 +0100 Subject: [PATCH 08/16] Add test without rounding before computing the kernel --- .../test_mesh_kernel.cpp | 87 +++++++++++++++++-- 1 file changed, 82 insertions(+), 5 deletions(-) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp index 0aedae81a73..20dd27b399d 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp @@ -13,6 +13,45 @@ using EMesh=CGAL::Surface_mesh; namespace PMP=CGAL::Polygon_mesh_processing; +template +void round_mesh(Mesh &m){ + auto exp = [](const double v) + { + int n; + frexp(v, &n); + return n; + }; + auto pow_2 = [](const int k) + { + return (k>=0)?std::pow(2,k):1./std::pow(2,-k); + }; + + CGAL::Bbox_3 bb; + for(auto vh : m.vertices()){ + bb += m.point(vh).bbox(); + } + size_t grid_size=26; + std::array max_abs{(std::max)(-bb.xmin(), bb.xmax()), + (std::max)(-bb.ymin(), bb.ymax()), + (std::max)(-bb.zmin(), bb.zmax())}; + // Compute scale so that the exponent of max absolute value are 52-1. + std::array scale{pow_2(grid_size - exp(max_abs[0]) - 1), + pow_2(grid_size - exp(max_abs[1]) - 1), + pow_2(grid_size - exp(max_abs[2]) - 1)}; + + // auto round=[&](Point_3 &p){ + auto round=[&](typename Mesh::Point &p){ + return typename Mesh::Point(std::ceil((CGAL::to_double(p.x()) * scale[0]) - 0.5) / scale[0], + std::ceil((CGAL::to_double(p.y()) * scale[1]) - 0.5) / scale[1], + std::ceil((CGAL::to_double(p.z()) * scale[2]) - 0.5) / scale[2]); + }; + + for(auto vh : m.vertices()){ + auto &p = m.point(vh); + p = round(p); + } +} + void test_traits() { Mesh m; @@ -44,6 +83,24 @@ void test_kernel(std::string fname) std::cout << "test_kernel done in " << timer.time() << "\n"; } +void test_kernel_with_rounding(std::string fname) +{ + CGAL::Real_timer timer; + Mesh m; + if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m)) + { + std::cerr << "ERROR: cannot read " << fname << "\n"; + exit(1); + } + round_mesh(m); + timer.start(); + Mesh kernel = PMP::experimental::kernel(m); + timer.stop(); + + std::ofstream("kernel.off") << kernel; + std::cout << "test_kernel done in " << timer.time() << "\n"; +} + void test_exact_kernel(std::string fname) { CGAL::Real_timer timer; @@ -61,6 +118,24 @@ void test_exact_kernel(std::string fname) std::cout << "test_exact_kernel done in " << timer.time() << "\n"; } +void test_exact_kernel_with_rounding(std::string fname) +{ + CGAL::Real_timer timer; + EMesh m; + if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m)) + { + std::cerr << "ERROR: cannot read " << fname << "\n"; + exit(1); + } + round_mesh(m); + timer.start(); + EMesh kernel = PMP::experimental::kernel(m); + timer.stop(); + + std::ofstream("ekernel.off") << kernel; + std::cout << "test_exact_kernel done in " << timer.time() << "\n"; +} + void test_kernel_with_chull(std::string fname) { CGAL::Real_timer timer; @@ -98,10 +173,12 @@ void test_kernel_with_chull_and_constructions(std::string fname) int main(int argc, char** argv) { - test_traits(); + // test_traits(); const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"); - test_kernel(filename); - test_exact_kernel(filename); - test_kernel_with_chull(filename); - test_kernel_with_chull_and_constructions(filename); + // test_kernel(filename); + // test_kernel_with_rounding(filename); + // test_exact_kernel(filename); + test_exact_kernel_with_rounding(filename); + // test_kernel_with_chull(filename); + // test_kernel_with_chull_and_constructions(filename); } From 82adcbcc8e10c960fd0682cbea72fb5336986219 Mon Sep 17 00:00:00 2001 From: lvalque Date: Wed, 19 Nov 2025 18:05:26 +0100 Subject: [PATCH 09/16] Modify 3 planes intersections to used Ring type instead of field type --- .../Plane_3_Plane_3_Plane_3_intersection.h | 66 +++++++++++-------- 1 file changed, 38 insertions(+), 28 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Plane_3_Plane_3_Plane_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Plane_3_Plane_3_Plane_3_intersection.h index 138d2c9bd45..4b66d399837 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Plane_3_Plane_3_Plane_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Plane_3_Plane_3_Plane_3_intersection.h @@ -28,49 +28,59 @@ namespace internal { // triple plane intersection template std::optional -intersection_point(const typename K::Plane_3& plane1, - const typename K::Plane_3& plane2, - const typename K::Plane_3& plane3, - const K&) +intersection_point_RT(const typename K::Plane_3& plane1, + const typename K::Plane_3& plane2, + const typename K::Plane_3& plane3, + const K&) { - typedef typename K::FT FT; + typedef typename K::RT RT; - const FT& m00 = plane1.a(); - const FT& m01 = plane1.b(); - const FT& m02 = plane1.c(); - const FT& b0 = - plane1.d(); - const FT& m10 = plane2.a(); - const FT& m11 = plane2.b(); - const FT& m12 = plane2.c(); - const FT& b1 = - plane2.d(); - const FT& m20 = plane3.a(); - const FT& m21 = plane3.b(); - const FT& m22 = plane3.c(); - const FT& b2 = - plane3.d(); + const RT& m00 = plane1.a(); + const RT& m01 = plane1.b(); + const RT& m02 = plane1.c(); + const RT& b0 = - plane1.d(); + const RT& m10 = plane2.a(); + const RT& m11 = plane2.b(); + const RT& m12 = plane2.c(); + const RT& b1 = - plane2.d(); + const RT& m20 = plane3.a(); + const RT& m21 = plane3.b(); + const RT& m22 = plane3.c(); + const RT& b2 = - plane3.d(); // Minors common to two determinants - const FT minor_0 = m00*m11 - m10*m01; - const FT minor_1 = m00*m21 - m20*m01; - const FT minor_2 = m10*m21 - m20*m11; + const RT minor_0 = m00*m11 - m10*m01; + const RT minor_1 = m00*m21 - m20*m01; + const RT minor_2 = m10*m21 - m20*m11; - const FT den = minor_0*m22 - minor_1*m12 + minor_2*m02; // determinant of M + const RT den = minor_0*m22 - minor_1*m12 + minor_2*m02; // determinant of M if(is_zero(den)){ return std::nullopt; } - const FT num3 = minor_0*b2 - minor_1*b1 + minor_2*b0; // determinant of M with M[x:2] swapped with [b0,b1,b2] + const RT num3 = minor_0*b2 - minor_1*b1 + minor_2*b0; // determinant of M with M[x:2] swapped with [b0,b1,b2] // Minors common to two determinants - const FT minor_3 = b0*m12 - b1*m02; - const FT minor_4 = b0*m22 - b2*m02; - const FT minor_5 = b1*m22 - b2*m12; + const RT minor_3 = b0*m12 - b1*m02; + const RT minor_4 = b0*m22 - b2*m02; + const RT minor_5 = b1*m22 - b2*m12; // num1 has opposite signs because b0 and M[:1] have been swapped - const FT num1 = - minor_3*m21 + minor_4*m11 - minor_5*m01; // determinant of M with M[x:0] swapped with [b0,b1,b2] - const FT num2 = minor_3*m20 - minor_4*m10 + minor_5*m00; // determinant of M with M[x:1] swapped with [b0,b1,b2] + const RT num1 = - minor_3*m21 + minor_4*m11 - minor_5*m01; // determinant of M with M[x:0] swapped with [b0,b1,b2] + const RT num2 = minor_3*m20 - minor_4*m10 + minor_5*m00; // determinant of M with M[x:1] swapped with [b0,b1,b2] - return std::make_optional(typename K::Point_3(num1/den, num2/den, num3/den)); + return typename K::Point_3(num1, num2, num3, den); +} + +template +std::optional +intersection_point(const typename K::Plane_3& plane1, + const typename K::Plane_3& plane2, + const typename K::Plane_3& plane3, + const K& k) +{ + return intersection_point_RT(plane1, plane2, plane3, k); } template From a06e595e2f80b50a7104201435e38b45fae844f3 Mon Sep 17 00:00:00 2001 From: lvalque Date: Fri, 21 Nov 2025 18:18:44 +0100 Subject: [PATCH 10/16] Working but not optimized trettner inspired kernel algorithm --- .../include/CGAL/Kernel/function_objects.h | 21 +- .../Polygon_mesh_processing/internal/kernel.h | 305 +++++++++++++++++- .../test_mesh_kernel.cpp | 60 +++- 3 files changed, 380 insertions(+), 6 deletions(-) diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index b71470a6262..c2df37bf337 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/Kernel_23/include/CGAL/Kernel/function_objects.h @@ -2097,18 +2097,24 @@ namespace CommonKernelFunctors { template class Construct_plane_line_intersection_point_3 { - typedef typename K::Plane_3 Plane; + typedef typename K::FT FT; typedef typename K::Line_3 Line; typedef typename K::Point_3 Point; - typename K::Construct_plane_3 construct_plane; - typename K::Construct_line_3 construct_line; + typedef typename K::Plane_3 Plane; + typedef typename K::Vector_3 Vector; + typename K::Construct_cross_product_vector_3 construct_cross_product_vector_3; + typename K::Compute_scalar_product_3 scalar_product; + + typename K::Construct_line_3 construct_line; + typename K::Construct_vector_3 vector; + typename K::Construct_plane_3 construct_plane; public: Point operator()(const Point& p1, const Point& p2, const Point& p3, const Point& l1, const Point& l2) const { - Plane plane = construct_plane(p1, p2, p3); + Plane plane = construct_plane( p1, p2, p3 ); Line line = construct_line( l1, l2 ); const auto res = typename K::Intersect_3()(plane,line); @@ -2116,6 +2122,13 @@ namespace CommonKernelFunctors { const Point* e_pt = std::get_if(&(*res)); CGAL_assertion(e_pt!=nullptr); return *e_pt; + + // Vector n = construct_cross_product_vector_3( vector(p1, p2), vector(p1, p3)); + // CGAL_assertion(n!=CGAL::NULL_VECTOR); + // Vector l = vector(l1, l2); + // CGAL_assertion(l!=CGAL::NULL_VECTOR); + // FT t = scalar_product(n, vector(p1,l1)) / scalar_product(n, l); + // return l1 + t * l; } Point diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h index e060a6cbbc4..b85c52b5507 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h @@ -18,6 +18,10 @@ #include #include +#include +#include +#include +#include namespace CGAL { namespace Polygon_mesh_processing { @@ -107,16 +111,47 @@ kernel(const TriangleMesh& pm, kernel); Three_point_cut_plane_traits kgt; +#ifdef CGAL_USE_CONCAVE_FACE_FIRST + auto compute_dihedral_angle=[](const TriangleMesh& mesh, const halfedge_descriptor h) { + // Calcul de l'angle dihédral entre les deux faces à partir des quatre points + return to_double(approximate_dihedral_angle(mesh.point(mesh.source(h)), mesh.point(mesh.target(h)), + mesh.point(mesh.target(mesh.next(h))), mesh.point(mesh.target(mesh.next(mesh.opposite(h)))))); + }; + + std::vector> faces_with_angles; + for(face_descriptor f : faces(pm)){ + double max_dihedral_angle = 0.0; + auto h = pm.halfedge(f); + auto hf_circ = pm.halfedges_around_face(h); + for (auto he : hf_circ) { + if (!pm.is_border(he)) { + max_dihedral_angle = (std::max)(max_dihedral_angle,compute_dihedral_angle(pm, he)); + } + faces_with_angles.push_back({f, max_dihedral_angle}); + } + } + std::sort(faces_with_angles.begin(), faces_with_angles.end(), + [](const std::pair& a, const std::pair& b) { + return a.second > b.second; + }); + for(auto pair : faces_with_angles) + { + auto h = halfedge(pair.first, pm); + auto plane = make_array( get(vpm,source(h, pm)), + get(vpm,target(h, pm)), + get(vpm,target(next(h, pm), pm)) ); +#else for (auto f : faces(pm)) { auto h = halfedge(f, pm); auto plane = make_array( get(vpm,source(h, pm)), get(vpm,target(h, pm)), get(vpm,target(next(h, pm), pm)) ); +#endif #ifdef CGAL_USE_OPTI_WITH_BBOX auto pred = kgt.oriented_side_3_object(); - auto gbox = bbox(kernel); + bb3 = bbox(kernel); std::array corners = CGAL::make_array(Point_3(bb3.xmax(),bb3.ymin(),bb3.zmin()), Point_3(bb3.xmax(),bb3.ymax(),bb3.zmin()), Point_3(bb3.xmin(),bb3.ymax(),bb3.zmin()), @@ -243,6 +278,274 @@ kernel_using_chull_and_constructions(const TriangleMesh& pm, +//__________________________________________________________________________________________________________________________________ + +template +struct Plane_based_traits +{ + using FT = typename Kernel::FT; + using RT = typename Kernel::RT; + using plane_descriptor = std::size_t; + using Plane_3 = std::pair; + + std::vector *m_planes; + + using Input_point_3 = typename Kernel::Point_3; + using Geometric_point_3 = typename Kernel::Point_3; + + struct Point_3 : public Geometric_point_3{ + std::array supports; + mutable std::set other_coplanar; + using Base = Geometric_point_3; + + Point_3(){} + Point_3(plane_descriptor a, plane_descriptor b, plane_descriptor c, const std::vector &planes): + Base(*CGAL::Intersections::internal::intersection_point(planes[a].first, planes[b].first, planes[c].first, Kernel())), + supports({a,b,c}){ + } + + Point_3(plane_descriptor a, plane_descriptor b, plane_descriptor c, Geometric_point_3 p): Base(p), supports({a,b,c}){} + }; + + Point_3 construct_point_3(plane_descriptor a, plane_descriptor b, plane_descriptor c){ + return Point_3(a, b, c, *m_planes); + } + + Point_3 construct_point_3(plane_descriptor a, plane_descriptor b, plane_descriptor c, Geometric_point_3 p){ + return Point_3(a, b, c, p); + } + + struct Does_not_support_CDT2{}; + + struct Oriented_side_3 + { + Oriented_side operator()(const Plane_3& plane, const Point_3& p) const + { + if((p.supports[0]==plane.second) || (p.supports[1]==plane.second) || (p.supports[2]==plane.second)) + return COPLANAR; + Oriented_side ori = plane.first.oriented_side(p); + if(ori==COPLANAR) + p.other_coplanar.emplace(plane.second); + return ori; + } + }; + + struct Construct_plane_line_intersection_point_3 + { + std::vector* m_planes; + Construct_plane_line_intersection_point_3(std::vector* planes) : m_planes(planes){} + Point_3 operator()(const Plane_3& plane, const Point_3& p, const Point_3& q) + { + auto get_common_supports=[](const Point_3& p, const Point_3 &q){ + plane_descriptor first=-1; + plane_descriptor second=-1; + for(short int i=0; i!=3; ++i) + for(short int j=0; j!=3; ++j) + if(p.supports[i]==q.supports[j]) + if(first==-1){ + first=p.supports[i]; + break; + } else { + second=p.supports[i]; + return std::make_pair(first, second); + } + + for(plane_descriptor plane: p.other_coplanar) + for(short int j=0; j!=3; ++j) + if(plane==q.supports[j]) + if(first==-1){ + first=plane; + break; + } else { + second=plane; + return std::make_pair(first, second); + } + + for(plane_descriptor plane: q.other_coplanar) + for(short int j=0; j!=3; ++j) + if(plane==p.supports[j]) + if(first==-1){ + first=plane; + break; + } else { + second=plane; + return std::make_pair(first, second); + } + + for(plane_descriptor plane: p.other_coplanar) + for(plane_descriptor qlane: q.other_coplanar) + if(plane==qlane) + if(first==-1){ + first=plane; + break; + } else { + second=plane; + return std::make_pair(first, second); + } + // The two points do not shair a common support + std::cout << p.supports[0] << " " << p.supports[1] << " " << p.supports[2] << std::endl; + std::cout << q.supports[0] << " " << q.supports[1] << " " << q.supports[2] << std::endl; + CGAL_assertion_code(std::cout << "The two points do not share two common supporting planes" << std::endl;) + CGAL_assertion(0); + }; + + auto orient=[](const typename Kernel::Plane_3& plane, const Geometric_point_3& p){ + return sign (plane.a()*p.hx() + plane.b()*p.hy() + plane.c()*p.hz() + plane.d() * p.hw()); + }; + std::pair line_supports=get_common_supports(p, q); + Point_3 res(plane.second, line_supports.first, line_supports.second, *m_planes); + return res; + } + }; + + template + Plane_based_traits(const Mesh &m){ + m_planes = new std::vector; + m_planes->reserve(faces(m).size()); + get_planes(m, std::back_inserter(*m_planes)); + } + + template + void get_planes(const Mesh &m, OutputIterator out){ + using face_descriptor = typename boost::graph_traits::face_descriptor; + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + size_t k=0; + auto to_int=[](const typename Mesh::Point &p){ + return Geometric_point_3(int(p.x()),int(p.y()),int(p.z())); + }; + auto to_int_plane=[&](face_descriptor f){ + auto pmap=boost::make_function_property_map([&](vertex_descriptor v){ + return to_int(m.point(v)); + }); + // auto pl = compute_face_normal(f, m); + auto pl = compute_face_normal(f, m, parameters::vertex_point_map(pmap)); + return typename Kernel::Vector_3(pl.x(),pl.y(),pl.z()); + }; + for(face_descriptor f : faces(m)){ + Plane_3 pl(typename Kernel::Plane_3(to_int(m.point(m.target(m.halfedge(f)))), + to_int_plane(f)), + k); + out++=pl; + k++; + } + + } + + Oriented_side_3 oriented_side_3_object() const + { + return Oriented_side_3(); + } + + Construct_plane_line_intersection_point_3 construct_plane_line_intersection_point_3_object() + { + return Construct_plane_line_intersection_point_3(m_planes); + } + + struct Construct_plane_3{ + std::vector* m_planes; + Construct_plane_3(std::vector* planes) : m_planes(planes){} + Plane_3 operator()(const typename Kernel::Plane_3 &pl){ + m_planes->emplace_back(pl, m_planes->size()); + return m_planes->back(); + } + + Plane_3 operator()(const RT &a, const RT &b, const RT &c, const RT &d){ + return (*this)(typename Kernel::Plane_3(a,b,c,d)); + } + }; + + Construct_plane_3 construct_plane_3_object() + { + return Construct_plane_3(m_planes); + } + + const std::vector& planes(){ + return *m_planes; + } + +#ifndef CGAL_PLANE_CLIP_DO_NOT_USE_BOX_INTERSECTION_D +// for does self-intersect + using Segment_3 = typename Kernel::Segment_3; + using Triangle_3 = typename Kernel::Triangle_3; + using Construct_segment_3 = typename Kernel::Construct_segment_3; + using Construct_triangle_3 =typename Kernel::Construct_triangle_3; + using Do_intersect_3 = typename Kernel::Do_intersect_3; + Construct_segment_3 construct_segment_3_object() const { return Construct_segment_3(); } + Construct_triangle_3 construct_triangle_3_object() const { return Construct_triangle_3(); } + Do_intersect_3 do_intersect_3_object() const { return Do_intersect_3(); } +#endif +}; + +template +TriangleMesh +trettner_kernel(const TriangleMesh& pm, + const NamedParameters& np = parameters::default_values()) +{ + // TODO: bench if with EPECK we can directly use Kernel::Plane_3 + // TODO: TriangleMesh as output is not correct since it is actually a PolygonMesh + using parameters::choose_parameter; + using parameters::get_parameter; + + // using int256 = boost::multiprecision::int256_t; + using int256 = Exact_integer; + using GT = Plane_based_traits>; + auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_const_property_map(vertex_point, pm)); + + using Point_3 = typename GT::Point_3; + using Plane_3 = typename GT::Plane_3; + using parameters::choose_parameter; + using parameters::get_parameter; + + GT gt(pm); + std::vector planes=gt.planes(); + + auto plane_3 = gt.construct_plane_3_object(); + + //TODO: what do we do with a mesh that is not closed? + //TODO: what do we do if the input is not a triangle mesh? + + if (vertices(pm).size() - edges(pm).size() + faces(pm).size() != 2) + return TriangleMesh(); + + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef Surface_mesh InternMesh; + + CGAL::Bbox_3 bb3 = bbox(pm, np); + InternMesh kernel; + Plane_3 xl=plane_3(1,0,0,int(-bb3.xmin())); + Plane_3 yl=plane_3(0,1,0,int(-bb3.ymin())); + Plane_3 zl=plane_3(0,0,1,int(-bb3.zmin())); + Plane_3 xr=plane_3(1,0,0,int(-bb3.xmax())); + Plane_3 yr=plane_3(0,1,0,int(-bb3.ymax())); + Plane_3 zr=plane_3(0,0,1,int(-bb3.zmax())); + CGAL::make_hexahedron(gt.construct_point_3(xl.second, yl.second, zl.second), + gt.construct_point_3(xl.second, yl.second, zr.second), + gt.construct_point_3(xl.second, yr.second, zr.second), + gt.construct_point_3(xl.second, yr.second, zl.second), + gt.construct_point_3(xr.second, yr.second, zl.second), + gt.construct_point_3(xr.second, yl.second, zl.second), + gt.construct_point_3(xr.second, yl.second, zr.second), + gt.construct_point_3(xr.second, yr.second, zr.second), + kernel); + + for (auto plane: planes) + { + clip(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(gt).do_not_triangulate_faces(true).used_for_kernel(true)); + if (is_empty(kernel)) break; + } + + // return kernel; + return pm; +} + + + } } } // end of CGAL::Polygon_mesh_processing::experimental diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp index 20dd27b399d..95bed6e7d01 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp @@ -52,6 +52,45 @@ void round_mesh(Mesh &m){ } } +template +void to_integer_mesh(Mesh &m){ + auto exp = [](const double v) + { + int n; + frexp(v, &n); + return n; + }; + auto pow_2 = [](const int k) + { + return (k>=0)?std::pow(2,k):1./std::pow(2,-k); + }; + + CGAL::Bbox_3 bb; + for(auto vh : m.vertices()){ + bb += m.point(vh).bbox(); + } + size_t grid_size=26; + std::array max_abs{(std::max)(-bb.xmin(), bb.xmax()), + (std::max)(-bb.ymin(), bb.ymax()), + (std::max)(-bb.zmin(), bb.zmax())}; + // Compute scale so that the exponent of max absolute value are 52-1. + std::array scale{pow_2(grid_size - exp(max_abs[0]) - 1), + pow_2(grid_size - exp(max_abs[1]) - 1), + pow_2(grid_size - exp(max_abs[2]) - 1)}; + + // auto round=[&](Point_3 &p){ + auto round=[&](typename Mesh::Point &p){ + return typename Mesh::Point(std::ceil((CGAL::to_double(p.x()) * scale[0]) - 0.5), + std::ceil((CGAL::to_double(p.y()) * scale[1]) - 0.5), + std::ceil((CGAL::to_double(p.z()) * scale[2]) - 0.5)); + }; + + for(auto vh : m.vertices()){ + auto &p = m.point(vh); + p = round(p); + } +} + void test_traits() { Mesh m; @@ -171,6 +210,24 @@ void test_kernel_with_chull_and_constructions(std::string fname) std::cout << "test_kernel_with_chull_and_constructions done in " << timer.time() << "\n"; } +void test_trettner_kernel(std::string fname) +{ + CGAL::Real_timer timer; + Mesh m; + if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m)) + { + std::cerr << "ERROR: cannot read " << fname << "\n"; + exit(1); + } + to_integer_mesh(m); + timer.start(); + Mesh kernel = PMP::experimental::trettner_kernel(m); + timer.stop(); + + std::ofstream("ekernel.off") << kernel; + std::cout << "test_trettner_kernel done in " << timer.time() << "\n"; +} + int main(int argc, char** argv) { // test_traits(); @@ -178,7 +235,8 @@ int main(int argc, char** argv) // test_kernel(filename); // test_kernel_with_rounding(filename); // test_exact_kernel(filename); - test_exact_kernel_with_rounding(filename); + // test_exact_kernel_with_rounding(filename); + test_trettner_kernel(filename); // test_kernel_with_chull(filename); // test_kernel_with_chull_and_constructions(filename); } From 66bcb882d03efa77bf2401265eb3ee1805f79571 Mon Sep 17 00:00:00 2001 From: lvalque Date: Mon, 24 Nov 2025 15:00:18 +0100 Subject: [PATCH 11/16] clean kernel.h and add assertions --- .../Polygon_mesh_processing/internal/kernel.h | 95 +++++++++++-------- .../test_mesh_kernel.cpp | 9 +- 2 files changed, 58 insertions(+), 46 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h index b85c52b5507..bb9fbfb6acc 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h @@ -280,19 +280,21 @@ kernel_using_chull_and_constructions(const TriangleMesh& pm, //__________________________________________________________________________________________________________________________________ +using int256 = boost::multiprecision::int256_t; +using int512 = boost::multiprecision::int512_t; +using intexact = Exact_integer; + template struct Plane_based_traits { using FT = typename Kernel::FT; using RT = typename Kernel::RT; + using Geometric_point_3 = typename Kernel::Point_3; + using plane_descriptor = std::size_t; using Plane_3 = std::pair; - std::vector *m_planes; - using Input_point_3 = typename Kernel::Point_3; - using Geometric_point_3 = typename Kernel::Point_3; - struct Point_3 : public Geometric_point_3{ std::array supports; mutable std::set other_coplanar; @@ -301,19 +303,23 @@ struct Plane_based_traits Point_3(){} Point_3(plane_descriptor a, plane_descriptor b, plane_descriptor c, const std::vector &planes): Base(*CGAL::Intersections::internal::intersection_point(planes[a].first, planes[b].first, planes[c].first, Kernel())), - supports({a,b,c}){ + supports({a,b,c}){ } Point_3(plane_descriptor a, plane_descriptor b, plane_descriptor c, Geometric_point_3 p): Base(p), supports({a,b,c}){} }; - Point_3 construct_point_3(plane_descriptor a, plane_descriptor b, plane_descriptor c){ - return Point_3(a, b, c, *m_planes); - } + struct Construct_point_3{ + std::vector* m_planes; + Construct_point_3(std::vector* planes) : m_planes(planes){} + Point_3 operator()(plane_descriptor a, plane_descriptor b, plane_descriptor c){ + return Point_3(a, b, c, *m_planes); + } - Point_3 construct_point_3(plane_descriptor a, plane_descriptor b, plane_descriptor c, Geometric_point_3 p){ - return Point_3(a, b, c, p); - } + Point_3 operator()(plane_descriptor a, plane_descriptor b, plane_descriptor c, Geometric_point_3 p){ + return Point_3(a, b, c, p); + } + }; struct Does_not_support_CDT2{}; @@ -394,6 +400,12 @@ struct Plane_based_traits }; std::pair line_supports=get_common_supports(p, q); Point_3 res(plane.second, line_supports.first, line_supports.second, *m_planes); + + int512 max=57896044618658097711785492504343953926634992332820282019728792003956564819968; //2^256 + assert(abs(p.hx())<=max); + assert(abs(p.hy())<=max); + assert(abs(p.hz())<=max); + assert(abs(p.hw())<=max); return res; } }; @@ -417,7 +429,6 @@ struct Plane_based_traits auto pmap=boost::make_function_property_map([&](vertex_descriptor v){ return to_int(m.point(v)); }); - // auto pl = compute_face_normal(f, m); auto pl = compute_face_normal(f, m, parameters::vertex_point_map(pmap)); return typename Kernel::Vector_3(pl.x(),pl.y(),pl.z()); }; @@ -442,10 +453,16 @@ struct Plane_based_traits } struct Construct_plane_3{ + std::vector* m_planes; Construct_plane_3(std::vector* planes) : m_planes(planes){} Plane_3 operator()(const typename Kernel::Plane_3 &pl){ m_planes->emplace_back(pl, m_planes->size()); + int512 max = 1208925819614629174706176; //2^80 + assert(abs(pl.a())<=max); + assert(abs(pl.b())<=max); + assert(abs(pl.c())<=max); + assert(abs(pl.d())<=max); return m_planes->back(); } @@ -459,6 +476,11 @@ struct Plane_based_traits return Construct_plane_3(m_planes); } + Construct_point_3 construct_point_3_object() + { + return Construct_point_3(m_planes); + } + const std::vector& planes(){ return *m_planes; } @@ -478,43 +500,33 @@ struct Plane_based_traits template -TriangleMesh +Surface_mesh>::Point_3> trettner_kernel(const TriangleMesh& pm, const NamedParameters& np = parameters::default_values()) { - // TODO: bench if with EPECK we can directly use Kernel::Plane_3 - // TODO: TriangleMesh as output is not correct since it is actually a PolygonMesh using parameters::choose_parameter; using parameters::get_parameter; - // using int256 = boost::multiprecision::int256_t; - using int256 = Exact_integer; - using GT = Plane_based_traits>; + using GT = Plane_based_traits>; auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), get_const_property_map(vertex_point, pm)); using Point_3 = typename GT::Point_3; using Plane_3 = typename GT::Plane_3; - using parameters::choose_parameter; - using parameters::get_parameter; + + using Construct_plane_3 = GT::Construct_plane_3; + using Construct_point_3 = GT::Construct_point_3; + + using InternMesh = Surface_mesh; GT gt(pm); std::vector planes=gt.planes(); - auto plane_3 = gt.construct_plane_3_object(); - - //TODO: what do we do with a mesh that is not closed? - //TODO: what do we do if the input is not a triangle mesh? + Construct_plane_3 plane_3 = gt.construct_plane_3_object(); + Construct_point_3 point_3 = gt.construct_point_3_object(); if (vertices(pm).size() - edges(pm).size() + faces(pm).size() != 2) - return TriangleMesh(); - - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::edge_descriptor edge_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; - - typedef Surface_mesh InternMesh; + return Surface_mesh(); CGAL::Bbox_3 bb3 = bbox(pm, np); InternMesh kernel; @@ -524,14 +536,14 @@ trettner_kernel(const TriangleMesh& pm, Plane_3 xr=plane_3(1,0,0,int(-bb3.xmax())); Plane_3 yr=plane_3(0,1,0,int(-bb3.ymax())); Plane_3 zr=plane_3(0,0,1,int(-bb3.zmax())); - CGAL::make_hexahedron(gt.construct_point_3(xl.second, yl.second, zl.second), - gt.construct_point_3(xl.second, yl.second, zr.second), - gt.construct_point_3(xl.second, yr.second, zr.second), - gt.construct_point_3(xl.second, yr.second, zl.second), - gt.construct_point_3(xr.second, yr.second, zl.second), - gt.construct_point_3(xr.second, yl.second, zl.second), - gt.construct_point_3(xr.second, yl.second, zr.second), - gt.construct_point_3(xr.second, yr.second, zr.second), + CGAL::make_hexahedron(point_3(xl.second, yl.second, zl.second), + point_3(xl.second, yl.second, zr.second), + point_3(xl.second, yr.second, zr.second), + point_3(xl.second, yr.second, zl.second), + point_3(xr.second, yr.second, zl.second), + point_3(xr.second, yl.second, zl.second), + point_3(xr.second, yl.second, zr.second), + point_3(xr.second, yr.second, zr.second), kernel); for (auto plane: planes) @@ -540,8 +552,7 @@ trettner_kernel(const TriangleMesh& pm, if (is_empty(kernel)) break; } - // return kernel; - return pm; + return kernel; } diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp index 95bed6e7d01..350a8dc2679 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp @@ -166,7 +166,8 @@ void test_exact_kernel_with_rounding(std::string fname) std::cerr << "ERROR: cannot read " << fname << "\n"; exit(1); } - round_mesh(m); + // round_mesh(m); + to_integer_mesh(m); timer.start(); EMesh kernel = PMP::experimental::kernel(m); timer.stop(); @@ -221,10 +222,10 @@ void test_trettner_kernel(std::string fname) } to_integer_mesh(m); timer.start(); - Mesh kernel = PMP::experimental::trettner_kernel(m); + auto kernel = PMP::experimental::trettner_kernel(m); timer.stop(); - std::ofstream("ekernel.off") << kernel; + std::ofstream("tkernel.off") << kernel; std::cout << "test_trettner_kernel done in " << timer.time() << "\n"; } @@ -235,7 +236,7 @@ int main(int argc, char** argv) // test_kernel(filename); // test_kernel_with_rounding(filename); // test_exact_kernel(filename); - // test_exact_kernel_with_rounding(filename); + test_exact_kernel_with_rounding(filename); test_trettner_kernel(filename); // test_kernel_with_chull(filename); // test_kernel_with_chull_and_constructions(filename); From ca66e33149de2f5c98802bfc89ac1c9b01aa7fb4 Mon Sep 17 00:00:00 2001 From: lvalque Date: Mon, 24 Nov 2025 18:01:19 +0100 Subject: [PATCH 12/16] Getting homogenous coordinate instead of common coordinates and assertion on the bound of coordinates --- .../Polygon_mesh_processing/internal/kernel.h | 36 +++++++++---------- .../test_mesh_kernel.cpp | 1 - 2 files changed, 16 insertions(+), 21 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h index bb9fbfb6acc..74a59ff7f58 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h @@ -401,11 +401,9 @@ struct Plane_based_traits std::pair line_supports=get_common_supports(p, q); Point_3 res(plane.second, line_supports.first, line_supports.second, *m_planes); - int512 max=57896044618658097711785492504343953926634992332820282019728792003956564819968; //2^256 - assert(abs(p.hx())<=max); - assert(abs(p.hy())<=max); - assert(abs(p.hz())<=max); - assert(abs(p.hw())<=max); + namespace mp = boost::multiprecision; + CGAL_assertion_code(int256 max2E200=mp::pow(int256(2),200);) + CGAL_assertion((mp::abs(p.hx())<=max2E200) && (mp::abs(p.hy())<=max2E200) && (mp::abs(p.hz())<=max2E200) && (mp::abs(p.hw())<=max2E200)); return res; } }; @@ -421,7 +419,10 @@ struct Plane_based_traits void get_planes(const Mesh &m, OutputIterator out){ using face_descriptor = typename boost::graph_traits::face_descriptor; using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + + Construct_plane_3 plane_3 = construct_plane_3_object(); size_t k=0; + namespace mp = boost::multiprecision; auto to_int=[](const typename Mesh::Point &p){ return Geometric_point_3(int(p.x()),int(p.y()),int(p.z())); }; @@ -430,15 +431,11 @@ struct Plane_based_traits return to_int(m.point(v)); }); auto pl = compute_face_normal(f, m, parameters::vertex_point_map(pmap)); - return typename Kernel::Vector_3(pl.x(),pl.y(),pl.z()); + return typename Kernel::Vector_3(pl.hx(),pl.hy(),pl.hz()); }; - for(face_descriptor f : faces(m)){ - Plane_3 pl(typename Kernel::Plane_3(to_int(m.point(m.target(m.halfedge(f)))), - to_int_plane(f)), - k); - out++=pl; - k++; - } + for(face_descriptor f : faces(m)) + plane_3(typename Kernel::Plane_3(to_int(m.point(m.target(m.halfedge(f)))), + to_int_plane(f))); } @@ -457,12 +454,11 @@ struct Plane_based_traits std::vector* m_planes; Construct_plane_3(std::vector* planes) : m_planes(planes){} Plane_3 operator()(const typename Kernel::Plane_3 &pl){ + namespace mp = boost::multiprecision; m_planes->emplace_back(pl, m_planes->size()); - int512 max = 1208925819614629174706176; //2^80 - assert(abs(pl.a())<=max); - assert(abs(pl.b())<=max); - assert(abs(pl.c())<=max); - assert(abs(pl.d())<=max); + CGAL_assertion_code(int256 max2E60=mp::pow(int256(2),60);) + CGAL_assertion_code(int256 max2E90=mp::pow(int256(2),90);) + CGAL_assertion((mp::abs(pl.a())<=max2E60) && (mp::abs(pl.b())<=max2E60) && (mp::abs(pl.c())<=max2E60) && (mp::abs(pl.d())<=max2E90)); return m_planes->back(); } @@ -500,14 +496,14 @@ struct Plane_based_traits template -Surface_mesh>::Point_3> +Surface_mesh>::Point_3> trettner_kernel(const TriangleMesh& pm, const NamedParameters& np = parameters::default_values()) { using parameters::choose_parameter; using parameters::get_parameter; - using GT = Plane_based_traits>; + using GT = Plane_based_traits>; auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), get_const_property_map(vertex_point, pm)); diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp index 350a8dc2679..88c58aaf678 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp @@ -39,7 +39,6 @@ void round_mesh(Mesh &m){ pow_2(grid_size - exp(max_abs[1]) - 1), pow_2(grid_size - exp(max_abs[2]) - 1)}; - // auto round=[&](Point_3 &p){ auto round=[&](typename Mesh::Point &p){ return typename Mesh::Point(std::ceil((CGAL::to_double(p.x()) * scale[0]) - 0.5) / scale[0], std::ceil((CGAL::to_double(p.y()) * scale[1]) - 0.5) / scale[1], From 5dbcc09821beb19a4d7a73f5be688eacf3a82b05 Mon Sep 17 00:00:00 2001 From: lvalque Date: Tue, 25 Nov 2025 15:47:27 +0100 Subject: [PATCH 13/16] Clean planes storage of Plane_based_traits and refine bound check of numbers --- .../Polygon_mesh_processing/internal/kernel.h | 101 ++++++++++-------- 1 file changed, 55 insertions(+), 46 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h index 74a59ff7f58..9a33d57caa4 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h @@ -293,7 +293,13 @@ struct Plane_based_traits using plane_descriptor = std::size_t; using Plane_3 = std::pair; - std::vector *m_planes; + +private: + using plane_range_pointer = std::shared_ptr>; + plane_range_pointer m_planes; +public: + + struct Point_3 : public Geometric_point_3{ std::array supports; @@ -310,10 +316,13 @@ struct Plane_based_traits }; struct Construct_point_3{ - std::vector* m_planes; - Construct_point_3(std::vector* planes) : m_planes(planes){} + plane_range_pointer m_planes; + Construct_point_3(plane_range_pointer planes) : m_planes(planes){} Point_3 operator()(plane_descriptor a, plane_descriptor b, plane_descriptor c){ - return Point_3(a, b, c, *m_planes); + const std::vector &planes = *m_planes; + auto res = CGAL::Intersections::internal::intersection_point(planes[a].first, planes[b].first, planes[c].first, Kernel()); + CGAL_assertion(res); + return Point_3(a, b, c, *res); } Point_3 operator()(plane_descriptor a, plane_descriptor b, plane_descriptor c, Geometric_point_3 p){ @@ -338,11 +347,14 @@ struct Plane_based_traits struct Construct_plane_line_intersection_point_3 { - std::vector* m_planes; - Construct_plane_line_intersection_point_3(std::vector* planes) : m_planes(planes){} + plane_range_pointer m_planes; + Construct_plane_line_intersection_point_3(plane_range_pointer planes) : m_planes(planes){} Point_3 operator()(const Plane_3& plane, const Point_3& p, const Point_3& q) { - auto get_common_supports=[](const Point_3& p, const Point_3 &q){ + Construct_point_3 point_3(m_planes); + const std::vector &planes=*m_planes; + + auto get_common_supports=[&](const Point_3& p, const Point_3 &q){ plane_descriptor first=-1; plane_descriptor second=-1; for(short int i=0; i!=3; ++i) @@ -356,36 +368,36 @@ struct Plane_based_traits return std::make_pair(first, second); } - for(plane_descriptor plane: p.other_coplanar) + for(plane_descriptor pd: p.other_coplanar) for(short int j=0; j!=3; ++j) - if(plane==q.supports[j]) + if(pd==q.supports[j]) if(first==-1){ - first=plane; + first=pd; break; - } else { - second=plane; + } else if(planes[first].first!=planes[pd].first){ + second=pd; return std::make_pair(first, second); } - for(plane_descriptor plane: q.other_coplanar) + for(plane_descriptor pd: q.other_coplanar) for(short int j=0; j!=3; ++j) - if(plane==p.supports[j]) + if(pd==p.supports[j]) if(first==-1){ - first=plane; + first=pd; break; - } else { - second=plane; + } else if(planes[first].first!=planes[pd].first){ + second=pd; return std::make_pair(first, second); } - for(plane_descriptor plane: p.other_coplanar) - for(plane_descriptor qlane: q.other_coplanar) - if(plane==qlane) + for(plane_descriptor pd: p.other_coplanar) + for(plane_descriptor qd: q.other_coplanar) + if(pd==qd) if(first==-1){ - first=plane; + first=pd; break; - } else { - second=plane; + } else if(planes[first].first!=planes[pd].first){ + second=pd; return std::make_pair(first, second); } // The two points do not shair a common support @@ -395,34 +407,25 @@ struct Plane_based_traits CGAL_assertion(0); }; - auto orient=[](const typename Kernel::Plane_3& plane, const Geometric_point_3& p){ - return sign (plane.a()*p.hx() + plane.b()*p.hy() + plane.c()*p.hz() + plane.d() * p.hw()); - }; std::pair line_supports=get_common_supports(p, q); - Point_3 res(plane.second, line_supports.first, line_supports.second, *m_planes); + assert((planes[line_supports.first].first != planes[line_supports.second].first) && + (planes[line_supports.first].first != planes[line_supports.second].first.opposite())); + Point_3 res=point_3(plane.second, line_supports.first, line_supports.second); namespace mp = boost::multiprecision; - CGAL_assertion_code(int256 max2E200=mp::pow(int256(2),200);) - CGAL_assertion((mp::abs(p.hx())<=max2E200) && (mp::abs(p.hy())<=max2E200) && (mp::abs(p.hz())<=max2E200) && (mp::abs(p.hw())<=max2E200)); + CGAL_assertion_code(int256 max2E195=mp::pow(int256(2),195);) + CGAL_assertion_code(int256 max2E169=mp::pow(int256(2),169);) + CGAL_assertion((mp::abs(p.hx())<=max2E195) && (mp::abs(p.hy())<=max2E195) && (mp::abs(p.hz())<=max2E195) && (mp::abs(p.hw())<=max2E169)); return res; } }; template Plane_based_traits(const Mesh &m){ - m_planes = new std::vector; - m_planes->reserve(faces(m).size()); - get_planes(m, std::back_inserter(*m_planes)); - } - - template - void get_planes(const Mesh &m, OutputIterator out){ + namespace mp = boost::multiprecision; using face_descriptor = typename boost::graph_traits::face_descriptor; using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; - Construct_plane_3 plane_3 = construct_plane_3_object(); - size_t k=0; - namespace mp = boost::multiprecision; auto to_int=[](const typename Mesh::Point &p){ return Geometric_point_3(int(p.x()),int(p.y()),int(p.z())); }; @@ -433,10 +436,16 @@ struct Plane_based_traits auto pl = compute_face_normal(f, m, parameters::vertex_point_map(pmap)); return typename Kernel::Vector_3(pl.hx(),pl.hy(),pl.hz()); }; + + m_planes = std::make_shared>(); + m_planes->reserve(faces(m).size()); + + // The ptr need to be initialized to get the functor + Construct_plane_3 plane_3 = construct_plane_3_object(); + for(face_descriptor f : faces(m)) plane_3(typename Kernel::Plane_3(to_int(m.point(m.target(m.halfedge(f)))), to_int_plane(f))); - } Oriented_side_3 oriented_side_3_object() const @@ -451,14 +460,14 @@ struct Plane_based_traits struct Construct_plane_3{ - std::vector* m_planes; - Construct_plane_3(std::vector* planes) : m_planes(planes){} + plane_range_pointer m_planes; + Construct_plane_3(plane_range_pointer planes) : m_planes(planes){} Plane_3 operator()(const typename Kernel::Plane_3 &pl){ namespace mp = boost::multiprecision; m_planes->emplace_back(pl, m_planes->size()); - CGAL_assertion_code(int256 max2E60=mp::pow(int256(2),60);) - CGAL_assertion_code(int256 max2E90=mp::pow(int256(2),90);) - CGAL_assertion((mp::abs(pl.a())<=max2E60) && (mp::abs(pl.b())<=max2E60) && (mp::abs(pl.c())<=max2E60) && (mp::abs(pl.d())<=max2E90)); + CGAL_assertion_code(int256 max2E55=mp::pow(int256(2),55);) + CGAL_assertion_code(int256 max2E82=mp::pow(int256(2),82);) + CGAL_assertion((mp::abs(pl.a())<=max2E55) && (mp::abs(pl.b())<=max2E55) && (mp::abs(pl.c())<=max2E82) && (mp::abs(pl.d())<=max2E82)); return m_planes->back(); } @@ -516,7 +525,7 @@ trettner_kernel(const TriangleMesh& pm, using InternMesh = Surface_mesh; GT gt(pm); - std::vector planes=gt.planes(); + const std::vector &planes=gt.planes(); Construct_plane_3 plane_3 = gt.construct_plane_3_object(); Construct_point_3 point_3 = gt.construct_point_3_object(); From f142aa040a5672df4fa837af22a767519a21da14 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Valque=20L=C3=A9o?= Date: Wed, 3 Dec 2025 14:08:08 +0100 Subject: [PATCH 14/16] Add clip_convex a specialization of clip for convex mesh --- .../Polygon_mesh_processing/clip_convex.h | 315 ++++++++++++++++++ .../Polygon_mesh_processing/internal/kernel.h | 104 +++++- .../test_mesh_kernel.cpp | 47 ++- 3 files changed, 457 insertions(+), 9 deletions(-) create mode 100644 Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h new file mode 100644 index 00000000000..3f6b461a55c --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h @@ -0,0 +1,315 @@ +// Copyright (c) 2024-2025 GeometryFactory (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) : Sébastien Loriot + + +#ifndef CGAL_POLYGON_MESH_PROCESSING_CLIP_CONVEX_H +#define CGAL_POLYGON_MESH_PROCESSING_CLIP_CONVEX_H + +#include + +#include +#include +#include +#include +#include +#ifndef CGAL_PLANE_CLIP_DO_NOT_USE_BOX_INTERSECTION_D +#include +#endif +#include + +namespace CGAL { +namespace Polygon_mesh_processing { + +template +void clip_convex(PolygonMesh& pm, + const typename GetGeomTraits::type::Plane_3& plane, + const NamedParameters& np = parameters::default_values()) +{ + using parameters::choose_parameter; + using parameters::get_parameter; + using parameters::get_parameter_reference; + using parameters::is_default_parameter; + + // graph typedefs + using BGT = boost::graph_traits; + using face_descriptor = typename BGT::face_descriptor; + using edge_descriptor = typename BGT::edge_descriptor; + using halfedge_descriptor = typename BGT::halfedge_descriptor; + using vertex_descriptor = typename BGT::vertex_descriptor; + + // np typedefs + using Default_ecm = Static_boolean_property_map; + using Default_visitor = Default_cut_visitor; + using Visitor_ref = typename internal_np::Lookup_named_param_def::reference; + using GT = typename GetGeomTraits::type; + GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits)); + + using FT = typename GT::FT; + using Point_3 = typename GT::Point_3; + using Vector_3 = typename GT::Vector_3; + + Default_visitor default_visitor; + Visitor_ref visitor = choose_parameter(get_parameter_reference(np, internal_np::visitor), default_visitor); + constexpr bool has_visitor = !std::is_same_v>>; + + auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_property_map(vertex_point, pm)); + + bool triangulate = !choose_parameter(get_parameter(np, internal_np::do_not_triangulate_faces), true); + if (triangulate && !is_triangle_mesh(pm)) + triangulate = false; + + bool throw_on_self_intersection = choose_parameter(get_parameter(np, internal_np::use_compact_clipper), false); + if (throw_on_self_intersection && !is_triangle_mesh(pm)) + throw_on_self_intersection = false; + + typedef typename internal_np::Lookup_named_param_def < + internal_np::concurrency_tag_t, + NamedParameters, + Sequential_tag + > ::type Concurrency_tag; + + // constexpr bool parallel_execution = std::is_same_v; + + auto oriented_side = traits.oriented_side_3_object(); + auto intersection_point = traits.construct_plane_line_intersection_point_3_object(); + auto scalar_product = traits.compute_scalar_product_3_object(); + auto point_on_plane = traits.construct_point_on_plane_3_object(); + // auto vector_3 = traits.construct_vector_3_object(); + + // TODO: the default is not thread-safe for example for Polyhedron + using V_os_tag = dynamic_vertex_property_t; + static constexpr bool use_default_vosm = + is_default_parameter::value; + + using Vertex_oriented_side_map = + std::conditional_t::type, + typename internal_np::Get_param::type>; + Vertex_oriented_side_map vertex_os; + + // ____________________ Find a crossing edge _____________________ + Vector_3 normal = traits.construct_orthogonal_vector_3_object()(plane); + + vertex_descriptor src=*vertices(pm).begin(); + Vector_3 vec(point_on_plane(plane), get(vpm,src)); + FT sp_src = scalar_product(vec, normal); // Not normalized distance to the plane + Sign direction_to_zero = sign(sp_src); + + vertex_descriptor trg; + FT sp_trg; + + bool is_local_max; + bool is_crossing_edge=false; + + if(direction_to_zero!=EQUAL){ + do{ + is_local_max=true; + for(auto v: vertices_around_target(src ,pm)){ + vec = Vector_3(point_on_plane(plane), get(vpm, v)); + sp_trg = scalar_product(vec, normal); + // TODO with EPICK, I need a predicate that compare scalar_product + // Check if v in the direction to the plane + if(compare(sp_src, sp_trg)==direction_to_zero){ + if(sign(sp_trg)!=direction_to_zero){ + // Fund a crossing edge + trg = v; + is_crossing_edge=true; + } else { + // Continue from v + sp_src = sp_trg; + src = v; + } + is_local_max=false; // repeat with the new vertex + break; + } + } + // No intersection with the plane, the is empty or full + if(is_local_max){ + if(direction_to_zero==POSITIVE) + pm.clear(); // The result is empty + return; + } + } while(!is_crossing_edge); + } else { + // src on the plane + trg=src; + sp_trg = sp_src; + } + + if(sign(sp_trg)==EQUAL && direction_to_zero!=POSITIVE){ + // Search a vertex around trg coming from positive side + bool no_positive_side = true; + bool no_negative_side = true; + for(auto v: vertices_around_target(trg ,pm)){ + Oriented_side side_v = oriented_side(plane, get(vpm, v)); + if(side_v==ON_POSITIVE_SIDE){ + src = v; + vec = Vector_3(point_on_plane(plane), get(vpm, src)); + sp_src = scalar_product(vec, normal); + no_positive_side = false; + break; + } + } + // Nothing to clip + if(no_positive_side) + return; + } else if(direction_to_zero==NEGATIVE){ + // Orient the edge from negative to positive + std::swap(src, trg); + } + + CGAL_assertion(oriented_side(plane, get(vpm, src)) == ON_POSITIVE_SIDE); + CGAL_assertion(oriented_side(plane, get(vpm, trg)) != ON_POSITIVE_SIDE); + + // Cut the convex along the plane by marching along crossing edges starting from the previous edge + std::vector boundaries; + std::set boundaries_vertices; + + halfedge_descriptor h = halfedge(src, trg, pm).first; + if(sign(sp_trg)!=EQUAL) + { + //split the first edge + auto pts = make_sorted_pair(get(vpm, src), get(vpm, trg)); + typename GT::Point_3 ip = intersection_point(plane, pts.first, pts.second); + //visitor.before_edge_split(h, pm); + h = CGAL::Euler::split_edge(h, pm); + put(vpm, target(h, pm), ip); + } + + vertex_descriptor v_start = target(h, pm); + halfedge_descriptor h_start=h; + do{ + halfedge_descriptor h_previous = h; + CGAL_assertion(oriented_side(plane, get(vpm, source(h,pm)))==ON_POSITIVE_SIDE); + CGAL_assertion(oriented_side(plane, get(vpm, target(h,pm)))==ON_ORIENTED_BOUNDARY); + + h = next(h, pm); + Oriented_side side_trg = oriented_side(plane, get(vpm, target(h,pm))); + + if(side_trg != ON_NEGATIVE_SIDE){ + // The face does not cross the plane + while(side_trg == ON_ORIENTED_BOUNDARY){ + // The edge is along the plane, add it to boundaries + boundaries.emplace_back(h); + boundaries_vertices.emplace(target(h, pm)); + pm.set_halfedge(target(h, pm), h); + + h = next(h, pm); + side_trg=oriented_side(plane, get(vpm, target(h,pm))); + } + // continue on next face + h = opposite(h, pm); + continue; + } + + // Search a crossing edge + h = next(h, pm); + side_trg=oriented_side(plane, get(vpm, target(h,pm))); + while(side_trg == ON_NEGATIVE_SIDE){ + h = next(h,pm); + side_trg=oriented_side(plane, get(vpm, target(h,pm))); + } + + if(side_trg != ON_ORIENTED_BOUNDARY){ + // Split the edge + auto pts = make_sorted_pair(get(vpm, source(h,pm)), get(vpm, target(h,pm))); + typename GT::Point_3 ip = intersection_point(plane, pts.first, pts.second); + // visitor.before_edge_split(h, pm); + h = CGAL::Euler::split_edge(h, pm); + put(vpm, target(h, pm), ip); + // visitor.new_vertex_added(vid, target(h,pm), pm); + // visitor.edge_split(h, pm); + // visitor.after_edge_split(); + + } + + // Split the face + halfedge_descriptor sh = CGAL::Euler::split_face(h_previous, h, pm); + boundaries.emplace_back(sh); + boundaries_vertices.emplace(target(sh, pm)); + pm.set_halfedge(target(sh, pm), sh); + + CGAL_assertion(target(sh, pm) == target(h, pm)); + h = opposite(next(sh,pm), pm); + } while(target(h, pm)!=v_start || (boundaries.size()==0 && h!=h_start)); + + CGAL_assertion(is_valid_polygon_mesh(pm)); + CGAL_assertion(boundaries.size()); + + // Remove the negative side + std::set vertices_to_remove; + std::set edges_to_remove; + std::set faces_to_remove; + + // Go through to find all elements to delete + face_descriptor start_face(face(halfedge(src, pm), pm)); + std::stack stack; + stack.push(start_face); + faces_to_remove.emplace(start_face); + + while (!stack.empty()) { + face_descriptor f = stack.top(); + stack.pop(); + + // Walk adjacent faces via halfedges + halfedge_descriptor h_start = halfedge(f, pm); + halfedge_descriptor h = h_start; + + do { + // if(boundaries.find(h)==boundaries.end()){ + if((std::find(boundaries.begin(), boundaries.end(), h)==boundaries.end()) && + (std::find(boundaries.begin(), boundaries.end(), opposite(h,pm))==boundaries.end())){ // TODO avoid this linear operation + edges_to_remove.emplace(edge(h, pm)); + if(boundaries_vertices.find(target(h, pm))==boundaries_vertices.end()) + vertices_to_remove.emplace(target(h, pm)); + CGAL_assertion(oriented_side(plane, get(vpm, target(h, pm)))!=ON_NEGATIVE_SIDE); + face_descriptor fn = face(opposite(h, pm), pm); + if (faces_to_remove.find(fn)==faces_to_remove.end()) { + faces_to_remove.emplace(fn); + stack.push(fn); + } + } + h = next(h, pm); + } while (h != h_start); + } + + for (vertex_descriptor v : vertices_to_remove){ + remove_vertex(v, pm); + CGAL_assertion(oriented_side(plane, get(vpm, v))==ON_POSITIVE_SIDE); + } + for (edge_descriptor e : edges_to_remove) + remove_edge(e, pm); + for (face_descriptor f : faces_to_remove) + remove_face(f, pm); + + // Reorder halfedges of the hole + for(size_t i=1; i #include +#include #include #include +#include + #include #include #include @@ -33,6 +36,7 @@ struct Three_point_cut_plane_traits using FT = typename Kernel::FT; using Plane_3 = std::array; using Point_3 = typename Kernel::Point_3; + using Vector_3 = typename Kernel::Vector_3; struct Does_not_support_CDT2{}; @@ -53,6 +57,20 @@ struct Three_point_cut_plane_traits } }; + struct Construct_orthogonal_vector_3{ + Vector_3 operator()(const Plane_3& plane) + { + return typename Kernel::Plane_3(plane[0], plane[1], plane[2]).orthogonal_vector(); + } + }; + + struct Construct_point_on_plane_3{ + Point_3 operator()(const Plane_3& plane) + { + return plane[0]; + } + }; + Oriented_side_3 oriented_side_3_object() const { return Oriented_side_3(); @@ -63,6 +81,19 @@ struct Three_point_cut_plane_traits return Construct_plane_line_intersection_point_3(); } + Construct_orthogonal_vector_3 construct_orthogonal_vector_3_object() const + { + return Construct_orthogonal_vector_3(); + } + + Construct_point_on_plane_3 construct_point_on_plane_3_object() const + { + return Construct_point_on_plane_3(); + } + + using Compute_scalar_product_3 = typename Kernel::Compute_scalar_product_3; + Compute_scalar_product_3 compute_scalar_product_3_object() const { return Compute_scalar_product_3(); } + #ifndef CGAL_PLANE_CLIP_DO_NOT_USE_BOX_INTERSECTION_D // for does self-intersect using Segment_3 = typename Kernel::Segment_3; @@ -186,8 +217,8 @@ kernel(const TriangleMesh& pm, } } #endif - - clip(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(kgt).do_not_triangulate_faces(true).used_for_kernel(true)); + clip_convex(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(kgt).do_not_triangulate_faces(true).used_for_kernel(true)); + // clip(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(kgt).do_not_triangulate_faces(true).used_for_kernel(true)); if (is_empty(kernel)) break; } @@ -290,6 +321,7 @@ struct Plane_based_traits using FT = typename Kernel::FT; using RT = typename Kernel::RT; using Geometric_point_3 = typename Kernel::Point_3; + using Vector_3 = typename Kernel::Vector_3; using plane_descriptor = std::size_t; using Plane_3 = std::pair; @@ -405,6 +437,7 @@ public: std::cout << q.supports[0] << " " << q.supports[1] << " " << q.supports[2] << std::endl; CGAL_assertion_code(std::cout << "The two points do not share two common supporting planes" << std::endl;) CGAL_assertion(0); + return std::make_pair(first, second); }; std::pair line_supports=get_common_supports(p, q); @@ -415,7 +448,7 @@ public: namespace mp = boost::multiprecision; CGAL_assertion_code(int256 max2E195=mp::pow(int256(2),195);) CGAL_assertion_code(int256 max2E169=mp::pow(int256(2),169);) - CGAL_assertion((mp::abs(p.hx())<=max2E195) && (mp::abs(p.hy())<=max2E195) && (mp::abs(p.hz())<=max2E195) && (mp::abs(p.hw())<=max2E169)); + // CGAL_assertion((mp::abs(p.hx())<=max2E195) && (mp::abs(p.hy())<=max2E195) && (mp::abs(p.hz())<=max2E195) && (mp::abs(p.hw())<=max2E169)); return res; } }; @@ -467,7 +500,7 @@ public: m_planes->emplace_back(pl, m_planes->size()); CGAL_assertion_code(int256 max2E55=mp::pow(int256(2),55);) CGAL_assertion_code(int256 max2E82=mp::pow(int256(2),82);) - CGAL_assertion((mp::abs(pl.a())<=max2E55) && (mp::abs(pl.b())<=max2E55) && (mp::abs(pl.c())<=max2E82) && (mp::abs(pl.d())<=max2E82)); + // CGAL_assertion((mp::abs(pl.a())<=max2E55) && (mp::abs(pl.b())<=max2E55) && (mp::abs(pl.c())<=max2E82) && (mp::abs(pl.d())<=max2E82)); return m_planes->back(); } @@ -553,7 +586,65 @@ trettner_kernel(const TriangleMesh& pm, for (auto plane: planes) { - clip(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(gt).do_not_triangulate_faces(true).used_for_kernel(true)); + // clip(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(gt).do_not_triangulate_faces(true).used_for_kernel(true)); + if (is_empty(kernel)) break; + } + + return kernel; +} + +template +Surface_mesh::Point_3> +trettner_epeck_kernel(const TriangleMesh& pm, + const NamedParameters& np = parameters::default_values()) +{ + using parameters::choose_parameter; + using parameters::get_parameter; + + // using GT = Plane_based_traits>; + using GT = Plane_based_traits; + auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_const_property_map(vertex_point, pm)); + + using Point_3 = typename GT::Point_3; + using Plane_3 = typename GT::Plane_3; + + using Construct_plane_3 = GT::Construct_plane_3; + using Construct_point_3 = GT::Construct_point_3; + + using InternMesh = Surface_mesh; + + GT gt(pm); + const std::vector &planes=gt.planes(); + + Construct_plane_3 plane_3 = gt.construct_plane_3_object(); + Construct_point_3 point_3 = gt.construct_point_3_object(); + + if (vertices(pm).size() - edges(pm).size() + faces(pm).size() != 2) + return Surface_mesh(); + + CGAL::Bbox_3 bb3 = bbox(pm, np); + InternMesh kernel; + Plane_3 xl=plane_3(1,0,0,int(-bb3.xmin())); + Plane_3 yl=plane_3(0,1,0,int(-bb3.ymin())); + Plane_3 zl=plane_3(0,0,1,int(-bb3.zmin())); + Plane_3 xr=plane_3(1,0,0,int(-bb3.xmax())); + Plane_3 yr=plane_3(0,1,0,int(-bb3.ymax())); + Plane_3 zr=plane_3(0,0,1,int(-bb3.zmax())); + CGAL::make_hexahedron(point_3(xl.second, yl.second, zl.second), + point_3(xl.second, yl.second, zr.second), + point_3(xl.second, yr.second, zr.second), + point_3(xl.second, yr.second, zl.second), + point_3(xr.second, yr.second, zl.second), + point_3(xr.second, yl.second, zl.second), + point_3(xr.second, yl.second, zr.second), + point_3(xr.second, yr.second, zr.second), + kernel); + + for (auto plane: planes) + { + // clip(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(gt).do_not_triangulate_faces(true).used_for_kernel(true)); if (is_empty(kernel)) break; } @@ -561,9 +652,6 @@ trettner_kernel(const TriangleMesh& pm, } - - - } } } // end of CGAL::Polygon_mesh_processing::experimental #endif // CGAL_POLYGON_MESH_PROCESSING_INTERNAL_KERNEL_H \ No newline at end of file diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp index 88c58aaf678..366734d4611 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp @@ -103,6 +103,31 @@ void test_traits() std::ofstream("clipped.off") << m; } +void elementary_test_kernel() +{ + using Point_3 = typename EK::Point_3; + CGAL::Real_timer timer; + EMesh m; + // make_hexahedron(Point_3(0, 1, 0), Point_3(0, 1, 1), Point_3(1, 0, 1), Point_3(1, 0, 0), + // Point_3(0.75, 0.25, 0.25), Point_3(0.25, 0.75, 0.25), Point_3(0.25, 0.75, 0.75), Point_3(0.75, 0.25, 0.75), + // m); + // make_tetrahedron(Point_3(1, 0, 0), Point_3(0, 1, 0), Point_3(0, 0, 1), Point_3(0.5, 0.5, 0.5), + // m); + // make_tetrahedron(Point_3(1, 0, 0), Point_3(1, 1, 0), Point_3(0, 1, 0), Point_3(0, 0, 0), + // m); + make_hexahedron(Point_3(1,0,0), Point_3(1,1,0), Point_3(0,1,0), Point_3(0,0,0), + Point_3(0,0,1), Point_3(1,0,1), Point_3(1,1,1), Point_3(0,1,1), + m); + + timer.start(); + EMesh kernel = PMP::experimental::kernel(m); + timer.stop(); + + std::ofstream("kernel.off") << kernel; + std::cout << "test_kernel done in " << timer.time() << "\n"; +} + + void test_kernel(std::string fname) { CGAL::Real_timer timer; @@ -228,15 +253,35 @@ void test_trettner_kernel(std::string fname) std::cout << "test_trettner_kernel done in " << timer.time() << "\n"; } +void test_trettner_epeck_kernel(std::string fname) +{ + CGAL::Real_timer timer; + Mesh m; + if (!CGAL::IO::read_polygon_mesh(fname, m)|| is_empty(m)) + { + std::cerr << "ERROR: cannot read " << fname << "\n"; + exit(1); + } + to_integer_mesh(m); + timer.start(); + auto kernel = PMP::experimental::trettner_epeck_kernel(m); + timer.stop(); + + std::ofstream("tkernel.off") << kernel; + std::cout << "test_trettner_kernel with epeck done in " << timer.time() << "\n"; +} + int main(int argc, char** argv) { // test_traits(); const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"); + // elementary_test_kernel(); // test_kernel(filename); // test_kernel_with_rounding(filename); // test_exact_kernel(filename); test_exact_kernel_with_rounding(filename); - test_trettner_kernel(filename); + // test_trettner_kernel(filename); + // test_trettner_epeck_kernel(filename); // test_kernel_with_chull(filename); // test_kernel_with_chull_and_constructions(filename); } From c380cdada36d79eb3884977efe1fb968ca761b6d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Valque=20L=C3=A9o?= Date: Wed, 3 Dec 2025 17:18:45 +0100 Subject: [PATCH 15/16] clean clip convex --- .../Polygon_mesh_processing/clip_convex.h | 66 ++++++-------- .../Polygon_mesh_processing/internal/kernel.h | 86 ++++++++++--------- 2 files changed, 73 insertions(+), 79 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h index 3f6b461a55c..55a6b5c9b1f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h @@ -1,4 +1,4 @@ -// Copyright (c) 2024-2025 GeometryFactory (France). +// Copyright (c) 2025 GeometryFactory (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -8,7 +8,7 @@ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // -// Author(s) : Sébastien Loriot +// Author(s) : Léo Valque #ifndef CGAL_POLYGON_MESH_PROCESSING_CLIP_CONVEX_H @@ -18,13 +18,6 @@ #include #include -#include -#include -#include -#ifndef CGAL_PLANE_CLIP_DO_NOT_USE_BOX_INTERSECTION_D -#include -#endif -#include namespace CGAL { namespace Polygon_mesh_processing { @@ -47,7 +40,7 @@ void clip_convex(PolygonMesh& pm, using vertex_descriptor = typename BGT::vertex_descriptor; // np typedefs - using Default_ecm = Static_boolean_property_map; + // using Default_ecm = Static_boolean_property_map; using Default_visitor = Default_cut_visitor; using Visitor_ref = typename internal_np::Lookup_named_param_def::reference; using GT = typename GetGeomTraits::type; @@ -72,18 +65,18 @@ void clip_convex(PolygonMesh& pm, if (throw_on_self_intersection && !is_triangle_mesh(pm)) throw_on_self_intersection = false; - typedef typename internal_np::Lookup_named_param_def < - internal_np::concurrency_tag_t, - NamedParameters, - Sequential_tag - > ::type Concurrency_tag; + // typedef typename internal_np::Lookup_named_param_def < + // internal_np::concurrency_tag_t, + // NamedParameters, + // Sequential_tag + // > ::type Concurrency_tag; // constexpr bool parallel_execution = std::is_same_v; auto oriented_side = traits.oriented_side_3_object(); auto intersection_point = traits.construct_plane_line_intersection_point_3_object(); - auto scalar_product = traits.compute_scalar_product_3_object(); - auto point_on_plane = traits.construct_point_on_plane_3_object(); + auto sq = traits.compute_squared_distance_3_object(); + // auto csq = traits.compare_squared_distance_3_object(); // auto vector_3 = traits.construct_vector_3_object(); // TODO: the default is not thread-safe for example for Polyhedron @@ -99,26 +92,22 @@ void clip_convex(PolygonMesh& pm, Vertex_oriented_side_map vertex_os; // ____________________ Find a crossing edge _____________________ - Vector_3 normal = traits.construct_orthogonal_vector_3_object()(plane); vertex_descriptor src=*vertices(pm).begin(); - Vector_3 vec(point_on_plane(plane), get(vpm,src)); - FT sp_src = scalar_product(vec, normal); // Not normalized distance to the plane + FT sp_src = sq(plane, get(vpm, src)); // Not normalized distance Sign direction_to_zero = sign(sp_src); vertex_descriptor trg; FT sp_trg; - bool is_local_max; bool is_crossing_edge=false; - if(direction_to_zero!=EQUAL){ do{ - is_local_max=true; + bool is_local_max=true; for(auto v: vertices_around_target(src ,pm)){ - vec = Vector_3(point_on_plane(plane), get(vpm, v)); - sp_trg = scalar_product(vec, normal); - // TODO with EPICK, I need a predicate that compare scalar_product + sp_trg = sq(plane, get(vpm, v)); + CGAL_assertion(sq(plane, get(vpm, v)) == sp_trg); + // TODO with EPICK, use compare_distance(plane, src, plane, trg) (But no possibility to memorize some computations for the next) // Check if v in the direction to the plane if(compare(sp_src, sp_trg)==direction_to_zero){ if(sign(sp_trg)!=direction_to_zero){ @@ -134,10 +123,10 @@ void clip_convex(PolygonMesh& pm, break; } } - // No intersection with the plane, the is empty or full + // No intersection with the plane, kernel is either empty or full if(is_local_max){ if(direction_to_zero==POSITIVE) - pm.clear(); // The result is empty + clear(pm); // The result is empty return; } } while(!is_crossing_edge); @@ -155,8 +144,8 @@ void clip_convex(PolygonMesh& pm, Oriented_side side_v = oriented_side(plane, get(vpm, v)); if(side_v==ON_POSITIVE_SIDE){ src = v; - vec = Vector_3(point_on_plane(plane), get(vpm, src)); - sp_src = scalar_product(vec, normal); + sp_src = sq(plane, get(vpm, src)); + // CGAL_assertion(sq(plane, get(vpm, src)) == sp_src); no_positive_side = false; break; } @@ -203,7 +192,7 @@ void clip_convex(PolygonMesh& pm, // The edge is along the plane, add it to boundaries boundaries.emplace_back(h); boundaries_vertices.emplace(target(h, pm)); - pm.set_halfedge(target(h, pm), h); + set_halfedge(target(h, pm), h, pm); h = next(h, pm); side_trg=oriented_side(plane, get(vpm, target(h,pm))); @@ -238,14 +227,14 @@ void clip_convex(PolygonMesh& pm, halfedge_descriptor sh = CGAL::Euler::split_face(h_previous, h, pm); boundaries.emplace_back(sh); boundaries_vertices.emplace(target(sh, pm)); - pm.set_halfedge(target(sh, pm), sh); + set_halfedge(target(sh, pm), sh, pm); CGAL_assertion(target(sh, pm) == target(h, pm)); h = opposite(next(sh,pm), pm); - } while(target(h, pm)!=v_start || (boundaries.size()==0 && h!=h_start)); + } while(target(h, pm)!=v_start || (boundaries.empty() && h!=h_start)); CGAL_assertion(is_valid_polygon_mesh(pm)); - CGAL_assertion(boundaries.size()); + CGAL_assertion(!boundaries.empty()); // Remove the negative side std::set vertices_to_remove; @@ -267,7 +256,6 @@ void clip_convex(PolygonMesh& pm, halfedge_descriptor h = h_start; do { - // if(boundaries.find(h)==boundaries.end()){ if((std::find(boundaries.begin(), boundaries.end(), h)==boundaries.end()) && (std::find(boundaries.begin(), boundaries.end(), opposite(h,pm))==boundaries.end())){ // TODO avoid this linear operation edges_to_remove.emplace(edge(h, pm)); @@ -295,15 +283,15 @@ void clip_convex(PolygonMesh& pm, // Reorder halfedges of the hole for(size_t i=1; i &planes = *m_planes; auto res = CGAL::Intersections::internal::intersection_point(planes[a].first, planes[b].first, planes[c].first, Kernel()); - CGAL_assertion(res); return Point_3(a, b, c, *res); } @@ -391,7 +386,7 @@ public: plane_descriptor second=-1; for(short int i=0; i!=3; ++i) for(short int j=0; j!=3; ++j) - if(p.supports[i]==q.supports[j]) + if(p.supports[i]==q.supports[j]){ if(first==-1){ first=p.supports[i]; break; @@ -399,10 +394,11 @@ public: second=p.supports[i]; return std::make_pair(first, second); } + } for(plane_descriptor pd: p.other_coplanar) for(short int j=0; j!=3; ++j) - if(pd==q.supports[j]) + if(pd==q.supports[j]){ if(first==-1){ first=pd; break; @@ -410,10 +406,11 @@ public: second=pd; return std::make_pair(first, second); } + } for(plane_descriptor pd: q.other_coplanar) for(short int j=0; j!=3; ++j) - if(pd==p.supports[j]) + if(pd==p.supports[j]){ if(first==-1){ first=pd; break; @@ -421,10 +418,11 @@ public: second=pd; return std::make_pair(first, second); } + } for(plane_descriptor pd: p.other_coplanar) for(plane_descriptor qd: q.other_coplanar) - if(pd==qd) + if(pd==qd){ if(first==-1){ first=pd; break; @@ -432,6 +430,7 @@ public: second=pd; return std::make_pair(first, second); } + } // The two points do not shair a common support std::cout << p.supports[0] << " " << p.supports[1] << " " << p.supports[2] << std::endl; std::cout << q.supports[0] << " " << q.supports[1] << " " << q.supports[2] << std::endl; @@ -453,6 +452,28 @@ public: } }; + struct Construct_orthogonal_vector_3 + { + Vector_3 operator()(const Plane_3& plane, const Point_3& p) const + { + if((p.supports[0]==plane.second) || (p.supports[1]==plane.second) || (p.supports[2]==plane.second)) + return COPLANAR; + Oriented_side ori = plane.first.oriented_side(p); + if(ori==COPLANAR) + p.other_coplanar.emplace(plane.second); + return ori; + } + }; + + struct Compute_squared_distance_3 + { + using Compute_scalar_product_3 = typename Kernel::Compute_scalar_product_3; + FT operator()(const Plane_3& plane, const Point_3& p) + { + return Compute_scalar_product_3()(Vector_3(ORIGIN, p), plane.first.orthogonal_vector())+plane.first.d(); + } + }; + template Plane_based_traits(const Mesh &m){ namespace mp = boost::multiprecision; @@ -481,16 +502,6 @@ public: to_int_plane(f))); } - Oriented_side_3 oriented_side_3_object() const - { - return Oriented_side_3(); - } - - Construct_plane_line_intersection_point_3 construct_plane_line_intersection_point_3_object() - { - return Construct_plane_line_intersection_point_3(m_planes); - } - struct Construct_plane_3{ plane_range_pointer m_planes; @@ -509,14 +520,13 @@ public: } }; - Construct_plane_3 construct_plane_3_object() - { - return Construct_plane_3(m_planes); - } - - Construct_point_3 construct_point_3_object() - { - return Construct_point_3(m_planes); + Oriented_side_3 oriented_side_3_object() const { return Oriented_side_3(); } + Construct_plane_3 construct_plane_3_object() { return Construct_plane_3(m_planes); } + Construct_point_3 construct_point_3_object() { return Construct_point_3(m_planes); } + Construct_orthogonal_vector_3 construct_orthogonal_vector_3_object() const { return Construct_orthogonal_vector_3(); } + Compute_squared_distance_3 compute_squared_distance_3_object() const { return Compute_squared_distance_3(); } + Construct_plane_line_intersection_point_3 construct_plane_line_intersection_point_3_object() const { + return Construct_plane_line_intersection_point_3(m_planes); } const std::vector& planes(){ @@ -546,8 +556,6 @@ trettner_kernel(const TriangleMesh& pm, using parameters::get_parameter; using GT = Plane_based_traits>; - auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), - get_const_property_map(vertex_point, pm)); using Point_3 = typename GT::Point_3; using Plane_3 = typename GT::Plane_3; @@ -586,7 +594,7 @@ trettner_kernel(const TriangleMesh& pm, for (auto plane: planes) { - // clip(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(gt).do_not_triangulate_faces(true).used_for_kernel(true)); + clip_convex(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(gt).do_not_triangulate_faces(true).used_for_kernel(true)); if (is_empty(kernel)) break; } @@ -604,8 +612,6 @@ trettner_epeck_kernel(const TriangleMesh& pm, // using GT = Plane_based_traits>; using GT = Plane_based_traits; - auto vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), - get_const_property_map(vertex_point, pm)); using Point_3 = typename GT::Point_3; using Plane_3 = typename GT::Plane_3; @@ -644,7 +650,7 @@ trettner_epeck_kernel(const TriangleMesh& pm, for (auto plane: planes) { - // clip(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(gt).do_not_triangulate_faces(true).used_for_kernel(true)); + clip_convex(kernel, plane, CGAL::parameters::clip_volume(true).geom_traits(gt).do_not_triangulate_faces(true).used_for_kernel(true)); if (is_empty(kernel)) break; } From 0fb60af4e0c4820718b73c937036b53efb938d7a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Valque=20L=C3=A9o?= Date: Wed, 3 Dec 2025 18:09:21 +0100 Subject: [PATCH 16/16] use vector and binary_search instead of set for boundaries vertices --- .../CGAL/Polygon_mesh_processing/clip_convex.h | 15 +++++++-------- .../Polygon_mesh_processing/internal/kernel.h | 4 ++-- .../Polygon_mesh_processing/test_mesh_kernel.cpp | 4 ++-- 3 files changed, 11 insertions(+), 12 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h index 55a6b5c9b1f..1f275a35c2d 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h @@ -48,7 +48,6 @@ void clip_convex(PolygonMesh& pm, using FT = typename GT::FT; using Point_3 = typename GT::Point_3; - using Vector_3 = typename GT::Vector_3; Default_visitor default_visitor; Visitor_ref visitor = choose_parameter(get_parameter_reference(np, internal_np::visitor), default_visitor); @@ -139,13 +138,11 @@ void clip_convex(PolygonMesh& pm, if(sign(sp_trg)==EQUAL && direction_to_zero!=POSITIVE){ // Search a vertex around trg coming from positive side bool no_positive_side = true; - bool no_negative_side = true; for(auto v: vertices_around_target(trg ,pm)){ Oriented_side side_v = oriented_side(plane, get(vpm, v)); if(side_v==ON_POSITIVE_SIDE){ src = v; sp_src = sq(plane, get(vpm, src)); - // CGAL_assertion(sq(plane, get(vpm, src)) == sp_src); no_positive_side = false; break; } @@ -163,7 +160,7 @@ void clip_convex(PolygonMesh& pm, // Cut the convex along the plane by marching along crossing edges starting from the previous edge std::vector boundaries; - std::set boundaries_vertices; + std::vector boundaries_vertices; halfedge_descriptor h = halfedge(src, trg, pm).first; if(sign(sp_trg)!=EQUAL) @@ -191,7 +188,7 @@ void clip_convex(PolygonMesh& pm, while(side_trg == ON_ORIENTED_BOUNDARY){ // The edge is along the plane, add it to boundaries boundaries.emplace_back(h); - boundaries_vertices.emplace(target(h, pm)); + boundaries_vertices.emplace_back(target(h, pm)); set_halfedge(target(h, pm), h, pm); h = next(h, pm); @@ -226,7 +223,7 @@ void clip_convex(PolygonMesh& pm, // Split the face halfedge_descriptor sh = CGAL::Euler::split_face(h_previous, h, pm); boundaries.emplace_back(sh); - boundaries_vertices.emplace(target(sh, pm)); + boundaries_vertices.emplace_back(target(sh, pm)); set_halfedge(target(sh, pm), sh, pm); CGAL_assertion(target(sh, pm) == target(h, pm)); @@ -236,11 +233,13 @@ void clip_convex(PolygonMesh& pm, CGAL_assertion(is_valid_polygon_mesh(pm)); CGAL_assertion(!boundaries.empty()); - // Remove the negative side + // ________________________ Remove the negative side _________________________ std::set vertices_to_remove; std::set edges_to_remove; std::set faces_to_remove; + std::sort(boundaries_vertices.begin(), boundaries_vertices.end()); + // Go through to find all elements to delete face_descriptor start_face(face(halfedge(src, pm), pm)); std::stack stack; @@ -259,7 +258,7 @@ void clip_convex(PolygonMesh& pm, if((std::find(boundaries.begin(), boundaries.end(), h)==boundaries.end()) && (std::find(boundaries.begin(), boundaries.end(), opposite(h,pm))==boundaries.end())){ // TODO avoid this linear operation edges_to_remove.emplace(edge(h, pm)); - if(boundaries_vertices.find(target(h, pm))==boundaries_vertices.end()) + if(!std::binary_search(boundaries_vertices.begin(), boundaries_vertices.end(), target(h, pm))) vertices_to_remove.emplace(target(h, pm)); CGAL_assertion(oriented_side(plane, get(vpm, target(h, pm)))!=ON_NEGATIVE_SIDE); face_descriptor fn = face(opposite(h, pm), pm); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h index e1c422d072b..28117776a20 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h @@ -348,8 +348,8 @@ public: Construct_point_3(plane_range_pointer planes) : m_planes(planes){} Point_3 operator()(plane_descriptor a, plane_descriptor b, plane_descriptor c){ const std::vector &planes = *m_planes; - auto res = CGAL::Intersections::internal::intersection_point(planes[a].first, planes[b].first, planes[c].first, Kernel()); - return Point_3(a, b, c, *res); + auto res = intersection(planes[a].first, planes[b].first, planes[c].first); + return Point_3(a, b, c, std::get(*res)); } Point_3 operator()(plane_descriptor a, plane_descriptor b, plane_descriptor c, Geometric_point_3 p){ diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp index 366734d4611..4cc977354b9 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp @@ -279,9 +279,9 @@ int main(int argc, char** argv) // test_kernel(filename); // test_kernel_with_rounding(filename); // test_exact_kernel(filename); - test_exact_kernel_with_rounding(filename); + // test_exact_kernel_with_rounding(filename); // test_trettner_kernel(filename); - // test_trettner_epeck_kernel(filename); + test_trettner_epeck_kernel(filename); // test_kernel_with_chull(filename); // test_kernel_with_chull_and_constructions(filename); }