Merge pull request #5435 from afabri/Intersect_3-fix_isocuboid_plane-GF

Intersect_3: Fix Iso_cuboid_3/Plane_3 and Tetrahedon_/Plane_3 intersection using PMP::clip internal function

# Conflicts:
#	Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h

Intersections_3/include/CGAL/Intersections_3/internal/Iso_cuboid_3_Plane_3_intersection.h

@@ -8,7 +8,7 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
 //
 //
-// Author(s)     : Maxime Gimeno
+// Author(s)     : Sebastien Loriot and Andreas Fabri
 //
 
 #ifndef CGAL_INTERNAL_INTERSECTIONS_3_ISO_CUBOID_3_PLANE_3_INTERSECTION_H
@@ -16,320 +16,195 @@
 
 #include <CGAL/kernel_basic.h>
 #include <CGAL/intersections.h>
+#include <CGAL/utility.h>
 
-#include <CGAL/Intersections_3/internal/tetrahedron_intersection_helpers.h>
+#include <vector>
-#include <CGAL/Intersections_3/Iso_cuboid_3_Segment_3.h>
+#include <array>
-#include <CGAL/Intersections_3/Plane_3_Plane_3.h>
-
-#include <set>
 
 namespace CGAL {
 namespace Intersections {
 namespace internal {
 
-template<typename Point>
+//Iso_cuboid_3 Plane_3
-void filter_points(std::vector<Point>& input,
-                   std::vector<Point>& output)
-{
-  std::sort(input.begin(), input.end());
-  auto last = std::unique(input.begin(), input.end());
-  input.erase(last, input.end());
-  output = std::move(input);
-}
-
-//Tetrahedron_3 Line_3
 template <class K>
 typename Intersection_traits<K, typename K::Iso_cuboid_3, typename K::Plane_3>::result_type
 intersection(
              const typename K::Iso_cuboid_3& cub,
-    const typename K::Plane_3 &pl,
+             const typename K::Plane_3& plane,
              const K& k)
 {
   typedef typename K::Point_3 Point_3;
-  typedef typename K::Segment_3 Segment_3;
+  typedef typename Intersection_traits<K, typename K::Iso_cuboid_3, typename K::Plane_3>::result_type result_type;
-  typedef typename K::Line_3 Line_3;
+  typename K::Oriented_side_3 oriented_side = k.oriented_side_3_object();
-  typedef typename K::Plane_3 Plane_3;
-  typedef std::vector<Point_3> Poly;
 
-  typedef typename Intersection_traits<K, CGAL::Iso_cuboid_3<K>,
+  std::vector<Point_3> corners(8);
-      CGAL::Plane_3<K> >::result_type result_type;
+  corners.reserve(14); // 8 corners + up to 6 polygon points
+  corners[0] = cub[0];
+  corners[1] = cub[3];
+  corners[2] = cub[2];
+  corners[3] = cub[1];
+  corners[4] = cub[5];
+  corners[5] = cub[4];
+  corners[6] = cub[7];
+  corners[7] = cub[6];
 
-  typedef typename Intersection_traits<K, CGAL::Segment_3<K>,
+  const std::array<CGAL::Oriented_side, 8> orientations  { {
-      CGAL::Plane_3<K> >::result_type Inter_type;
+      oriented_side(plane, corners[0]),
+      oriented_side(plane, corners[1]),
+      oriented_side(plane, corners[2]),
+      oriented_side(plane, corners[3]),
+      oriented_side(plane, corners[4]),
+      oriented_side(plane, corners[5]),
+      oriented_side(plane, corners[6]),
+      oriented_side(plane, corners[7])
+    } };
 
-  std::vector<Segment_3> edges;
+  // description of faces of the bbox
-  edges.reserve(12);
+  constexpr std::array<int, 24> face_indices
+    { { 0, 1, 2, 3,
+        2, 1, 5, 6,
+        3, 2, 6, 7,
+        1, 0, 4, 5,
+        4, 0, 3, 7,
+        6, 5, 4, 7 } };
 
-  //get all edges of cub
+  constexpr std::array<int, 24> edge_indices
-  for(int i=0; i< 4; ++i)
+    { { 0,  1,  2, 3,
-    edges.emplace_back(cub.vertex(i), cub.vertex((i+1)%4));
+        1,  4,  5, 6,
-  for(int i=0; i < 4; ++i)
+        2,  6,  7, 8,
-    edges.emplace_back(cub.vertex(i+4), cub.vertex((i+1)%4+4));
+        0,  9, 10, 4,
-  for(int i=0; i < 4; ++i)
+        9,  3,  8, 11,
-    edges.emplace_back(cub.vertex(i), cub.vertex((i+1)%4+4));
+        5, 10, 11, 7 } };
 
-  //get all intersections between pl and cub edges
+  std::array<int, 12> edge_ipt_id;
-  std::vector<Segment_3> segments;
+  edge_ipt_id.fill(-1);
-  std::vector<Point_3> points;
 
-  for(int i=0; i < 12; ++i)
+  auto inter_pt_index =
+    [&plane, &corners, &edge_ipt_id](int i, int j, int edge_id)
   {
-    // Intersect_3 checks the orientation of the segment's extremities to avoid actually computing
+    if (edge_ipt_id[edge_id]==-1)
-    // the intersection if possible
-    Inter_type inter = typename K::Intersect_3()(pl, edges[i]);
-    if(inter)
     {
-      if(const Segment_3* seg = boost::get<Segment_3>(&*inter))
+      edge_ipt_id[edge_id] = static_cast<int> (corners.size());
-        segments.push_back(*seg);
+      corners.push_back(typename K::Construct_plane_line_intersection_point_3()
-      else if(const Point_3* p = boost::get<Point_3>(&*inter))
+                      (plane, corners[i], corners[j]));
-        points.push_back(*p);
-    }
     }
 
-  if(segments.empty() && points.empty())
+    return edge_ipt_id[edge_id];
-    return result_type();
+  };
 
-  switch(segments.size())
+  bool all_in = true;
-  {
+  bool all_out = true;
-    case 0:
-    //points dealt with later
-    break;
-    case 1:
-    {
-      //adj to an edge
-      if(points.size() == 4)
-      {
-        return result_type(std::forward<Segment_3>(segments.front()));
-      }
-      //plane intersecting through an edge (not 2)
-      else
-      {
-        Poly res(4);
-        const Segment_3& entry_seg(segments.front());
-        Point_3 p1, p2;
-        bool p1_found(false),
-             p2_found(false);
 
-        for(const Point_3& p : points)
+  std::vector<std::array<int,2>> neighbor_ids(14, {-1,-1});
+
+  auto add_neighbor = [&neighbor_ids](int i, int j)
   {
-          if(!k.equal_3_object()(entry_seg.source(), p)
+    if (neighbor_ids[i][0] == -1 ) {
-             && ! k.equal_3_object()(entry_seg.target(), p))
+      neighbor_ids[i][0] = j;
-          {
-            if(!p1_found)
-            {
-              p1 = p;
-              p1_found = true;
     }
     else {
-              p2 = p;
+      if (neighbor_ids[i][0]!=j && neighbor_ids[i][1]==-1)
-              p2_found = true;
+      {
+        neighbor_ids[i][1] = j;
+      }
+    }
+  };
+
+  int start_id = -1;
+  int solo_id = -1;
+  // for each face of the bbox, we look for intersection of the plane with its edges
+  std::vector<int> ids;
+  for (int i = 0; i < 6; ++i)
+  {
+    ids.clear();
+    for (int k = 0; k < 4; ++k)
+    {
+
+      int current_id = face_indices[4 * i + k];
+      int next_id = face_indices[4 * i + (k + 1) % 4];
+      int edge_id = edge_indices[4 * i + k];
+
+      switch (orientations[current_id])
+      {
+        case ON_NEGATIVE_SIDE:
+        {
+          all_out = false;
+          // check for intersection of the edge
+          if (orientations[next_id] == ON_POSITIVE_SIDE)
+          {
+            ids.push_back(
+                          inter_pt_index(current_id, next_id, edge_id));
+          }
           break;
         }
-          }
+        case ON_POSITIVE_SIDE:
-        }
-        CGAL_USE(p2_found);
-        CGAL_assertion(p1_found && p2_found);
-        res[0] = entry_seg.target();
-        res[1] = p2;
-        res[2] = p1;
-        res[3] = entry_seg.source();
-
-        typename K::Coplanar_orientation_3 coplanar_orientation =
-          k.coplanar_orientation_3_object();
-
-        if( coplanar_orientation(res[0], res[1], res[2])
-            != coplanar_orientation(res[0], res[1], res[3]))
         {
-          std::swap(res[1], res[2]);
+          all_in = false;
-        }
+          // check for intersection of the edge
-
+          if (orientations[next_id] == ON_NEGATIVE_SIDE)
-        return result_type(std::forward<Poly>(res));
+          {
-      }
+            ids.push_back(inter_pt_index(current_id, next_id, edge_id));
           }
           break;
-
+        }
-    case 2: //intersects diagonally
+        case ON_ORIENTED_BOUNDARY:
         {
-      Poly res(4);
+          all_in = all_out = false;
-      Segment_3 &front(segments.front()),
+          ids.push_back(current_id);
-          &back(segments.back());
+        }
-      res[0] = front.target();
+      }
-      res[1] = back.target();
-      res[2] = back.source();
-      res[3] = front.source();
-      typename K::Coplanar_orientation_3 coplanar_orientation =
-        k.coplanar_orientation_3_object();
-
-      if( coplanar_orientation(res[0], res[1], res[2])
-          != coplanar_orientation(res[0], res[1], res[3]))
-      {
-        std::swap(res[1], res[2]);
     }
 
-      return result_type(std::forward<Poly>(res));
+    switch (ids.size())
+    {
+      case 4:
+      {
+        std::vector<Point_3> res({ corners[ids[0]],
+                                   corners[ids[1]],
+                                   corners[ids[2]],
+                                   corners[ids[3]] });
+        return result_type(res);
       }
+      case 2:
+      {
+        if (start_id == -1) start_id = ids[0];
+        add_neighbor(ids[0], ids[1]);
+        add_neighbor(ids[1], ids[0]);
         break;
-  case 4: // intersects a face
-  {
-    Poly res;
-    res.reserve(4);
-    typename K::Has_on_3 has_on;
-    if(has_on(pl, cub.vertex(0))
-       && has_on(pl, cub.vertex(5))
-       && has_on(pl, cub.vertex(4)))
-    {
-      res.push_back(cub.vertex(0));
-      res.push_back(cub.vertex(5));
-      res.push_back(cub.vertex(4));
-      res.push_back(cub.vertex(3));
       }
-    else if(has_on(pl, cub.vertex(0))
+      case 1:
-            && has_on(pl, cub.vertex(1))
+        solo_id = ids[0];
-            && has_on(pl, cub.vertex(6)))
-    {
-      res.push_back(cub.vertex(0));
-      res.push_back(cub.vertex(1));
-      res.push_back(cub.vertex(6));
-      res.push_back(cub.vertex(5));
-
-    }
-    else if(has_on(pl, cub.vertex(1))
-            && has_on(pl, cub.vertex(2))
-            && has_on(pl, cub.vertex(7)))
-    {
-      res.push_back(cub.vertex(1));
-      res.push_back(cub.vertex(2));
-      res.push_back(cub.vertex(7));
-      res.push_back(cub.vertex(6));
-    }
-    else if(has_on(pl, cub.vertex(2))
-            && has_on(pl, cub.vertex(3))
-            && has_on(pl, cub.vertex(4)))
-    {
-      res.push_back(cub.vertex(2));
-      res.push_back(cub.vertex(3));
-      res.push_back(cub.vertex(4));
-      res.push_back(cub.vertex(7));
-    }
-    else if(has_on(pl, cub.vertex(6))
-            && has_on(pl, cub.vertex(7))
-            && has_on(pl, cub.vertex(4)))
-    {
-      res.push_back(cub.vertex(6));
-      res.push_back(cub.vertex(7));
-      res.push_back(cub.vertex(4));
-      res.push_back(cub.vertex(5));
-    }
-    else if(has_on(pl, cub.vertex(2))
-            && has_on(pl, cub.vertex(1))
-            && has_on(pl, cub.vertex(0)))
-    {
-      res.push_back(cub.vertex(2));
-      res.push_back(cub.vertex(1));
-      res.push_back(cub.vertex(0));
-      res.push_back(cub.vertex(3));
-    }
-    return result_type(std::forward<Poly>(res));
-  }
-    break;
       default:
         break;
     }
-
-  Poly filtered_points;
-  filter_points(points, filtered_points);
-
-  //adjacent to a vertex
-  if(filtered_points.size() == 1)
-  {
-    return result_type(std::forward<Point_3>(filtered_points.front()));
   }
 
-  //get intersections between pl and each face -> line. Foreach line, creates segment with points. Then use helper_function to recover polygon.
+  if (all_in || all_out) return boost::none;
-  typedef typename Intersection_traits<K,
+  if (start_id == -1) return { result_type(corners[solo_id]) };
-      CGAL::Plane_3<K>,
-      CGAL::Plane_3<K> >::result_type Pl_pl_type;
 
-  std::vector<Line_3> plane_intersections;
+  int prv_id = -1;
-  Pl_pl_type pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(0),
+  int cur_id = start_id;
-                                                       cub.vertex(1),
+  std::vector<Point_3> res;
-                                                       cub.vertex(5)));
+  res.reserve(6);
-  if(const Line_3* line = boost::get<Line_3>(&*pl_inter)){
+  do {
-    plane_intersections.push_back(*line);
+    res.push_back(corners[cur_id]);
+    int nxt_id = neighbor_ids[cur_id][0] == prv_id
+               ? neighbor_ids[cur_id][1]
+               : neighbor_ids[cur_id][0];
+    if (nxt_id == -1 || nxt_id == start_id){
+      if(res.size() == 2){
+        typename K::Segment_3 seg(res[0], res[1]);
+        return result_type(seg);
       }
-  pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(0),
-                                            cub.vertex(3),
-                                            cub.vertex(4)));
-  if(const Line_3* line = boost::get<Line_3>(&*pl_inter)){
-    plane_intersections.push_back(*line);
-  }
-  pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(0),
-                                            cub.vertex(1),
-                                            cub.vertex(3)));
-  if(const Line_3* line = boost::get<Line_3>(&*pl_inter)){
-    plane_intersections.push_back(*line);
-  }
-  pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(7),
-                                            cub.vertex(6),
-                                            cub.vertex(1)));
-  if(const Line_3* line = boost::get<Line_3>(&*pl_inter)){
-    plane_intersections.push_back(*line);
-  }
-  pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(7),
-                                            cub.vertex(4),
-                                            cub.vertex(3)));
-  if(const Line_3* line = boost::get<Line_3>(&*pl_inter)){
-    plane_intersections.push_back(*line);
-  }
-  pl_inter = CGAL::intersection(pl, Plane_3(cub.vertex(7),
-                                            cub.vertex(6),
-                                            cub.vertex(4)));
-  if(const Line_3* line = boost::get<Line_3>(&*pl_inter)){
-    plane_intersections.push_back(*line);
-  }
-
-  std::list<Segment_3> tmp_segs;
-  for(const auto& line : plane_intersections)
-  {
-    bool first_found = false;
-    Point_3 first_p;
-    typename K::Has_on_3 has_on;
-    for(const auto  &p : filtered_points)
-    {
-      if(has_on(line, p))
-      {
-        if(!first_found)
-        {
-          first_found = true;
-          first_p = p;
-        }
-        else
-        {
-          tmp_segs.emplace_back(first_p, p);
-          break;
-        }
-      }
-    }
-  }
-
-  if(tmp_segs.size() < 3)
-    return result_type();
-
-  std::list<Point_3> tmp_pts;
-  fill_points_list(tmp_segs,tmp_pts);
-
-  Poly res;
-  for(const auto& p : tmp_pts)
-    res.push_back(p);
-
       if(res.size() == 3){
         typename K::Triangle_3 tr(res[0], res[1], res[2]);
-    return result_type(std::forward<typename K::Triangle_3>(tr));
+        return result_type(tr);
       }
