Merge branch 'Spatial_searching-misc-glisse-old' into Spatial_searching-misc-glisse

Installation/changes.html

@@ -149,6 +149,12 @@ and <code>src/</code> directories).
 <!-- Surface Reconstruction -->
 <!-- Geometry Processing -->
 <!-- Spatial Searching and Sorting -->
+<h3>Spatial Searching</h3>
+  <ul>
+    <li>
+      Add function <code>Kd_tree::remove(Point)</code>.
+    </li>
+  </ul>
 <!-- Geometric Optimization -->
 <!-- Interpolation -->
 <!-- Kinetic Data Structures -->

Spatial_searching/doc/Spatial_searching/CGAL/Kd_tree.h

@@ -96,12 +96,20 @@ template <class InputIterator> Kd_tree(InputIterator first, InputIterator beyond
 The constructor does not build the internal data structure, and it
 is also not updated after calls to `insert()`.  
 The method `build()` is called implicitly
-at the first call to a query member function. You can call
+at the first call to a query or removal member function. You can call
 `build()` explicitly to ensure that the next call to
 query functions will not trigger the reconstruction of the
 data structure.
 */
 void build();
+
+/*!
+This clears the internal data structure, which then gets rebuilt either by an
+explicit call to `build()` or implicitly by the next query or removal. The only
+reason to call this function explicitly is to rebalance the tree after some
+number of removals.
+*/
+void invalidate_built();
 /// @}
 
 /// \name Operations
@@ -109,6 +117,9 @@ void build();
 
 /*!
 Inserts the point `p` in the `k-d` tree.
+\note Insertions do not dynamically update the internal data structure. The
+next query, or a call to `build()`, automatically triggers a rebuild of the
+whole structure.
 */
 void insert(Point_d p);
 
@@ -118,6 +129,23 @@ The value type of the `InputIterator` must be `Point_d`.
 */
 template <class InputIterator> void insert(InputIterator first, InputIterator beyond);
 
+/*!
+Removes the point `p` from the `k-d` tree. It uses `equal_to_p` to identify
+the point after locating it, which can matter in particular when 2 points are
+in the same place. `Identify_point` is a unary functor that takes a `Point_d`
+and returns a `bool`.  This is a limited and naive implementation that does not
+rebalance the tree. On the other hand, the tree remains valid and ready for
+queries. If the internal data structure is not already built, for instance
+because the last operation was an insertion, it first calls `build()`.
+*/
+template<class Identify_point>
+void remove(Point_d p, Identify_point equal_to_p);
+
+/*!
+Removes point `p`, calling the 2-argument function `remove()` with a functor
+that simply compares coordinates.
+*/
+void remove(Point_d p);
 
 /*
 Pre-allocates memory in order to store at least 'size' points.

Spatial_searching/include/CGAL/Kd_tree.h

@@ -106,6 +106,7 @@ private:
   mutable CGAL_MUTEX building_mutex;//mutex used to protect const calls inducing build()
   #endif
   bool built_;
+  bool removed_;
 
   // protected copy constructor
   Kd_tree(const Tree& tree)
@@ -239,13 +240,13 @@ private:
 public:
 
   Kd_tree(Splitter s = Splitter(),const SearchTraits traits=SearchTraits())
-    : traits_(traits),split(s), built_(false)
+    : traits_(traits),split(s), built_(false), removed_(false)
   {}
 
   template <class InputIterator>
   Kd_tree(InputIterator first, InputIterator beyond,
 	  Splitter s = Splitter(),const SearchTraits traits=SearchTraits())
-    : traits_(traits),split(s), built_(false)
+    : traits_(traits),split(s), built_(false), removed_(false)
   {
     pts.insert(pts.end(), first, beyond);
   }
@@ -257,6 +258,10 @@ public:
   void
   build()
   {
+    // This function is not ready to be called when a tree already exists, one
+    // must call invalidate_built() first.
+    CGAL_assertion(!is_built());
+    CGAL_assertion(!removed_);
     const Point_d& p = *pts.begin();
     typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits_.construct_cartesian_const_iterator_d_object();
     int dim = static_cast<int>(std::distance(ccci(p), ccci(p,0)));
@@ -309,6 +314,16 @@ public:
 
   void invalidate_built()
   {
+    if(removed_){
+      // Walk the tree to collect the remaining points.
+      // Writing directly to pts would likely work, but better be safe.
+      std::vector<Point_d> ptstmp;
+      //ptstmp.resize(root()->num_items());
+      root()->tree_items(std::back_inserter(ptstmp));
+      pts.swap(ptstmp);
+      removed_=false;
+      CGAL_assertion(is_built()); // the rest of the cleanup must happen
+    }
     if(is_built()){
       internal_nodes.clear();
       leaf_nodes.clear();
@@ -322,6 +337,7 @@ public:
   {
     invalidate_built();
     pts.clear();
+    removed_ = false;
   }
 
   void
@@ -339,6 +355,98 @@ public:
     pts.insert(pts.end(),first, beyond);
   }
 
+private:
+  struct Equal_by_coordinates {
+    SearchTraits const* traits;
+    Point_d const* pp;
+    bool operator()(Point_d const&q) const {
+      typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits->construct_cartesian_const_iterator_d_object();
+      return std::equal(ccci(*pp), ccci(*pp,0), ccci(q));
+    }
+  };
+  Equal_by_coordinates equal_by_coordinates(Point_d const&p){
+    Equal_by_coordinates ret = { &traits(), &p };
+    return ret;
+  }
+
+public:
+  void
+  remove(const Point_d& p)
+  {
+    remove(p, equal_by_coordinates(p));
+  }
+
+  template<class Equal>
+  void
+  remove(const Point_d& p, Equal const& equal_to_p)
+  {
+#if 0
+    // This code could have quadratic runtime.
+    if (!is_built()) {
+      std::vector<Point_d>::iterator pi = std::find(pts.begin(), pts.end(), p);
+      // Precondition: the point must be there.
+      CGAL_assertion (pi != pts.end());
+      pts.erase(pi);
+      return;
+    }
+#endif
+    bool success = remove_(p, 0, false, 0, false, root(), equal_to_p);
+    CGAL_assertion(success);
+
+    // Do not set the flag is the tree has been cleared.
+    if(is_built())
+      removed_ |= success;
+  }
+private:
+  template<class Equal>
+  bool remove_(const Point_d& p,
+      Internal_node_handle grandparent, bool parent_islower,
+      Internal_node_handle parent, bool islower,
+      Node_handle node, Equal const& equal_to_p) {
+    // Recurse to locate the point
+    if (!node->is_leaf()) {
+      Internal_node_handle newparent = static_cast<Internal_node_handle>(node);
+      // FIXME: This should be if(x<y) remove low; else remove up;
+      if (traits().construct_cartesian_const_iterator_d_object()(p)[newparent->cutting_dimension()] <= newparent->cutting_value()) {
+	if (remove_(p, parent, islower, newparent, true, newparent->lower(), equal_to_p))
+	  return true;
+      }
+      //if (traits().construct_cartesian_const_iterator_d_object()(p)[newparent->cutting_dimension()] >= newparent->cutting_value())
+	return remove_(p, parent, islower, newparent, false, newparent->upper(), equal_to_p);
+
+      CGAL_assertion(false); // Point was not found
+    }
+
+    // Actual removal
+    Leaf_node_handle lnode = static_cast<Leaf_node_handle>(node);
+    if (lnode->size() > 1) {
+      iterator pi = std::find_if(lnode->begin(), lnode->end(), equal_to_p);
+      // FIXME: we should ensure this never happens
+      if (pi == lnode->end()) return false;
+      iterator lasti = lnode->end() - 1;
+      if (pi != lasti) {
+	// Hack to get a non-const iterator
+	std::iter_swap(pts.begin()+(pi-pts.begin()), pts.begin()+(lasti-pts.begin()));
+      }
+      lnode->drop_last_point();
+    } else if (!equal_to_p(*lnode->begin())) {
+      // FIXME: we should ensure this never happens
+      return false;
+    } else if (grandparent) {
+      Node_handle brother = islower ? parent->upper() : parent->lower();
+      if (parent_islower)
+	grandparent->set_lower(brother);
+      else
+	grandparent->set_upper(brother);
+    } else if (parent) {
+      tree_root = islower ? parent->upper() : parent->lower();
+    } else {
+      clear();
+    }
+    return true;
+  }
+
+public:
   //For efficiency; reserve the size of the points vectors in advance (if the number of points is already known).
   void reserve(size_t size)
   {

Spatial_searching/include/CGAL/Kd_tree_node.h

@@ -206,8 +206,7 @@ namespace CGAL {
           static_cast<Internal_node_const_handle>(this);
 	// after splitting b denotes the lower part of b
 	Kd_tree_rectangle<FT,D> b_upper(b);
-	b.split(b_upper, node->cutting_dimension(),
+	node->split_bbox(b, b_upper);
-		node->cutting_value());
 
 	if (q.outer_range_contains(b)) 	
 	  it=node->lower()->tree_items(it);
@@ -236,15 +235,14 @@ namespace CGAL {
 	if (node->size()>0) 
 	  for (iterator i=node->begin(); i != node->end(); i++) 
 	    if (q.contains(*i)) 
-	      { result = boost::make_optional(*i);}
+	      { result = *i; break; }
       }
       else {
          Internal_node_const_handle node = 
           static_cast<Internal_node_const_handle>(this);
 	// after splitting b denotes the lower part of b
 	Kd_tree_rectangle<FT,D> b_upper(b);
-	b.split(b_upper, node->cutting_dimension(),
+	node->split_bbox(b, b_upper);
-		node->cutting_value());
                              
 	if (q.outer_range_contains(b)){ 	
           result = node->lower()->any_tree_item();
@@ -323,6 +321,13 @@ namespace CGAL {
       return data + n;
     }
  
+    inline
+    void
+    drop_last_point()
+    {
+      --n;
+    }
+
   }; //leaf node
 
 
@@ -338,6 +343,7 @@ namespace CGAL {
 
     typedef typename TreeTraits::FT FT;
     typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Separator Separator;
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::D D;
 
   private:
     
@@ -395,6 +401,20 @@ namespace CGAL {
       return upper_ch; 
     }
   	
+    inline
+    void
+    set_lower(Node_handle nh)
+    {
+      lower_ch = nh;
+    }
+
+    inline
+    void
+    set_upper(Node_handle nh)
+    {
+      upper_ch = nh;
+    }
+
     // inline Separator& separator() {return sep; }
     // use instead
     inline
@@ -452,7 +472,12 @@ namespace CGAL {
       return Separator(cutting_dimension,cutting_value);
     }*/
 
-
+    void split_bbox(Kd_tree_rectangle<FT,D>& l, Kd_tree_rectangle<FT,D>& u) const {
+      l.lower()[cut_dim]=lower_low_val;
+      l.upper()[cut_dim]=lower_high_val;
+      u.lower()[cut_dim]=upper_low_val;
+      u.upper()[cut_dim]=upper_high_val;
+    }
   };//internal node
 
  template < class TreeTraits, class Splitter> 
