Ajout de Grid_cell.h

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

Former-commit-id: e1b4ccdcb482c992c444adc91d67531bfa3e90dc

6 files changed, 231 insertions, 28 deletions

diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h
index 17412e6c..2bcc6a61 100644
--- a/src/Bottleneck/include/gudhi/Graph_matching.h
+++ b/src/Bottleneck/include/gudhi/Graph_matching.h
@@ -74,13 +74,13 @@ private:
     void update(std::deque<int> & path);
 };
 
-Graph_matching::Graph_matching()
+inline Graph_matching::Graph_matching()
     : r(0.), v_to_u(G::size(), null_point_index()), unmatched_in_u() {
     for (int u_point_index = 0; u_point_index < G::size(); ++u_point_index)
         unmatched_in_u.emplace_back(u_point_index);
 }
 
-Graph_matching& Graph_matching::operator=(const Graph_matching& m) {
+inline Graph_matching& Graph_matching::operator=(const Graph_matching& m) {
     r = m.r;
     v_to_u = m.v_to_u;
     unmatched_in_u = m.unmatched_in_u;
diff --git a/src/Bottleneck/include/gudhi/Grid_cell.h b/src/Bottleneck/include/gudhi/Grid_cell.h
new file mode 100644
index 00000000..eee938cb
--- /dev/null
+++ b/src/Bottleneck/include/gudhi/Grid_cell.h
@@ -0,0 +1,188 @@
+/*    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_GRID_CELL_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRID_CELL_H_
+
+#include <list>
+#include <set>
+#include <map>
+
+#include "Persistence_diagrams_graph.h"
+
+namespace Gudhi {
+
+namespace bottleneck {
+
+/** \internal \brief TODO
+ *
+ * \ingroup bottleneck_distance
+ */
+class Grid_cell {
+public:
+    Grid_cell(double r);
+    void add(int v_point_index);
+    void remove(int v_point_index);
+    bool contains(int v_point_index) const;
+    int pull_center();
+    int pull_xi(int u_point_index);
+    int pull_xd(int u_point_index);
+    int pull_yi(int u_point_index);
+    int pull_yd(int u_point_index);
+    int pull_xi_yi(int u_point_index);
+    int pull_xi_yd(int u_point_index);
+    int pull_xd_yi(int u_point_index);
+    int pull_xd_yd(int u_point_index);
+
+private:
+    double r;
+    std::set<int, G::Compare_x> xi_order;
+    std::set<int, G::Compare_y> yi_order;
+    struct Hidden_points_tree{int point; std::list<Hidden_points_tree> hidden;};
+    typedef std::map<int, std::list<Hidden_points_tree>, G::Compare_x> Corner_tree;
+    Corner_tree xi_yi_order;
+    Corner_tree xi_yd_order;
+    Corner_tree xd_yi_order;
+    Corner_tree xd_yd_order;
+    void remove_aux(Corner_tree t, int v_point_index);
+    void build_xi_yi();
+    void build_xi_yd();
+    void build_xd_yi();
+    void build_xd_yd();
+};
+
+
+inline Grid_cell::Grid_cell(double r)
+    : r(r), xi_order(G::Compare_x(r)), yi_order(G::Compare_y(r)), xi_yi_order(G::Compare_x(r)),
+      xi_yd_order(G::Compare_x(r)), xd_yi_order(G::Compare_x(r)), xd_yd_order(G::Compare_x(r)) {}
+
+inline void Grid_cell::add(int v_point_index){
+    xi_order.emplace(v_point_index);
+}
+
+inline bool Grid_cell::contains(int v_point_index) const{
+    return (xi_order.count(v_point_index) > 0);
+}
+
+inline void Grid_cell::remove_aux(Corner_tree t, int v_point_index){
+    if(t.empty())
+        return;
+    std::list<Hidden_points_tree> hidden_points = t.at(v_point_index);
+    t.erase(v_point_index);
+    for(auto it = hidden_points.begin(); it != hidden_points.end(); ++it)
+        t.emplace(it->point,it->hidden);
+
+}
+
+inline void Grid_cell::remove(int v_point_index){
+    xi_order.erase(v_point_index);
+    yi_order.erase(v_point_index);
+    remove_aux(xi_yi_order,v_point_index);
+    remove_aux(xi_yd_order,v_point_index);
+    remove_aux(xd_yi_order,v_point_index);
+    remove_aux(xd_yd_order,v_point_index);
+}
+
+//factorization needed \/ \/ \/
+
+inline int Grid_cell::pull_center(){
+    if(xi_order.empty())
+        return null_point_index();
+    int v_point_index = *xi_order.begin();
+    remove(v_point_index);
+    return v_point_index;
+}
+
+inline int Grid_cell::pull_xi(int u_point_index){
+    if(xi_order.empty())
+        return null_point_index();
+    int v_point_index = *xi_order.begin(); //!
+    if(G::distance(u_point_index,v_point_index)<=r){
+        remove(v_point_index);
+        return v_point_index;
+    }
+    return null_point_index();
+}
+
+inline int Grid_cell::pull_xd(int u_point_index){
+    if(xi_order.empty())
+        return null_point_index();
+    int v_point_index = *xi_order.rbegin(); //!
+    if(G::distance(u_point_index,v_point_index)<=r){
+        remove(v_point_index);
+        return v_point_index;
+    }
+    return null_point_index();
+}
+
+inline int Grid_cell::pull_yi(int u_point_index){
+    if(xi_order.empty())
+        return null_point_index();
+    if(yi_order.empty())
+        for(auto it = xi_order.begin(); it!= xi_order.end(); ++it)
+            yi_order.emplace(*it);
+    int v_point_index = *yi_order.begin(); //!
+    if(G::distance(u_point_index,v_point_index)<=r){
+        remove(v_point_index);
+        return v_point_index;
+    }
+    return null_point_index();
+}
+
+inline int Grid_cell::pull_yd(int u_point_index){
+    if(xi_order.empty())
+        return null_point_index();
+    if(yi_order.empty())
+        for(auto it = xi_order.begin(); it!= xi_order.end(); ++it)
+            yi_order.emplace(*it);
+    int v_point_index = *yi_order.rbegin(); //!
+    if(G::distance(u_point_index,v_point_index)<=r){
+        remove(v_point_index);
+        return v_point_index;
+    }
+    return null_point_index();
+}
+
+inline int Grid_cell::pull_xi_yi(int u_point_index){
+    if(xi_order.empty())
+        return null_point_index();
+    if(xi_yi_order.empty())
+        build_xi_yi();
+    auto it = xi_yi_order.upper_bound(u_point_index+G::size());
+    if(it==xi_yi_order.cbegin()) //!
+        return null_point_index();
+    it--; //!
+    int v_point_index = it->first;
+    for(auto it2=it->second.begin();it2!=it->second.end();it2++)
+        xi_yi_order.emplace(it2->point,it2->hidden);
+    if(G::distance(u_point_index,v_point_index)<=r){
+        remove(v_point_index);
+        return v_point_index;
+    }
+    return null_point_index();
+}
+
+}  // namespace bottleneck
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_GRID_CELL_H_
diff --git a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
index e6b7f30d..58805b86 100644
--- a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
+++ b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
@@ -54,7 +54,7 @@ private:
     std::vector<Neighbors_finder> neighbors_finder;
 };
 
