save

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

Former-commit-id: 4be2335267f7ad5253cde7f0573199013034c00e

5 files changed, 268 insertions, 91 deletions

diff --git a/src/Bipartite_graphs_matching/include/gudhi/Graph_matching.h b/src/Bipartite_graphs_matching/include/gudhi/Graph_matching.h
index a2754333..228de584 100644
--- a/src/Bipartite_graphs_matching/include/gudhi/Graph_matching.h
+++ b/src/Bipartite_graphs_matching/include/gudhi/Graph_matching.h
@@ -65,7 +65,7 @@ private:
     std::list<int> unmatched_in_u;
 
     /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */
-    std::unique_ptr<Layered_neighbors_finder> layering() const;
+    std::shared_ptr<Layered_neighbors_finder> layering() const;
     /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */
     bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth);
     /** \internal \brief Update the matching with the simple augmenting path given as parameter. */
@@ -138,20 +138,20 @@ inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u
     return true;
 }
 
-inline std::unique_ptr<Layered_neighbors_finder> Graph_matching::layering() const {
+inline std::shared_ptr<Layered_neighbors_finder> Graph_matching::layering() const {
     std::list<int> u_vertices(unmatched_in_u);
     std::list<int> v_vertices;
     Neighbors_finder nf(r);
     for (int v_point_index = 0; v_point_index < G::size(); ++v_point_index)
         nf.add(v_point_index);
-    std::unique_ptr<Layered_neighbors_finder> layered_nf(new Layered_neighbors_finder(r));
+    std::shared_ptr<Layered_neighbors_finder> layered_nf(new Layered_neighbors_finder(r));
     for(int layer = 0; !u_vertices.empty(); layer++) {
         // one layer is one step in the BFS
-        for (auto it = u_vertices.cbegin(); it != u_vertices.cend(); ++it) {
-            std::unique_ptr< std::list<int> > u_succ = std::move(nf.pull_all_near(*it));
-            for (auto it = u_succ->cbegin(); it != u_succ->cend(); ++it) {
-                layered_nf->add(*it, layer);
-                v_vertices.emplace_back(*it);
+        for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) {
+            std::shared_ptr<std::list<int>> u_succ(nf.pull_all_near(*it1));
+            for (auto it2 = u_succ->begin(); it2 != u_succ->end(); ++it2) {
+                layered_nf->add(*it2, layer);
+                v_vertices.emplace_back(*it2);
             }
         }
         // When the above for finishes, we have progress of one half-step (from U to V) in the BFS
@@ -183,7 +183,7 @@ inline void Graph_matching::update(std::deque<int>& path) {
 template<typename Persistence_diagram1, typename Persistence_diagram2>
 double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) {
     G::initialize(diag1, diag2, e);
-    std::unique_ptr< std::vector<double> > sd = std::move(G::sorted_distances());
+    std::shared_ptr< std::vector<double> > sd(G::sorted_distances());
     int idmin = 0;
     int idmax = sd->size() - 1;
     // alpha can be modified, this will change the complexity
diff --git a/src/Bipartite_graphs_matching/include/gudhi/Internal_point.h b/src/Bipartite_graphs_matching/include/gudhi/Internal_point.h
new file mode 100644
index 00000000..633a3767
--- /dev/null
+++ b/src/Bipartite_graphs_matching/include/gudhi/Internal_point.h
@@ -0,0 +1,134 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
+ *
+ *    This program is free software: 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.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_INTERNAL_POINT_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_INTERNAL_POINT_H_
+
+//namespace Gudhi {
+
+//namespace bipartite_graph_matching {
+
+/** \internal \brief Returns the used index for encoding none of the points */
+int null_point_index();
+
+/** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */
+struct Internal_point {
+    double vec[2];
+    int point_index;
+    Internal_point() { vec[0]= vec[1] = 0.; point_index = null_point_index(); }
+    Internal_point(double x, double y) { vec[0]=x; vec[1]=y; point_index = null_point_index(); }
+    double x() const { return vec[ 0 ]; }
+    double y() const { return vec[ 1 ]; }
+    double& x() { return vec[ 0 ]; }
+    double& y() { return vec[ 1 ]; }
+    bool operator==(const Internal_point& p) const
+    {
+        return (x() == p.x()) && (y() == p.y());
+    }
+    bool  operator!=(const Internal_point& p) const { return ! (*this == p); }
+};
+
+namespace CGAL {
+template <>
+struct Kernel_traits<Internal_point> {
+    struct Kernel {
+        typedef double FT;
+        typedef double RT;
+    };
+};
+}
+
+struct Construct_coord_iterator {
+    typedef  const double* result_type;
+    const double* operator()(const Internal_point& p) const
+    { return static_cast<const double*>(p.vec); }
+    const double* operator()(const Internal_point& p, int)  const
+    { return static_cast<const double*>(p.vec+2); }
+};
+/*
+struct Distance {
+    typedef Internal_point Query_item;
+    typedef Internal_point Point_d;
+    typedef double FT;
+    typedef CGAL::Dimension_tag<2> D;
+
+    FT transformed_distance(const Query_item& q, const Point_d& p) const {
+        FT d0= std::abs(q.x()-p.x());
+        FT d1= std::abs(q.y()-p.y());
+        return std::max(d0,d1);
+    }
+   FT min_distance_to_rectangle(const Query_item& q,
+                                     const CGAL::Kd_tree_rectangle<FT,D>& b) const {
+        FT d0(0.), d1(0.);
+        if (q.x() < b.min_coord(0))
+            d0 += (b.min_coord(0)-q.x());
+        if (q.x() > b.max_coord(0))
+            d0 += (q.x()-b.max_coord(0));
+        if (q.y() < b.min_coord(1))
+            d1 += (b.min_coord(1)-q.y());
+        if (q.y() > b.max_coord(1))
+            d1 += (q.y()-b.max_coord(1));
+        return std::max(d0,d1);
+    }
+    FT min_distance_to_rectangle(const Query_item& q, const CGAL::Kd_tree_rectangle<FT,D>& b,std::vector<FT>& dists){
+        dists[0] = dists[1] = 0.;
+        if (q.x() < b.min_coord(0))
+            dists[0] = (b.min_coord(0)- q.x());
+        if (q.x() > b.max_coord(0))
+            dists[0] = (q.x()-b.max_coord(0));
+        if (q.y() < b.min_coord(1))
+            dists[1] = (b.min_coord(1)-q.y());
+        if (q.y() > b.max_coord(1))
+            dists[1] = (q.y()-b.max_coord(1));
+        return std::max(dists[0],dists[1]);
+    }
+    FT max_distance_to_rectangle(const Query_item& q, const CGAL::Kd_tree_rectangle<FT,D>& b) const {
+        FT d0 = (q.x() >= (b.min_coord(0)+b.max_coord(0))/2.) ?
+                    (q.x()-b.min_coord(0)) : (b.max_coord(0)-q.x());
+        FT d1 = (q.y() >= (b.min_coord(1)+b.max_coord(1))/2.) ?
+                    (q.y()-b.min_coord(1)) : (b.max_coord(1)-q.y());
+        return std::max(d0, d1);
+    }
+    FT max_distance_to_rectangle(const Query_item& q, const CGAL::Kd_tree_rectangle<FT,D>& b,std::vector<FT>& dists){
+        dists[0] = (q.x() >= (b.min_coord(0)+b.max_coord(0))/2.0) ?
+                    (q.x()-b.min_coord(0)) : (b.max_coord(0)-q.x());
+        dists[1] = (q.y() >= (b.min_coord(1)+b.max_coord(1))/2.0) ?
+                    (q.y()-b.min_coord(1)) : (b.max_coord(1)-q.y());
+        return std::max(dists[0], dists[1]);
+    }
+    FT new_distance(FT& dist, FT old_off, FT new_off,
+                        int)  const {
+        return dist + new_off - old_off; //works ?
+    }
+    FT transformed_distance(FT d) const { return d; }
+    FT inverse_of_transformed_distance(FT d) { return d; }
+}; // end of struct Distance
+*/
+inline int null_point_index() {
+    return -1;
+}
+
+//}  // namespace bipartite_graph_matching
+
+//}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_
diff --git a/src/Bipartite_graphs_matching/include/gudhi/Neighbors_finder.h b/src/Bipartite_graphs_matching/include/gudhi/Neighbors_finder.h
index 78b9debc..2ef7f768 100644
--- a/src/Bipartite_graphs_matching/include/gudhi/Neighbors_finder.h
+++ b/src/Bipartite_graphs_matching/include/gudhi/Neighbors_finder.h
@@ -47,7 +47,7 @@ public:
     /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */
     int pull_near(int u_point_index);
     /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */
-    std::unique_ptr< std::list<int> > pull_all_near(int u_point_index);
+    std::shared_ptr< std::list<int> > pull_all_near(int u_point_index);
 
 private:
     const double r;
@@ -77,7 +77,7 @@ public:
 
 private:
     const double r;
-    std::vector<Neighbors_finder> neighbors_finder;
+    std::vector<std::shared_ptr<Neighbors_finder>> neighbors_finder;
 };
 
 inline Neighbors_finder::Neighbors_finder(double r) :
@@ -119,8 +119,8 @@ inline int Neighbors_finder::pull_near(int u_point_index) {
     return tmp;
 }
 
-inline std::unique_ptr< std::list<int> > Neighbors_finder::pull_all_near(int u_point_index) {
-    std::unique_ptr< std::list<int> > all_pull = std::move(planar_neighbors_f.pull_all_near(u_point_index));
+inline std::shared_ptr< std::list<int> > Neighbors_finder::pull_all_near(int u_point_index) {
+    std::shared_ptr< std::list<int> > all_pull(planar_neighbors_f.pull_all_near(u_point_index));
     int last_pull = pull_near(u_point_index);
     while (last_pull != null_point_index()) {
         all_pull->emplace_back(last_pull);
@@ -134,14 +134,14 @@ inline Layered_neighbors_finder::Layered_neighbors_finder(double r) :
 
 inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) {
     for (int l = neighbors_finder.size(); l <= vlayer; l++)
-        neighbors_finder.emplace_back(Neighbors_finder(r));
-    neighbors_finder.at(vlayer).add(v_point_index);
+        neighbors_finder.emplace_back(std::shared_ptr<Neighbors_finder>(new Neighbors_finder(r)));
+    neighbors_finder.at(vlayer)->add(v_point_index);
 }
 
 inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) {
     if (static_cast<int> (neighbors_finder.size()) <= vlayer)
         return null_point_index();
-    return neighbors_finder.at(vlayer).pull_near(u_point_index);
+    return neighbors_finder.at(vlayer)->pull_near(u_point_index);
 }
 
 inline int Layered_neighbors_finder::vlayers_number() const {
diff --git a/src/Bipartite_graphs_matching/include/gudhi/Persistence_diagrams_graph.h b/src/Bipartite_graphs_matching/include/gudhi/Persistence_diagrams_graph.h
index 6420b772..2ba960a7 100644
--- a/src/Bipartite_graphs_matching/include/gudhi/Persistence_diagrams_graph.h
+++ b/src/Bipartite_graphs_matching/include/gudhi/Persistence_diagrams_graph.h
@@ -30,13 +30,12 @@
 #include <algorithm>
 #include <math.h>
 #include <memory>
+#include "Internal_point.h"
 
 namespace Gudhi {
 
 namespace bipartite_graph_matching {
 
-/** \internal \brief Returns the used index for encoding none of the points */
-int null_point_index();
 
 /** \internal \brief Structure representing an euclidean bipartite graph containing
  *  the points from the two persistence diagrams (including the projections).
@@ -61,30 +60,24 @@ public:
     /** \internal \brief Returns size = |U| = |V|. */
     static int size();
     /** \internal \brief Returns the O(n^2) sorted distances between the points. */
-    static std::unique_ptr< std::vector<double> > sorted_distances();
+    static std::shared_ptr< std::vector<double> > sorted_distances();
 
 private:
-    /** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */
-    typedef std::pair<double, double> Internal_point;
     static std::vector<Internal_point> u;
     static std::vector<Internal_point> v;
     static Internal_point get_u_point(int u_point_index);
     static Internal_point get_v_point(int v_point_index);
 
     friend class Naive_pnf;
-    friend class Upper_envelope_tree;
+    friend class Cgal_pnf;
 };
 
 /** \internal \typedef \brief Shorter alias */
 typedef Persistence_diagrams_graph G;
 
-// static initialization, seems to work but strange
-std::vector<G::Internal_point> G::u = [] {return std::vector<G::Internal_point>();}();
-std::vector<G::Internal_point> G::v = [] {return std::vector<G::Internal_point>();}();
-
-inline int null_point_index() {
-    return -1;
-}
+// static initialization
+std::vector<Internal_point> G::u = [] {return std::vector<Internal_point>();}();
+std::vector<Internal_point> G::v = [] {return std::vector<Internal_point>();}();
 
 template<typename Persistence_diagram1, typename Persistence_diagram2>
 inline void G::initialize(const Persistence_diagram1 &diag1,
@@ -92,10 +85,10 @@ inline void G::initialize(const Persistence_diagram1 &diag1,
     u.clear();
     v.clear();
     for (auto it = diag1.cbegin(); it != diag1.cend(); ++it)
-        if (it->second - it->first > e)
+        if (it->y() - it->x() > e)
             u.emplace_back(*it);
     for (auto it = diag2.cbegin(); it != diag2.cend(); ++it)
-        if (it->second - it->first > e)
+        if (it->y() - it->x() > e)
             v.emplace_back(*it);
     if (u.size() < v.size())
         swap(u, v);
@@ -124,37 +117,37 @@ inline double G::distance(int u_point_index, int v_point_index) {
         return 0;
     Internal_point p_u = get_u_point(u_point_index);
     Internal_point p_v = get_v_point(v_point_index);
-    return std::max(std::fabs(p_u.first - p_v.first), std::fabs(p_u.second - p_v.second));
+    return std::max(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y()));
 }
 
 inline int G::size() {
     return static_cast<int> (u.size() + v.size());
 }
 
-inline std::unique_ptr< std::vector<double> > G::sorted_distances() {
+inline std::shared_ptr< std::vector<double> > G::sorted_distances() {
     // could be optimized
     std::set<double> sorted_distances;
     for (int u_point_index = 0; u_point_index < size(); ++u_point_index)
         for (int v_point_index = 0; v_point_index < size(); ++v_point_index)
             sorted_distances.emplace(distance(u_point_index, v_point_index));
-    std::unique_ptr< std::vector<double> > sd_up(new std::vector<double>(sorted_distances.cbegin(), sorted_distances.cend()));
+    std::shared_ptr< std::vector<double> > sd_up(new std::vector<double>(sorted_distances.cbegin(), sorted_distances.cend()));
     return sd_up;
 }
 
-inline G::Internal_point G::get_u_point(int u_point_index) {
+inline Internal_point G::get_u_point(int u_point_index) {
     if (!on_the_u_diagonal(u_point_index))
         return u.at(u_point_index);
     Internal_point projector = v.at(corresponding_point_in_v(u_point_index));
-    double x = (projector.first + projector.second) / 2;
-    return Internal_point(x, x);
+    double m = (projector.x() + projector.y()) / 2;
+    return Internal_point(m, m);
 }
 
-inline G::Internal_point G::get_v_point(int v_point_index) {
+inline Internal_point G::get_v_point(int v_point_index) {
     if (!on_the_v_diagonal(v_point_index))
         return v.at(v_point_index);
     Internal_point projector = u.at(corresponding_point_in_u(v_point_index));
-    double x = (projector.first + projector.second) / 2;
-    return Internal_point(x, x);
+    double m = (projector.x() + projector.y()) / 2;
+    return Internal_point(m, m);
 }
 
 }  // namespace bipartite_graph_matching
diff --git a/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h b/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h
index 3564dcf9..47c1a8fc 100644
--- a/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h
+++ b/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h
@@ -25,9 +25,12 @@
 
 #include <list>
 #include <map>
-
+#include <CGAL/Search_traits.h>
+#include <CGAL/Orthogonal_incremental_neighbor_search.h>
+#include <CGAL/Weighted_Minkowski_distance.h>
 #include "Persistence_diagrams_graph.h"
 
+
 namespace Gudhi {
 
 namespace bipartite_graph_matching {
@@ -39,34 +42,39 @@ namespace bipartite_graph_matching {
  *
  * \ingroup bottleneck_distance
  */
-class Abstract_planar_neighbors_finder {
+class Naive_pnf {
 public:
-    /** \internal \brief Constructor TODO. */
-    Abstract_planar_neighbors_finder(double r);
-    virtual ~Abstract_planar_neighbors_finder() = 0;
+    /** \internal \brief Constructor taking the near distance definition as parameter. */
+    Naive_pnf(double r_);
     /** \internal \brief A point added will be possibly pulled. */
-    virtual void add(int v_point_index) = 0;
+    void add(int v_point_index);
     /** \internal \brief A point manually removed will no longer be possibly pulled. */
-    virtual void remove(int v_point_index) = 0;
+    void remove(int v_point_index);
     /** \internal \brief Can the point given as parameter be returned ? */
-    virtual bool contains(int v_point_index) const = 0;
+    bool contains(int v_point_index) const;
     /** \internal \brief Provide and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */
-    virtual int pull_near(int u_point_index) = 0;
+    int pull_near(int u_point_index);
     /** \internal \brief Provide and remove all the V points near to the U point given as parameter. */
-    virtual std::unique_ptr< std::list<int> > pull_all_near(int u_point_index);
+    virtual std::shared_ptr< std::list<int> > pull_all_near(int u_point_index);
 
-protected:
-    const double r;
+private:
+    double r;
+    std::pair<int,int> get_v_key(int v_point_index) const;
+    std::multimap<std::pair<int,int>,int> grid;
 };
 
-/** \internal \brief Naive_pnf is an naïve Abstract_planar_neighbors_finder implementation
- *
- * \ingroup bottleneck_distance
- */
-class Naive_pnf : public Abstract_planar_neighbors_finder {
+class Cgal_pnf {
+
+    typedef CGAL::Dimension_tag<2> D;
+    typedef CGAL::Search_traits<double, Internal_point, const double*, Construct_coord_iterator, D> Cgal_traits;
+    typedef CGAL::Weighted_Minkowski_distance<Cgal_traits> Distance;
+    typedef CGAL::Orthogonal_incremental_neighbor_search<Cgal_traits, Distance> K_neighbor_search;
+    typedef K_neighbor_search::Tree Kd_tree;
+
+
 public:
     /** \internal \brief Constructor taking the near distance definition as parameter. */
-    Naive_pnf(double r);
+    Cgal_pnf(double r_);
     /** \internal \brief A point added will be possibly pulled. */
     void add(int v_point_index);
     /** \internal \brief A point manually removed will no longer be possibly pulled. */
@@ -76,39 +84,24 @@ public:
     /** \internal \brief Provide and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */
     int pull_near(int u_point_index);
     /** \internal \brief Provide and remove all the V points near to the U point given as parameter. */
-    virtual std::unique_ptr< std::list<int> > pull_all_near(int u_point_index);
+    virtual std::shared_ptr< std::list<int> > pull_all_near(int u_point_index);
 
 private:
-    std::pair<int,int> get_v_key(int v_point_index) const;
-    std::multimap<std::pair<int,int>,int> grid;
+    double r;
+    std::set<int> contents;
+    Kd_tree kd_t;
 };
 
-/** \internal \typedef \brief Planar_neighbors_finder is the used Abstract_planar_neighbors_finder's implementation. */
-typedef Naive_pnf Planar_neighbors_finder;
-
-
-inline Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(double r) :
-    r(r) { }
-
-inline Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() {}
-
-inline std::unique_ptr< std::list<int> > Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) {
-    std::unique_ptr< std::list<int> > all_pull(new std::list<int>);
-    int last_pull = pull_near(u_point_index);
-    while (last_pull != null_point_index()) {
-        all_pull->emplace_back(last_pull);
-        last_pull = pull_near(u_point_index);
-    }
-    return all_pull;
-}
+/** \internal \typedef \brief Planar_neighbors_finder is the used implementation. */
+typedef Cgal_pnf Planar_neighbors_finder;
 
-inline Naive_pnf::Naive_pnf(double r) :
-    Abstract_planar_neighbors_finder(r), grid() { }
+inline Naive_pnf::Naive_pnf(double r_) :
+    r(r_), grid() { }
 
 
 inline std::pair<int,int> Naive_pnf::get_v_key(int v_point_index) const{
-    G::Internal_point v_point = G::get_v_point(v_point_index);
-    return std::make_pair(static_cast<int>(v_point.first/r), static_cast<int>(v_point.second/r));
+    Internal_point v_point = G::get_v_point(v_point_index);
+    return std::make_pair(static_cast<int>(v_point.x()/r), static_cast<int>(v_point.y()/r));
 }
 
 inline void Naive_pnf::add(int v_point_index) {
@@ -133,9 +126,9 @@ inline bool Naive_pnf::contains(int v_point_index) const {
 }
 
 inline int Naive_pnf::pull_near(int u_point_index) {
-    G::Internal_point u_point = G::get_u_point(u_point_index);
-    int i0 = static_cast<int>(u_point.first/r);
-    int j0 = static_cast<int>(u_point.second/r);
+    Internal_point u_point = G::get_u_point(u_point_index);
+    int i0 = static_cast<int>(u_point.x()/r);
+    int j0 = static_cast<int>(u_point.y()/r);
     for(int i = 1; i<= 3; i++)
         for(int j = 1; j<= 3; j++)
             for(auto it = grid.find(std::make_pair(i0 +(i%3)-1, j0+(j%3)-1)); it!=grid.end(); it++)
@@ -147,11 +140,11 @@ inline int Naive_pnf::pull_near(int u_point_index) {
     return null_point_index();
 }
 
-inline std::unique_ptr< std::list<int> > Naive_pnf::pull_all_near(int u_point_index) {
-    std::unique_ptr< std::list<int> > all_pull(new std::list<int>);
-    G::Internal_point u_point = G::get_u_point(u_point_index);
-    int i0 = static_cast<int>(u_point.first/r);
-    int j0 = static_cast<int>(u_point.second/r);
+inline std::shared_ptr< std::list<int> > Naive_pnf::pull_all_near(int u_point_index) {
+    std::shared_ptr< std::list<int> > all_pull(new std::list<int>);
+    Internal_point u_point = G::get_u_point(u_point_index);
+    int i0 = static_cast<int>(u_point.x()/r);
+    int j0 = static_cast<int>(u_point.y()/r);
     for(int i = 1; i<= 3; i++)
         for(int j = 1; j<= 3; j++)
             for(auto it = grid.find(std::make_pair(i0 +(i%3)-1, j0+(j%3)-1)); it!=grid.end(); it++)
@@ -163,6 +156,63 @@ inline std::unique_ptr< std::list<int> > Naive_pnf::pull_all_near(int u_point_in
     return all_pull;
 }
 
+
+/** \internal \brief Constructor taking the near distance definition as parameter. */
+inline Cgal_pnf::Cgal_pnf(double r_)
+    : r(r_), contents(), kd_t() {}
+
+
+/** \internal \brief A point added will be possibly pulled. */
+inline void Cgal_pnf::add(int v_point_index){
+    Internal_point v_point = G::get_v_point(v_point_index);
+    v_point.point_index = v_point_index;
+    kd_t.insert(v_point);
+    contents.insert(v_point_index);
+}
+
+/** \internal \brief A point manually removed will no longer be possibly pulled. */
+inline void Cgal_pnf::remove(int v_point_index){
+    contents.erase(v_point_index);
+}
+
+/** \internal \brief Can the point given as parameter be returned ? */
+inline bool Cgal_pnf::contains(int v_point_index) const{
+    return contents.count(v_point_index)>0;
+}
+
+/** \internal \brief Provide and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */
+inline int Cgal_pnf::pull_near(int u_point_index){
+    Internal_point u_point = G::get_u_point(u_point_index);
+    K_neighbor_search search(kd_t, u_point, 0., true, Distance(0.));
+    for (auto it = search.begin(); it != search.end(); it++)
+        if(contents.count(it->first.point_index)==0)
+            kd_t.remove(it->first);
+        else if(G::distance(u_point_index, it->first.point_index) > r){
+            for(auto itc=contents.cbegin(); itc != contents.cend(); itc++)
+                if(G::distance(u_point_index, *itc) <= r)
+                    std::cout << G::distance(u_point_index, *itc) << " ! > " << r << std::endl;
+            return null_point_index();
+        }
+        else
+        {
+            kd_t.remove(it->first);
+            contents.erase(it->first.point_index);
+            return it->first.point_index;
+        }
+    return null_point_index();
+}
+
+inline std::shared_ptr< std::list<int> > Cgal_pnf::pull_all_near(int u_point_index) {
+    std::shared_ptr< std::list<int> > all_pull(new std::list<int>);
+    int last_pull = pull_near(u_point_index);
+    while (last_pull != null_point_index()) {
+        all_pull->emplace_back(last_pull);
+        last_pull = pull_near(u_point_index);
+    }
+    return all_pull;
+}
+
+
 }  // namespace bipartite_graph_matching
 
 }  // namespace Gudhi