@@ -466,6 +491,7 @@ namespace CGAL {
 
     typedef typename TreeTraits::FT FT;
     typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Separator Separator;
+    typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::D D;
 
   private:
     
@@ -516,6 +542,20 @@ namespace CGAL {
       return upper_ch; 
     }
   	
+    inline
+    void
+    set_lower(Node_handle nh)
+    {
+      lower_ch = nh;
+    }
+
+    inline
+    void
+    set_upper(Node_handle nh)
+    {
+      upper_ch = nh;
+    }
+
     // inline Separator& separator() {return sep; }
     // use instead
 
@@ -539,29 +579,16 @@ namespace CGAL {
       return cut_dim;
     }
 
-    // members for extended internal node only
-    inline 
-    FT
-    low_value() const 
-    { 
-      return this->low_val;
-    }
-    
-    inline 
-    FT
-    high_value() const 
-    {
-      return this->high_val;
-    }
-       
-
    /* Separator& 
     separator() 
     {
       return Separator(cutting_dimension,cutting_value);
     }*/
 
-
+    void split_bbox(Kd_tree_rectangle<FT,D>& l, Kd_tree_rectangle<FT,D>& u) const {
+      l.upper()[cut_dim]=cut_val;
+      u.lower()[cut_dim]=cut_val;
+    }
   };//internal node
 
 

