Merge last trunk modifications

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

Former-commit-id: 29b4678c6e1ecaec8c43fc68a36aa988ddcbf298

22 files changed, 198 insertions, 142 deletions

diff --git a/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp b/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp
index fd9f0858..fd164b22 100644
--- a/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp
+++ b/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp
@@ -74,7 +74,7 @@ int main(int argc, char * argv[]) {
   // --------------------------------------------
   // Rips persistence
   // --------------------------------------------
-  Rips_complex rips_complex(off_reader.get_point_cloud(), threshold, Euclidean_distance());
+  Rips_complex rips_complex(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance());
 
   // Construct the Rips complex in a Simplex Tree
   Simplex_tree rips_stree;
diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h
index e1708c5b..f51e22e9 100644
--- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h
+++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h
@@ -26,7 +26,8 @@
 #include <gudhi/Neighbors_finder.h>
 
 #include <vector>
-#include <list>
+#include <unordered_set>
+#include <algorithm>
 
 namespace Gudhi {
 
@@ -40,8 +41,6 @@ class Graph_matching {
  public:
   /** \internal \brief Constructor constructing an empty matching. */
   explicit Graph_matching(Persistence_graph &g);
-  /** \internal \brief Copy operator. */
-  Graph_matching& operator=(const Graph_matching& m);
   /** \internal \brief Is the matching perfect ? */
   bool perfect() const;
   /** \internal \brief Augments the matching with a maximal set of edge-disjoint shortest augmenting paths. */
@@ -50,12 +49,12 @@ class Graph_matching {
   void set_r(double r);
 
  private:
-  Persistence_graph& g;
+  Persistence_graph* gp;
   double r;
   /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */
   std::vector<int> v_to_u;
   /** \internal \brief All the unmatched points in U. */
-  std::list<int> unmatched_in_u;
+  std::unordered_set<int> unmatched_in_u;
 
   /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */
   Layered_neighbors_finder layering() const;
@@ -66,17 +65,9 @@ class Graph_matching {
 };
 
 inline Graph_matching::Graph_matching(Persistence_graph& g)
-    : g(g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u() {
+    : gp(&g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u(g.size()) {
   for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index)
-    unmatched_in_u.emplace_back(u_point_index);
-}
-
-inline Graph_matching& Graph_matching::operator=(const Graph_matching& m) {
-  g = m.g;
-  r = m.r;
-  v_to_u = m.v_to_u;
-  unmatched_in_u = m.unmatched_in_u;
-  return *this;
+    unmatched_in_u.insert(u_point_index);
 }
 
 inline bool Graph_matching::perfect() const {
@@ -88,12 +79,12 @@ inline bool Graph_matching::multi_augment() {
     return false;
   Layered_neighbors_finder layered_nf(layering());
   int max_depth = layered_nf.vlayers_number()*2 - 1;
-  double rn = sqrt(g.size());
+  double rn = sqrt(gp->size());
   // verification of a necessary criterion in order to shortcut if possible
   if (max_depth < 0 || (unmatched_in_u.size() > rn && max_depth >= rn))
     return false;
   bool successful = false;
-  std::list<int> tries(unmatched_in_u);
+  std::vector<int> tries(unmatched_in_u.cbegin(), unmatched_in_u.cend());
   for (auto it = tries.cbegin(); it != tries.cend(); it++)
     // 'augment' has side-effects which have to be always executed, don't change order
     successful = augment(layered_nf, *it, max_depth) || successful;
@@ -133,12 +124,12 @@ inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u
 }
 
 inline Layered_neighbors_finder Graph_matching::layering() const {
-  std::list<int> u_vertices(unmatched_in_u);
-  std::list<int> v_vertices;
-  Neighbors_finder nf(g, r);
-  for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index)
+  std::vector<int> u_vertices(unmatched_in_u.cbegin(), unmatched_in_u.cend());
+  std::vector<int> v_vertices;
+  Neighbors_finder nf(*gp, r);
+  for (int v_point_index = 0; v_point_index < gp->size(); ++v_point_index)
     nf.add(v_point_index);
-  Layered_neighbors_finder layered_nf(g, r);
+  Layered_neighbors_finder layered_nf(*gp, r);
   for (int layer = 0; !u_vertices.empty(); layer++) {
     // one layer is one step in the BFS
     for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) {
@@ -166,7 +157,8 @@ inline Layered_neighbors_finder Graph_matching::layering() const {
 }
 
 inline void Graph_matching::update(std::vector<int>& path) {
-  unmatched_in_u.remove(path.front());
+  // Must return 1.
+  unmatched_in_u.erase(path.front());
   for (auto it = path.cbegin(); it != path.cend(); ++it) {
     // Be careful, the iterator is incremented twice each time
     int tmp = *it;
diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h
index cd5486f8..bdc47578 100644
--- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h
+++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h
@@ -25,9 +25,6 @@
 
 // Inclusion order is important for CGAL patch
 #include <CGAL/Kd_tree.h>
-#include <CGAL/Kd_tree_node.h>
-#include <CGAL/Orthogonal_k_neighbor_search.h>
-#include <CGAL/Weighted_Minkowski_distance.h>
 #include <CGAL/Search_traits.h>
 
 #include <gudhi/Persistence_graph.h>
@@ -40,6 +37,33 @@ namespace Gudhi {
 
 namespace persistence_diagram {
 
+/** \internal \brief Variant of CGAL::Fuzzy_iso_box to ensure that the box ic closed. It isn't clear how necessary that is.
+ */
+struct Square_query {
+  typedef CGAL::Dimension_tag<2> D;
+  typedef Internal_point Point_d;
+  typedef double FT;
+  bool contains(Point_d p) const {
+    return std::abs(p.x()-c.x())<=size && std::abs(p.y()-c.y())<=size;
+  }
+  bool inner_range_intersects(CGAL::Kd_tree_rectangle<FT,D> const&r) const {
+    return
+      r.max_coord(0) >= c.x() - size &&
+      r.min_coord(0) <= c.x() + size &&
+      r.max_coord(1) >= c.y() - size &&
+      r.min_coord(1) <= c.y() + size;
+  }
+  bool outer_range_contains(CGAL::Kd_tree_rectangle<FT,D> const&r) const {
+    return
+      r.min_coord(0) >= c.x() - size &&
+      r.max_coord(0) <= c.x() + size &&
+      r.min_coord(1) >= c.y() - size &&
+      r.max_coord(1) <= c.y() + size;
+  }
+  Point_d c;
+  FT size;
+};
+
 /** \internal \brief data structure used to find any point (including projections) in V near to a query point from U
  * (which can be a projection).
  *
@@ -51,9 +75,7 @@ namespace persistence_diagram {
 class Neighbors_finder {
   typedef CGAL::Dimension_tag<2> D;
   typedef CGAL::Search_traits<double, Internal_point, const double*, Construct_coord_iterator, D> Traits;
-  typedef CGAL::Weighted_Minkowski_distance<Traits> Distance;
-  typedef CGAL::Orthogonal_k_neighbor_search<Traits, Distance> K_neighbor_search;
-  typedef K_neighbor_search::Tree Kd_tree;
+  typedef CGAL::Kd_tree<Traits> Kd_tree;
 
  public:
   /** \internal \brief Constructor taking the near distance definition as parameter. */
@@ -123,15 +145,13 @@ inline int Neighbors_finder::pull_near(int u_point_index) {
   } else {
     // Is the query point near to a V point in the plane ?
     Internal_point u_point = g.get_u_point(u_point_index);
-    std::array<double, 2> w = {
-      {1., 1.}
-    };
-    K_neighbor_search search(kd_t, u_point, 1, 0., true, Distance(0, 2, w.begin(), w.end()));
-    auto it = search.begin();
-    if (it == search.end() || g.distance(u_point_index, it->first.point_index) > r)
+    auto neighbor = kd_t.search_any_point(Square_query{u_point, r});
+    if(!neighbor)
       return null_point_index();
-    tmp = it->first.point_index;
-    kd_t.remove(g.get_v_point(tmp));
+    tmp = neighbor->point_index;
+    auto point = g.get_v_point(tmp);
+    int idx = point.point_index;
+    kd_t.remove(point, [idx](Internal_point const&p){return p.point_index == idx;});
   }
   return tmp;
 }
diff --git a/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp b/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp
index ba752999..252e8aef 100644
--- a/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp
+++ b/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp
@@ -83,7 +83,7 @@ int main(int argc, char * argv[]) {
 
   // Compute the proximity graph of the points
   start = std::chrono::system_clock::now();
-  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Euclidean_distance());
+  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance());
   end = std::chrono::system_clock::now();
   elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
   std::cout << "Compute Rips graph in " << elapsed_sec << " ms.\n";
diff --git a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp
index 7674b5a5..dae36ed2 100644
--- a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp
+++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp
@@ -62,7 +62,7 @@ int main(int argc, char * argv[]) {
   program_options(argc, argv, off_file_points, filediag, threshold, dim_max, min_p, max_p, min_persistence);
 
   Points_off_reader off_reader(off_file_points);
-  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Euclidean_distance());
+  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance());
 
   // Construct the Rips complex in a Simplex Tree
   Simplex_tree simplex_tree;
diff --git a/src/Persistent_cohomology/example/rips_persistence.cpp b/src/Persistent_cohomology/example/rips_persistence.cpp
index c6378de7..d504798b 100644
--- a/src/Persistent_cohomology/example/rips_persistence.cpp
+++ b/src/Persistent_cohomology/example/rips_persistence.cpp
@@ -60,7 +60,7 @@ int main(int argc, char * argv[]) {
   program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence);
 
   Points_off_reader off_reader(off_file_points);
-  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Euclidean_distance());
+  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance());
 
   // Construct the Rips complex in a Simplex Tree
   Simplex_tree simplex_tree;
diff --git a/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp
index b159c62e..554eeba6 100644
--- a/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp
+++ b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp
@@ -82,7 +82,7 @@ int main(int argc, char * argv[]) {
 
   // Compute the proximity graph of the points
   Graph_t prox_graph = compute_proximity_graph(off_reader.get_point_cloud(), threshold
-                                               , Euclidean_distance());
+                                               , Gudhi::Euclidean_distance());
 
   // Construct the Rips complex in a Simplex Tree
   Simplex_tree st;
diff --git a/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp
index 63da9847..9618f278 100644
--- a/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp
+++ b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp
@@ -70,7 +70,7 @@ int main(int argc, char * argv[]) {
   program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence);
 
   Points_off_reader off_reader(off_file_points);
-  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Euclidean_distance());
+  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance());
 
   // Construct the Rips complex in a Simplex Tree
   Simplex_tree& st = *new Simplex_tree;
diff --git a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp
index 0ed3fddf..da418034 100644
--- a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp
+++ b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp
@@ -9,8 +9,8 @@
 #define BOOST_TEST_MODULE "betti_numbers"
 #include <boost/test/unit_test.hpp>
 
-#include "gudhi/Simplex_tree.h"
-#include "gudhi/Persistent_cohomology.h"
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
 
 struct MiniSTOptions : Gudhi::Simplex_tree_options_full_featured {
   // Implicitly use 0 as filtration value for all simplices
diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp
index 6efd749e..e2e0bc71 100644
--- a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp
+++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp
@@ -10,10 +10,10 @@
 #define BOOST_TEST_MODULE "persistent_cohomology"
 #include <boost/test/unit_test.hpp>
 
-#include "gudhi/graph_simplicial_complex.h"
-#include "gudhi/reader_utils.h"
-#include "gudhi/Simplex_tree.h"
-#include "gudhi/Persistent_cohomology.h"
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/reader_utils.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
 
 using namespace Gudhi;
 using namespace Gudhi::persistent_cohomology;
diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp
index 1a6e3296..6a7ecbec 100644
--- a/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp
+++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp
@@ -9,11 +9,11 @@
 #define BOOST_TEST_MODULE "persistent_cohomology_multi_field"
 #include <boost/test/unit_test.hpp>
 
-#include "gudhi/graph_simplicial_complex.h"
-#include "gudhi/reader_utils.h"
-#include "gudhi/Simplex_tree.h"
-#include "gudhi/Persistent_cohomology.h"
-#include "gudhi/Persistent_cohomology/Multi_field.h"
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/reader_utils.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Persistent_cohomology/Multi_field.h>
 
 using namespace Gudhi;
 using namespace Gudhi::persistent_cohomology;
diff --git a/src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp b/src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp
index 3fd69ebc..a1db8910 100644
--- a/src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp
+++ b/src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp
@@ -27,7 +27,7 @@ int main() {
   // Init of a Rips complex from points
   // ----------------------------------------------------------------------------
   double threshold = 12.0;
-  Rips_complex rips_complex_from_points(points, threshold, Euclidean_distance());
+  Rips_complex rips_complex_from_points(points, threshold, Gudhi::Euclidean_distance());
 
   Simplex_tree stree;
   rips_complex_from_points.create_complex(stree, 1);
diff --git a/src/Rips_complex/example/example_rips_complex_from_off_file.cpp b/src/Rips_complex/example/example_rips_complex_from_off_file.cpp
index a1e4e255..de2e4ea4 100644
--- a/src/Rips_complex/example/example_rips_complex_from_off_file.cpp
+++ b/src/Rips_complex/example/example_rips_complex_from_off_file.cpp
@@ -32,7 +32,7 @@ int main(int argc, char **argv) {
   // Init of a Rips complex from an OFF file
   // ----------------------------------------------------------------------------
   Gudhi::Points_off_reader<Point> off_reader(off_file_name);
-  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Euclidean_distance());
+  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance());
 
   std::streambuf* streambufffer;
   std::ofstream ouput_file_stream;
diff --git a/src/Rips_complex/test/test_rips_complex.cpp b/src/Rips_complex/test/test_rips_complex.cpp
index ae68ba0d..fc2179f2 100644
--- a/src/Rips_complex/test/test_rips_complex.cpp
+++ b/src/Rips_complex/test/test_rips_complex.cpp
@@ -60,7 +60,7 @@ BOOST_AUTO_TEST_CASE(RIPS_DOC_OFF_file) {
       rips_threshold << "==========" << std::endl;
 
   Gudhi::Points_off_reader<Point> off_reader(off_file_name);
-  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), rips_threshold, Euclidean_distance());
+  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), rips_threshold, Gudhi::Euclidean_distance());
 
   const int DIMENSION_1 = 1;
   Simplex_tree st;
@@ -89,10 +89,10 @@ BOOST_AUTO_TEST_CASE(RIPS_DOC_OFF_file) {
         std::cout << vertex << ",";
         vp.push_back(off_reader.get_point_cloud().at(vertex));
       }
-      std::cout << ") - distance =" << Euclidean_distance()(vp.at(0), vp.at(1)) <<
+      std::cout << ") - distance =" << Gudhi::Euclidean_distance()(vp.at(0), vp.at(1)) <<
           " - filtration =" << st.filtration(f_simplex) << std::endl;
       BOOST_CHECK(vp.size() == 2);
-      BOOST_CHECK(are_almost_the_same(st.filtration(f_simplex), Euclidean_distance()(vp.at(0), vp.at(1))));
+      BOOST_CHECK(are_almost_the_same(st.filtration(f_simplex), Gudhi::Euclidean_distance()(vp.at(0), vp.at(1))));
     }
   }
 
@@ -341,7 +341,7 @@ BOOST_AUTO_TEST_CASE(Rips_create_complex_throw) {
       rips_threshold << "==========" << std::endl;
 
   Gudhi::Points_off_reader<Point> off_reader(off_file_name);
-  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), rips_threshold, Euclidean_distance());
+  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), rips_threshold, Gudhi::Euclidean_distance());
 
   Simplex_tree stree;
   std::vector<int> simplex = {0, 1, 2};