-Layered_neighbors_finder::Layered_neighbors_finder(double r) :
+inline Layered_neighbors_finder::Layered_neighbors_finder(double r) :
     r(r), neighbors_finder() { }
 
 inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) {
diff --git a/src/Bottleneck/include/gudhi/Neighbors_finder.h b/src/Bottleneck/include/gudhi/Neighbors_finder.h
index ac3f8600..be81877a 100644
--- a/src/Bottleneck/include/gudhi/Neighbors_finder.h
+++ b/src/Bottleneck/include/gudhi/Neighbors_finder.h
@@ -58,7 +58,7 @@ private:
     bool contains(int v_point_index);
 };
 
-Neighbors_finder::Neighbors_finder(double r) :
+inline Neighbors_finder::Neighbors_finder(double r) :
     r(r), planar_neighbors_f(r), projections_f() { }
 
 inline void Neighbors_finder::add(int v_point_index) {
diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
index 66e43145..aed328e2 100644
--- a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
+++ b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
@@ -28,6 +28,7 @@
 #include <cmath>
 #include <utility>
 #include <algorithm>
+#include <math.h>
 
 namespace Gudhi {
 
@@ -60,10 +61,10 @@ public:
     static int size();
     /** \internal \brief Returns the O(n^2) sorted distances between the points. */
     static std::unique_ptr< std::vector<double> > sorted_distances();
-    /** \internal \brief Compare points regarding x coordinate. Use v_point_index for V points and u_point_index + G::size() for U points.  */
-    struct Compare_x{bool operator() (const int point_index_1, const int point_index_2) const;};
-    /** \internal \brief Compare points regarding y coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */
-    struct Compare_y{bool operator() (const int point_index_1, const int point_index_2) const;};
+    /** \internal \brief Compare points regarding x%r coordinate. Use v_point_index for V points and u_point_index + G::size() for U points.  */
+    struct Compare_x{double r; Compare_x(double r); bool operator()(const int point_index_1, const int point_index_2) const;};
+    /** \internal \brief Compare points regarding y%r coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */
+    struct Compare_y{double r; Compare_y(double r);bool operator()(const int point_index_1, const int point_index_2) const;};
 
 private:
     /** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */
@@ -74,20 +75,20 @@ private:
     static Internal_point get_v_point(int v_point_index);
 };
 
-/** \internal \typedef \brief Alias used outside. */
+/** \internal \typedef \brief Shorter alias */
 typedef Persistence_diagrams_graph G;
 
 // static initialization, seems to work but strange
-std::vector<G::Internal_point> Persistence_diagrams_graph::u = [] {return std::vector<G::Internal_point>();}();
-std::vector<G::Internal_point> Persistence_diagrams_graph::v = [] {return std::vector<G::Internal_point>();}();
+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;
 }
 
 template<typename Persistence_diagram1, typename Persistence_diagram2>