-  else
+      return result_type(res);
-  {
-    return result_type(std::forward<Poly>(res));
     }
+    prv_id = cur_id;
+    cur_id = nxt_id;
+  } while (true);
 }
 
 template <class K>

Intersections_3/include/CGAL/Intersections_3/internal/Tetrahedron_3_Plane_3_intersection.h

@@ -32,101 +32,170 @@ template <class K>
 typename Intersection_traits<K, typename K::Tetrahedron_3, typename K::Plane_3>::result_type
 intersection(
     const typename K::Tetrahedron_3 &tet,
-    const typename K::Plane_3 &pl,
+    const typename K::Plane_3 &plane,
-    const K&)
+    const K& k)
 {
-  typedef typename Intersection_traits<K,
+ typedef typename K::Point_3 Point_3;
-      typename K::Tetrahedron_3,
+  typedef typename K::Triangle_3 Triangle_3;
-      typename K::Plane_3>::result_type result_type;
+ typedef typename Intersection_traits<K, typename K::Tetrahedron_3, typename K::Plane_3>::result_type result_type;
+  typename K::Oriented_side_3 oriented_side = k.oriented_side_3_object();
 
-  typedef typename Intersection_traits<K,
+  std::vector<Point_3> corners(4);
-      typename K::Triangle_3,
+  corners.reserve(8); // 4 corners + up to 4 polygon points
-      typename K::Plane_3>::result_type Inter_type;
+  corners[0] = tet[0];
+  corners[1] = tet[1];
+  corners[2] = tet[2];
+  corners[3] = tet[3];
 
-  typedef typename K::Segment_3 Segment_3;
+  const std::array<CGAL::Oriented_side, 4> orientations  { {
-  Inter_type intersections[4];
+      oriented_side(plane, corners[0]),
-  int p_id = -1;
+      oriented_side(plane, corners[1]),
-  std::vector<typename K::Point_3> points;
+      oriented_side(plane, corners[2]),
-  std::vector<Segment_3> segments;
+      oriented_side(plane, corners[3])
+    } };
+
+  // description of faces of the bbox
+  constexpr std::array<int, 12> face_indices
+    { { 0, 1, 2,
+        0, 1, 3,
+        1, 2, 3,
+        2, 0, 3 } };
+
+  constexpr std::array<int, 12> edge_indices
+    { { 0,  1,  2,
+        0,  3,  5,
+        1,  4,  3,
+        2,  5,  4 } };
+
+  std::array<int, 12> edge_ipt_id;
+  edge_ipt_id.fill(-1);
+
+  auto inter_pt_index =
+    [&plane, &corners, &edge_ipt_id](int i, int j, int edge_id)
+  {
+    if (edge_ipt_id[edge_id]==-1)
+    {
+      edge_ipt_id[edge_id] = static_cast<int> (corners.size());
+      corners.push_back(typename K::Construct_plane_line_intersection_point_3()
+                      (plane, corners[i], corners[j]));
+    }
+
+    return edge_ipt_id[edge_id];
+  };
+
+  bool all_in = true;
+  bool all_out = true;
+
+  std::vector<std::array<int,2>> neighbor_ids(8, {-1,-1});
+
+  auto add_neighbor = [&neighbor_ids](int i, int j)
+  {
+    if (neighbor_ids[i][0] == -1 ) {
+      neighbor_ids[i][0] = j;
+    }
+    else {
+      if (neighbor_ids[i][0]!=j && neighbor_ids[i][1]==-1)
+      {
+        neighbor_ids[i][1] = j;
+      }
+    }
+  };
+
+  int start_id = -1;
+  int solo_id = -1;
+  // for each face of the bbox, we look for intersection of the plane with its edges
+  std::vector<int> ids;
   for (int i = 0; i < 4; ++i)
   {
-    const typename K::Triangle_3 triangle(tet.vertex((i+1)%4),
+    ids.clear();
-                                          tet.vertex((i+2)%4),
+    for (int k = 0; k < 3; ++k)
-                                          tet.vertex((i+3)%4));
-    intersections[i] = typename K::Intersect_3()(pl, triangle);
-    if(intersections[i]){
-      if(const typename K::Triangle_3* tr = boost::get<typename K::Triangle_3>(&*intersections[i]))
     {
-        typename K::Triangle_3 res = *tr;
-        return result_type(std::forward<typename K::Triangle_3>(res));
-      }
-      else if( const Segment_3* s
-               = boost::get<Segment_3>(&*intersections[i]))
-      {
-        segments.push_back(*s);
-      }
-      else if( const typename K::Point_3* p
-               = boost::get<typename K::Point_3>(&*intersections[i]))
-      {
-        points.push_back(*p);
-        p_id = i;
-      }
-    }
-  }
-  CGAL_assertion(segments.size() != 1);
 
-  switch(segments.size())
+      int current_id = face_indices[3 * i + k];
-  {
+      int next_id = face_indices[3 * i + (k + 1) % 3];
-  case 0:
+      int edge_id = edge_indices[3 * i + k];
-  {
-    if(p_id == -1)
-      return result_type();
-    else
-    {
-      typename K::Point_3 p
-          = *boost::get<typename K::Point_3>(&*intersections[p_id]);
 
-      return result_type(std::forward<typename  K::Point_3>(p));
+      switch (orientations[current_id])
-    }
+      {
+        case ON_NEGATIVE_SIDE:
+        {
+          all_out = false;
+          // check for intersection of the edge
+          if (orientations[next_id] == ON_POSITIVE_SIDE)
+          {
+            ids.push_back(
+                          inter_pt_index(current_id, next_id, edge_id));
           }
           break;
+        }
+        case ON_POSITIVE_SIDE:
+        {
+          all_in = false;
+          // check for intersection of the edge
+          if (orientations[next_id] == ON_NEGATIVE_SIDE)
+          {
+            ids.push_back(inter_pt_index(current_id, next_id, edge_id));
+          }
+          break;
+        }
+        case ON_ORIENTED_BOUNDARY:
+        {
+          all_in = all_out = false;
+          ids.push_back(current_id);
+        }
+      }
+    }
+
+    switch (ids.size())
+    {
+      case 3:
+      {
+        Triangle_3 res(corners[ids[0]],
+                       corners[ids[1]],
+                       corners[ids[2]]);
+        return result_type(res);
+      }
       case 2:
       {
-    return result_type(std::forward<typename  K::Segment_3>(segments.back()));
+        if (start_id == -1) start_id = ids[0];
-  }
+        add_neighbor(ids[0], ids[1]);
+        add_neighbor(ids[1], ids[0]);
         break;
+      }
+      case 1:
+        solo_id = ids[0];
       default:
-  {
-    std::set<typename K::Point_3> all_points;
-    for (const auto& s : segments)
-    {
-      all_points.insert(s.source());
-      all_points.insert(s.target());
-    }
-    if(all_points.size() == 3)
-    {
-      auto p_it = all_points.begin();
-      ++p_it;
-      typename K::Point_3 mid_point = *p_it;
-      ++p_it;
-      typename K::Triangle_3 result(*all_points.begin(), mid_point, *p_it );
-      return result_type(std::forward<typename K::Triangle_3>(result));
-    }
-    else //size = 4
-    {
-      std::list<Segment_3> segs(segments.begin(), segments.end());
-      std::list<typename K::Point_3> tmp;
-      fill_points_list(segs, tmp);
-      std::vector<typename K::Point_3> res;
-      for( const auto& p : tmp)
-        res.push_back(p);
-      return result_type(std::forward<std::vector<typename K::Point_3> >(res));
-    }
-  }
         break;
     }
-  CGAL_assertion(false);
+  }
-  return result_type();
+
+  if (all_in || all_out) return boost::none;
+  if (start_id == -1) return { result_type(corners[solo_id]) };
+
+  int prv_id = -1;
+  int cur_id = start_id;
+  std::vector<Point_3> res;
+  res.reserve(4);
+  do {
+    res.push_back(corners[cur_id]);
+    int nxt_id = neighbor_ids[cur_id][0] == prv_id
+               ? neighbor_ids[cur_id][1]
+               : neighbor_ids[cur_id][0];
+    if (nxt_id == -1 || nxt_id == start_id){
+      if(res.size() == 2){
+        typename K::Segment_3 seg(res[0], res[1]);
+        return result_type(seg);
+      }
+      if(res.size() == 3){
+        typename K::Triangle_3 tr(res[0], res[1], res[2]);
+        return result_type(tr);
+      }
+      return result_type(res);
+    }
+    prv_id = cur_id;
+    cur_id = nxt_id;
+  } while (true);
+
 }
 
 template <class K>