diff --git a/src/Subsampling/example/example_choose_n_farthest_points.cpp b/src/Subsampling/example/example_choose_n_farthest_points.cpp
index 533aba74..ebf631fc 100644
--- a/src/Subsampling/example/example_choose_n_farthest_points.cpp
+++ b/src/Subsampling/example/example_choose_n_farthest_points.cpp
@@ -19,7 +19,9 @@ int main(void) {
 
   K k;
   std::vector<Point_d> results;
-  Gudhi::subsampling::choose_n_farthest_points(k, points, 100, std::back_inserter(results));
+  Gudhi::subsampling::choose_n_farthest_points(k, points, 100,
+                                               Gudhi::subsampling::random_starting_point,
+                                               std::back_inserter(results));
   std::cout << "Before sparsification: " << points.size() << " points.\n";
   std::cout << "After  sparsification: " << results.size() << " points.\n";
 
diff --git a/src/Subsampling/example/example_custom_kernel.cpp b/src/Subsampling/example/example_custom_kernel.cpp
index 25b5bf6c..2d42bdde 100644
--- a/src/Subsampling/example/example_custom_kernel.cpp
+++ b/src/Subsampling/example/example_custom_kernel.cpp
@@ -54,7 +54,9 @@ int main(void) {
   std::vector<Point_d> points = {0, 1, 2, 3};
   std::vector<Point_d> results;
 
-  Gudhi::subsampling::choose_n_farthest_points(k, points, 2, std::back_inserter(results));
+  Gudhi::subsampling::choose_n_farthest_points(k, points, 2,
+                                               Gudhi::subsampling::random_starting_point,
+                                               std::back_inserter(results));
   std::cout << "Before sparsification: " << points.size() << " points.\n";
   std::cout << "After  sparsification: " << results.size() << " points.\n";
   std::cout << "Result table: {" << results[0] << "," << results[1] << "}\n";
diff --git a/src/Subsampling/include/gudhi/choose_n_farthest_points.h b/src/Subsampling/include/gudhi/choose_n_farthest_points.h
index 5e908090..86500b28 100644
--- a/src/Subsampling/include/gudhi/choose_n_farthest_points.h
+++ b/src/Subsampling/include/gudhi/choose_n_farthest_points.h
@@ -25,16 +25,9 @@
 
 #include <boost/range.hpp>
 
-#include <gudhi/Kd_tree_search.h>
-
-#include <gudhi/Clock.h>
-
-#include <CGAL/Search_traits.h>
-#include <CGAL/Search_traits_adapter.h>
-#include <CGAL/Fuzzy_sphere.h>
+#include <gudhi/Null_output_iterator.h>
 
 #include <iterator>
-#include <algorithm>  // for sort
 #include <vector>
 #include <random>
 #include <limits>  // for numeric_limits<>
@@ -43,36 +36,51 @@ namespace Gudhi {
 
 namespace subsampling {
 
+/**
+ *  \ingroup subsampling
+ */
+enum : std::size_t {
+/**
+ *  Argument for `choose_n_farthest_points` to indicate that the starting point should be picked randomly.
+ */
+  random_starting_point = std::size_t(-1)
+};
+
 /** 
  *  \ingroup subsampling
  *  \brief Subsample by a greedy strategy of iteratively adding the farthest point from the
  *  current chosen point set to the subsampling. 
- *  The iteration starts with the landmark `starting point`.
+ *  The iteration starts with the landmark `starting point` or, if `starting point==random_starting_point`, with a random landmark.
  *  \tparam Kernel must provide a type Kernel::Squared_distance_d which is a model of the 
  *          concept <a target="_blank"
- *   href="http://doc.cgal.org/latest/Kernel_d/classKernel__d_1_1Squared__distance__d.html">Kernel_d::Squared_distance_d</a>
- *   concept.
- *  It must also contain a public member 'squared_distance_d_object' of this type.
+ *   href="http://doc.cgal.org/latest/Kernel_d/classKernel__d_1_1Squared__distance__d.html">Kernel_d::Squared_distance_d</a> (despite the name, taken from CGAL, this can be any kind of metric or proximity measure).
+ *  It must also contain a public member `squared_distance_d_object()` that returns an object of this type.
  *  \tparam Point_range Range whose value type is Kernel::Point_d.  It must provide random-access 
  *         via `operator[]` and the points should be stored contiguously in memory.
- *  \tparam OutputIterator Output iterator whose value type is Kernel::Point_d.
- *  \details It chooses `final_size` points from a random access range `input_pts` and
- *  outputs it in the output iterator `output_it`.
+ *  \tparam PointOutputIterator Output iterator whose value type is Kernel::Point_d.
+ *  \tparam DistanceOutputIterator Output iterator for distances.
+ *  \details It chooses `final_size` points from a random access range
+ *  `input_pts` and outputs them in the output iterator `output_it`. It also
+ *  outputs the distance from each of those points to the set of previous
+ *  points in `dist_it`.
  * @param[in] k A kernel object.
  * @param[in] input_pts Const reference to the input points.
  * @param[in] final_size The size of the subsample to compute.
  * @param[in] starting_point The seed in the farthest point algorithm.
- * @param[out] output_it The output iterator.
+ * @param[out] output_it The output iterator for points.
+ * @param[out] dist_it The optional output iterator for distances.
  *  
  */
 template < typename Kernel,
 typename Point_range,
-typename OutputIterator>
+typename PointOutputIterator,
+typename DistanceOutputIterator = Null_output_iterator>
 void choose_n_farthest_points(Kernel const &k,
                               Point_range const &input_pts,
                               std::size_t final_size,
                               std::size_t starting_point,
-                              OutputIterator output_it) {
+                              PointOutputIterator output_it,
+                              DistanceOutputIterator dist_it = {}) {
   std::size_t nb_points = boost::size(input_pts);
   if (final_size > nb_points)
     final_size = nb_points;
@@ -81,6 +89,14 @@ void choose_n_farthest_points(Kernel const &k,
   if (final_size < 1)
     return;
 
+  if (starting_point == random_starting_point) {
+    // Choose randomly the first landmark
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_int_distribution<std::size_t> dis(0, (input_pts.size() - 1));
+    starting_point = dis(gen);
+  }
+
   typename Kernel::Squared_distance_d sqdist = k.squared_distance_d_object();
 
   std::size_t current_number_of_landmarks = 0;  // counter for landmarks
@@ -92,6 +108,7 @@ void choose_n_farthest_points(Kernel const &k,
   for (current_number_of_landmarks = 0; current_number_of_landmarks != final_size; current_number_of_landmarks++) {
     // curr_max_w at this point is the next landmark
     *output_it++ = input_pts[curr_max_w];
+    *dist_it++ = dist_to_L[curr_max_w];
     std::size_t i = 0;
     for (auto& p : input_pts) {
       double curr_dist = sqdist(p, *(std::begin(input_pts) + curr_max_w));
@@ -109,47 +126,6 @@ void choose_n_farthest_points(Kernel const &k,
   }
 }
 
-/** 
- *  \ingroup subsampling
- *  \brief Subsample by a greedy strategy of iteratively adding the farthest point from the
- *  current chosen point set to the subsampling. 
- *  The iteration starts with a random landmark.
- *  \tparam Kernel must provide a type Kernel::Squared_distance_d which is a model of the 
- *          concept <a target="_blank"
- *   href="http://doc.cgal.org/latest/Kernel_d/classKernel__d_1_1Squared__distance__d.html">Kernel_d::Squared_distance_d</a>
- *   concept.
- *  It must also contain a public member 'squared_distance_d_object' of this type.
- *  \tparam Point_range Range whose value type is Kernel::Point_d.  It must provide random-access 
- *         via `operator[]` and the points should be stored contiguously in memory.
- *  \tparam OutputIterator Output iterator whose value type is Kernel::Point_d.
- *  \details It chooses `final_size` points from a random access range `input_pts` and
- *  outputs it in the output iterator `output_it`.
- * @param[in] k A kernel object.
- * @param[in] input_pts Const reference to the input points.
- * @param[in] final_size The size of the subsample to compute.
- * @param[out] output_it The output iterator.
- *  
- */
-template < typename Kernel,
-typename Point_range,
-typename OutputIterator>
-void choose_n_farthest_points(Kernel const& k,
-                              Point_range const &input_pts,
-                              unsigned final_size,
-                              OutputIterator output_it) {
-  // Tests to the limit
-  if ((final_size < 1) || (input_pts.size() == 0))
-    return;
-
-  // Choose randomly the first landmark
-  std::random_device rd;
-  std::mt19937 gen(rd());
-  std::uniform_int_distribution<> dis(0, (input_pts.size() - 1));
-  std::size_t starting_point = dis(gen);
-
-  choose_n_farthest_points(k, input_pts, final_size, starting_point, output_it);
-}
-
 }  // namespace subsampling
 
 }  // namespace Gudhi
diff --git a/src/Subsampling/test/test_choose_n_farthest_points.cpp b/src/Subsampling/test/test_choose_n_farthest_points.cpp
index 0bc0dff4..ee9d4c77 100644
--- a/src/Subsampling/test/test_choose_n_farthest_points.cpp
+++ b/src/Subsampling/test/test_choose_n_farthest_points.cpp
@@ -25,7 +25,7 @@
 // #endif
 
 #define BOOST_TEST_DYN_LINK
-#define BOOST_TEST_MODULE "witness_complex_points"
+#define BOOST_TEST_MODULE Subsampling - test choose_n_farthest_points
 #include <boost/test/unit_test.hpp>
 #include <boost/mpl/list.hpp>
 
@@ -56,7 +56,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_choose_farthest_point, Kernel, list_of_tested
 
   landmarks.clear();
   Kernel k;
-  Gudhi::subsampling::choose_n_farthest_points(k, points, 100, std::back_inserter(landmarks));
+  Gudhi::subsampling::choose_n_farthest_points(k, points, 100, Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks));
 
   BOOST_CHECK(landmarks.size() == 100);
   for (auto landmark : landmarks)
@@ -70,34 +70,45 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_choose_farthest_point_limits, Kernel, list_of
   typedef typename Kernel::FT FT;
   typedef typename Kernel::Point_d Point_d;
   std::vector< Point_d > points, landmarks;
+  std::vector< FT > distances;
   landmarks.clear();
   Kernel k;
   // Choose -1 farthest points in an empty point cloud
-  Gudhi::subsampling::choose_n_farthest_points(k, points, -1, std::back_inserter(landmarks));
+  Gudhi::subsampling::choose_n_farthest_points(k, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances));
   BOOST_CHECK(landmarks.size() == 0);
-  landmarks.clear();
+  landmarks.clear(); distances.clear();
   // Choose 0 farthest points in an empty point cloud
-  Gudhi::subsampling::choose_n_farthest_points(k, points, 0, std::back_inserter(landmarks));
+  Gudhi::subsampling::choose_n_farthest_points(k, points, 0, -1, std::back_inserter(landmarks), std::back_inserter(distances));
   BOOST_CHECK(landmarks.size() == 0);
-  landmarks.clear();
+  landmarks.clear(); distances.clear();
   // Choose 1 farthest points in an empty point cloud
-  Gudhi::subsampling::choose_n_farthest_points(k, points, 1, std::back_inserter(landmarks));
+  Gudhi::subsampling::choose_n_farthest_points(k, points, 1, -1, std::back_inserter(landmarks), std::back_inserter(distances));
   BOOST_CHECK(landmarks.size() == 0);
-  landmarks.clear();
+  landmarks.clear(); distances.clear();
 
   std::vector<FT> point({0.0, 0.0, 0.0, 0.0});
   points.push_back(Point_d(point.begin(), point.end()));
-  // Choose -1 farthest points in an empty point cloud
-  Gudhi::subsampling::choose_n_farthest_points(k, points, -1, std::back_inserter(landmarks));
-  BOOST_CHECK(landmarks.size() == 1);
-  landmarks.clear();
+  // Choose -1 farthest points in a one point cloud
+  Gudhi::subsampling::choose_n_farthest_points(k, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances));
+  BOOST_CHECK(landmarks.size() == 1 && distances.size() == 1);
+  BOOST_CHECK(distances[0] == std::numeric_limits<FT>::infinity());
+  landmarks.clear(); distances.clear();
   // Choose 0 farthest points in a one point cloud
-  Gudhi::subsampling::choose_n_farthest_points(k, points, 0, std::back_inserter(landmarks));
-  BOOST_CHECK(landmarks.size() == 0);
-  landmarks.clear();
+  Gudhi::subsampling::choose_n_farthest_points(k, points, 0, -1, std::back_inserter(landmarks), std::back_inserter(distances));
+  BOOST_CHECK(landmarks.size() == 0 && distances.size() == 0);
+  landmarks.clear(); distances.clear();
   // Choose 1 farthest points in a one point cloud
-  Gudhi::subsampling::choose_n_farthest_points(k, points, 1, std::back_inserter(landmarks));
-  BOOST_CHECK(landmarks.size() == 1);
-  landmarks.clear();
+  Gudhi::subsampling::choose_n_farthest_points(k, points, 1, -1, std::back_inserter(landmarks), std::back_inserter(distances));
+  BOOST_CHECK(landmarks.size() == 1 && distances.size() == 1);
+  BOOST_CHECK(distances[0] == std::numeric_limits<FT>::infinity());
+  landmarks.clear(); distances.clear();
 
+  std::vector<FT> point2({1.0, 0.0, 0.0, 0.0});
+  points.push_back(Point_d(point2.begin(), point2.end()));
+  // Choose all farthest points in a one point cloud
+  Gudhi::subsampling::choose_n_farthest_points(k, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances));
+  BOOST_CHECK(landmarks.size() == 2 && distances.size() == 2);
+  BOOST_CHECK(distances[0] == std::numeric_limits<FT>::infinity());
+  BOOST_CHECK(distances[1] == 1);
+  landmarks.clear(); distances.clear();
 }
diff --git a/src/common/include/gudhi/Null_output_iterator.h b/src/common/include/gudhi/Null_output_iterator.h
new file mode 100644
index 00000000..42e6e449
--- /dev/null
+++ b/src/common/include/gudhi/Null_output_iterator.h
@@ -0,0 +1,48 @@
+/*    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):       Marc Glisse
+ *
+ *    Copyright (C) 2017  INRIA
+ *
+ *    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 NULL_OUTPUT_ITERATOR_H_
+#define NULL_OUTPUT_ITERATOR_H_
+
+#include <iterator>
+
+namespace Gudhi {
+
+/** An output iterator that ignores whatever it is given. */
+struct Null_output_iterator {
+  typedef std::output_iterator_tag iterator_category;
+  typedef void                     value_type;
+  typedef void                     difference_type;
+  typedef void                     pointer;
+  typedef void                     reference;
+
+  Null_output_iterator& operator++() {return *this;}
+  Null_output_iterator operator++(int) {return *this;}
+  struct proxy {
+    template<class T>
+    proxy& operator=(T&&){return *this;}
+  };
+  proxy operator*()const{return {};}
+};
+}  // namespace Gudhi
+
+#endif  // NULL_OUTPUT_ITERATOR_H_
diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h
index 22747637..f6e2ab5a 100644
--- a/src/common/include/gudhi/distance_functions.h
+++ b/src/common/include/gudhi/distance_functions.h
@@ -27,6 +27,8 @@
 #include <type_traits>  // for std::decay
 #include <iterator>  // for std::begin, std::end
 
