Build K3_tree from seperate vertex, edge, and facet lists

Nef_3/include/CGAL/Nef_3/K3_tree.h

@@ -41,8 +41,7 @@ template <class Traits>
 class K3_tree
 {
 
-template <typename Kernel, typename Object,
+template <typename Kernel, typename Vertex, typename Coordinate>
-          typename Vertex, typename Coordinate>
 class Smaller_than
 {
 public:
@@ -59,26 +58,13 @@ public:
     return false;
   }
 
-  bool operator()( const Object& o1, const Object& o2) {
-    Vertex v1,v2;
-    CGAL::assign(v1,o1);
-    CGAL::assign(v2,o2);
-    switch(coord) {
-    case 0: return CGAL::compare_x(v1->point(), v2->point()) == SMALLER;
-    case 1: return CGAL::compare_y(v1->point(), v2->point()) == SMALLER;
-    case 2: return CGAL::compare_z(v1->point(), v2->point()) == SMALLER;
-    default: CGAL_error();
-    }
-    return false;
-  }
 private:
   Coordinate coord;
 };
 
 
-template <typename Object, typename Vertex,
+template <typename Vertex, typename Coordinate, typename EK>
-          typename Coordinate, typename EK>
+  class Smaller_than<CGAL::Lazy_kernel<EK>, Vertex, Coordinate>
-  class Smaller_than<CGAL::Lazy_kernel<EK>, Object, Vertex, Coordinate>
 {
 public:
   Smaller_than(Coordinate c) : coord(c) {
@@ -97,21 +83,6 @@ public:
     return false;
   }
 
-  bool operator()( const Object& o1, const Object& o2) {
-    Vertex v1,v2;
-    CGAL::assign(v1,o1);
-    CGAL::assign(v2,o2);
-    switch(coord) {
-    case 0: return CGAL::to_interval(v1->point().x()).second <
-                           CGAL::to_interval(v2->point().x()).first;
-    case 1: return CGAL::to_interval(v1->point().y()).second <
-                           CGAL::to_interval(v2->point().y()).first;
-    case 2: return CGAL::to_interval(v1->point().z()).second <
-                           CGAL::to_interval(v2->point().z()).first;
-    default: CGAL_error();
-    }
-    return false;
-  }
 private:
   Coordinate coord;
 };
@@ -126,13 +97,17 @@ public:
 typedef typename Traits::SNC_decorator SNC_decorator;
 typedef typename Traits::Infimaximal_box Infimaximal_box;
 typedef typename Traits::Vertex_handle Vertex_handle;
+typedef typename Traits::Vertex_list Vertex_list;
+typedef typename Vertex_list::iterator Vertex_iterator;
+typedef typename Vertex_list::const_iterator Vertex_const_iterator;
 typedef typename Traits::Halfedge_handle Halfedge_handle;
 typedef typename Traits::Halffacet_handle Halffacet_handle;
 typedef typename Traits::Object_handle Object_handle;
-typedef std::vector<Object_handle> Object_list;
+typedef typename Traits::Object_list Object_list;
-typedef typename Object_list::const_iterator Object_const_iterator;
+typedef typename Traits::Halfedge_list Halfedge_list;
-typedef typename Object_list::iterator Object_iterator;
+typedef typename Halfedge_list::const_iterator Halfedge_const_iterator;
-typedef typename Object_list::size_type size_type;
+typedef typename Traits::Halffacet_list Halffacet_list;
+typedef typename Halffacet_list::const_iterator Halffacet_const_iterator;
 
 typedef typename Traits::Point_3 Point_3;
 typedef typename Traits::Segment_3 Segment_3;
@@ -144,7 +119,6 @@ typedef typename Traits::Aff_transformation_3 Aff_transformation_3;
 
 typedef typename Traits::Bounding_box_3 Bounding_box_3;
 typedef typename Traits::Side_of_plane Side_of_plane;
-typedef typename Traits::Objects_bbox Objects_bbox;
 
 typedef typename Traits::Kernel Kernel;
 typedef typename Kernel::RT RT;
@@ -152,29 +126,36 @@ typedef typename Kernel::FT FT;
 
 typedef Smaller_than<
   Kernel,
-  Object_handle,
   Vertex_handle,
-  int> Smaller_;
+  int> Smaller;
 
   class Node  {
   friend class K3_tree<Traits>;
 public:
     typedef Node* Node_handle;
-
+  Node(const Vertex_list& V, const Halfedge_list& E, const Halffacet_list& F) :
-  Node( Node_handle p, Node_handle l, Node_handle r, Plane_3 pl, const Object_list& O) :
+    left_node(nullptr), right_node(nullptr)
-    parent_node(p), left_node(l), right_node(r), splitting_plane(pl),
+  {
-        object_list(O) {
+      object_list.reserve(V.size()+E.size()+F.size());
-    if(l == nullptr)
+      for(Vertex_const_iterator vi=V.begin(); vi!=V.end(); ++vi)
-      point_on_plane = Point_3();
+          object_list.push_back(make_object(*vi));
-    else
+      for(Halfedge_const_iterator ei=E.begin(); ei!=E.end(); ++ei)
-      point_on_plane = pl.point();
+          object_list.push_back(make_object(*ei));
+      for(Halffacet_const_iterator fi=F.begin(); fi!=F.end(); ++fi)
+          object_list.push_back(make_object(*fi));
   }
+
+  Node(Node_handle l, Node_handle r, const Plane_3& pl) :
+    left_node(l), right_node(r), splitting_plane(pl)
+  {
+  }
+
   bool is_leaf() const {
     CGAL_assertion( (left_node != nullptr && right_node != nullptr) ||
                     (left_node == nullptr && right_node == nullptr));
     return (left_node == nullptr && right_node == nullptr);
   }
-  Node_handle parent() const { return parent_node; }
+
   Node_handle left() const { return left_node; }
   Node_handle right() const { return right_node; }
   const Plane_3& plane() const { return splitting_plane; }
@@ -189,45 +170,42 @@ public:
     }
   }
 
