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

Former-commit-id: 93dd0d2099fce5d2b5b05f2ddc22cfebecac14fc

14 files changed, 429 insertions, 427 deletions

diff --git a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h
index 9d7883e0..3a6c6c85 100644
--- a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h
+++ b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h
@@ -33,6 +33,8 @@ namespace graph_induced_complex {
  * 
  * @{
  * 
+ * Visualizations of the simplicial complexes require neato, python and firefox!!
+ *
  * \section covers Covers
  *
  * Nerves and Graph Induced Complexes require a cover C of the input point cloud P,
@@ -63,7 +65,7 @@ namespace graph_induced_complex {
  *
  * When launching:
  *
- * \code $> ./Nerve ../../../data/points/human.off 2 10 0.3
+ * \code $> ./Nerve ../../../data/points/human.off 2 10 0.3 --v
  * \endcode
  *
  * the program output is:
@@ -90,7 +92,7 @@ namespace graph_induced_complex {
  *
  * \image html "gic_complex.png" "GIC of a point cloud."
  *
- * \subsection gicexample Example
+ * \subsection gicexample Example with cover from function
  *
  * This example builds the GIC of a point cloud sampled on a 3D human shape (human.off).
  * The cover C comes from the preimages of intervals (with length 0.075 and gain 0)
@@ -103,13 +105,32 @@ namespace graph_induced_complex {
  *
  * When launching:
  *
- * \code $> ./GIC ../../../data/points/human.off 0.075 2 0.075 0
+ * \code $> ./GIC ../../../data/points/human.off 0.075 2 0.075 0 --v
  * \endcode
  *
  * the program output is:
  *
  * \include Nerve_GIC/GIC.txt
  *
+ * \subsection gicexamplevor Example with cover from Voronoï
+ *
+ * This example builds the GIC of a point cloud sampled on a 3D human shape (human.off).
+ * We randomly subsampled 100 points in the point cloud, which act as seeds of
+ * a geodesic Voronoï diagram. Each cell of the diagram is then an element of C.
+ * The graph G (used to compute both the geodesics for Voronoï and the GIC)
+ * comes from the triangulation of the human shape.
+ *
+ * \include Nerve_GIC/GICvoronoi.cpp
+ *
+ * When launching:
+ *
+ * \code $> ./GICvoronoi ../../../data/points/human.off 100 --v
+ * \endcode
+ *
+ * the program output is:
+ *
+ * \include Nerve_GIC/GICvoronoi.txt
+ *
  * \subsection mapperdeltadefinition Mapper Delta
  *
  * If one restricts to the cliques in G whose nodes all belong to preimages of consecutive
@@ -132,7 +153,7 @@ namespace graph_induced_complex {
  *
  * When launching:
  *
- * \code $> ./MapperDeltaCoord ../../../data/points/human.off 2
+ * \code $> ./MapperDeltaCoord ../../../data/points/human.off 2 --v
  * \endcode
  *
  * the program output is:
@@ -147,21 +168,13 @@ namespace graph_induced_complex {
  *
  * When launching:
  *
- * \code $> ./MapperDeltaFunc ../../../data/points/COIL_database/lucky_cat.off ../../../data/points/COIL_database/lucky_cat_PCA1
+ * \code $> ./MapperDeltaFunc ../../../data/points/COIL_database/lucky_cat.off ../../../data/points/COIL_database/lucky_cat_PCA1 --v
  * \endcode
  *
  * the program output is:
  *
  * \include MapperDeltaFunc.txt
  *
- * If you have python and firefox, all the previous .txt files can then be displayed in a browser.
- * We provide a .py file called visu.py that comes from the
- * <a target="_blank" href="https://github.com/MLWave/kepler-mapper"> Kepler-Mapper </a> library.
- * One can visualize data by launching:
- *
- * \code python visu.py && firefox SC.html
- * \endcode
- *
  * \copyright GNU General Public License v3.                         
  * \verbatim  Contact: gudhi-users@lists.gforge.inria.fr \endverbatim
  */
diff --git a/src/Nerve_GIC/example/CMakeLists.txt b/src/Nerve_GIC/example/CMakeLists.txt
index a2e4efe5..b2c501c3 100644
--- a/src/Nerve_GIC/example/CMakeLists.txt
+++ b/src/Nerve_GIC/example/CMakeLists.txt
@@ -13,6 +13,9 @@ target_link_libraries(MapperDeltaCoord  ${Boost_SYSTEM_LIBRARY})
 add_executable ( MapperDeltaFunc MapperDeltaFunc.cpp )
 target_link_libraries(MapperDeltaFunc  ${Boost_SYSTEM_LIBRARY})
 
+add_executable ( GICvoronoi GICvoronoi.cpp )
+target_link_libraries(MapperDeltaFunc  ${Boost_SYSTEM_LIBRARY})
+
 file(COPY visu.py km.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
 
 if (TBB_FOUND)
diff --git a/src/Nerve_GIC/example/GIC.cpp b/src/Nerve_GIC/example/GIC.cpp
index cae0d55d..1889cb33 100644
--- a/src/Nerve_GIC/example/GIC.cpp
+++ b/src/Nerve_GIC/example/GIC.cpp
@@ -2,8 +2,8 @@
 
 void usage(int nbArgs, char * const progName) {
   std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n";
-  std::cerr << "Usage: " << progName << " filename.off threshold coordinate resolution gain\n";
-  std::cerr << "       i.e.: " << progName << " ../../../../data/points/human.off 0.075 2 0.075 0 \n";
+  std::cerr << "Usage: " << progName << " filename.off threshold coordinate resolution gain [--v] \n";
+  std::cerr << "       i.e.: " << progName << " ../../../../data/points/human.off 0.075 2 0.075 0 --v \n";
   exit(-1);  // ----- >>
 }
 
diff --git a/src/Nerve_GIC/example/GIC.txt b/src/Nerve_GIC/example/GIC.txt
index 6d9f07f5..79f61e92 100644
--- a/src/Nerve_GIC/example/GIC.txt
+++ b/src/Nerve_GIC/example/GIC.txt
@@ -1,92 +1,89 @@
-../../../../data/points/human.off
-coordinate 2
-coordinate 2
-0.075 0
-44 43
-0 -0.954369 220
-1 -0.954377 220
-2 -0.869559 45
-3 -0.869405 45
-4 -0.794255 25
-5 -0.795222 26
-6 -0.723847 48
-7 -0.722679 48
-8 -0.633796 48
-9 -0.634232 48
-10 -0.564455 48
-11 -0.564428 48
-12 -0.491837 98
-13 -0.491864 98
-14 -0.422733 88
-15 -0.422739 88
-16 -0.343376 74
-17 -0.343367 74
-18 -0.271635 96
-19 -0.190912 96
-20 -0.158869 3
-21 -0.158905 3
-22 -0.114958 97
-23 -0.108561 96
-24 -0.115241 100
-25 -0.0317845 185
-26 -0.041866 50
-27 -0.0429161 51
-28 0.0353033 158
-29 0.0275106 31
-30 0.0252492 30
-31 0.107097 170
-32 0.110226 51
-33 0.110952 50
-34 0.186188 57
-35 0.176564 97
-36 0.186215 57
-37 0.256898 232
-38 0.338811 268
-39 0.406299 273
-40 0.481925 220
-41 0.555454 241
-42 0.632415 198
-43 0.744473 407
-0 3
-1 2
-2 4
-3 5
-4 7
-5 6
-6 9
-7 8
-8 10
-9 11
-10 13
-11 12
-12 14
-13 15
-14 16
-15 17
-16 18
-17 18
-18 19
-19 23
-20 22
-21 24
-22 27
-23 25
-24 26
-25 28
-26 29
-27 30
-28 31
-29 33
-30 32
-31 35
-32 34
-33 36
-34 37
-35 37
-36 37
-37 38
-38 39
-39 40
-40 41
-41 42
-42 43
+Graph induced complex is of dimension 1 - 87 simplices - 44 vertices.
+Iterator on graph induced complex simplices
+0 
+1 
+2 
+2 1 
+3 
+3 0 
+4 
+4 2 
+5 
+5 3 
+6 
+6 5 
+7 
+7 4 
+8 
+8 7 
+9 
+9 6 
+10 
+10 8 
+11 
+11 9 
+12 
+12 11 
+13 
+13 10 
+14 
+14 12 
+15 
+15 13 
+16 
+16 14 
+17 
+17 15 
+18 
+18 16 
+18 17 
+19 
+19 18 
+20 
+21 
+22 
+22 20 
+23 
+23 19 
+24 
+24 21 
+25 
+25 23 
+26 
+26 24 
+27 
+27 22 
+28 
+28 25 
+29 
+29 26 
+30 
+30 27 
+31 
+31 28 
+32 
+32 30 
+33 
+33 29 
+34 
+34 32 
+35 
+35 31 
+36 
+36 33 
+37 
+37 34 
+37 35 
+37 36 
+38 
+38 37 
+39 
+39 38 
+40 
+40 39 
+41 
+41 40 
+42 
+42 41 
+43 
+43 42 
diff --git a/src/Nerve_GIC/example/GICvoronoi.cpp b/src/Nerve_GIC/example/GICvoronoi.cpp
index 141268a5..9bf3de5e 100644
--- a/src/Nerve_GIC/example/GICvoronoi.cpp
+++ b/src/Nerve_GIC/example/GICvoronoi.cpp
@@ -2,18 +2,17 @@
 
 void usage(int nbArgs, char * const progName) {
   std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n";
-  std::cerr << "Usage: " << progName << " filename.off threshold N\n";
-  std::cerr << "       i.e.: " << progName << " ../../../data/points/human.off 0.075 100 \n";
+  std::cerr << "Usage: " << progName << " filename.off N [--v] \n";
+  std::cerr << "       i.e.: " << progName << " ../../../../data/points/human.off 100 --v \n";
   exit(-1);  // ----- >>
 }
 
 int main(int argc, char **argv) {
-  if ((argc != 4) && (argc != 5)) usage(argc, (argv[0] - 1));
+  if ((argc != 3) && (argc != 4)) usage(argc, (argv[0] - 1));
 
   std::string off_file_name(argv[1]);
-  double threshold = atof(argv[2]);
-  int m = atoi(argv[3]);
-  bool verb = 0; if(argc == 5)  verb = 1;
+  int m = atoi(argv[2]);
+  bool verb = 0; if(argc == 4)  verb = 1;
 
   // ----------------------------------------------------------------------------
   // Init of a graph induced complex from an OFF file
diff --git a/src/Nerve_GIC/example/MapperDeltaCoord.cpp b/src/Nerve_GIC/example/MapperDeltaCoord.cpp
index c4d86caf..9445e988 100644
--- a/src/Nerve_GIC/example/MapperDeltaCoord.cpp
+++ b/src/Nerve_GIC/example/MapperDeltaCoord.cpp
@@ -2,7 +2,7 @@
 
 void usage(int nbArgs, char * const progName) {
   std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n";
-  std::cerr << "Usage: " << progName << " filename.off coordinate \n";
+  std::cerr << "Usage: " << progName << " filename.off coordinate [--v] \n";
   std::cerr << "       i.e.: " << progName << " ../../../../data/points/human.off 2 --v \n";
   exit(-1);  // ----- >>
 }
diff --git a/src/Nerve_GIC/example/MapperDeltaCoord.txt b/src/Nerve_GIC/example/MapperDeltaCoord.txt
index c1ead814..faf832c7 100644
--- a/src/Nerve_GIC/example/MapperDeltaCoord.txt
+++ b/src/Nerve_GIC/example/MapperDeltaCoord.txt
@@ -1,134 +1,143 @@
-../../../../data/points/human.off
-coordinate 2
-coordinate 2
-0.072934 0.3
-65 64
-0 -0.954816 218
-1 -0.954824 218
-2 -0.903705 57
-3 -0.9035 57
-4 -0.84796 27
-5 -0.847937 27
-6 -0.794255 25
-7 -0.791001 29
-8 -0.73369 36
-9 -0.733071 35
-10 -0.683795 40
-11 -0.682735 40
-12 -0.633796 48
-13 -0.634232 48
-14 -0.585137 64
-15 -0.585197 64
-16 -0.527997 66
-17 -0.528025 66
-18 -0.478646 93
-19 -0.478673 93
-20 -0.43674 99
-21 -0.436761 99
-22 -0.380976 70
-23 -0.380955 70
-24 -0.33528 72
-25 -0.335279 72
-26 -0.277496 112
-27 -0.232679 96
-28 -0.1842 77
-29 -0.149795 16
-30 -0.150304 16
-31 -0.119698 90
-32 -0.122291 74
-33 -0.119626 94
-34 -0.0849586 74
-35 -0.0681724 107
-36 -0.0856915 73
-37 -0.0263688 211
-38 -0.025358 49
-39 -0.025088 50
-40 0.0287585 146
-41 0.022349 34
-42 0.0213867 34
-43 0.0768417 175
-44 0.081987 34
-45 0.0819086 34
-46 0.128047 150
-47 0.129348 53
-48 0.128561 54
-49 0.173705 104
-50 0.180884 53
-51 0.180911 53
-52 0.229553 68
-53 0.233083 77
-54 0.23311 77
-55 0.283572 215
-56 0.338574 266
-57 0.381249 299
-58 0.433627 244
-59 0.485309 214
-60 0.535849 237
-61 0.588559 193
-62 0.638771 196
-63 0.692993 222
-64 0.763344 307
-0 3
-1 2
-2 4
-3 5
-4 6
-5 7
-6 9
-7 8
-8 10
-9 11
-10 13
-11 12
-12 14
-13 15
-14 17
-15 16
-16 18
-17 19
-18 20
-19 21
-20 22
-21 23
-22 24
-23 25
-24 26
-25 26
-26 27
-27 28
-28 32
-29 31
-30 33
-31 34
-32 35
-33 36
-34 39
-35 37
-36 38
-37 40
-38 41
-39 42
-40 43
-41 45
-42 44
-43 46
-44 47
-45 48
-46 49
-47 50
-48 51
-49 52
-50 53
-51 54
-52 55
-53 55
-54 55
-55 56
-56 57
-57 58
-58 59
-59 60
-60 61
-61 62
-62 63
-63 64
+Mapper Delta is of dimension 1 - 141 simplices - 71 vertices.
+Iterator on Mapper Delta simplices
+0 
+1 
+2 
+2 1 
+3 
+3 0 
+4 
+4 2 
+5 
+5 3 
+6 
+6 4 
+7 
+7 5 
+8 
+8 6 
+9 
+9 7 
+10 
+10 9 
+11 
+11 8 
+12 
+12 11 
+13 
+13 10 
+14 
+14 12 
+15 
+15 13 
+16 
+16 14 
+17 
+17 15 
+18 
+18 17 
+19 
+19 16 
+20 
+20 18 
+21 
+21 19 
+22 
+22 20 
+23 
+23 21 
+24 
+24 22 
+25 
+25 23 
+26 
+26 24 
+27 
+27 25 
+28 
+28 26 
+28 27 
+29 
+29 28 
+30 
+30 29 
+31 
+32 
+33 
+33 30 
+34 
+34 31 
+35 
+35 32 
+36 
+36 34 
+37 
+37 33 
+38 
+38 35 
+39 
+39 37 
+40 
+40 38 
+41 
+41 36 
+42 
+42 39 
+43 
+43 40 
+44 
+44 41 
+45 
+45 42 
+46 
+46 44 
+47 
+47 43 
+48 
+48 45 
+49 
+49 46 
+50 
+50 47 
+51 
+51 48 
+52 
+52 49 
+53 
+53 50 
+54 
+54 51 
+55 
+55 52 
+56 
+56 53 
+57 
+57 54 
+58 
+58 55 
+59 
+59 56 
+60 
+60 57 
+60 58 
+60 59 
+61 
+61 60 
+62 
+62 61 
+63 
+63 62 
+64 
+64 63 
+65 
+65 64 
+66 
+66 65 
+67 
+67 66 
+68 
+68 67 
+69 
+69 68 
+70 
+70 69 
diff --git a/src/Nerve_GIC/example/MapperDeltaFunc.cpp b/src/Nerve_GIC/example/MapperDeltaFunc.cpp
index 1c3a77a5..448323b1 100644
--- a/src/Nerve_GIC/example/MapperDeltaFunc.cpp
+++ b/src/Nerve_GIC/example/MapperDeltaFunc.cpp
@@ -2,7 +2,7 @@
 
 void usage(int nbArgs, char * const progName) {
   std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n";
-  std::cerr << "Usage: " << progName << " filename.off coordinate \n";
+  std::cerr << "Usage: " << progName << " filename.off function [--v] \n";
   std::cerr << "       i.e.: " << progName << " ../../../../data/points/COIL_database/lucky_cat.off ../../../../data/points/COIL_database/lucky_cat_PCA1 --v \n";
   exit(-1);  // ----- >>
 }
diff --git a/src/Nerve_GIC/example/MapperDeltaFunc.txt b/src/Nerve_GIC/example/MapperDeltaFunc.txt
index 01b0fbea..d13d0695 100644
--- a/src/Nerve_GIC/example/MapperDeltaFunc.txt
+++ b/src/Nerve_GIC/example/MapperDeltaFunc.txt
@@ -1,17 +1,10 @@
-../../../../data/points/COIL_database/lucky_cat.off
-../../../../data/points/COIL_database/lucky_cat_PCA1
-../../../../data/points/COIL_database/lucky_cat_PCA1
-2592.99 0.3
-6 6
-0 -1382.24 20
-1 385.162 3
-2 246.455 12
-3 2204.66 4
-4 2238.94 5
-5 5346.1 39
-0 1
-0 2
-1 3
-2 4
-3 5
-4 5
+Mapper Delta is of dimension 1 - 8 simplices - 4 vertices.
+Iterator on Mapper Delta simplices
+0 
+1 
+1 0 
+2 
+2 0 
+3 
+3 1 
+3 2 
diff --git a/src/Nerve_GIC/example/Nerve.cpp b/src/Nerve_GIC/example/Nerve.cpp
index a549f544..84f74625 100644
--- a/src/Nerve_GIC/example/Nerve.cpp
+++ b/src/Nerve_GIC/example/Nerve.cpp
@@ -2,7 +2,7 @@
 
 void usage(int nbArgs, char * const progName) {
   std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n";
-  std::cerr << "Usage: " << progName << " filename.off coordinate resolution gain --v \n";
+  std::cerr << "Usage: " << progName << " filename.off coordinate resolution gain [--v] \n";
   std::cerr << "       i.e.: " << progName << " ../../../../data/human.off 2 10 0.3 --v \n";
   exit(-1);  // ----- >>
 }
diff --git a/src/Nerve_GIC/example/Nerve.txt b/src/Nerve_GIC/example/Nerve.txt
index a6fad665..2a861c5f 100644
--- a/src/Nerve_GIC/example/Nerve.txt
+++ b/src/Nerve_GIC/example/Nerve.txt
@@ -1,46 +1,43 @@
-../../../../data/points/human.off
-coordinate 2
-coordinate 2
-0 0.3
-21 20
-0 -0.927412 290
-1 -0.927011 291
-2 -0.708477 127
-3 -0.709946 128
-4 -0.513775 250
-5 -0.514319 251
-6 -0.376304 247
-7 -0.376315 247
-8 -0.104886 127
-9 -0.168391 165
-10 -0.16871 165
-11 -0.105984 128
-12 -0.0174106 183
-13 0.0197303 542
-14 -0.0161482 181
-15 0.159925 388
-16 0.204922 208
-17 0.204892 208
-18 0.368798 854
-19 0.542323 764
-20 0.709288 597
-0 2
-1 3
-2 5
-3 4
-4 6
-5 7
-6 9
-7 10
-8 12
-9 13
-10 13
-11 14
-12 16
-13 15
-14 17
-15 18
-16 18
-17 18
-18 19
-19 20
+Nerve is of dimension 1 - 41 simplices - 21 vertices.
+Iterator on Nerve simplices
+0 
+1 
+2 
+2 0 
+3 
+3 1 
+4 
+4 3 
+5 
+5 2 
+6 
+6 4 
+7 
+7 5 
+8 
+9 
+9 6 
+10 
+10 7 
+11 
+12 
+12 8 
+13 
+13 9 
+13 10 
+14 
+14 11 
+15 
+15 13 
+16 
+16 12 
+17 
+17 14 
+18 
+18 15 
+18 16 
+18 17 
+19 
+19 18 
+20 
+20 19 
diff --git a/src/Nerve_GIC/example/km.py b/src/Nerve_GIC/example/km.py
index 881273ad..a1a08156 100755
--- a/src/Nerve_GIC/example/km.py
+++ b/src/Nerve_GIC/example/km.py
@@ -248,9 +248,9 @@ class KeplerMapper(object):
       if custom_tooltips is not None:
         tooltip_s = "<h2>Cluster %s</h2>"%k + " ".join(str(custom_tooltips[complex["nodes"][k][0]]).split(" "))
         if maximum == minimum:
-          tooltip_i = int(30*(custom_tooltips[complex["nodes"][k][0]]-minimum)/(maximum-minimum))
-        else:
           tooltip_i = 0
+        else:
+          tooltip_i = int(30*(custom_tooltips[complex["nodes"][k][0]]-minimum)/(maximum-minimum))
         json_s["nodes"].append({"name": str(k), "tooltip": tooltip_s, "group": 2 * int(np.log(complex["nodes"][k][2])), "color": tooltip_i})
       else:
         tooltip_s = "<h2>Cluster %s</h2>Contains %s members."%(k,len(complex["nodes"][k]))
diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h
index 7be71df3..72e08513 100644
--- a/src/Nerve_GIC/include/gudhi/GIC.h
+++ b/src/Nerve_GIC/include/gudhi/GIC.h
@@ -53,21 +53,13 @@
 #define ETA 0.001
 #define MASK 0
 
-namespace gss = Gudhi::spatial_searching;
-
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Filtration_value = Simplex_tree::Filtration_value;
 using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
 using Point = std::vector<float>;
 
-typedef CGAL::Epick_d<CGAL::Dimension_tag<4> > K;
-typedef typename K::Point_d Pointd;
-typedef std::vector<Pointd> Pointsd;
-typedef gss::Kd_tree_search<K, Pointsd> Points_ds;
-
 std::map<int, double> func;
 std::map<int, double> func_color;
-Gudhi::Points_off_reader<Point> off_reader("tmp");
 
 namespace Gudhi {
 
@@ -99,19 +91,22 @@ class Graph_induced_complex {
 
  private:
    bool verbose; // whether to display information.
+   std::vector<Point> point_cloud;
    typedef int Cover_t; // elements of cover C are indexed by integers.
    std::vector<std::vector<Cover_t> > simplices;
    std::map<int, std::vector<Cover_t> > cover;
    int maximal_dim; // maximal dimension of output simplicial complex.
    int data_dimension; // dimension of input data.
+   int n; // number of points.
    std::map<Cover_t,int> cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not.
    std::map<Cover_t,std::pair<int,double> > cover_color; // size and coloring of the vertices of the output simplicial complex.
    Simplex_tree<> st;
    std::map<int,std::vector<int> > adjacency_matrix;
+   std::vector<std::vector<double> > distances;
    int resolution_int;
    double resolution_double;
    double gain;
-   std::vector<int> subsamples;
+   std::vector<int> voronoi_subsamples;
    std::string cover_name;
    std::string point_cloud_name;
    std::string color_name;
@@ -167,8 +162,11 @@ class Graph_induced_complex {
 
  public:
    void read_point_cloud(const std::string& off_file_name){
-     off_reader = Points_off_reader<Point>(off_file_name);
+     Gudhi::Points_off_reader<Point> off_reader = Points_off_reader<Point>(off_file_name);
+     point_cloud = off_reader.get_point_cloud();
      point_cloud_name = off_file_name;
+     n = point_cloud.size();
+     data_dimension = point_cloud[0].size();
    }
 
  // *******************************************************************************************************************
@@ -207,7 +205,7 @@ class Graph_induced_complex {
     *
     */
    void set_graph_from_OFF(const std::string& off_file_name){
-     int n, numedges, numfaces, i; std::vector<int> edge(2); double x; int num; std::vector<int> simplex;
+     int numedges, numfaces, i; std::vector<int> edge(2); double x; int num; std::vector<int> simplex;
      std::ifstream input(off_file_name); std::string line; getline(input, line);
      input >> n; input >> numfaces; input >> numedges;
      i = 0;  while(i < n){input >> x; input >> x; input >> x; i++;}
@@ -242,10 +240,10 @@ class Graph_induced_complex {
     *
     */
    void set_graph_from_rips(const double& threshold){
-     Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), threshold, Euclidean_distance());
-     rips_complex_from_points.create_complex(st, 1); data_dimension = off_reader.get_point_cloud()[0].size();
+     Rips_complex rips_complex_from_points(point_cloud, threshold, Euclidean_distance());
+     rips_complex_from_points.create_complex(st, 1);
 
-     std::vector<int> empty; int n = off_reader.get_point_cloud().size();
+     std::vector<int> empty;
      for(int i = 0; i < n; i++)  adjacency_matrix.insert(std::pair<int,std::vector<int> >(i,empty));
      for (auto simplex : st.complex_simplex_range()) {
        if(st.dimension(simplex) == 1){
@@ -257,58 +255,59 @@ class Graph_induced_complex {
 
    }
 
- public: // Automatic tuning of Rips complex.
-   /** \brief Creates the graph G from a Rips complex whose threshold value is automatically tuned with subsampling.
-    *
-    * @param[in] off_file_name name of the input .OFF file.
-    * @param[in] N number of subsampling iteration (default value 100).
-    *
+ public: // Pairwise distances.
+   /** \brief Computes all pairwise distances (Euclidean norm).
     */
-   void set_graph_from_automatic_rips(int N = 100){
-
-     int n = off_reader.get_point_cloud().size();
-     int m = floor(n/pow(log(n)/log(CONSTANT),1+ETA)); m = std::min(m,n-1);
-     std::vector<int> samples(m); double delta = 0; int dim = off_reader.get_point_cloud()[0].size(); data_dimension = dim;
+   void compute_pairwise_distances(){
 
-     if(verbose)  std::cout << n << " points in R^" << dim << std::endl;
-     if(verbose)  std::cout << "Subsampling " << m << " points" << std::endl;
-
-     std::vector<std::vector<double> > dist(n); std::vector<double> dumb(n);
-     for(int i = 0; i < n; i++)  dist[i] = dumb;
-     double d;
-
-     char distances[100]; sprintf(distances,"%s_dist",(char*) point_cloud_name.c_str());
-     std::ifstream input(distances, std::ios::out | std::ios::binary);
+     double d; std::vector<double> dumb(n); for(int i = 0; i < n; i++)  distances.push_back(dumb);
+     char distance[100]; sprintf(distance,"%s_dist",(char*) point_cloud_name.c_str());
+     std::ifstream input(distance, std::ios::out | std::ios::binary);
 
      if(input.good()){
        if(verbose)  std::cout << "Reading distances..." << std::endl;
        for(int i = 0; i < n; i++){
-         for (int j = 0; j < n; j++){
-           input.read((char*) &d,8); dist[i][j] = d;
+         for (int j = i; j < n; j++){
+           input.read((char*) &d,8); distances[i][j] = d; distances[j][i] = d;
          }
        }
        input.close();
      }
+
      else{
        if(verbose)  std::cout << "Computing distances..." << std::endl;
-       input.close(); std::ofstream output(distances, std::ios::out | std::ios::binary);
+       input.close(); std::ofstream output(distance, std::ios::out | std::ios::binary);
        for(int i = 0; i < n; i++){
          if( (int) floor( 100*(i*1.0+1)/(n*1.0) ) %10 == 0  )
            if(verbose)  std::cout << "\r" << floor( 100*(i*1.0+1)/(n*1.0) ) << "%" << std::flush;
-         for (int j = 0; j < n; j++){
-           double dis = 0; for(int k = 0; k < dim; k++)  dis += pow(off_reader.get_point_cloud()[i][k]-off_reader.get_point_cloud()[j][k],2);
-           dis = std::sqrt(dis); dist[i][j] = dis; output.write((char*) &dis,8);
+         for (int j = i; j < n; j++){
+           double dis = 0; for(int k = 0; k < data_dimension; k++)
+               dis += pow(point_cloud[i][k]-point_cloud[j][k],2);
+           dis = std::sqrt(dis); distances[i][j] = dis; distances[j][i] = dis;
+           output.write((char*) &dis,8);
          }
        }
        output.close(); if(verbose)  std::cout << std::endl;
      }
 
-     for(int i = 0; i < n; i++){
-       for(int j = i; j < n; j++){
-         double dis = 0; for(int k = 0; k < dim; k++)  dis += pow(off_reader.get_point_cloud()[i][k]-off_reader.get_point_cloud()[j][k],2);
-         dist[i][j] = std::sqrt(dis); dist[j][i] = std::sqrt(dis);
-       }
-     }
+   }
+
+ public: // Automatic tuning of Rips complex.
+   /** \brief Creates the graph G from a Rips complex whose threshold value is automatically tuned with subsampling.
+    *
+    * @param[in] off_file_name name of the input .OFF file.
+    * @param[in] N number of subsampling iteration (default value 100).
+    *
+    */
+   void set_graph_from_automatic_rips(int N = 100){
+
+     int m = floor(n/pow(log(n)/log(CONSTANT),1+ETA)); m = std::min(m,n-1);
+     std::vector<int> samples(m); double delta = 0;
+
+     if(verbose)  std::cout << n << " points in R^" << data_dimension << std::endl;
+     if(verbose)  std::cout << "Subsampling " << m << " points" << std::endl;
+
+     if(distances.size() == 0)  compute_pairwise_distances();
 
      //#pragma omp parallel for
      for (int i = 0; i < N; i++){
@@ -316,7 +315,7 @@ class Graph_induced_complex {
        SampleWithoutReplacement(n,m,samples);
        double hausdorff_dist = 0;
        for (int j = 0; j < n; j++){
-         double mj = dist[j][samples[0]]; for (int k = 1; k < m; k++)  mj = std::min(mj, dist[j][samples[k]]);
+         double mj = distances[j][samples[0]]; for (int k = 1; k < m; k++)  mj = std::min(mj, distances[j][samples[k]]);
          hausdorff_dist = std::max(hausdorff_dist, mj);
        }
        delta += hausdorff_dist/N;
@@ -324,7 +323,7 @@ class Graph_induced_complex {
      }
 
      if(verbose)  std::cout << "delta = " << delta << std::endl;
-     Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), delta, Euclidean_distance());
+     Rips_complex rips_complex_from_points(point_cloud, delta, Euclidean_distance());
      rips_complex_from_points.create_complex(st, 1);
 
      std::vector<int> empty;
@@ -363,12 +362,10 @@ class Graph_induced_complex {
    /** \brief Creates the function f from the k-th coordinate of the point cloud P.
     *
     * @param[in] k coordinate to use (start at 0).
-    * @param[in] off_file_name name of the input .OFF file.
     *
     */
    void set_function_from_coordinate(const int& k){
-     int n = off_reader.get_point_cloud().size();
-     for(int i = 0; i < n; i++)  func.insert(std::pair<int,double>(i,off_reader.get_point_cloud()[i][k]));
+     for(int i = 0; i < n; i++)  func.insert(std::pair<int,double>(i,point_cloud[i][k]));
      char coordinate[100]; sprintf(coordinate, "coordinate %d", k);
      cover_name = coordinate;
    }
@@ -380,7 +377,6 @@ class Graph_induced_complex {
     *
     */
    void set_function_from_vector(const std::vector<double>& function){
-     int n = function.size();
      for(int i = 0; i < n; i++)  func.insert(std::pair<int,double>(i, function[i]));
    }
 
@@ -413,48 +409,56 @@ class Graph_induced_complex {
      cover_name = cover_file_name;
    }
 
- /* TODO: complete method with nearest geodesic neighbor
-
  public: // Set cover from Voronoi
    /** \brief Creates the cover C from the Voronoï cells of a subsampling of the point cloud.
     *
     * @param[in] m number of points in the subsample.
-    * @param[in] off_file_name name of the input .OFF file.
     *
-
+    */
    void set_cover_from_Voronoi(const int& m){
 
-     int n = off_reader.get_point_cloud().size(); data_dimension = off_reader.get_point_cloud()[0].size();
-     Pointsd pointsd(m+1); std::vector<int> samples(m); SampleWithoutReplacement(n,m,samples);
-     double* coord = new double[data_dimension]; subsamples = samples;
-
-     for(int i = 1; i <= m; i++){
-       for(int j = 0; j < data_dimension; j++)  coord[j] = off_reader.get_point_cloud()[samples[i-1]][j];
-       pointsd[i] = Pointd(data_dimension, coord+0, coord + data_dimension); cover_fct[i-1] = i-1;
-     } std::cout << std::endl;
-
-     int curr_subsample = 0;
-     for(int i = 0; i < n; i++){
-       for(int j = 0; j < data_dimension; j++)  coord[j] = off_reader.get_point_cloud()[i][j];
-       pointsd[0] = Pointd(data_dimension, coord+0, coord + data_dimension); Points_ds points_ds(pointsd);
-       auto knn_range = points_ds.query_k_nearest_neighbors(pointsd[0], 2, true);
-       Cover_t cluster; // = nearest geodesic neighbor.
-       if(cluster >= 0){ // Case where i is not a subsample point.
-         cover[i].push_back(cluster); cover_color[cluster].second += func_color[i]; cover_color[cluster].first++;
-       }
-       else{ // Case where i is a subsample point.
-         cover[i].push_back(curr_subsample); cover_color[curr_subsample].second += func_color[i];
-         cover_color[curr_subsample].first++; curr_subsample++;
+     voronoi_subsamples.resize(m); SampleWithoutReplacement(n,m,voronoi_subsamples);
+     if(distances.size() == 0)  compute_pairwise_distances();
+     std::vector<double> mindist(n); for(int j = 0; j < n; j++)  mindist[j] = std::numeric_limits<double>::max();
+
+     // Compute the geodesic distances to subsamples with Dijkstra
+     for(int i = 0; i < m; i++){
+       if(verbose)  std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl;
+       int seed = voronoi_subsamples[i];
+       std::vector<double> dist(n); std::vector<int> process(n);
+       for(int j = 0; j < n; j++){ dist[j] = std::numeric_limits<double>::max(); process[j] = j; }
+       dist[seed] = 0; int curr_size = process.size(); int min_point, min_index; double min_dist;
+       std::vector<int> neighbors; int num_neighbors;
+
+       while(curr_size > 0){
+         min_dist = std::numeric_limits<double>::max(); min_index = -1; min_point = -1;
+         for(int j = 0; j < curr_size; j++){
+           if(dist[process[j]] < min_dist){
+             min_point = process[j]; min_dist = dist[process[j]]; min_index = j;
+           }
+         }
+         assert(min_index != -1); process.erase(process.begin() + min_index);
+         assert(min_point != -1); neighbors = adjacency_matrix[min_point]; num_neighbors = neighbors.size();
+         for(int j = 0; j < num_neighbors; j++){
+           double d = dist[min_point] + distances[min_point][neighbors[j]];
+           dist[neighbors[j]] = std::min(dist[neighbors[j]], d);
+         }
+         curr_size = process.size();
        }
+
+       for(int j = 0; j < n; j++)
+         if(mindist[j] > dist[j]){
+           mindist[j] = dist[j];
+           if(cover[j].size() == 0) cover[j].push_back(i);
+           else cover[j][0] = i;
+         }
      }
 
+     for(int i = 0; i < n; i++){ cover_color[cover[i][0]].second += func_color[i]; cover_color[cover[i][0]].first++; }
      for(int i = 0; i < m; i++)  cover_color[i].second /= cover_color[i].first;
-
-     delete [] coord;
      maximal_dim = m-1; cover_name = "Voronoi";
 
    }
-   */
 
  public: // Automatic tuning of resolution for Mapper Delta.
    /** \brief Computes the optimal length of intervals for a Mapper Delta.
@@ -653,8 +657,7 @@ class Graph_induced_complex {
     *
     */
    void set_color_from_coordinate(int k = 0){
-     int n = off_reader.get_point_cloud().size();
-     for(int i = 0; i < n; i++)  func_color[i] = off_reader.get_point_cloud()[i][k];
+     for(int i = 0; i < n; i++)  func_color.insert(std::pair<int,double>(i, point_cloud[i][k]));
      char coordinate[100]; sprintf(coordinate, "coordinate %d", k);
      color_name = coordinate;
    }
@@ -670,15 +673,15 @@ class Graph_induced_complex {
    }
 
  public: // Create a .dot file that can be compiled with neato to produce a .pdf file.
-   /** \brief Creates a .dot file for neato once the simplicial complex is computed to get a .pdf output.
-    * For Mapper Delta only.
+   /** \brief Creates a .dot file for neato once the simplicial complex is computed to get a nice visualization
+    * of its 1-skeleton in a .pdf file.
     */
    void plot_with_neato(){
      char mapp[11] = "SC.dot";   std::ofstream graphic(mapp); graphic << "graph Mapper {" << std::endl;
      double maxv, minv; maxv = std::numeric_limits<double>::min(); minv = std::numeric_limits<double>::max();
      for (std::map<Cover_t,std::pair<int,double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++){
        maxv = std::max(maxv, iit->second.second);  minv = std::min(minv, iit->second.second);
-     } //std::cout << minv << " " << maxv << std::endl;
+     }
      int k = 0; std::vector<int> nodes; nodes.clear();
      for (std::map<Cover_t,std::pair<int,double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++){
        if(iit->second.first > MASK){
@@ -701,8 +704,8 @@ class Graph_induced_complex {
    }
 
  public: // Create a .txt file that can be compiled with KeplerMapper to produce a .html file.
-   /** \brief Creates a .html file for KeplerMapper once the simplicial complex is computed to get a nice visualization
-    * of its 1_skeleton in browser.
+   /** \brief Creates a .txt file for KeplerMapper once the simplicial complex is computed to get a nice visualization
+    * of its 1-skeleton in browser.
     */
    void plot_with_KeplerMapper(){
 
@@ -732,16 +735,16 @@ class Graph_induced_complex {
      if(systemRet == -1)  std::cout << "Visualization failed. Do you have python and firefox?" << std::endl;
    }
 
- /*
+
  public: // Create a .off file that can be visualized with Geomview.
    /** \brief Creates a .off file for visualization with Geomview.
     * For GIC computed with Voronoi only.
-    *
+    */
    void plot_with_Geomview(){
 
      assert(data_dimension <= 3);
-     char mapp[11] = "SC.off";  std::ofstream graphic(mapp);
-     graphic << "OFF" << std::endl; int m = subsamples.size(); int numedges = 0; int numfaces = 0;
+     char gic[11] = "SC.off";  std::ofstream graphic(gic);
+     graphic << "OFF" << std::endl; int m = voronoi_subsamples.size(); int numedges = 0; int numfaces = 0;
      std::vector<std::vector<int> > edges, faces; int numsimplices = simplices.size();
      for (int i = 0; i < numsimplices; i++) {
        if(simplices[i].size() == 2){ numedges++;
@@ -752,26 +755,18 @@ class Graph_induced_complex {
        }
      }
      graphic << m << " " << numedges + numfaces << std::endl;
-     for(int i = 0; i < m; i++)  graphic << off_reader.get_point_cloud()[subsamples[i]][0] << " " \
-                                         << off_reader.get_point_cloud()[subsamples[i]][1] << " " \
-                                         << off_reader.get_point_cloud()[subsamples[i]][2] << std::endl;
+     for(int i = 0; i < m; i++)  graphic << point_cloud[voronoi_subsamples[i]][0] << " " \
+                                         << point_cloud[voronoi_subsamples[i]][1] << " " \
+                                         << point_cloud[voronoi_subsamples[i]][2] << std::endl;
      for(int i = 0; i < numedges; i++)  graphic << 2 << " " << edges[i][0] << " " << edges[i][1] << std::endl;
      for(int i = 0; i < numfaces; i++)  graphic << 3 << " " << faces[i][0] << " " << faces[i][1] << " " << faces[i][2] << std::endl;
      graphic.close();
 
-     char cl[11] = "cloud.off";  std::ofstream cloud(cl);
-     cloud << "COFF" << std::endl << 4706 << " " << 9408 << std::endl;
-     for(int i = 0; i < off_reader.get_point_cloud().size(); i++)
-       cloud << off_reader.get_point_cloud()[i][0] << " " << off_reader.get_point_cloud()[i][1] << " " << off_reader.get_point_cloud()[i][2] << " " <<\
-                                                      cover[i][0]*1.0/(maximal_dim+1) << " " <<\
-                                                      0 << " " <<\
-                                                      cover[i][0]*1.0/(maximal_dim+1) << " " << 0.5 << std::endl;
-
      char command[100]; sprintf(command, "geomview SC.off");
      int systemRet = system(command);
      if(systemRet == -1)  std::cout << "Visualization failed. Do you have geomview?" << std::endl;
 
-   }*/
+   }
 
  // *******************************************************************************************************************
  // *******************************************************************************************************************
diff --git a/src/Nerve_GIC/test/test_GIC.cpp b/src/Nerve_GIC/test/test_GIC.cpp
index baf494e3..e8051cf0 100644
--- a/src/Nerve_GIC/test/test_GIC.cpp
+++ b/src/Nerve_GIC/test/test_GIC.cpp
@@ -53,11 +53,10 @@ BOOST_AUTO_TEST_CASE(check_nerve) {
 
 BOOST_AUTO_TEST_CASE(check_GICMAP) {
 
-  Gudhi::graph_induced_complex::Graph_induced_complex GIC; GIC.set_verbose(1);
-  std::string cloud_file_name("data/cloud"); GIC.read_point_cloud(cloud_file_name);
+  Gudhi::graph_induced_complex::Graph_induced_complex GIC;
+  std::string cloud_file_name("data/cloud"); GIC.read_point_cloud(cloud_file_name); GIC.set_color_from_coordinate();
   std::string graph_file_name("data/graph"); GIC.set_graph_from_file(graph_file_name);
   std::string cover_file_name("data/cover"); GIC.set_cover_from_file(cover_file_name);
-  GIC.set_color_from_coordinate();
   GIC.find_GICMAP_simplices_with_functional_minimal_cover(); Simplex_tree stree; GIC.create_complex(stree);
 
   BOOST_CHECK(stree.num_vertices() == 3);
@@ -68,10 +67,9 @@ BOOST_AUTO_TEST_CASE(check_GICMAP) {
 BOOST_AUTO_TEST_CASE(check_GICcover) {
 
   Gudhi::graph_induced_complex::Graph_induced_complex GIC;
-  std::string cloud_file_name("data/cloud"); GIC.read_point_cloud(cloud_file_name);
+  std::string cloud_file_name("data/cloud"); GIC.read_point_cloud(cloud_file_name); GIC.set_color_from_coordinate();
   std::string graph_file_name("data/graph"); GIC.set_graph_from_file(graph_file_name);
   std::string cover_file_name("data/cover"); GIC.set_cover_from_file(cover_file_name);
-  GIC.set_color_from_coordinate();
   GIC.find_GIC_simplices(); Simplex_tree stree; GIC.create_complex(stree);
 
   BOOST_CHECK(stree.num_vertices() == 3);
@@ -79,12 +77,11 @@ BOOST_AUTO_TEST_CASE(check_GICcover) {
   BOOST_CHECK(stree.dimension() == 2);
 }
 
-/*
+
 BOOST_AUTO_TEST_CASE(check_GICvoronoi) {
 
   Gudhi::graph_induced_complex::Graph_induced_complex GIC;
-  std::string cloud_file_name("data/cloud"); GIC.read_point_cloud(cloud_file_name);
-  GIC.set_color_from_coordinate();
+  std::string cloud_file_name("data/cloud"); GIC.read_point_cloud(cloud_file_name); GIC.set_color_from_coordinate();
   std::string graph_file_name("data/graph"); GIC.set_graph_from_file(graph_file_name);
   GIC.set_cover_from_Voronoi(2);
   GIC.find_GIC_simplices(); Simplex_tree stree; GIC.create_complex(stree);
@@ -93,5 +90,4 @@ BOOST_AUTO_TEST_CASE(check_GICvoronoi) {
   BOOST_CHECK((stree.num_simplices()-stree.num_vertices()) == 1);
   BOOST_CHECK(stree.dimension() == 1);
 }
-*/