-void Persistence_diagrams_graph::initialize(const Persistence_diagram1 &diag1,
-                                            const Persistence_diagram2 &diag2, double e){
+inline void G::initialize(const Persistence_diagram1 &diag1,
+                          const Persistence_diagram2 &diag2, double e){
     u.clear();
     v.clear();
     for (auto it = diag1.cbegin(); it != diag1.cend(); ++it)
@@ -100,25 +101,25 @@ void Persistence_diagrams_graph::initialize(const Persistence_diagram1 &diag1,
         swap(u, v);
 }
 
-inline bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) {
+inline bool G::on_the_u_diagonal(int u_point_index) {
     return u_point_index >= static_cast<int> (u.size());
 }
 
-inline bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) {
+inline bool G::on_the_v_diagonal(int v_point_index) {
     return v_point_index >= static_cast<int> (v.size());
 }
 
-inline int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) {
+inline int G::corresponding_point_in_u(int v_point_index) {
     return on_the_v_diagonal(v_point_index) ?
                 v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size());
 }
 
-inline int Persistence_diagrams_graph::corresponding_point_in_v(int u_point_index) {
+inline int G::corresponding_point_in_v(int u_point_index) {
     return on_the_u_diagonal(u_point_index) ?
                 u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size());
 }
 
-inline double Persistence_diagrams_graph::distance(int u_point_index, int v_point_index) {
+inline double G::distance(int u_point_index, int v_point_index) {
     if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index))
         return 0;
     Internal_point p_u = get_u_point(u_point_index);
@@ -126,11 +127,11 @@ inline double Persistence_diagrams_graph::distance(int u_point_index, int v_poin
     return std::max(std::fabs(p_u.first - p_v.first), std::fabs(p_u.second - p_v.second));
 }
 
-inline int Persistence_diagrams_graph::size() {
+inline int G::size() {
     return static_cast<int> (u.size() + v.size());
 }
 
-inline std::unique_ptr< std::vector<double> > Persistence_diagrams_graph::sorted_distances() {
+inline std::unique_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)
@@ -140,7 +141,7 @@ inline std::unique_ptr< std::vector<double> > Persistence_diagrams_graph::sorted
     return sd_up;
 }
 
-inline Persistence_diagrams_graph::Internal_point Persistence_diagrams_graph::get_u_point(int u_point_index) {
+inline G::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));
@@ -148,7 +149,7 @@ inline Persistence_diagrams_graph::Internal_point Persistence_diagrams_graph::ge
     return Internal_point(x, x);
 }
 