Intersections_3/test/Intersections_3/issue_5428.cpp

@@ -0,0 +1,24 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK;
+
+typedef EPICK::Point_3 Point;
+
+int main()
+{
+  EPICK::Plane_3 pl(Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0));
+
+  std::vector<Point> pts;
+  pts.push_back(Point(-10,-10,-10));
+  pts.push_back(Point(10,10,10));
+  CGAL::Bbox_3 cub = CGAL::bbox_3(pts.begin(), pts.end());
+
+  auto result = intersection(cub,pl);
+
+  const std::vector<Point>* res = boost::get <std::vector<Point> >(&*result);
+  for (Point p : *res) {
+    std::cout << p << std::endl;
+  }
+
+  return 0;
+}

Intersections_3/test/Intersections_3/test_intersections_3.cpp

@@ -672,7 +672,8 @@ struct Test {
                           P(0,1,0));
 
       check_intersection (tet, Pl(P(0,0.5,0), P(1,0.5,-5), P(0.5,0.5,0.5)),
-                          Tr(P(0.5,0.5,0), P(0,0.5,0), P(0,0.5,0.5)));
+                          Tr(P(0, 0.5, 0), P(0.5,0.5,0),  P(0,0.5,0.5)));
+
       Pl pl(P(0,0.9,0), P(0.9,0,0), P(0.9,0.01,0.06));
 
       typedef typename CGAL::Intersection_traits<K,