-  template<typename Depth>
+  void add_facet(Halffacet_handle f, int depth) {
-  void add_facet(Halffacet_handle f, Depth depth) {
     if(left_node == nullptr) {
       object_list.push_back(make_object(f));
       return;
     }
 
-    Side_of_plane sop;
+    Side_of_plane sop(splitting_plane.point(), depth);
-    Oriented_side side = sop(splitting_plane.point(), f, depth);
+    Oriented_side side = sop(f);
     if( side == ON_NEGATIVE_SIDE || side == ON_ORIENTED_BOUNDARY)
       left_node->add_facet(f, depth+1);
     if( side == ON_POSITIVE_SIDE || side == ON_ORIENTED_BOUNDARY)
       right_node->add_facet(f, depth+1);
   }
 
-  template<typename Depth>
+  void add_edge(Halfedge_handle e, int depth) {
-  void add_edge(Halfedge_handle e, Depth depth) {
     if(left_node == nullptr) {
       object_list.push_back(make_object(e));
       return;
     }
 
-    Side_of_plane sop;
+    Side_of_plane sop(splitting_plane.point(), depth);
-    Oriented_side side = sop(splitting_plane.point(), e, depth);
+    Oriented_side side = sop(e);
     if( side == ON_NEGATIVE_SIDE || side == ON_ORIENTED_BOUNDARY)
       left_node->add_edge(e, depth+1);
     if( side == ON_POSITIVE_SIDE || side == ON_ORIENTED_BOUNDARY)
       right_node->add_edge(e, depth+1);
   }
 
-  template<typename Depth>
+  void add_vertex(Vertex_handle v, int depth) {
-  void add_vertex(Vertex_handle v, Depth depth) {
     if(left_node == nullptr) {
       object_list.push_back(make_object(v));
       return;
     }
 
-    Side_of_plane sop;
+    Side_of_plane sop(splitting_plane.point(), depth);
-    Oriented_side side = sop(splitting_plane.point(), v, depth);
+    Oriented_side side = sop(v);
     if( side == ON_NEGATIVE_SIDE || side == ON_ORIENTED_BOUNDARY)
       left_node->add_vertex(v, depth+1);
     if( side == ON_POSITIVE_SIDE || side == ON_ORIENTED_BOUNDARY)
@@ -252,27 +230,11 @@ friend std::ostream& operator<<
   return os;
 }
 
-    /*
-~Node() noexcept(!CGAL_ASSERTIONS_ENABLED)
-{
-  CGAL_NEF_TRACEN("~Node: deleting node...");
-  CGAL_destructor_assertion_catch(
-    if( !is_leaf()) {
-      delete left_node;
-      delete right_node;
-    }
-  );
-}
-    */
 private:
 
-
-
-  Node_handle parent_node;
   Node_handle left_node;
   Node_handle right_node;
   Plane_3 splitting_plane;
-  Point_3 point_on_plane;
   Object_list object_list;
 };
 
@@ -437,23 +399,15 @@ void divide_segment_by_plane( Segment_3 s, Plane_3 pl,
       // First of all, we need to find out wheather we are working over an extended kernel or on a standard kernel. As precondition we have that ray is oriented in the minus x axis direction.  When having an extended kernel, the ray can be subtituted by a segment with the endpoint on the 'intersection' between the ray and the bounding infimaximal box.  In the presence of a standard kernel, the intersection is computed with the bounding box with the vertices of the Nef polyhedron.
       Point_3 p(r.source()), q;
       Bounding_box_3 b = k.bounding_box;
-      int c = (CGAL::abs(vec[0]) > CGAL::abs(vec[1]) ? 0 : 1);
+      typename Kernel::Non_zero_dimension_3 non_zero_dimension_3;
-      c = (CGAL::abs(vec[2]) > CGAL::abs(vec[c]) ? 2 : c);
+      int c = non_zero_dimension_3(vec);
 
       Point_3 pt_on_minus_x_plane = vec[c] < 0 ?
         Point_3(FT(b.min_coord(0)), FT(b.min_coord(1)),FT(b.min_coord(2))) :
         Point_3(FT(b.max_coord(0)), FT(b.max_coord(1)),FT(b.max_coord(2)));
       // We compute the intersection between a plane with normal vector in
       // the minus x direction and located at the minimum point of the bounding box, and the input ray.  When the ray does not intersect the bounding volume, there won't be any object hit, so it is safe to construct a segment that simply lay in the unbounded side of the bounding box.  This approach is taken instead of somehow (efficiently) report that there was no hit object, in order to mantain a clear interface with the Iterator class.
-      Plane_3 pl_on_minus_x;
+      Plane_3 pl_on_minus_x = K3_tree::construct_splitting_plane(pt_on_minus_x_plane, c);
-      if(c==0)
-        pl_on_minus_x = Plane_3(pt_on_minus_x_plane, Vector_3( 1, 0, 0));
-      else if(c==1)
-        pl_on_minus_x = Plane_3(pt_on_minus_x_plane, Vector_3( 0, 1, 0));
-      else {
-        CGAL_assertion_msg(c==2, "wrong value");
-         pl_on_minus_x = Plane_3(pt_on_minus_x_plane, Vector_3( 0, 0, 1));
-      }
       Object o = traits.intersect_object()( pl_on_minus_x, r);
       if( !CGAL::assign( q, o) || pl_on_minus_x.has_on(p))
         q = r.source() + vec;
@@ -471,7 +425,7 @@ private:
 
 
   Node_handle root;
-  boost::container::deque<Node> nodes;
+  Node_range nodes;
 
   int max_depth;
   Bounding_box_3 bounding_box;
@@ -488,31 +442,32 @@ public:
     typedef typename SNC_structure::Halffacet_iterator Halffacet_iterator;
 
     CGAL_assertion( W != nullptr);
-    Object_list objects;
+
+    Vertex_list vertices;
+    vertices.reserve(W->number_of_vertices());
     Vertex_iterator v;
-    Halfedge_iterator e;
-    Halffacet_iterator f;
     CGAL_forall_vertices(v, *W)
-      objects.push_back(make_object(Vertex_handle(v)));
+      vertices.push_back(v);
-    typename Object_list::difference_type v_end = objects.size();
+
+    Halfedge_list edges;
+    edges.reserve(W->number_of_halfedges());
+    Halfedge_iterator e;
     CGAL_forall_edges(e, *W)
-      objects.push_back(make_object(Halfedge_handle(e)));
+      edges.push_back(e);
-    CGAL_forall_facets( f, *W) {
-      objects.push_back(make_object(Halffacet_handle(f)));
-    }
-    Object_iterator oli=objects.begin()+v_end;
-    root = build_kdtree( objects, oli, 0);
-  }
 
-  K3_tree(Object_list& objects, Object_iterator& v_end) {
+    Halffacet_list facets;
+    facets.reserve(W->number_of_halffacets());
+    Halffacet_iterator f;
+    CGAL_forall_facets(f, *W)
+      facets.push_back(f);
 
-typename Object_list::difference_type n_vertices = std::distance(objects.begin(),v_end);
+    CGAL_NEF_TRACEN("K3_tree(): n_vertices = " << vertices.size());
- CGAL_NEF_TRACEN("K3_tree(): n_vertices = " << std::distance(objects.begin(),v_end));
+    std::frexp( (double) vertices.size(), &max_depth);
- std::frexp( (double) n_vertices, &max_depth);
 
     // TODO: in the presence of a infimaximal bounding box, the bounding box does not have to be computed
- Objects_bbox objects_bbox = traits.objects_bbox_object();
+    bounding_box = Bounding_box_3();
- bounding_box = objects_bbox(objects);
+    for(typename Vertex_list::iterator vi=vertices.begin(); vi!=vertices.end(); ++vi)
+        bounding_box.extend((*vi)->point());
     //CGAL_NEF_TRACEN("bounding box:"<<objects_bbox);
 
 #ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
@@ -520,7 +475,7 @@ typename Object_list::difference_type n_vertices = std::distance(objects.begin()
     reference_counted = (&(p1.hx()) == &(p2.hx()));
     CGAL_NEF_TRACEN("reference counted " << reference_counted);
 #endif
-    root = build_kdtree( objects, v_end, 0);
+    root = build_kdtree(vertices, edges, facets, 0);
   }
   const Object_list& objects_around_point( const Point_3& p) const {
     return locate( p, root);
@@ -664,144 +619,114 @@ std::string dump_object_list( const Object_list& O, int level = 0) {
   return os.str();
  }
 
-  /*
-~K3_tree() noexcept(!CGAL_ASSERTIONS_ENABLED)
-{
-  CGAL_NEF_TRACEN("~K3_tree: deleting root...");
-  CGAL_destructor_assertion_catch(
-    delete root;
-  );
-}
-  */
 private:
 
-template <typename Depth>
+Node_handle build_kdtree(Vertex_list& V, Halfedge_list& E, Halffacet_list& F,
-Node_handle build_kdtree(Object_list& O, Object_iterator v_end,
+                         int depth, int non_efective_splits=0) {
-                   Depth depth, Node_handle parent=nullptr, int non_efective_splits=0) {
   CGAL_precondition( depth >= 0);
-  CGAL_NEF_TRACEN( "build_kdtree: "<<O.size()<<" objects, "<<"depth "<<depth);
+
-  CGAL_NEF_TRACEN( "build_kdtree: "<<dump_object_list(O,1));
+  if( !can_set_be_divided(depth, V.size())) {
-  if( !can_set_be_divided(O.begin(), v_end, depth)) {
     CGAL_NEF_TRACEN("build_kdtree: set cannot be divided");
-    nodes.push_back(Node( parent, nullptr, nullptr, Plane_3(), O));
+    nodes.push_back(Node(V, E, F));
     return &(nodes.back());
   }
-  Object_iterator median;
-  Plane_3 partition_plane = construct_splitting_plane(O.begin(), v_end, median, depth);
-  CGAL_NEF_TRACEN("build_kdtree: plane: "<<partition_plane<< " " << partition_plane.point());
 
-  Object_list O1, O2;
+  int coord = depth%3;
-  Vertex_handle vm,vx;
+  Vertex_iterator median;
-  CGAL::assign(vm,*median);
+  Point_3 point_on_plane = find_median_point(V, median, coord);
-  Smaller_ smaller(depth%3);
+  CGAL_NEF_TRACEN("build_kdtree: plane: "<<partition_plane<< " " << point_on_plane);
-  for(Object_iterator oi=O.begin();oi!=median;++oi) {
-    O1.push_back(*oi);
-    CGAL::assign(vx,*oi);
-    if(!smaller(vx, vm))
-      O2.push_back(*oi);
-  }
-
-  O1.push_back(*median);
-  O2.push_back(*median);
-
-  for(Object_iterator oi=median+1;oi!=v_end;++oi) {
-    O2.push_back(*oi);
-    CGAL::assign(vx,*oi);
-    if(!smaller(vm, vx))
-      O1.push_back(*oi);
-  }
 
 #ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
-  Side_of_plane sop(reference_counted);
+  Side_of_plane sop(point_on_plane, coord, reference_counted);
 #else
-  Side_of_plane sop;
+  Side_of_plane sop(point_on_plane, coord);
 #endif
-  typename Object_list::size_type v_end1 = O1.size();
-  typename Object_list::size_type v_end2 = O2.size();
-  bool splitted = classify_objects( v_end, O.end(), partition_plane, sop,
-                                    std::back_inserter(O1),
-                                    std::back_inserter(O2), depth);
 
-  bool non_efective_split = false;
+  Vertex_list V1,V2;
+  V1.reserve(std::distance(V.begin(), median));
+  V2.reserve(std::distance(median, V.end()));
+  classify_objects(V, sop, V1, V2);
+
+  Halfedge_list E1,E2;
+  bool splitted = classify_objects(E, sop, E1, E2);
+
+  Halffacet_list F1,F2;
+  splitted = splitted && classify_objects(F, sop, F1, F2);
+
   if(!splitted) {
     CGAL_NEF_TRACEN("build_kdtree: splitting plane not found");
-    //    if(depth > max_depth)
+    nodes.push_back(Node(V, E, F));
-    nodes.push_back(Node( parent, nullptr, nullptr, Plane_3(), O));
     return &(nodes.back());
-  } else {
-    CGAL_NEF_TRACEN("Sizes " << O1.size() << ", " << O2.size() << ", " << O.size());
-    CGAL_assertion( O1.size() <= O.size() && O2.size() <= O.size());
-    CGAL_assertion( O1.size() + O2.size() >= O.size());
-    non_efective_split = ((O1.size() == O.size()) || (O2.size() == O.size()));
   }
-  if( non_efective_split)
+  auto O_size = V.size() + E.size() + F.size();
+  auto O1_size = V1.size() + E1.size() + F1.size();
+  auto O2_size = V2.size() + E2.size() + F2.size();
+  CGAL_NEF_TRACEN("Sizes " << O1_size << ", " << O2_size << ", " << O_size);
+  CGAL_assertion( O1_size <= O_size && O2_size <= O_size);
+  CGAL_assertion( O1_size + O2_size >= O_size);
+  if((O1_size == O_size) || (O2_size == O_size))
     non_efective_splits++;
   else
     non_efective_splits = 0;
+
   if(non_efective_splits > 2) {
     CGAL_NEF_TRACEN("build_kdtree: non efective splits reached maximum");
-    nodes.push_back(Node( parent, nullptr, nullptr, Plane_3(), O));
+    nodes.push_back(Node(V, E, F));
     return &(nodes.back());
   }
-  nodes.push_back(Node( parent, nullptr, nullptr, partition_plane, Object_list()));
+
-  Node_handle node = &(nodes.back());
+  Node_handle left_node = build_kdtree(V1, E1, F1, depth + 1, non_efective_splits);
-  node->left_node = build_kdtree( O1, O1.begin()+v_end1, depth + 1, node, non_efective_splits);
+  Node_handle right_node = build_kdtree(V2, E2, F2, depth + 1, non_efective_splits);
-  node->right_node = build_kdtree( O2, O2.begin()+v_end2, depth + 1, node, non_efective_splits);
+  nodes.push_back(Node(left_node, right_node, construct_splitting_plane(point_on_plane, coord)));
-  return node;
+  return &(nodes.back());
 }
 
-template <typename Depth>
+bool can_set_be_divided(int depth, typename Vertex_list::size_type size) {
-bool can_set_be_divided(Object_iterator start, Object_iterator end, Depth depth) {
   if(depth >= max_depth)
     return false;
-  if(std::distance(start,end)<2)
+  if(size <= 2)
     return false;
   return true;
 }
 
-template <typename OutputIterator, typename Depth>
+template <typename List>
-bool classify_objects(Object_iterator start, Object_iterator end,
+static bool classify_objects(const List& l,Side_of_plane& sop,
-                      Plane_3 partition_plane, Side_of_plane& sop,
+                      List& l1, List& l2) {
-                      OutputIterator o1, OutputIterator o2, Depth depth) {
+  typename List::size_type on_oriented_boundary = 0;
-  typename Object_list::difference_type on_oriented_boundary = 0;
+  for(typename List::const_iterator i=l.begin(); i!=l.end(); ++i) {
-  typename Object_list::const_iterator o;
+    Oriented_side side = sop(*i);
-
+    switch(side) {
-  Point_3 point_on_plane(partition_plane.point());
+      case ON_NEGATIVE_SIDE:
-
+        l1.push_back(*i);
-  for( o = start; o != end; ++o) {
+        break;
-    Oriented_side side = sop( point_on_plane, *o, depth);
+      case ON_POSITIVE_SIDE:
-    if( side == ON_NEGATIVE_SIDE || side == ON_ORIENTED_BOUNDARY) {
+        l2.push_back(*i);
-      *o1 = *o;
+        break;
-      ++o1;
+      case ON_ORIENTED_BOUNDARY:
-    }
-    if( side == ON_POSITIVE_SIDE || side == ON_ORIENTED_BOUNDARY) {
-      *o2 = *o;
-      ++o2;
-    }
-    if( side == ON_ORIENTED_BOUNDARY)
         ++on_oriented_boundary;
+        l1.push_back(*i);
+        l2.push_back(*i);
     }
-  return (on_oriented_boundary != std::distance(start,end));
+  }
+  return (on_oriented_boundary != l.size());
 }
 
+static Point_3 find_median_point(Vertex_list& V,
+                          Vertex_iterator& median, int coord) {
+  CGAL_assertion(V.size() > 1);
 
-template <typename Depth>
+  median = V.begin() + V.size()/2;
-Plane_3 construct_splitting_plane(Object_iterator start, Object_iterator end,
+  std::nth_element(V.begin(), median, V.end(), Smaller(coord));
-                                  Object_iterator& median, Depth depth) {
-  CGAL_precondition( depth >= 0);
-  typename Object_list::difference_type n=std::distance(start,end);
-  CGAL_assertion(n>1);
 
-  std::nth_element(start, start+n/2, end,
+  return (*median)->point();
-                     Smaller_(depth%3));
+}
 
-  Vertex_handle v;
+static Plane_3 construct_splitting_plane(const Point_3& pt, int coord)
-  median = start+n/2;
+{
-  CGAL::assign(v,*median);
+  switch(coord) {
-  switch( depth % 3) {
+  case 0: return Plane_3(1, 0, 0, -pt.x());
-  case 0: return Plane_3( v->point(), Vector_3( 1, 0, 0)); break;
+  case 1: return Plane_3(0, 1, 0, -pt.y());
-  case 1: return Plane_3( v->point(), Vector_3( 0, 1, 0)); break;
+  case 2: return Plane_3(0, 0, 1, -pt.z());
-  case 2: return Plane_3( v->point(), Vector_3( 0, 0, 1)); break;
   }
 
   CGAL_error_msg( "never reached");

Nef_3/include/CGAL/Nef_3/SNC_k3_tree_traits.h

@@ -18,6 +18,7 @@
 
 #include <CGAL/Nef_3/Bounding_box_3.h>
 #include <CGAL/Lazy_kernel.h>
+#include <CGAL/Unique_hash_map.h>
 #include <list>
 
 #undef CGAL_NEF_DEBUG
@@ -28,11 +29,11 @@ namespace CGAL {
 
 
 template <typename Kernel, typename Coordinate>
-class ComparePoints {
+class Compare_points {
 
   typedef typename Kernel::Point_3  Point_3;
  public:
-  ComparePoints(Coordinate c) : coord(c) {
+  Compare_points(Coordinate c) : coord(c) {
     CGAL_assertion( c >= 0 && c <=2);
   }
   CGAL::Comparison_result operator()(const Point_3& p1, const Point_3& p2) {
@@ -56,15 +57,15 @@ private:
 
 
 template <typename Coordinate, typename EK>
-class ComparePoints<CGAL::Lazy_kernel<EK>, Coordinate> {
+class Compare_points<CGAL::Lazy_kernel<EK>, Coordinate> {
 
   typedef CGAL::Lazy_kernel<EK>     Kernel;
   typedef typename Kernel::Point_3  Point_3;
  public:
-  ComparePoints(Coordinate c) : coord(c) {
+  Compare_points(Coordinate c) : coord(c) {
     CGAL_assertion( c >= 0 && c <=2);
   }
-  CGAL::Comparison_result operator()( const Point_3 p1, const Point_3 p2) {
+  CGAL::Comparison_result operator()( const Point_3& p1, const Point_3& p2) {
     switch(coord) {
     case 0:
       if(CGAL::to_interval(p1.x()).second <
@@ -100,16 +101,13 @@ private:
 
 template <class SNC_decorator>
 class Side_of_plane {
-public:
+
-  typedef typename SNC_decorator::SNC_structure SNC_structure;
   typedef typename SNC_decorator::Decorator_traits Decorator_traits;
 
   typedef typename Decorator_traits::Vertex_handle Vertex_handle;
   typedef typename Decorator_traits::Halfedge_handle Halfedge_handle;
   typedef typename Decorator_traits::Halffacet_handle Halffacet_handle;
 
-  typedef typename SNC_structure::Object_handle Object_handle;
-
   typedef typename Decorator_traits::Halffacet_cycle_iterator
     Halffacet_cycle_iterator;
   typedef typename Decorator_traits::SHalfedge_around_facet_circulator
@@ -118,154 +116,34 @@ public:
 
   typedef typename SNC_decorator::Kernel Kernel;
   typedef typename Kernel::Point_3 Point_3;
-  typedef typename Kernel::Segment_3 Segment_3;
+#ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
-  typedef typename Kernel::Plane_3 Plane_3;
-  typedef typename Kernel::Triangle_3 Triangle_3;
-  typedef typename Kernel::Vector_3 Vector_3;
   typedef typename Kernel::RT RT;
+#endif
+  typedef Compare_points<Kernel, int> Compare;
+  static constexpr Oriented_side unknown_side = static_cast<Oriented_side>(-2);
 
-  typedef ComparePoints<Kernel, int> ComparePoints_;
+public:
-
 #ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
-  Side_of_plane(bool rc = false) : reference_counted(rc) {}
+  Side_of_plane(const Point_3& p, int c, bool rc = false) : reference_counted(rc), coord(c), pop(p) {}
 #else
-    Side_of_plane() {}
+  Side_of_plane(const Point_3& p, int c) : OnSideMap(unknown_side), coord(c), pop(p) {}
 #endif
 
-
+  Oriented_side operator()(Vertex_handle v);
-  template<typename Depth> Oriented_side operator()
+  Oriented_side operator()(Halfedge_handle e);
-    ( const Point_3& pop, Object_handle o, Depth depth);
+  Oriented_side operator()(Halffacet_handle f);
-  template<typename Depth> Oriented_side operator()
+private:
-    ( const Point_3& pop, Vertex_handle v, Depth depth);
-  template<typename Depth> Oriented_side operator()
-    ( const Point_3& pop, Halfedge_handle e, Depth depth);
-  template<typename Depth> Oriented_side operator()
-    ( const Point_3& pop, Halffacet_handle f, Depth depth);
-#ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
-  bool reference_counted;
-#endif
   Unique_hash_map<Vertex_handle,Oriented_side> OnSideMap;
 #ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
   Unique_hash_map<const RT*, Oriented_side> OnSideMapRC;
+  bool reference_counted;
 #endif
+  int coord;
+  const Point_3& pop;
 };
 
 template <class SNC_decorator>
-class Objects_bbox {
+constexpr Oriented_side Side_of_plane<SNC_decorator>::unknown_side;
-public:
-  typedef typename SNC_decorator::SNC_structure SNC_structure;
-  typedef typename SNC_decorator::Decorator_traits Decorator_traits;
-
-  typedef typename Decorator_traits::Vertex_handle Vertex_handle;
-  typedef typename Decorator_traits::Halfedge_handle Halfedge_handle;
-  typedef typename Decorator_traits::Halffacet_handle Halffacet_handle;
-
-  typedef typename SNC_structure::Halffacet_triangle_handle Halffacet_triangle_handle;
-  typedef typename SNC_structure::Object_handle Object_handle;
-  typedef std::vector<Object_handle> Object_list;
-
-  typedef typename Decorator_traits::Halffacet_cycle_iterator
-    Halffacet_cycle_iterator;
-  typedef typename Decorator_traits::SHalfedge_around_facet_circulator
-    SHalfedge_around_facet_circulator;
-  typedef typename Decorator_traits::SHalfedge_handle SHalfedge_handle;
-
-  typedef typename SNC_decorator::Kernel Kernel;
-  typedef typename Kernel::Plane_3 Plane_3;
-  typedef typename Kernel::Segment_3 Segment_3;
-  typedef typename Kernel::Point_3 Point_3;
-  typedef typename Kernel::Triangle_3 Triangle_3;
-
-  typedef typename Kernel::RT RT;
-  typedef typename Kernel::FT FT;
-  typedef typename Kernel::Kernel_tag Kernel_tag;
-  typedef CGAL::Bounding_box_3<Tag_true, Kernel>
-    Bounding_box_3;
-  //  typedef CGAL::Bounding_box_3
-  //    <typename Is_extended_kernel<Kernel>::value_type, Kernel>
-  //    Bounding_box_3;
-
-  Bounding_box_3 operator()( const Object_list& O) const {
-    Bounding_box_3 b;
-    typename Object_list::const_iterator o = O.begin();
-    Vertex_handle v;
-    while(o != O.end() && !CGAL::assign(v, *o)) ++o;
-    if(o != O.end()) {
-      FT q[3];
-      q[0] = v->point().x();
-      q[1] = v->point().y();
-      q[2] = v->point().z();
-      Bounding_box_3 b(q);
-      for(++o; o != O.end(); ++o) {
-        if( CGAL::assign( v, *o)) {
-          b.extend(v->point());
-        }
-      }
-      return b;
-    }
-    FT q[3];
-    q[0] = q[1] = q[2] = 0;
-    return Bounding_box_3(q);
-  }
-
-  /*
-  Bounding_box_3 operator()(Object_handle o) const {
-    Vertex_handle v;
-    Halfedge_handle e;
-    Halffacet_handle f;
-    if( CGAL::assign( v, o))
-      return operator()(v);
-    else if( CGAL::assign( e, o))
-      return operator()(e);
-    else if( CGAL::assign( f, o))
-      return operator()(f);
-    else {
-      Halffacet_triangle_handle t;
-      if( CGAL::assign( t, o))
-        return operator()(t);
-      else
-        CGAL_error_msg( "wrong handle");
-    }
-    return Bounding_box_3(); // never reached
-  }
-
-  Bounding_box_3 operator()(Vertex_handle v) const {
-    Bounding_box_3 b;
-    b.extend(v->point());
-    return b;
-  }
-
-  Bounding_box_3 operator()(Halfedge_handle e) const {
-    Bounding_box_3 b;
-    b.extend(e->source()->point());
-    b.extend(e->twin()->source()->point());
-    return b;
-  }
-
-  Bounding_box_3 operator()(Halffacet_triangle_handle t) const {
-    Bounding_box_3 bbox;
-    Triangle_3 tr(t.get_triangle());
-    for( int i = 0; i < 3; ++i)
-      bbox.extend(tr[i]);
-    return bbox;
-  }
-
-  Bounding_box_3 operator()(Halffacet_handle f) const {
-    CGAL_assertion( f->facet_cycles_begin() !=
-                    Halffacet_cycle_iterator());
-    Halffacet_cycle_iterator fc(f->facet_cycles_begin());
-    SHalfedge_handle e;
-    CGAL_assertion(fc.is_shalfedge());
-    e = SHalfedge_handle(fc);
-    SHalfedge_around_facet_circulator sc(e), send(sc);
-    CGAL_assertion( !is_empty_range( sc, send));
-    Bounding_box_3 b;
-    CGAL_For_all( sc, send)
-      b.extend(sc->source()->source()->point());
-    return b;
-  }
-  */
-};
 
 template <class Decorator>
 class SNC_k3_tree_traits {
@@ -283,6 +161,9 @@ public:
 
   typedef typename SNC_structure::Object_handle Object_handle;
   typedef std::vector<Object_handle> Object_list;
+  typedef std::vector<Vertex_handle> Vertex_list;
+  typedef std::vector<Halfedge_handle> Halfedge_list;
+  typedef std::vector<Halffacet_handle> Halffacet_list;
 
   typedef typename Kernel::Point_3 Point_3;
   typedef typename Kernel::Segment_3 Segment_3;
@@ -294,55 +175,24 @@ public:
 
   typedef typename Kernel::RT RT;
   typedef typename Kernel::Kernel_tag Kernel_tag;
-  //  typedef CGAL::Bounding_box_3
+  typedef CGAL::Bounding_box_3<Tag_true, Kernel> Bounding_box_3;
-  //    <typename Is_extended_kernel<Kernel>::value_type, Kernel>
-  //    Bounding_box_3;
-  typedef CGAL::Bounding_box_3<Tag_true, Kernel>
-    Bounding_box_3;
-
   typedef typename Kernel::Intersect_3 Intersect;
-  typedef CGAL::Objects_bbox<SNC_decorator> Objects_bbox;
   typedef CGAL::Side_of_plane<SNC_decorator> Side_of_plane;
 
   Intersect intersect_object() const {
     return Intersect();
   }
 
-  Objects_bbox objects_bbox_object() const {
-    return Objects_bbox();
-  }
 };
 
 template <class SNC_decorator>
-template <typename Depth>
 Oriented_side
-Side_of_plane<SNC_decorator>::operator()
+Side_of_plane<SNC_decorator>::operator()(Vertex_handle v) {
-  (const Point_3& pop, Object_handle o, Depth depth) {
-  Vertex_handle v;
-  Halfedge_handle e;
-  Halffacet_handle f;
-  if( CGAL::assign( v, o))
-    return (*this)(pop, v, depth);
-  else if( CGAL::assign( e, o))
-    return (*this)(pop, e, depth);
-  else if( CGAL::assign( f, o))
-    return (*this)(pop, f, depth);
-  else
-    CGAL_error_msg( "wrong handle");
-
-  return Oriented_side(); // never reached
-}
-
-template <class SNC_decorator>
-template <typename Depth>
-Oriented_side
-Side_of_plane<SNC_decorator>::operator()
-( const Point_3& pop, Vertex_handle v, Depth depth) {
 Comparison_result cr;
 #ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
   if(reference_counted) {
     if(!OnSideMapRC.is_defined(&(v->point().hw())))
-      switch(depth%3) {
+      switch(coord) {
       case 0:
         cr = CGAL::compare_x(v->point(), pop);
         OnSideMapRC[&(v->point().hw())] = cr == LARGER ? ON_POSITIVE_SIDE :
@@ -363,13 +213,14 @@ Side_of_plane<SNC_decorator>::operator()
     return OnSideMapRC[&(v->point().hw())];
   } else {
 #endif
-    ComparePoints_ compare(depth%3);
+  Oriented_side& side = OnSideMap[v];
-    if(!OnSideMap.is_defined(v)) {
+  if(side == unknown_side) {
+    Compare compare(coord);
     cr = compare(v->point(), pop);
-      OnSideMap[v] = cr == LARGER ? ON_POSITIVE_SIDE :
+    side = cr == LARGER ? ON_POSITIVE_SIDE :
            cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
   }
-    return OnSideMap[v];
+  return side;
 #ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
   }
 #endif
@@ -383,57 +234,13 @@ Side_of_plane<SNC_decorator>::operator()
  */
 
 template <class SNC_decorator>
-template <typename Depth>
 Oriented_side
-Side_of_plane<SNC_decorator>::operator()
+Side_of_plane<SNC_decorator>::operator()(Halfedge_handle e) {
-( const Point_3& pop, Halfedge_handle e, Depth depth) {
   Vertex_handle v = e->source();
   Vertex_handle vt = e->twin()->source();
-  /*
+
-  Comparison_result cr;
+  Oriented_side src_side = (*this) (v);
-  if(!OnSideMap.is_defined(v))
+  Oriented_side tgt_side = (*this) (vt);
-    switch(depth%3) {
-    case 0:
-      cr = CGAL::compare_x(v->point(), pop);
-      OnSideMap[v] = cr == LARGER ? ON_POSITIVE_SIDE :
-                       cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
-      break;
-    case 1:
-      cr = CGAL::compare_y(v->point(), pop);
-      OnSideMap[v] = cr == LARGER ? ON_POSITIVE_SIDE :
-                       cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
-      break;
-    case 2:
-      cr = CGAL::compare_z(v->point(), pop);
-      OnSideMap[v] = cr == LARGER ? ON_POSITIVE_SIDE :
-                       cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
-      break;
-    default: CGAL_error_msg( "wrong value");
-    }
-  if(!OnSideMap.is_defined(vt))
-    switch(depth%3) {
-    case 0:
-      cr = CGAL::compare_x(vt->point(), pop);
-      OnSideMap[vt] = cr == LARGER ? ON_POSITIVE_SIDE :
-                       cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
-      break;
-    case 1:
-      cr = CGAL::compare_y(vt->point(), pop);
-      OnSideMap[vt] = cr == LARGER ? ON_POSITIVE_SIDE :
-                       cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
-      break;
-    case 2:
-      cr = CGAL::compare_z(vt->point(), pop);
-      OnSideMap[vt] = cr == LARGER ? ON_POSITIVE_SIDE :
-                       cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
-      break;
-    default: CGAL_error_msg( "wrong value");
-    }
-  Oriented_side src_side = OnSideMap[v];
-  Oriented_side tgt_side = OnSideMap[vt];
-*/
-  Oriented_side src_side = (*this) (pop, v, depth);
-  Oriented_side tgt_side = (*this) (pop, vt, depth);
   if( src_side == tgt_side)
     return src_side;
   if( src_side == ON_ORIENTED_BOUNDARY)
@@ -457,37 +264,10 @@ Side_of_plane<SNC_decorator>::operator()
 
 
 template <class SNC_decorator>
-template <typename Depth>
 Oriented_side
-Side_of_plane<SNC_decorator>::operator()
+Side_of_plane<SNC_decorator>::operator()(Halffacet_handle f) {
-  (const Point_3& pop, Halffacet_handle f, Depth depth) {
     CGAL_assertion( std::distance( f->facet_cycles_begin(), f->facet_cycles_end()) > 0);
-    /*
+
-#ifdef CGAL_NEF3_FACET_WITH_BOX
-    switch(depth%3) {
-    case 0:
-      if(f->b.min_coord(0) > pop.x())
-        return ON_POSITIVE_SIDE;
-      if(f->b.max_coord(0) < pop.x())
-        return ON_NEGATIVE_SIDE;
-      break;
-    case 1:
-      if(f->b.min_coord(1) > pop.y())
-        return ON_POSITIVE_SIDE;
-      if(f->b.max_coord(1) < pop.y())
-        return ON_NEGATIVE_SIDE;
-      break;
-    case 2:
-      if(f->b.min_coord(2) > pop.z())
-        return ON_POSITIVE_SIDE;
-      if(f->b.max_coord(2) < pop.z())
-        return ON_NEGATIVE_SIDE;
-      break;
-    default: CGAL_error_msg( "wrong value");
-    }
-    return ON_ORIENTED_BOUNDARY;
-#else
-    */
   Halffacet_cycle_iterator fc(f->facet_cycles_begin());
   SHalfedge_handle e;
   CGAL_assertion(fc.is_shalfedge());
@@ -499,7 +279,7 @@ Side_of_plane<SNC_decorator>::operator()
   Vertex_handle v;
   do {
     v = sc->source()->center_vertex();
-    facet_side = (*this) (pop, v, depth);
+    facet_side = (*this) (v);
     ++sc;
   }
   while( facet_side == ON_ORIENTED_BOUNDARY && sc != send);
@@ -509,7 +289,7 @@ Side_of_plane<SNC_decorator>::operator()
   Oriented_side point_side;
   while( sc != send) {
     v = sc->source()->center_vertex();
-    point_side = (*this) (pop, v, depth);
+    point_side = (*this) (v);
     ++sc;
     if( point_side == ON_ORIENTED_BOUNDARY)
       continue;

Nef_3/include/CGAL/Nef_3/SNC_point_locator.h

@@ -190,37 +190,13 @@ public:
 
 
   virtual void initialize(SNC_structure* W) {
-#ifdef CGAL_NEF_LIST_OF_TRIANGLES
+
-    this->set_snc(*W);
-    candidate_provider = new SNC_candidate_provider(W);
-#else // CGAL_NEF_LIST_OF_TRIANGLES
-    CGAL_NEF_TIMER(ct_t.start());
-    this->version_ = std::string("Point Locator by Spatial Subdivision (tm)");
-    CGAL_NEF_CLOG(version());
-    CGAL_assertion( W != nullptr);
-//    (Base) *this = SNC_decorator(*W);
-    this->set_snc(*W);
-    Object_list objects;
-    Vertex_iterator v;
-    Halfedge_iterator e;
-    Halffacet_iterator f;
-    CGAL_forall_vertices( v, *this->sncp())
-      objects.push_back(make_object(Vertex_handle(v)));
-    typename Object_list::difference_type v_end = objects.size();
-    CGAL_forall_edges( e, *this->sncp())
-      objects.push_back(make_object(Halfedge_handle(e)));
-    CGAL_forall_facets( f, *this->sncp()) {
-      objects.push_back(make_object(Halffacet_handle(f)));
-    }
     if(initialized)
       delete candidate_provider;
 
-    Object_list_iterator oli=objects.begin()+v_end;
+    this->set_snc(*W);
-    candidate_provider = new SNC_candidate_provider(objects,oli);
+    candidate_provider = new SNC_candidate_provider(W);
 
-    // CGAL_NEF_TRACEN(*candidate_provider);
-    CGAL_NEF_TIMER(ct_t.stop());
-#endif // CGAL_NEF_LIST_OF_TRIANGLES
     initialized = true;
   }