-inline Persistence_diagrams_graph::Internal_point Persistence_diagrams_graph::get_v_point(int v_point_index) {
+inline G::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));
@@ -156,16 +157,30 @@ inline Persistence_diagrams_graph::Internal_point Persistence_diagrams_graph::ge
     return Internal_point(x, x);
 }
 
-inline bool Persistence_diagrams_graph::Compare_x::operator() (const int point_index_1, const int point_index_2) const{
+G::Compare_x::Compare_x(double r)
+    : r(r){ }
+
+G::Compare_y::Compare_y(double r)
+    : r(r){ }
+
+inline bool G::Compare_x::operator()(const int point_index_1, const int point_index_2) const{
     G::Internal_point p1 = point_index_1 < G::size() ? G::get_v_point(point_index_1) : G::get_u_point(point_index_1 - G::size());
     G::Internal_point p2 = point_index_2 < G::size() ? G::get_v_point(point_index_2) : G::get_u_point(point_index_2 - G::size());
-    return p1.first < p2.first;
+    double x1 = fmod(p1.first,r);
+    double x2 = fmod(p2.first,r);
+    if(x1 == x2)
+        return point_index_1 > point_index_2;
+    return x1 < x2;
 }
 
-inline bool Persistence_diagrams_graph::Compare_y::operator() (const int point_index_1, const int point_index_2) const{
+inline bool G::Compare_y::operator()(const int point_index_1, const int point_index_2) const{
     G::Internal_point p1 = point_index_1 < G::size() ? G::get_v_point(point_index_1) : G::get_u_point(point_index_1 - G::size());
     G::Internal_point p2 = point_index_2 < G::size() ? G::get_v_point(point_index_2) : G::get_u_point(point_index_2 - G::size());
-    return p1.second < p2.second;
+    double y1 = fmod(p1.second,r);
+    double y2 = fmod(p2.second,r);
+    if(y1 == y2)
+        return  point_index_1 > point_index_2;
+    return y1 < y2;
 }
 
 }  // namespace bottleneck
diff --git a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
index 4d820c7f..e403735c 100644
--- a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
+++ b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
@@ -85,7 +85,7 @@ private:
 typedef Naive_pnf Planar_neighbors_finder;
 
 
-Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(double r) :
+inline Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(double r) :
     r(r) { }
 
 inline Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() {}
@@ -100,7 +100,7 @@ inline std::unique_ptr< std::list<int> > Abstract_planar_neighbors_finder::pull_
     return all_pull;
 }
 
-Naive_pnf::Naive_pnf(double r) :
+inline Naive_pnf::Naive_pnf(double r) :
     Abstract_planar_neighbors_finder(r), candidates() { }
 
 inline void Naive_pnf::add(int v_point_index) {