Merge last trunk modifications

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/distance_matrix_in_rips_module@1938 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 8af39a00b9134dbcc77883412c29d7184479979b

6 files changed, 1876 insertions, 6 deletions

diff --git a/src/common/include/gudhi/random_point_generators.h b/src/common/include/gudhi/random_point_generators.h
index 3050b7ea..2ec465ef 100644
--- a/src/common/include/gudhi/random_point_generators.h
+++ b/src/common/include/gudhi/random_point_generators.h
@@ -192,10 +192,9 @@ static void generate_uniform_points_on_torus_d(const Kernel &k, int dim, std::si
                                                double radius_noise_percentage = 0.,
                                                std::vector<typename Kernel::FT> current_point = std::vector<typename Kernel::FT>()) {
   CGAL::Random rng;
-  if (current_point.size() == 2 * dim) {
-    *out++ = k.construct_point_d_object()(
-                                          static_cast<int> (current_point.size()),
-                                          current_point.begin(), current_point.end());
+  int point_size = static_cast<int>(current_point.size());
+  if (point_size == 2 * dim) {
+    *out++ = k.construct_point_d_object()(point_size, current_point.begin(), current_point.end());
   } else {
     for (std::size_t slice_idx = 0; slice_idx < num_slices; ++slice_idx) {
       double radius_noise_ratio = 1.;
@@ -338,8 +337,6 @@ std::vector<typename Kernel::Point_d> generate_points_on_3sphere_and_circle(std:
   std::vector<Point> points;
   points.reserve(num_points);
 
-  typename Kernel::Translated_point_d k_transl =
-      k.translated_point_d_object();
   typename Kernel::Compute_coordinate_d k_coord =
       k.compute_coordinate_d_object();
   for (std::size_t i = 0; i < num_points;) {
diff --git a/src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt b/src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt
new file mode 100644
index 00000000..a588d113
--- /dev/null
+++ b/src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt
@@ -0,0 +1,3 @@
+CGAL/Kd_tree.h
+CGAL/Kd_tree_node.h
+CGAL/Orthogonal_incremental_neighbor_search.h
diff --git a/src/common/include/gudhi_patches/CGAL/Kd_tree.h b/src/common/include/gudhi_patches/CGAL/Kd_tree.h
new file mode 100644
index 00000000..f085b0da
--- /dev/null
+++ b/src/common/include/gudhi_patches/CGAL/Kd_tree.h
@@ -0,0 +1,582 @@
+// Copyright (c) 2002,2011,2014 Utrecht University (The Netherlands), Max-Planck-Institute Saarbruecken (Germany).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Hans Tangelder (<hanst@cs.uu.nl>),
+//               : Waqar Khan <wkhan@mpi-inf.mpg.de>
+
+#ifndef CGAL_KD_TREE_H
+#define CGAL_KD_TREE_H
+
+#include "Kd_tree_node.h"
+
+#include <CGAL/basic.h>
+#include <CGAL/assertions.h>
+#include <vector>
+
+#include <CGAL/algorithm.h>
+#include <CGAL/internal/Get_dimension_tag.h>
+#include <CGAL/Search_traits.h>
+
+
+#include <deque>
+#include <boost/container/deque.hpp>
+#include <boost/optional.hpp>
+
+#ifdef CGAL_HAS_THREADS
+#include <CGAL/mutex.h>
+#endif
+
+namespace CGAL {
+
+//template <class SearchTraits, class Splitter_=Median_of_rectangle<SearchTraits>, class UseExtendedNode = Tag_true >
+template <class SearchTraits, class Splitter_=Sliding_midpoint<SearchTraits>, class UseExtendedNode = Tag_true >
+class Kd_tree {
+
+public:
+  typedef SearchTraits Traits;
+  typedef Splitter_ Splitter;
+  typedef typename SearchTraits::Point_d Point_d;
+  typedef typename Splitter::Container Point_container;
+
+  typedef typename SearchTraits::FT FT;
+  typedef Kd_tree_node<SearchTraits, Splitter, UseExtendedNode > Node;
+  typedef Kd_tree_leaf_node<SearchTraits, Splitter, UseExtendedNode > Leaf_node;
+  typedef Kd_tree_internal_node<SearchTraits, Splitter, UseExtendedNode > Internal_node;
+  typedef Kd_tree<SearchTraits, Splitter> Tree;
+  typedef Kd_tree<SearchTraits, Splitter,UseExtendedNode> Self;
+
+  typedef Node* Node_handle;
+  typedef const Node* Node_const_handle;
+  typedef Leaf_node* Leaf_node_handle;
+  typedef const Leaf_node* Leaf_node_const_handle;
+  typedef Internal_node* Internal_node_handle;
+  typedef const Internal_node* Internal_node_const_handle;
+  typedef typename std::vector<const Point_d*>::const_iterator Point_d_iterator;
+  typedef typename std::vector<const Point_d*>::const_iterator Point_d_const_iterator;
+  typedef typename Splitter::Separator Separator;
+  typedef typename std::vector<Point_d>::const_iterator iterator;
+  typedef typename std::vector<Point_d>::const_iterator const_iterator;
+
+  typedef typename std::vector<Point_d>::size_type size_type;
+
+  typedef typename internal::Get_dimension_tag<SearchTraits>::Dimension D;
+
+private:
+  SearchTraits traits_;
+  Splitter split;
+
+
+  // wokaround for https://svn.boost.org/trac/boost/ticket/9332
+#if   (_MSC_VER == 1800) && (BOOST_VERSION == 105500)
+  std::deque<Internal_node> internal_nodes;
+  std::deque<Leaf_node> leaf_nodes;
+#else
+  boost::container::deque<Internal_node> internal_nodes;
+  boost::container::deque<Leaf_node> leaf_nodes;
+#endif
+
+  Node_handle tree_root;
+
+  Kd_tree_rectangle<FT,D>* bbox;
+  std::vector<Point_d> pts;
+
+  // Instead of storing the points in arrays in the Kd_tree_node
+  // we put all the data in a vector in the Kd_tree.
+  // and we only store an iterator range in the Kd_tree_node.
+  //
+  std::vector<const Point_d*> data;
+
+
+  #ifdef CGAL_HAS_THREADS
+  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)
+    : traits_(tree.traits_),built_(tree.built_)
+  {};
+
+
+  // Instead of the recursive construction of the tree in the class Kd_tree_node
+  // we do this in the tree class. The advantage is that we then can optimize
+  // the allocation of the nodes.
+
+  // The leaf node
+  Node_handle
+  create_leaf_node(Point_container& c)
+  {
+    Leaf_node node(true , static_cast<unsigned int>(c.size()));
+    std::ptrdiff_t tmp = c.begin() - data.begin();
+    node.data = pts.begin() + tmp;
+
+    leaf_nodes.push_back(node);
+    Leaf_node_handle nh = &leaf_nodes.back();
+
+
+    return nh;
+  }
+
+
+  // The internal node
+
+  Node_handle
+  create_internal_node(Point_container& c, const Tag_true&)
+  {
+    return create_internal_node_use_extension(c);
+  }
+
+  Node_handle
+  create_internal_node(Point_container& c, const Tag_false&)
+  {
+    return create_internal_node(c);
+  }
+
+
+
+  // TODO: Similiar to the leaf_init function above, a part of the code should be
+  //       moved to a the class Kd_tree_node.
+  //       It is not proper yet, but the goal was to see if there is
+  //       a potential performance gain through the Compact_container
+  Node_handle
+  create_internal_node_use_extension(Point_container& c)
+  {
+    Internal_node node(false);
+    internal_nodes.push_back(node);
+    Internal_node_handle nh = &internal_nodes.back();
+
+    Separator sep;
+    Point_container c_low(c.dimension(),traits_);
+    split(sep, c, c_low);
+    nh->set_separator(sep);
+
+    int cd  = nh->cutting_dimension();
+    if(!c_low.empty()){
+      nh->lower_low_val = c_low.tight_bounding_box().min_coord(cd);
+      nh->lower_high_val = c_low.tight_bounding_box().max_coord(cd);
+    }
+    else{
+      nh->lower_low_val = nh->cutting_value();
+      nh->lower_high_val = nh->cutting_value();
+    }
+    if(!c.empty()){
+      nh->upper_low_val = c.tight_bounding_box().min_coord(cd);
+      nh->upper_high_val = c.tight_bounding_box().max_coord(cd);
+    }
+    else{
+      nh->upper_low_val = nh->cutting_value();
+      nh->upper_high_val = nh->cutting_value();
+    }
+
+    CGAL_assertion(nh->cutting_value() >= nh->lower_low_val);
+    CGAL_assertion(nh->cutting_value() <= nh->upper_high_val);
+
+    if (c_low.size() > split.bucket_size()){
+      nh->lower_ch = create_internal_node_use_extension(c_low);
+    }else{
+      nh->lower_ch = create_leaf_node(c_low);
+    }
+    if (c.size() > split.bucket_size()){
+      nh->upper_ch = create_internal_node_use_extension(c);
+    }else{
+      nh->upper_ch = create_leaf_node(c);
+    }
+
+
+
+
+    return nh;
+  }
+
+
+  // Note also that I duplicated the code to get rid if the if's for
+  // the boolean use_extension which was constant over the construction
+  Node_handle
+  create_internal_node(Point_container& c)
+  {
+    Internal_node node(false);
+    internal_nodes.push_back(node);
+    Internal_node_handle nh = &internal_nodes.back();
+    Separator sep;
+
+    Point_container c_low(c.dimension(),traits_);
+    split(sep, c, c_low);
+    nh->set_separator(sep);
+
+    if (c_low.size() > split.bucket_size()){
+      nh->lower_ch = create_internal_node(c_low);
+    }else{
+      nh->lower_ch = create_leaf_node(c_low);
+    }
+    if (c.size() > split.bucket_size()){
+      nh->upper_ch = create_internal_node(c);
+    }else{
+      nh->upper_ch = create_leaf_node(c);
+    }
+
+
+
+    return nh;
+  }
+
+
+
+public:
+
+  Kd_tree(Splitter s = Splitter(),const SearchTraits traits=SearchTraits())
+    : 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), removed_(false)
+  {
+    pts.insert(pts.end(), first, beyond);
+  }
+
+  bool empty() const {
+    return pts.empty();
+  }
+
+  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)));
+
+    data.reserve(pts.size());
+    for(unsigned int i = 0; i < pts.size(); i++){
+      data.push_back(&pts[i]);
+    }
+    Point_container c(dim, data.begin(), data.end(),traits_);
+    bbox = new Kd_tree_rectangle<FT,D>(c.bounding_box());
+    if (c.size() <= split.bucket_size()){
+      tree_root = create_leaf_node(c);
+    }else {
+      tree_root = create_internal_node(c, UseExtendedNode());
+    }
+
+    //Reorder vector for spatial locality
+    std::vector<Point_d> ptstmp;
+    ptstmp.resize(pts.size());
+    for (std::size_t i = 0; i < pts.size(); ++i){
+      ptstmp[i] = *data[i];
+    }
+    for(std::size_t i = 0; i < leaf_nodes.size(); ++i){
+      std::ptrdiff_t tmp = leaf_nodes[i].begin() - pts.begin();
+      leaf_nodes[i].data = ptstmp.begin() + tmp;
+    }
+    pts.swap(ptstmp);
+
+    data.clear();
+
+    built_ = true;
+  }
+
+private:
+  //any call to this function is for the moment not threadsafe
+  void const_build() const {
+    #ifdef CGAL_HAS_THREADS
+    //this ensure that build() will be called once
+    CGAL_SCOPED_LOCK(building_mutex);
+    if(!is_built())
+    #endif
+      const_cast<Self*>(this)->build(); //THIS IS NOT THREADSAFE
+  }
+public:
+
+  bool is_built() const
+  {
+    return built_;
+  }
+
+  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();
+      data.clear();
+      delete bbox;
+      built_ = false;
+    }
+  }
+
+  void clear()
+  {
+    invalidate_built();
+    pts.clear();
+    removed_ = false;
+  }
+
+  void
+  insert(const Point_d& p)
+  {
+    invalidate_built();
+    pts.push_back(p);
+  }
+
+  template <class InputIterator>
+  void
+  insert(InputIterator first, InputIterator beyond)
+  {
+    invalidate_built();
+    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)
+  {
+    pts.reserve(size);
+  }
+
+  //Get the capacity of the underlying points vector.
+  size_t capacity()
+  {
+    return pts.capacity();
+  }
+
+
+  template <class OutputIterator, class FuzzyQueryItem>
+  OutputIterator
+  search(OutputIterator it, const FuzzyQueryItem& q) const
+  {
+    if(! pts.empty()){
+
+      if(! is_built()){
+        const_build();
+      }
+      Kd_tree_rectangle<FT,D> b(*bbox);
+      return tree_root->search(it,q,b);
+    }
+    return it;
+  }
+
+
+  template <class FuzzyQueryItem>
+  boost::optional<Point_d>
+  search_any_point(const FuzzyQueryItem& q) const
+  {
+    if(! pts.empty()){
+
+      if(! is_built()){
+        const_build();
+      }
+      Kd_tree_rectangle<FT,D> b(*bbox);
+      return tree_root->search_any_point(q,b);
+    }
+    return boost::none;
+  }
+
+
+  ~Kd_tree() {
+    if(is_built()){
+      delete bbox;
+    }
+  }
+
+
+  const SearchTraits&
+  traits() const
+  {
+    return traits_;
+  }
+
+  Node_const_handle
+  root() const
+  {
+    if(! is_built()){
+      const_build();
+    }
+    return tree_root;
+  }
+
+  Node_handle
+  root()
+  {
+    if(! is_built()){
+      build();
+    }
+    return tree_root;
+  }
+
+  void
+  print() const
+  {
+    if(! is_built()){
+      const_build();
+    }
+    root()->print();
+  }
+
+  const Kd_tree_rectangle<FT,D>&
+  bounding_box() const
+  {
+    if(! is_built()){
+      const_build();
+    }
+    return *bbox;
+  }
+
+  const_iterator
+  begin() const
+  {
+    return pts.begin();
+  }
+
+  const_iterator
+  end() const
+  {
+    return pts.end();
+  }
+
+  size_type
+  size() const
+  {
+    return pts.size();
+  }
+
+  // Print statistics of the tree.
+  std::ostream&
+  statistics(std::ostream& s) const
+  {
+    if(! is_built()){
+      const_build();
+    }
+    s << "Tree statistics:" << std::endl;
+    s << "Number of items stored: "
+      << root()->num_items() << std::endl;
+    s << "Number of nodes: "
+      << root()->num_nodes() << std::endl;
+    s << " Tree depth: " << root()->depth() << std::endl;
+    return s;
+  }
+
+
+};
+
+} // namespace CGAL
+
+#endif // CGAL_KD_TREE_H
diff --git a/src/common/include/gudhi_patches/CGAL/Kd_tree_node.h b/src/common/include/gudhi_patches/CGAL/Kd_tree_node.h
new file mode 100644
index 00000000..909ee260
--- /dev/null
+++ b/src/common/include/gudhi_patches/CGAL/Kd_tree_node.h
@@ -0,0 +1,586 @@
+// Copyright (c) 2002,2011  Utrecht University (The Netherlands).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// 
+//
+// Authors       : Hans Tangelder (<hanst@cs.uu.nl>)
+
+#ifndef CGAL_KD_TREE_NODE_H
+#define CGAL_KD_TREE_NODE_H
+
+#include "CGAL/Splitters.h"
+
+#include <CGAL/Compact_container.h>
+#include <boost/cstdint.hpp>
+
+namespace CGAL {
+
+  template <class SearchTraits, class Splitter, class UseExtendedNode>
+  class Kd_tree;
+
+  template < class TreeTraits, class Splitter, class UseExtendedNode >
+  class Kd_tree_node {
+
+     friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>;
+
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_handle Node_handle;
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_const_handle Node_const_handle;
+     typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Internal_node_handle Internal_node_handle;
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Internal_node_const_handle Internal_node_const_handle;
+     typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Leaf_node_handle Leaf_node_handle;
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Leaf_node_const_handle Leaf_node_const_handle;
+    typedef typename TreeTraits::Point_d Point_d;
+
+    typedef typename TreeTraits::FT FT;
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Separator Separator;
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Point_d_iterator Point_d_iterator;
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::iterator iterator;
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::D D;
+
+    bool leaf;
+
+  public :
+    Kd_tree_node(bool leaf_)
+      :leaf(leaf_){}
+
+    bool is_leaf() const{
+      return leaf;
+    }
+
+    std::size_t
+    num_items() const
+    {
+      if (is_leaf()){
+        Leaf_node_const_handle node =
+          static_cast<Leaf_node_const_handle>(this);
+        return node->size();
+      }
+      else {
+        Internal_node_const_handle node =
+          static_cast<Internal_node_const_handle>(this);
+        return node->lower()->num_items() + node->upper()->num_items();
+      }
+    }
+
+    std::size_t
+    num_nodes() const
+    {
+      if (is_leaf()) return 1;
+      else {
+        Internal_node_const_handle node =
+          static_cast<Internal_node_const_handle>(this);
+        return node->lower()->num_nodes() + node->upper()->num_nodes();
+      }
+    }
+
+    int
+    depth(const int current_max_depth) const
+    {
+      if (is_leaf()){
+        return current_max_depth;
+      }
+      else {
+        Internal_node_const_handle node =
+          static_cast<Internal_node_const_handle>(this);
+        return
+             (std::max)( node->lower()->depth(current_max_depth + 1),
+                         node->upper()->depth(current_max_depth + 1));
+      }
+    }
+
+    int
+    depth() const
+    {
+      return depth(1);
+    }
+
+    template <class OutputIterator>
+    OutputIterator
+    tree_items(OutputIterator it) const {
+      if (is_leaf()) {
+         Leaf_node_const_handle node =
+          static_cast<Leaf_node_const_handle>(this);
+         if (node->size()>0)
+            for (iterator i=node->begin(); i != node->end(); i++)
+              {*it=*i; ++it;}
+        }
+      else {
+         Internal_node_const_handle node =
+          static_cast<Internal_node_const_handle>(this);
+          it=node->lower()->tree_items(it);
+          it=node->upper()->tree_items(it);
+      }
+      return it;
+    }
+
+
+    boost::optional<Point_d>
+    any_tree_item() const {
+      boost::optional<Point_d> result = boost::none;
+      if (is_leaf()) {
+         Leaf_node_const_handle node =
+          static_cast<Leaf_node_const_handle>(this);
+         if (node->size()>0){
+           return boost::make_optional(*(node->begin()));
+         }
+        }
+      else {
+         Internal_node_const_handle node =
+          static_cast<Internal_node_const_handle>(this);
+          result = node->lower()->any_tree_item();
+          if(! result){
+            result = node->upper()->any_tree_item();
+          }
+      }
+      return result;
+    }
+
+
+     void
+    indent(int d) const
+    {
+      for(int i = 0; i < d; i++){
+        std::cout << " ";
+      }
+    }
+
+
+    void
+    print(int d = 0) const
+    {
+      if (is_leaf()) {
+        Leaf_node_const_handle node =
+          static_cast<Leaf_node_const_handle>(this);
+        indent(d);
+        std::cout << "leaf" << std::endl;
+        if (node->size()>0)
+          for (iterator i=node->begin(); i != node->end(); i++)
+          {indent(d);std::cout << *i << std::endl;}
+      }
+      else {
+        Internal_node_const_handle node =
+          static_cast<Internal_node_const_handle>(this);
+        indent(d);
+        std::cout << "lower tree" << std::endl;
+        node->lower()->print(d+1);
+        indent(d);
+        std::cout << "separator: dim = " << node->cutting_dimension() << "  val = " << node->cutting_value() << std::endl;
+        indent(d);
+        std::cout << "upper tree" << std::endl;
+        node->upper()->print(d+1);
+      }
+    }
+
+
+    template <class OutputIterator, class FuzzyQueryItem>
+    OutputIterator
+    search(OutputIterator it, const FuzzyQueryItem& q,
+           Kd_tree_rectangle<FT,D>& b) const
+    {
+      if (is_leaf()) {
+        Leaf_node_const_handle node =
+          static_cast<Leaf_node_const_handle>(this);
+        if (node->size()>0)
+          for (iterator i=node->begin(); i != node->end(); i++)
+            if (q.contains(*i))
+              {*it++=*i;}
+      }
+      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->cutting_value());
+
+        if (q.outer_range_contains(b))
+          it=node->lower()->tree_items(it);
+        else
+          if (q.inner_range_intersects(b))
+            it=node->lower()->search(it,q,b);
+        if  (q.outer_range_contains(b_upper))
+          it=node->upper()->tree_items(it);
+        else
+          if (q.inner_range_intersects(b_upper))
+            it=node->upper()->search(it,q,b_upper);
+      };
+      return it;
+    }
+
+
+    template <class FuzzyQueryItem>
+    boost::optional<Point_d>
+    search_any_point(const FuzzyQueryItem& q,
+                     Kd_tree_rectangle<FT,D>& b) const
+    {
+      boost::optional<Point_d> result = boost::none;
+      if (is_leaf()) {
+        Leaf_node_const_handle node =
+          static_cast<Leaf_node_const_handle>(this);
+        if (node->size()>0)
+          for (iterator i=node->begin(); i != node->end(); i++)
+            if (q.contains(*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->cutting_value());
+
+        if (q.outer_range_contains(b)){
+          result = node->lower()->any_tree_item();
+        }else{
+          if (q.inner_range_intersects(b)){
+            result = node->lower()->search_any_point(q,b);
+          }
+        }
+        if(result){
+          return result;
+        }
+        if  (q.outer_range_contains(b_upper)){
+          result = node->upper()->any_tree_item();
+        }else{
+          if (q.inner_range_intersects(b_upper))
+            result = node->upper()->search_any_point(q,b_upper);
+        }
+      }
+      return result;
+    }
+
+  };
+
+
+  template < class TreeTraits, class Splitter, class UseExtendedNode >
+  class Kd_tree_leaf_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{
+
+    friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>;
+
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::iterator iterator;
+    typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base;
+    typedef typename TreeTraits::Point_d Point_d;
+
+  private:
+
+    // private variables for leaf nodes
+    boost::int32_t n; // denotes number of items in a leaf node
+    iterator data; // iterator to data in leaf node
+
+
+  public:
+
+    // default constructor
+    Kd_tree_leaf_node()
+    {}
+
+    Kd_tree_leaf_node(bool leaf_ )
+      : Base(leaf_)
+    {}
+
+    Kd_tree_leaf_node(bool leaf_,unsigned int n_ )
+      : Base(leaf_), n(n_)
+    {}
+
+    // members for all nodes
+
+    // members for leaf nodes only
+    inline
+    unsigned int
+    size() const
+    {
+      return n;
+    }
+
+    inline
+    iterator
+    begin() const
+    {
+      return data;
+    }
+
+    inline
+    iterator
+    end() const
+    {
+      return data + n;
+    }
+
+    inline
+    void
+    drop_last_point()
+    {
+      --n;
+    }
+
+  }; //leaf node
+
+
+
+  template < class TreeTraits, class Splitter, class UseExtendedNode>
+  class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{
+
+    friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>;
+
+    typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base;
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_handle Node_handle;
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_const_handle Node_const_handle;
+
+    typedef typename TreeTraits::FT FT;
+    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Separator Separator;
+
+  private:
+
+       // private variables for internal nodes
+    boost::int32_t cut_dim;
+    FT cut_val;
+    Node_handle lower_ch, upper_ch;
+
+
+    // private variables for extended internal nodes
+    FT upper_low_val;
+    FT upper_high_val;
+    FT lower_low_val;
+    FT lower_high_val;
+
+
+  public:
+
+    // default constructor
+    Kd_tree_internal_node()
+    {}
+
+    Kd_tree_internal_node(bool leaf_)
+      : Base(leaf_)
+    {}
+
+
+    // members for internal node and extended internal node
+
+    inline
+    Node_const_handle
+    lower() const
+    {
+      return lower_ch;
+    }
+
+    inline
+    Node_const_handle
+    upper() const
+    {
+      return upper_ch;
+    }
+
+    inline
+    Node_handle
+    lower()
+    {
+      return lower_ch;
+    }
+
+    inline
+    Node_handle
+    upper()
+    {
+      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
+    void set_separator(Separator& sep){
+      cut_dim = sep.cutting_dimension();
+      cut_val = sep.cutting_value();
+    }
+
+    inline
+    FT
+    cutting_value() const
+    {
+      return cut_val;
+    }
+
+    inline
+    int
+    cutting_dimension() const
+    {
+      return cut_dim;
+    }
+
+    // members for extended internal node only
+    inline
+    FT
+    upper_low_value() const
+    {
+      return upper_low_val;
+    }
+
+    inline
+    FT
+    upper_high_value() const
+    {
+      return upper_high_val;
+    }
+
+    inline
+    FT
+    lower_low_value() const
+    {
+      return lower_low_val;
+    }
+
+    inline
+    FT
+    lower_high_value() const
+    {
+      return lower_high_val;
+    }
+
+    /*Separator&
+    separator()
+    {
+      return Separator(cutting_dimension,cutting_value);
+    }*/
+
+
+  };//internal node
+
+ template < class TreeTraits, class Splitter>
+ class Kd_tree_internal_node<TreeTraits,Splitter,Tag_false> : public Kd_tree_node< TreeTraits, Splitter, Tag_false >{
+
+    friend class Kd_tree<TreeTraits,Splitter,Tag_false>;
+
+    typedef Kd_tree_node< TreeTraits, Splitter, Tag_false> Base;
+    typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Node_handle Node_handle;
+    typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Node_const_handle Node_const_handle;
+
+    typedef typename TreeTraits::FT FT;
+    typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Separator Separator;
+
+  private:
+
+       // private variables for internal nodes
+    boost::uint8_t cut_dim;
+    FT cut_val;
+
+    Node_handle lower_ch, upper_ch;
+
+  public:
+
+    // default constructor
+    Kd_tree_internal_node()
+    {}
+
+    Kd_tree_internal_node(bool leaf_)
+      : Base(leaf_)
+    {}
+
+
+    // members for internal node and extended internal node
+
+    inline
+    Node_const_handle
+    lower() const
+    {
+      return lower_ch;
+    }
+
+    inline
+    Node_const_handle
+    upper() const
+    {
+      return upper_ch;
+    }
+
+    inline
+    Node_handle
+    lower()
+    {
+      return lower_ch;
+    }
+
+    inline
+    Node_handle
+    upper()
+    {
+      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
+    void set_separator(Separator& sep){
+      cut_dim = sep.cutting_dimension();
+      cut_val = sep.cutting_value();
+    }
+
+    inline
+    FT
+    cutting_value() const
+    {
+      return cut_val;
+    }
+
+    inline
+    int
+    cutting_dimension() const
+    {
+      return cut_dim;
+    }
+
+   /* Separator&
+    separator()
+    {
+      return Separator(cutting_dimension,cutting_value);
+    }*/
+
+
+  };//internal node
+
+
+
+} // namespace CGAL
+#endif // CGAL_KDTREE_NODE_H
diff --git a/src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h b/src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h
new file mode 100644
index 00000000..e29ce14f
--- /dev/null
+++ b/src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h
@@ -0,0 +1,620 @@
+// Copyright (c) 2002,2011 Utrecht University (The Netherlands).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// 
+//
+// Author(s)     : Hans Tangelder (<hanst@cs.uu.nl>)
+
+#ifndef CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH
+#define CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH
+
+#include <CGAL/Kd_tree.h>
+#include <cstring>
+#include <list>
+#include <queue>
+#include <memory>
+#include <CGAL/Euclidean_distance.h>
+#include <CGAL/tuple.h>
+
+namespace CGAL {
+
+  template <class SearchTraits,
+            class Distance_= typename internal::Spatial_searching_default_distance<SearchTraits>::type,
+            class Splitter_ = Sliding_midpoint<SearchTraits>,
+            class Tree_= Kd_tree<SearchTraits, Splitter_, Tag_true> >
+  class Orthogonal_incremental_neighbor_search {
+
+  public:
+    typedef Splitter_ Splitter;
+    typedef Tree_  Tree;
+    typedef Distance_ Distance;
+    typedef typename SearchTraits::Point_d Point_d;
+    typedef typename Distance::Query_item Query_item;
+    typedef typename SearchTraits::FT FT;
+    typedef typename Tree::Point_d_iterator Point_d_iterator;
+    typedef typename Tree::Node_const_handle Node_const_handle;
+
+    typedef std::pair<Point_d,FT> Point_with_transformed_distance;
+    typedef CGAL::cpp11::tuple<Node_const_handle,FT,std::vector<FT> > Node_with_distance;
+    typedef std::vector<Node_with_distance*> Node_with_distance_vector;
+    typedef std::vector<Point_with_transformed_distance*> Point_with_transformed_distance_vector;
+
+    template<class T>
+    struct Object_wrapper
+    {
+      T object;
+      Object_wrapper(const T& t):object(t){}
+      const T& operator* () const { return object; }
+      const T* operator-> () const { return &object; }
+    };
+
+    class Iterator_implementation {
+      SearchTraits traits;
+    public:
+
+      int number_of_neighbours_computed;
+      int number_of_internal_nodes_visited;
+      int number_of_leaf_nodes_visited;
+      int number_of_items_visited;
+
+    private:
+
+      typedef std::vector<FT> Distance_vector;
+
+      Distance_vector dists;
+
+      Distance Orthogonal_distance_instance;
+
+      FT multiplication_factor;
+
+      Query_item query_point;
+
+      FT distance_to_root;
+
+      bool search_nearest_neighbour;
+
+      FT rd;
+
+
+      class Priority_higher {
+      public:
+
+        bool search_nearest;
+
+        Priority_higher(bool search_the_nearest_neighbour)
+          : search_nearest(search_the_nearest_neighbour)
+        {}
+
+        //highest priority is smallest distance
+        bool
+        operator() (Node_with_distance* n1, Node_with_distance* n2) const
+        {
+          return (search_nearest) ? (CGAL::cpp11::get<1>(*n1) > CGAL::cpp11::get<1>(*n2)) : (CGAL::cpp11::get<1>(*n2) > CGAL::cpp11::get<1>(*n1));
+        }
+      };
+
+      class Distance_smaller {
+
+      public:
+
+        bool search_nearest;
+
+        Distance_smaller(bool search_the_nearest_neighbour)
+          : search_nearest(search_the_nearest_neighbour)
+        {}
+
+        //highest priority is smallest distance
+        bool operator() (Point_with_transformed_distance* p1, Point_with_transformed_distance* p2) const
+        {
+          return (search_nearest) ? (p1->second > p2->second) : (p2->second > p1->second);
+        }
+      };
+
+
+      std::priority_queue<Node_with_distance*, Node_with_distance_vector,
+                          Priority_higher> PriorityQueue;
+
+    public:
+      std::priority_queue<Point_with_transformed_distance*, Point_with_transformed_distance_vector,
+                          Distance_smaller> Item_PriorityQueue;
+
+
+    public:
+
+      int reference_count;
+
+
+
+      // constructor
+      Iterator_implementation(const Tree& tree,const Query_item& q, const Distance& tr,
+                              FT Eps=FT(0.0), bool search_nearest=true)
+        : traits(tree.traits()),number_of_neighbours_computed(0), number_of_internal_nodes_visited(0),
+        number_of_leaf_nodes_visited(0), number_of_items_visited(0),
+        Orthogonal_distance_instance(tr), multiplication_factor(Orthogonal_distance_instance.transformed_distance(FT(1.0)+Eps)),
+        query_point(q), search_nearest_neighbour(search_nearest),
+        PriorityQueue(Priority_higher(search_nearest)), Item_PriorityQueue(Distance_smaller(search_nearest)),
+        reference_count(1)
+
+
+      {
+        if (tree.empty()) return;
+
+        typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits.construct_cartesian_const_iterator_d_object();
+        int dim = static_cast<int>(std::distance(ccci(q), ccci(q,0)));
+
+        dists.resize(dim);
+        for(int i=0 ; i<dim ; ++i){
+          dists[i] = 0;
+        }
+
+        if (search_nearest){
+          distance_to_root=
+            Orthogonal_distance_instance.min_distance_to_rectangle(q, tree.bounding_box(),dists);
+          Node_with_distance *The_Root = new Node_with_distance(tree.root(),
+                                                                distance_to_root, dists);
+          PriorityQueue.push(The_Root);
+
+          // rd is the distance of the top of the priority queue to q
+          rd=CGAL::cpp11::get<1>(*The_Root);
+          Compute_the_next_nearest_neighbour();
+        }
+         else{
+           distance_to_root=
+         Orthogonal_distance_instance.max_distance_to_rectangle(q,
+                                                tree.bounding_box(), dists);
+           Node_with_distance *The_Root = new Node_with_distance(tree.root(),
+                                                                 distance_to_root, dists);
+        PriorityQueue.push(The_Root);
+
+        // rd is the distance of the top of the priority queue to q
+        rd=CGAL::cpp11::get<1>(*The_Root);
+        Compute_the_next_furthest_neighbour();
+         }
+
+
+      }
+
+      // * operator
+      const Point_with_transformed_distance&
+      operator* () const
+      {
+        return *(Item_PriorityQueue.top());
+      }
+
+      // prefix operator
+      Iterator_implementation&
+      operator++()
+      {
+        Delete_the_current_item_top();
+        if(search_nearest_neighbour)
+          Compute_the_next_nearest_neighbour();
+        else
+          Compute_the_next_furthest_neighbour();
+        return *this;
+      }
+
+      // postfix operator
+      Object_wrapper<Point_with_transformed_distance>
+      operator++(int)
+      {
+        Object_wrapper<Point_with_transformed_distance> result( *(Item_PriorityQueue.top()) );
+        ++*this;
+        return result;
+      }
+
+      // Print statistics of the general priority search process.
+      std::ostream&
+      statistics (std::ostream& s) const {
+        s << "Orthogonal priority search statistics:"
+          << std::endl;
+        s << "Number of internal nodes visited:"
+          << number_of_internal_nodes_visited << std::endl;
+        s << "Number of leaf nodes visited:"
+          << number_of_leaf_nodes_visited << std::endl;
+        s << "Number of items visited:"
+          << number_of_items_visited << std::endl;
+        s << "Number of neighbours computed:"
+          << number_of_neighbours_computed << std::endl;
+        return s;
+      }
+
+
+      //destructor
+      ~Iterator_implementation()
+      {
+        while (!PriorityQueue.empty()) {
+          Node_with_distance* The_top=PriorityQueue.top();
+          PriorityQueue.pop();
+          delete The_top;
+        }
+        while (!Item_PriorityQueue.empty()) {
+          Point_with_transformed_distance* The_top=Item_PriorityQueue.top();
+          Item_PriorityQueue.pop();
+          delete The_top;
+        }
+      }
+
+    private:
+
+      void
+      Delete_the_current_item_top()
+      {
+        Point_with_transformed_distance* The_item_top=Item_PriorityQueue.top();
+        Item_PriorityQueue.pop();
+        delete The_item_top;
+      }
+
+      void
+      Compute_the_next_nearest_neighbour()
+      {
+        // compute the next item
+        bool next_neighbour_found=false;
+        if (!(Item_PriorityQueue.empty())) {
+            next_neighbour_found=
+              (multiplication_factor*rd > Item_PriorityQueue.top()->second);
+        }
+        typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
+        typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point);
+        // otherwise browse the tree further
+        while ((!next_neighbour_found) && (!PriorityQueue.empty())) {
+          Node_with_distance* The_node_top=PriorityQueue.top();
+          Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top);
+          dists = CGAL::cpp11::get<2>(*The_node_top);
+          PriorityQueue.pop();
+          delete The_node_top;
+          FT copy_rd=rd;
+          while (!(N->is_leaf())) { // compute new distance
+            typename Tree::Internal_node_const_handle node =
+            static_cast<typename Tree::Internal_node_const_handle>(N);
+            number_of_internal_nodes_visited++;
+            int new_cut_dim=node->cutting_dimension();
+            FT new_rd,dst = dists[new_cut_dim];
+            FT val = *(query_point_it + new_cut_dim);
+            FT diff1 = val - node->upper_low_value();
+            FT diff2 = val - node->lower_high_value();
+            if (diff1 + diff2 < FT(0.0)) {
+              new_rd=
+                Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim);
+
+              CGAL_assertion(new_rd >= copy_rd);
+                dists[new_cut_dim] = diff1;
+                Node_with_distance *Upper_Child =
+                  new Node_with_distance(node->upper(), new_rd, dists);
+                PriorityQueue.push(Upper_Child);
+                dists[new_cut_dim] = dst;
+                N=node->lower();
+
+            }
+            else { // compute new distance
+              new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim);
+              CGAL_assertion(new_rd >= copy_rd);
+                dists[new_cut_dim] = diff2;
+                Node_with_distance *Lower_Child =
+                  new Node_with_distance(node->lower(), new_rd, dists);
+                PriorityQueue.push(Lower_Child);
+                dists[new_cut_dim] = dst;
+                N=node->upper();
+            }
+          }
+          // n is a leaf
+          typename Tree::Leaf_node_const_handle node =
+            static_cast<typename Tree::Leaf_node_const_handle>(N);
+          number_of_leaf_nodes_visited++;
+          if (node->size() > 0) {
+            for (typename Tree::iterator it=node->begin(); it != node->end(); it++) {
+              number_of_items_visited++;
+              FT distance_to_query_point=
+                Orthogonal_distance_instance.transformed_distance(query_point,*it);
+              Point_with_transformed_distance *NN_Candidate=
+                new Point_with_transformed_distance(*it,distance_to_query_point);
+              Item_PriorityQueue.push(NN_Candidate);
+            }
+            // old top of PriorityQueue has been processed,
+            // hence update rd
+
+            if (!(PriorityQueue.empty()))  {
+              rd = CGAL::cpp11::get<1>(*PriorityQueue.top());
+                next_neighbour_found =
+                  (multiplication_factor*rd >
+                   Item_PriorityQueue.top()->second);
+            }
+            else // priority queue empty => last neighbour found
+              {
+                next_neighbour_found=true;
+              }
+
+            number_of_neighbours_computed++;
+          }
+        }   // next_neighbour_found or priority queue is empty
+        // in the latter case also the item priority quee is empty
+      }
+
+
+      void
+      Compute_the_next_furthest_neighbour()
+      {
+        // compute the next item
+        bool next_neighbour_found=false;
+        if (!(Item_PriorityQueue.empty())) {
+            next_neighbour_found=
+              (rd < multiplication_factor*Item_PriorityQueue.top()->second);
+        }
+        typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
+        typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point);
+        // otherwise browse the tree further
+        while ((!next_neighbour_found) && (!PriorityQueue.empty())) {
+          Node_with_distance* The_node_top=PriorityQueue.top();
+          Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top);
+          dists = CGAL::cpp11::get<2>(*The_node_top);
+          PriorityQueue.pop();
+          delete The_node_top;
+          FT copy_rd=rd;
+          while (!(N->is_leaf())) { // compute new distance
+            typename Tree::Internal_node_const_handle node =
+              static_cast<typename Tree::Internal_node_const_handle>(N);
+            number_of_internal_nodes_visited++;
+            int new_cut_dim=node->cutting_dimension();
+            FT new_rd,dst = dists[new_cut_dim];
+            FT val = *(query_point_it + new_cut_dim);
+            FT diff1 = val - node->upper_low_value();
+            FT diff2 = val - node->lower_high_value();
+            if (diff1 + diff2 < FT(0.0)) {
+              diff1 = val - node->upper_high_value();
+              new_rd=
+                Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim);
+                Node_with_distance *Lower_Child =
+                  new Node_with_distance(node->lower(), copy_rd, dists);
+                PriorityQueue.push(Lower_Child);
+                N=node->upper();
+                dists[new_cut_dim] = diff1;
+                copy_rd=new_rd;
+
+            }
+            else { // compute new distance
+              diff2 = val - node->lower_low_value();
+              new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim);
+                Node_with_distance *Upper_Child =
+                  new Node_with_distance(node->upper(), copy_rd, dists);
+                PriorityQueue.push(Upper_Child);
+                N=node->lower();
+                dists[new_cut_dim] = diff2;
+                copy_rd=new_rd;
+            }
+          }
+          // n is a leaf
+          typename Tree::Leaf_node_const_handle node =
+            static_cast<typename Tree::Leaf_node_const_handle>(N);
+          number_of_leaf_nodes_visited++;
+          if (node->size() > 0) {
+            for (typename Tree::iterator it=node->begin(); it != node->end(); it++) {
+              number_of_items_visited++;
+              FT distance_to_query_point=
+                Orthogonal_distance_instance.transformed_distance(query_point,*it);
+              Point_with_transformed_distance *NN_Candidate=
+                new Point_with_transformed_distance(*it,distance_to_query_point);
+              Item_PriorityQueue.push(NN_Candidate);
+            }
+            // old top of PriorityQueue has been processed,
+            // hence update rd
+
+            if (!(PriorityQueue.empty()))  {
+              rd = CGAL::cpp11::get<1>(*PriorityQueue.top());
+                next_neighbour_found =
+                  (multiplication_factor*rd <
+                   Item_PriorityQueue.top()->second);
+            }
+            else // priority queue empty => last neighbour found
+              {
+                next_neighbour_found=true;
+              }
+
+            number_of_neighbours_computed++;
+          }
+        }   // next_neighbour_found or priority queue is empty
+        // in the latter case also the item priority quee is empty
+      }
+    }; // class Iterator_implementaion
+
+
+
+
+
+
+
+
+
+  public:
+    class iterator;
+    typedef iterator const_iterator;
+
+    // constructor
+    Orthogonal_incremental_neighbor_search(const Tree& tree,
+                                           const Query_item& q, FT Eps = FT(0.0),
+                                           bool search_nearest=true, const Distance& tr=Distance())
+      : m_tree(tree),m_query(q),m_dist(tr),m_Eps(Eps),m_search_nearest(search_nearest)
+    {}
+
+    iterator
+    begin() const
+    {
+      return iterator(m_tree,m_query,m_dist,m_Eps,m_search_nearest);
+    }
+
+    iterator
+    end() const
+    {
+      return iterator();
+    }
+
+    std::ostream&
+    statistics(std::ostream& s)
+    {
+      begin()->statistics(s);
+      return s;
+    }
+
+
+
+
+    class iterator {
+
+    public:
+
+      typedef std::input_iterator_tag iterator_category;
+      typedef Point_with_transformed_distance       value_type;
+      typedef Point_with_transformed_distance*      pointer;
+      typedef const Point_with_transformed_distance&      reference;
+      typedef std::size_t               size_type;
+      typedef std::ptrdiff_t            difference_type;
+      typedef int distance_type;
+
+      //class Iterator_implementation;
+      Iterator_implementation *Ptr_implementation;
+
+
+    public:
+
+      // default constructor
+      iterator()
+      : Ptr_implementation(0)
+      {}
+
+      int
+      the_number_of_items_visited()
+      {
+        return Ptr_implementation->number_of_items_visited;
+      }
+
+      // constructor
+      iterator(const Tree& tree,const Query_item& q, const Distance& tr=Distance(), FT eps=FT(0.0),
+               bool search_nearest=true)
+        : Ptr_implementation(new Iterator_implementation(tree, q, tr, eps, search_nearest))
+        {}
+
+      // copy constructor
+      iterator(const iterator& Iter)
+      {
+        Ptr_implementation = Iter.Ptr_implementation;
+        if (Ptr_implementation != 0) Ptr_implementation->reference_count++;
+      }
+
+      iterator& operator=(const iterator& Iter)
+      {
+        if (Ptr_implementation != Iter.Ptr_implementation){
+          if (Ptr_implementation != 0 && --(Ptr_implementation->reference_count)==0) {
+              delete Ptr_implementation;
+          }
+          Ptr_implementation = Iter.Ptr_implementation;
+          if (Ptr_implementation != 0) Ptr_implementation->reference_count++;
+        }
+        return *this;
+      }
+
+
+      const Point_with_transformed_distance&
+      operator* () const
+      {
+        return *(*Ptr_implementation);
+      }
+
+      // -> operator
+      const Point_with_transformed_distance*
+      operator-> () const
+      {
+        return &*(*Ptr_implementation);
+      }
+
+      // prefix operator
+      iterator&
+      operator++()
+      {
+        ++(*Ptr_implementation);
+        return *this;
+      }
+
+      // postfix operator
+      Object_wrapper<Point_with_transformed_distance>
+      operator++(int)
+      {
+        return (*Ptr_implementation)++;
+      }
+
+
+      bool
+      operator==(const iterator& It) const
+      {
+        if (
+            ((Ptr_implementation == 0) ||
+             Ptr_implementation->Item_PriorityQueue.empty()) &&
+            ((It.Ptr_implementation == 0) ||
+             It.Ptr_implementation->Item_PriorityQueue.empty())
+            )
+          return true;
+        // else
+        return (Ptr_implementation == It.Ptr_implementation);
+      }
+
+      bool
+      operator!=(const iterator& It) const
+      {
+        return !(*this == It);
+      }
+
+      std::ostream&
+      statistics (std::ostream& s)
+      {
+        Ptr_implementation->statistics(s);
+        return s;
+      }
+
+      ~iterator()
+      {
+        if (Ptr_implementation != 0) {
+          Ptr_implementation->reference_count--;
+          if (Ptr_implementation->reference_count==0) {
+            delete Ptr_implementation;
+            Ptr_implementation = 0;
+          }
+        }
+      }
+
+
+    }; // class iterator
+
+    //data members
+    const Tree& m_tree;
+    Query_item m_query;
+    Distance m_dist;
+    FT m_Eps;
+    bool m_search_nearest;
+  }; // class
+
+  template <class Traits, class Query_item, class Distance>
+  void swap (typename Orthogonal_incremental_neighbor_search<Traits,
+             Query_item, Distance>::iterator& x,
+             typename Orthogonal_incremental_neighbor_search<Traits,
+             Query_item, Distance>::iterator& y)
+  {
+    typename Orthogonal_incremental_neighbor_search<Traits,
+      Query_item, Distance>::iterator::Iterator_implementation
+      *tmp = x.Ptr_implementation;
+    x.Ptr_implementation  = y.Ptr_implementation;
+    y.Ptr_implementation = tmp;
+  }
+
+} // namespace CGAL
+
+#endif // CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH_H
diff --git a/src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt b/src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt
new file mode 100644
index 00000000..5b9581a0
--- /dev/null
+++ b/src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt
@@ -0,0 +1,82 @@
+CGAL/Regular_triangulation_traits_adapter.h
+CGAL/Triangulation_ds_vertex.h
+CGAL/Triangulation_data_structure.h
+CGAL/transforming_pair_iterator.h
+CGAL/NewKernel_d/static_int.h
+CGAL/NewKernel_d/Cartesian_LA_functors.h
+CGAL/NewKernel_d/Cartesian_change_FT.h
+CGAL/NewKernel_d/Wrapper/Vector_d.h
+CGAL/NewKernel_d/Wrapper/Hyperplane_d.h
+CGAL/NewKernel_d/Wrapper/Ref_count_obj.h
+CGAL/NewKernel_d/Wrapper/Cartesian_wrap.h
+CGAL/NewKernel_d/Wrapper/Point_d.h
+CGAL/NewKernel_d/Wrapper/Segment_d.h
+CGAL/NewKernel_d/Wrapper/Weighted_point_d.h
+CGAL/NewKernel_d/Wrapper/Sphere_d.h
+CGAL/NewKernel_d/Cartesian_per_dimension.h
+CGAL/NewKernel_d/Kernel_object_converter.h
+CGAL/NewKernel_d/KernelD_converter.h
+CGAL/NewKernel_d/Vector/sse2.h
+CGAL/NewKernel_d/Vector/avx4.h
+CGAL/NewKernel_d/Vector/determinant_of_vectors_small_dim_internal.h
+CGAL/NewKernel_d/Vector/determinant_of_iterator_to_points_from_points.h
+CGAL/NewKernel_d/Vector/determinant_of_points_from_vectors.h
+CGAL/NewKernel_d/Vector/array.h
+CGAL/NewKernel_d/Vector/determinant_of_iterator_to_points_from_iterator_to_vectors.h
+CGAL/NewKernel_d/Vector/determinant_of_iterator_to_vectors_from_vectors.h
+CGAL/NewKernel_d/Vector/determinant_of_vectors_small_dim.h
+CGAL/NewKernel_d/Vector/vector.h
+CGAL/NewKernel_d/Vector/v2int.h
+CGAL/NewKernel_d/Vector/mix.h
+CGAL/NewKernel_d/Cartesian_static_filters.h
+CGAL/NewKernel_d/Cartesian_LA_base.h
+CGAL/NewKernel_d/Lazy_cartesian.h
+CGAL/NewKernel_d/Coaffine.h
+CGAL/NewKernel_d/store_kernel.h
+CGAL/NewKernel_d/Dimension_base.h
+CGAL/NewKernel_d/Kernel_3_interface.h
+CGAL/NewKernel_d/Cartesian_complete.h
+CGAL/NewKernel_d/Cartesian_base.h
+CGAL/NewKernel_d/Cartesian_filter_K.h
+CGAL/NewKernel_d/functor_tags.h
+CGAL/NewKernel_d/Filtered_predicate2.h
+CGAL/NewKernel_d/functor_properties.h
+CGAL/NewKernel_d/Define_kernel_types.h
+CGAL/NewKernel_d/LA_eigen/LA.h
+CGAL/NewKernel_d/LA_eigen/constructors.h
+CGAL/NewKernel_d/Types/Aff_transformation.h
+CGAL/NewKernel_d/Types/Sphere.h
+CGAL/NewKernel_d/Types/Hyperplane.h
+CGAL/NewKernel_d/Types/Line.h
+CGAL/NewKernel_d/Types/Ray.h
+CGAL/NewKernel_d/Types/Iso_box.h
+CGAL/NewKernel_d/Types/Weighted_point.h
+CGAL/NewKernel_d/Types/Segment.h
+CGAL/NewKernel_d/Kernel_d_interface.h
+CGAL/NewKernel_d/utils.h
+CGAL/NewKernel_d/Kernel_2_interface.h
+CGAL/NewKernel_d/Cartesian_filter_NT.h
+CGAL/NewKernel_d/function_objects_cartesian.h
+CGAL/Convex_hull.h
+CGAL/Triangulation_ds_full_cell.h
+CGAL/Regular_triangulation.h
+CGAL/Epick_d.h
+CGAL/transforming_iterator.h
+CGAL/iterator_from_indices.h
+CGAL/Delaunay_triangulation.h
+CGAL/IO/Triangulation_off_ostream.h
+CGAL/typeset.h
+CGAL/Triangulation_full_cell.h
+CGAL/Triangulation.h
+CGAL/internal/Static_or_dynamic_array.h
+CGAL/internal/Combination_enumerator.h
+CGAL/internal/Triangulation/utilities.h
+CGAL/internal/Triangulation/Triangulation_ds_iterators.h
+CGAL/internal/Triangulation/Dummy_TDS.h
+CGAL/argument_swaps.h
+CGAL/Epeck_d.h
+CGAL/determinant_of_vectors.h
+CGAL/TDS_full_cell_default_storage_policy.h
+CGAL/TDS_full_cell_mirror_storage_policy.h
+CGAL/Triangulation_face.h
+CGAL/Triangulation_vertex.h