@@ -907,7 +908,7 @@ struct Test {
 
       //edge
       check_intersection (cub, Pl(P(1,1,1), P(1,2,1), P(1.5,0,0)),
-                          S(P(1,2,1), P(1,1,1)));
+                          S(P(1,1,1), P(1,2,1)));
 
 
       //face
@@ -924,6 +925,7 @@ struct Test {
         assert(p.x() == 1);
       }
       res = CGAL::intersection(cub, Pl(P(1,1,1), P(1,2,1), P(2,2,2)));
+
       poly = boost::get<std::vector<P> >(&*res);
       assert(poly != nullptr);
       assert(poly->size() == 4);
@@ -940,8 +942,8 @@ struct Test {
       check_intersection (cub, Pl(P(2, 1.66, 2),
                                   P(1.66,2,2),
                                   P(2,2,1.66)),
-                          Tr(P(2, 2, 1.66),
+                          Tr(P(1.66,2,2),
-                             P(1.66,2,2),
+                             P(2, 2, 1.66),
                              P(2,1.66,2)));
 
       //other edge
@@ -1100,7 +1102,7 @@ struct Test {
 
       //edge
       check_intersection (cub, Pl(P(1,1,1), P(1,2,1), P(1.5,0,0)),
-                          S(P(1,2,1), P(1,1,1)));
+                          S(P(1,1,1), P(1,2,1)));
       //face
       typedef typename CGAL::Intersection_traits<K,
           typename K::Plane_3,
@@ -1126,8 +1128,8 @@ struct Test {
       check_intersection (cub, Pl(P(2, 1.66, 2),
                                   P(1.66,2,2),
                                   P(2,2,1.66)),
-                          Tr(P(2, 2, 1.66),
+                          Tr(P(1.66,2,2),
-                             P(1.66,2,2),
+                             P(2, 2, 1.66),
                              P(2,1.66,2)));
 
       //random

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h

@@ -41,6 +41,7 @@
 #include <unordered_map>
 #include <utility>
 #include <vector>
+#include <bitset>
 
 namespace CGAL{
 namespace Polygon_mesh_processing {
@@ -48,24 +49,6 @@ namespace Polygon_mesh_processing {
 namespace internal
 {
 
-template <class Geom_traits, class Plane_3, class Point_3>
-int
-inter_pt_index(int i, int j,
-               const Plane_3& plane,
-               std::vector<Point_3>& points,
-               std::map<std::pair<int,int>, int>& id_map)
-{
-  std::pair<std::map<std::pair<int,int>, int>::iterator, bool> res =
-    id_map.insert(std::make_pair(make_sorted_pair(i,j),
-                  static_cast<int> (points.size())));
-  if(res.second)
-    points.push_back(
-      typename Geom_traits::Construct_plane_line_intersection_point_3()
-        (plane, points[i], points[j]));
-
-  return res.first->second;
-}
-
 template <class Plane_3,
           class TriangleMesh,
           class NamedParameters>
@@ -106,7 +89,7 @@ clip_to_bbox(const Plane_3& plane,
   }};
 
   // description of faces of the bbox
-  std::array<int, 24> face_indices =
+  constexpr std::array<int, 24> face_indices
     { { 0, 1, 2, 3,
         2, 1, 5, 6,
         3, 2, 6, 7,
@@ -114,11 +97,34 @@ clip_to_bbox(const Plane_3& plane,
         4, 0, 3, 7,
         6, 5, 4, 7 } };
 
-  std::map<std::pair<int,int>, int> id_map;
+  constexpr std::array<int, 24> edge_indices
+    { { 0,  1,  2, 3,
+        1,  4,  5, 6,
+        2,  6,  7, 8,
+        0,  9, 10, 4,
+        9,  3,  8, 11,
+        5, 10, 11, 7 } };
+
+  std::array<int, 12> edge_ipt_id;
+  edge_ipt_id.fill(-1);
+
+  auto inter_pt_index =
+    [&plane, &corners, &edge_ipt_id](int i, int j, int edge_id)
+  {
+    if (edge_ipt_id[edge_id]==-1)
+    {
+      edge_ipt_id[edge_id] = static_cast<int> (corners.size());
+      corners.push_back(typename Geom_traits::Construct_plane_line_intersection_point_3()
+                      (plane, corners[i], corners[j]));
+    }
+
+    return edge_ipt_id[edge_id];
+  };
+
   std::vector< std::vector<int> > output_faces(6);
   bool all_in = true;
   bool all_out = true;
-  std::set<int> in_point_ids; // to collect the set of points in the clipped bbox
+  std::bitset<14> in_point_bits; // to collect the set of points in the clipped bbox
 
   // for each face of the bbox, we look for intersection of the plane with its edges
   for(int i=0; i<6; ++i)
@@ -127,6 +133,7 @@ clip_to_bbox(const Plane_3& plane,
     {
       int current_id = face_indices[4*i + k];
       int next_id = face_indices[4*i + (k+1)%4];
+      int edge_id = edge_indices[4 * i + k];
 
       switch(orientations[ current_id ])
       {
@@ -135,13 +142,13 @@ clip_to_bbox(const Plane_3& plane,
           all_out=false;
           // point on or on the negative side
           output_faces[i].push_back(current_id);
-          in_point_ids.insert(output_faces[i].back());
+          in_point_bits.set(output_faces[i].back());
           // check for intersection of the edge
           if(orientations[ next_id ] == ON_POSITIVE_SIDE)
           {
             output_faces[i].push_back(
-              inter_pt_index<Geom_traits>(current_id, next_id, plane, corners, id_map));
+              inter_pt_index(current_id, next_id, edge_id));
-            in_point_ids.insert(output_faces[i].back());
+            in_point_bits.set(output_faces[i].back());
           }
           break;
         }
@@ -152,15 +159,15 @@ clip_to_bbox(const Plane_3& plane,
           if(orientations[ next_id ] == ON_NEGATIVE_SIDE)
           {
             output_faces[i].push_back(
-              inter_pt_index<Geom_traits>(current_id, next_id, plane, corners, id_map));
+              inter_pt_index(current_id, next_id, edge_id));
-            in_point_ids.insert(output_faces[i].back());
+            in_point_bits.set(output_faces[i].back());
           }
           break;
         }
         case ON_ORIENTED_BOUNDARY:
         {
           output_faces[i].push_back(current_id);
-          in_point_ids.insert(output_faces[i].back());
+          in_point_bits.set(output_faces[i].back());
         }
       }
     }
@@ -182,12 +189,15 @@ clip_to_bbox(const Plane_3& plane,
   typedef typename graph_traits::face_descriptor face_descriptor;
 
   std::map<int, vertex_descriptor> out_vertices;
-  for(int i : in_point_ids)
+  for(int i=0; i<14;++i)
+  {
+    if (in_point_bits.test(i))
     {
       vertex_descriptor v = add_vertex(tm_out);
       out_vertices.insert(std::make_pair(i, v));
       put(vpm_out, v, corners[i]);
     }
+  }
 
   std::map< std::pair<int,int>, halfedge_descriptor> hedge_map;
   const halfedge_descriptor null_hedge = graph_traits::null_halfedge();