Spatial_searching/include/CGAL/Kd_tree_rectangle.h

@@ -117,7 +117,7 @@ namespace CGAL {
     
     
     explicit 
-    Kd_tree_rectangle(const Kd_tree_rectangle<FT,D>& r) 
+    Kd_tree_rectangle(const Kd_tree_rectangle& r)
     : max_span_coord_(r.max_span_coord_) 
     {
       lower_ = r.lower_;
@@ -216,11 +216,13 @@ namespace CGAL {
       return D::value;
     }
 
-    const T* lower() const {return lower_;}
+    T* lower() {return lower_.data();}
-    const T* upper() const {return upper_;}    
+    T* upper() {return upper_.data();}
+    const T* lower() const {return lower_.data();}
+    const T* upper() const {return upper_.data();}
  
-    Kd_tree_rectangle<FT,D>& 
+    Kd_tree_rectangle&
-    operator=(const Kd_tree_rectangle<FT,D>& r) 
+    operator=(const Kd_tree_rectangle& r)
     {
       CGAL_assertion(dimension() == r.dimension());
       if (this != &r) {
@@ -247,9 +249,8 @@ namespace CGAL {
     
   private:
     
+    T* coords_;
     int dim;
-    T* lower_;
-    T* upper_;
     int max_span_coord_;
     
   public:
@@ -258,8 +259,8 @@ namespace CGAL {
     set_upper_bound(int i, const FT& x) 
     {
       CGAL_assertion(i >= 0 && i < dim);
-      CGAL_assertion(x >= lower_[i]);
+      CGAL_assertion(x >= lower()[i]);
-      upper_[i] = x;
+      upper()[i] = x;
       set_max_span();
     }
 
@@ -267,18 +268,18 @@ namespace CGAL {
     set_lower_bound(int i, const FT& x) 
     {
       CGAL_assertion(i >= 0 && i < dim);
-      CGAL_assertion(x <= upper_[i]);
+      CGAL_assertion(x <= upper()[i]);
-      lower_[i] = x;
+      lower()[i] = x;
       set_max_span();
     }
 
     inline void 
     set_max_span() 
     {
-      FT span = upper_[0]-lower_[0];
+      FT span = upper()[0]-lower()[0];
       max_span_coord_ = 0;
       for (int i = 1; i < dim; ++i) {
-	FT tmp = upper_[i] - lower_[i];
+	FT tmp = upper()[i] - lower()[i];
 	if (span < tmp) {
 	  span = tmp;
 	  max_span_coord_ = i;
@@ -287,25 +288,23 @@ namespace CGAL {
     }
     
     Kd_tree_rectangle(int d) 
-      : dim(d), lower_(new FT[d]), upper_(new FT[d]), max_span_coord_(0)
+      : coords_(new FT[2*d]), dim(d), max_span_coord_(0)
     {
-      std::fill(lower_, lower_ + dim, FT(0));
+      std::fill(coords_, coords_ + 2*dim, FT(0));
-      std::fill(upper_, upper_ + dim, FT(0));
     }
 
     Kd_tree_rectangle() 
-      : dim(0), lower_(0), upper_(0) 
+      : coords_(0), dim(0)
     {
 }
     
     
     explicit 
-    Kd_tree_rectangle(const Kd_tree_rectangle<FT,Dynamic_dimension_tag>& r) 
+    Kd_tree_rectangle(const Kd_tree_rectangle& r)
-      : dim(r.dim), lower_(new FT[dim]), upper_(new FT[dim]), 
+      : coords_(new FT[2*r.dim]), dim(r.dim),
 	max_span_coord_(r.max_span_coord_) 
     {
-      std::copy(r.lower_, r.lower_+dim, lower_);
+      std::copy(r.coords_, r.coords_+2*dim, lower());
-      std::copy(r.upper_, r.upper_+dim, upper_);
     }
 
     template <class Construct_cartesian_const_iterator_d,class PointPointerIter>
@@ -320,17 +319,17 @@ namespace CGAL {
       typename Construct_cartesian_const_iterator_d::result_type bit = construct_it(**begin);
       
       for (int i=0; i < dim; ++i, ++bit) {
-        lower_[i]= *bit; upper_[i]=lower_[i];
+        lower()[i]= *bit; upper()[i]=lower()[i];
       }
       begin++;
       typedef typename std::iterator_traits<PointPointerIter>::value_type P;
-      std::for_each(begin, end,set_bounds_from_pointer<Construct_cartesian_const_iterator_d,P,T>(dim, lower_, upper_,construct_it));
+      std::for_each(begin, end,set_bounds_from_pointer<Construct_cartesian_const_iterator_d,P,T>(dim, lower(), upper(),construct_it));
       set_max_span();
     }
     
     template <class Construct_cartesian_const_iterator_d,class PointPointerIter> // was PointIter
     Kd_tree_rectangle(int d,  PointPointerIter begin,  PointPointerIter end,const Construct_cartesian_const_iterator_d& construct_it)
-      : dim(d), lower_(new FT[d]), upper_(new FT[d])
+      : coords_(new FT[2*d]), dim(d)
     {
       update_from_point_pointers<Construct_cartesian_const_iterator_d>(begin,end,construct_it);
     }
@@ -344,20 +343,20 @@ namespace CGAL {
     inline FT 
     max_span() const 
     {
-      return upper_[max_span_coord_] - lower_[max_span_coord_];
+      return upper()[max_span_coord_] - lower()[max_span_coord_];
     }
 
     inline FT 
     min_coord(int i) const 
     {
-      CGAL_assertion(lower_ != NULL);
+      CGAL_assertion(coords_ != NULL);
-      return lower_[i];
+      return lower()[i];
     }
     
     inline FT 
     max_coord(int i) const 
     {
-      return upper_[i];
+      return upper()[i];
     }
     
     std::ostream& 
@@ -366,13 +365,13 @@ namespace CGAL {
       s << "Rectangle dimension = " << dim;
       s << "\n lower: ";
       for (int i=0; i < dim; ++i)
-	s << lower_[i] << " ";
+	s << lower()[i] << " ";
-      // std::copy(lower_, lower_ + dim,
+      // std::copy(lower(), lower() + dim,
       // 	      std::ostream_iterator<FT>(s," "));
       s << "\n upper: ";
       for (int j=0; j < dim; ++j)
-	s << upper_[j] << " ";
+	s << upper()[j] << " ";
-      // std::copy(upper_, upper_ + dim,
+      // std::copy(upper(), upper() + dim,
       //	      std::ostream_iterator<FT>(s," "));
       s << "\n maximum span " << max_span() <<
 	" at coordinate " << max_span_coord() << std::endl;
@@ -386,11 +385,11 @@ namespace CGAL {
     split(Kd_tree_rectangle& r, int d, FT value)
     {
       CGAL_assertion(d >= 0 && d < dim);
-      CGAL_assertion(lower_[d] <= value && value <= upper_[d]);
+      CGAL_assertion(lower()[d] <= value && value <= upper()[d]);
       
       //Kd_tree_rectangle* r = new Kd_tree_rectangle(*this);
-      upper_[d]=value;
+      upper()[d]=value;
-      r.lower_[d]=value;
+      r.lower()[d]=value;
       //return r;
     }
     
@@ -398,8 +397,7 @@ namespace CGAL {
     ~Kd_tree_rectangle() 
     {
       if (dim) {
-	if (lower_) delete [] lower_;
+	if (coords_) delete [] coords_;
-	if (upper_) delete [] upper_;
       }
     }
     
@@ -409,16 +407,17 @@ namespace CGAL {
       return dim;
     }
 
-    const T* lower() const {return lower_;}
+    T* lower() {return coords_;}
-    const T* upper() const {return upper_;}    
+    T* upper() {return coords_ + dim;}
+    const T* lower() const {return coords_;}
+    const T* upper() const {return coords_ + dim;}
  
-    Kd_tree_rectangle<FT,Dynamic_dimension_tag>& 
+    Kd_tree_rectangle&
-    operator=(const Kd_tree_rectangle<FT,Dynamic_dimension_tag>& r) 
+    operator=(const Kd_tree_rectangle& r)
     {
       CGAL_assertion(dimension() == r.dimension());
       if (this != &r) {
-        std::copy(r.lower_, r.lower_+dim, lower_);
+        std::copy(r.coords_, r.coords_+2*dim, coords_);
-	std::copy(r.upper_, r.upper_+dim, upper_);
 	set_max_span();
       }
       return *this;

Spatial_searching/include/CGAL/Weighted_Minkowski_distance.h

@@ -30,11 +30,24 @@
 #include <cmath>
 #include <vector>
 
+#include <CGAL/array.h>
 #include <CGAL/number_utils.h>
 #include <CGAL/Kd_tree_rectangle.h>
 #include <CGAL/internal/Get_dimension_tag.h>
 
 namespace CGAL {
+  namespace internal {
+    template<class T, class Dim>
+      struct Array_or_vector_selector {
+	typedef std::vector<T> type;
+	static void resize(type&v, std::size_t d) { v.resize(d); }
+      };
+    template<class T, int D>
+      struct Array_or_vector_selector<T, Dimension_tag<D> > {
+	typedef cpp11::array<T,D> type;
+	static void resize(type&, std::size_t CGAL_assertion_code(d)) { CGAL_assertion(d==D); }
+      };
+  }
 
   template <class SearchTraits>
   class Weighted_Minkowski_distance {
@@ -44,8 +57,9 @@ namespace CGAL {
     typedef typename SearchTraits::Point_d Point_d;
     typedef Point_d                        Query_item;
     typedef typename SearchTraits::FT      FT;
-    typedef std::vector<FT>                Weight_vector;
     typedef typename internal::Get_dimension_tag<SearchTraits>::Dimension Dimension;
+    typedef internal::Array_or_vector_selector<FT,Dimension> Weight_vector_traits;
+    typedef typename Weight_vector_traits::type Weight_vector;
 
     private:
 
@@ -71,28 +85,15 @@ namespace CGAL {
     //default copy constructor and destructor
     
 
-    Weighted_Minkowski_distance (FT pow, int dim,
-				 const Weight_vector& weights,
-                                 const SearchTraits& traits_=SearchTraits()) 
-      : traits(traits_),power(pow)
-    {
-      CGAL_assertion(power >= FT(0));
-      CGAL_assertion(dim==weights.size());
-      for (unsigned int i = 0; i < weights.size(); ++i)
-	CGAL_assertion(weights[i]>=FT(0));
-      the_weights.resize(weights.size());
-      the_weights = weights;
-    }
-
     template <class InputIterator>
     Weighted_Minkowski_distance (FT pow, int dim,
-				 InputIterator begin, InputIterator end,
+				 InputIterator begin,
+				 InputIterator CGAL_assertion_code(end),
                                  const SearchTraits& traits_=SearchTraits()) 
       : traits(traits_),power(pow)
     {
       CGAL_assertion(power >= FT(0));
-      the_weights.resize(dim);
+      Weight_vector_traits::resize(the_weights, dim);
-      std::copy(begin, end, the_weights.begin());
       for (int i = 0; i < dim; ++i){
 	the_weights[i] = *begin;
 	++begin;

Spatial_searching/test/Spatial_searching/remove.cpp

@@ -0,0 +1,29 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Search_traits_2.h>
+#include <CGAL/Kd_tree.h>
+
+typedef CGAL::Simple_cartesian<double> K;
+typedef K::Point_2 Point;
+typedef CGAL::Search_traits_2<K> Traits;
+typedef CGAL::Kd_tree<Traits> Tree;
+
+int main()
+{
+  Tree t;
+  t.insert(Point(0,0));
+  t.insert(Point(1,2));
+  t.insert(Point(2,0));
+  t.remove(Point(1,2));
+  t.insert(Point(3,4));
+  t.remove(Point(0,0));
+  t.remove(Point(3,4));
+  t.insert(Point(5,5));
+  t.build();
+  assert(t.size()==2);
+  t.clear();
+  for(int i=0;i<1000;++i)
+    t.insert(Point(i,-i));
+  for(int i=0;i<1000;++i)
+    t.remove(Point(i,-i));
+  t.print();
+}