+namespace Gudhi {
+
 /** @file
  * @brief Global distance functions
  */
@@ -48,4 +50,6 @@ class Euclidean_distance {
   }
 };
 
+}  // namespace Gudhi
+
 #endif  // DISTANCE_FUNCTIONS_H_
diff --git a/src/cython/include/Rips_complex_interface.h b/src/cython/include/Rips_complex_interface.h
index 1879bd74..6d813f4a 100644
--- a/src/cython/include/Rips_complex_interface.h
+++ b/src/cython/include/Rips_complex_interface.h
@@ -49,7 +49,7 @@ class Rips_complex_interface {
     if (euclidean) {
       // Rips construction where values is a vector of points
       rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(values, threshold,
-                                                                                   Euclidean_distance());
+                                                                                   Gudhi::Euclidean_distance());
     } else {
       // Rips construction where values is a distance matrix
       rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(values, threshold);
@@ -61,7 +61,8 @@ class Rips_complex_interface {
       // Rips construction where file_name is an OFF file
       Gudhi::Points_off_reader<Point_d> off_reader(file_name);
       rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(off_reader.get_point_cloud(),
-                                                                                   threshold, Euclidean_distance());
+                                                                                   threshold,
+                                                                                   Gudhi::Euclidean_distance());
     } else {
       // Rips construction where values is a distance matrix
       Distance_matrix distances =
diff --git a/src/cython/include/Subsampling_interface.h b/src/cython/include/Subsampling_interface.h
index 1c6032c0..b0f4a50a 100644
--- a/src/cython/include/Subsampling_interface.h
+++ b/src/cython/include/Subsampling_interface.h
@@ -46,7 +46,7 @@ std::vector<std::vector<double>> subsampling_n_farthest_points(const std::vector
                                                                unsigned nb_points) {
   std::vector<std::vector<double>> landmarks;
   Subsampling_dynamic_kernel k;
-  choose_n_farthest_points(k, points, nb_points, std::back_inserter(landmarks));
+  choose_n_farthest_points(k, points, nb_points, random_starting_point, std::back_inserter(landmarks));
 
   return landmarks;
 }