use vector and binary_search instead of set for boundaries vertices

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 <CGAL/license/Polygon_mesh_processing/corefinement.h>
+
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+namespace CGAL {
+namespace Polygon_mesh_processing {
+
+template <class PolygonMesh, class NamedParameters =  parameters::Default_named_parameters>
+void clip_convex(PolygonMesh& pm,
+                 const typename GetGeomTraits<PolygonMesh, NamedParameters>::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<PolygonMesh>;
+  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<edge_descriptor, false>;
+  using Default_visitor = Default_cut_visitor<PolygonMesh>;
+  using Visitor_ref = typename internal_np::Lookup_named_param_def<internal_np::visitor_t, NamedParameters, Default_visitor>::reference;
+  using GT = typename GetGeomTraits<PolygonMesh, NamedParameters>::type;
+  GT traits = choose_parameter<GT>(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<Default_visitor, std::remove_cv_t<std::remove_reference_t<Visitor_ref>>>;
+
+  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<Parallel_tag, Concurrency_tag>;
+
+  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<Oriented_side>;
+  static constexpr bool use_default_vosm =
+    is_default_parameter<NamedParameters, internal_np::vertex_oriented_side_map_t>::value;
+
+  using Vertex_oriented_side_map =
+    std::conditional_t<use_default_vosm,
+                       typename boost::property_map<PolygonMesh, V_os_tag>::type,
+                       typename internal_np::Get_param<typename NamedParameters::base,
+                                                       internal_np::vertex_oriented_side_map_t>::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<halfedge_descriptor> boundaries;
+  std::vector<vertex_descriptor> 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<vertex_descriptor> vertices_to_remove;
+  std::set<edge_descriptor> edges_to_remove;
+  std::set<face_descriptor> 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<face_descriptor> 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<boundaries.size(); ++i)
+    set_next(boundaries[i-1], boundaries[i], pm);
+  set_next(boundaries.back(), boundaries[0], pm);
+
+  // Fill the hole
+  face_descriptor f=pm.add_face();
+  for(auto h: boundaries){
+    set_face(h, f, pm);
+  }
+  set_halfedge(f, boundaries[0], pm);
+
+  // std::ofstream("clip.off") << pm;
+  CGAL_assertion(is_valid_polygon_mesh(pm));
+}
+
+} } // CGAL::Polygon_mesh_processing
+
+
+#endif // CGAL_POLYGON_MESH_PROCESSING_CLIP_CONVEX_H

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/kernel.h

@@ -15,9 +15,12 @@
 
 #include <CGAL/license/Polygon_mesh_processing/corefinement.h>
 #include <CGAL/Polygon_mesh_processing/clip.h>
+#include <CGAL/Polygon_mesh_processing/clip_convex.h>
 #include <CGAL/Convex_hull_3/dual/halfspace_intersection_3.h>
 #include <CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h>
 
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
 #include <CGAL/Homogeneous.h>
 #include <CGAL/Exact_integer.h>
 #include <CGAL/Surface_mesh.h>
@@ -33,6 +36,7 @@ struct Three_point_cut_plane_traits
   using FT = typename Kernel::FT;
   using Plane_3 = std::array<typename Kernel::Point_3, 3>;
   using Point_3 = typename Kernel::Point_3;
+  using Vector_3 = typename Kernel::Vector_3;
 
   struct Does_not_support_CDT2{};
 
@@ -53,6 +57,23 @@ 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 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();
@@ -63,6 +84,13 @@ 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();
+  }
+
+  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;
@@ -186,8 +214,7 @@ kernel(const TriangleMesh& pm,
       }
     }
 #endif
-
+    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 +317,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<typename Kernel::Plane_3, plane_descriptor>;
@@ -320,9 +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<Plane_3> &planes = *m_planes;
-      auto res = CGAL::Intersections::internal::intersection_point(planes[a].first, planes[b].first, planes[c].first, Kernel());
+      auto res = intersection(planes[a].first, planes[b].first, planes[c].first);
-      CGAL_assertion(res);
+      return Point_3(a, b, c, std::get<Geometric_point_3>(*res));
-      return Point_3(a, b, c, *res);
     }
 
     Point_3 operator()(plane_descriptor a, plane_descriptor b, plane_descriptor c, Geometric_point_3 p){
@@ -359,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;
@@ -367,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;
@@ -378,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;
@@ -389,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;
@@ -400,11 +430,13 @@ 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;
         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<plane_descriptor, plane_descriptor> line_supports=get_common_supports(p, q);
@@ -415,11 +447,33 @@ 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;
     }
   };
 
+  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<class Mesh>
   Plane_based_traits(const Mesh &m){
     namespace mp = boost::multiprecision;
@@ -448,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;
@@ -467,7 +511,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();
     }
 
@@ -476,14 +520,13 @@ public:
     }
   };
 
-  Construct_plane_3 construct_plane_3_object()
+  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); }
-    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_point_3 construct_point_3_object()
+  Construct_plane_line_intersection_point_3 construct_plane_line_intersection_point_3_object() const {
-  {
+    return Construct_plane_line_intersection_point_3(m_planes);
-    return Construct_point_3(m_planes);
   }
 
   const std::vector<Plane_3>& planes(){
@@ -513,8 +556,6 @@ trettner_kernel(const TriangleMesh& pm,
   using parameters::get_parameter;
 
   using GT = Plane_based_traits<Homogeneous<int256>>;
-  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;
@@ -553,7 +594,63 @@ 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;
+  }
+
+  return kernel;
+}
+
+template <class TriangleMesh,
+          class NamedParameters = parameters::Default_named_parameters>
+Surface_mesh<Plane_based_traits<Exact_predicates_exact_constructions_kernel>::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<Homogeneous<int256>>;
+  using GT = Plane_based_traits<Exact_predicates_exact_constructions_kernel>;
+
+  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<Point_3>;
+
+  GT gt(pm);
+  const std::vector<Plane_3> &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<Point_3>();
+
+  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;
   }
 
@@ -561,9 +658,6 @@ trettner_kernel(const TriangleMesh& pm,
 }
 
 
-
-
-
 } } } // end of CGAL::Polygon_mesh_processing::experimental
 
 #endif // CGAL_POLYGON_MESH_PROCESSING_INTERNAL_KERNEL_H

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_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);
 }