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 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/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/clip_convex.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h new file mode 100644 index 00000000000..1f275a35c2d --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip_convex.h @@ -0,0 +1,302 @@ +// Copyright (c) 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) : Léo Valque + + +#ifndef CGAL_POLYGON_MESH_PROCESSING_CLIP_CONVEX_H +#define CGAL_POLYGON_MESH_PROCESSING_CLIP_CONVEX_H + +#include + +#include +#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; + + 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 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 + 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 _____________________ + + vertex_descriptor src=*vertices(pm).begin(); + 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_crossing_edge=false; + if(direction_to_zero!=EQUAL){ + do{ + bool is_local_max=true; + for(auto v: vertices_around_target(src ,pm)){ + 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){ + // 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, kernel is either empty or full + if(is_local_max){ + if(direction_to_zero==POSITIVE) + clear(pm); // 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; + 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)); + 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::vector 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_back(target(h, pm)); + set_halfedge(target(h, pm), h, pm); + + 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_back(target(sh, pm)); + 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.empty() && h!=h_start)); + + CGAL_assertion(is_valid_polygon_mesh(pm)); + CGAL_assertion(!boundaries.empty()); + + // ________________________ 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; + 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((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(!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); + 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 +#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; + using Vector_3 = typename Kernel::Vector_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; + return construction(plane[0], plane[1], plane[2], p, q); + } + }; + + 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 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) + { + typename Kernel::Plane_3 pl(plane[0], plane[1], plane[2]); + return Compute_scalar_product_3()(Vector_3(ORIGIN, p), pl.orthogonal_vector())+pl.d(); + } + }; + + 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(); + } + + 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(); } + +#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; + + //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); + 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; +#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(); + 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()), + 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_convex(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)) ); + } + + auto origin_opt = CGAL::halfspace_intersection_interior_point_3(planes.begin(), planes.end()); + + TriangleMesh kernel; + + if (origin_opt.has_value()) + CGAL::halfspace_intersection_3(planes.begin(), + planes.end(), + kernel, + origin_opt); + + 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)) ); + } + + auto origin_opt = CGAL::halfspace_intersection_interior_point_3(planes.begin(), planes.end()); + + TriangleMesh kernel; + + if (origin_opt.has_value()) + CGAL::halfspace_intersection_with_constructions_3(planes.begin(), + planes.end(), + kernel, + origin_opt); + + return kernel; +} + + + +//__________________________________________________________________________________________________________________________________ + +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 Vector_3 = typename Kernel::Vector_3; + + using plane_descriptor = std::size_t; + using Plane_3 = std::pair; + +private: + using plane_range_pointer = std::shared_ptr>; + plane_range_pointer m_planes; +public: + + + + 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}){} + }; + + struct Construct_point_3{ + 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){ + const std::vector &planes = *m_planes; + 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){ + 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 + { + 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) + { + 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) + 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 pd: p.other_coplanar) + for(short int j=0; j!=3; ++j) + if(pd==q.supports[j]){ + if(first==-1){ + first=pd; + break; + } else if(planes[first].first!=planes[pd].first){ + 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(first==-1){ + first=pd; + break; + } else if(planes[first].first!=planes[pd].first){ + 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(first==-1){ + first=pd; + break; + } 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 + 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); + return std::make_pair(first, second); + }; + + std::pair line_supports=get_common_supports(p, q); + 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 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; + } + }; + + 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; + using face_descriptor = typename boost::graph_traits::face_descriptor; + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + + 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, 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))); + } + + struct Construct_plane_3{ + + 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 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(); + } + + 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)); + } + }; + + 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(){ + 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 +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 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_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; + } + + 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; + + 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_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; + } + + 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..4cc977354b9 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_kernel.cpp @@ -0,0 +1,287 @@ +#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; + +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=[&](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); + } +} + +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; + 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 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; + Mesh 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); + timer.stop(); + + std::ofstream("kernel.off") << kernel; + 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; + EMesh 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(); + + std::ofstream("ekernel.off") << kernel; + 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); + to_integer_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; + Mesh 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); + 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(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); + } + + 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"; +} + +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(); + auto kernel = PMP::experimental::trettner_kernel(m); + timer.stop(); + + std::ofstream("tkernel.off") << kernel; + 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_epeck_kernel(filename); + // test_kernel_with_chull(filename); + // test_kernel_with_chull_and_constructions(filename); +} 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)