diff options
author | mcarrier <mcarrier@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2017-05-15 14:33:05 +0000 |
---|---|---|
committer | mcarrier <mcarrier@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2017-05-15 14:33:05 +0000 |
commit | 7d6b227a4529c0b6f8be899f613b1299d73160b5 (patch) | |
tree | e0aa9be90e9267dc8d99970fb6f835824faa96bf /src/Nerve_GIC/include | |
parent | 494907f1b452974625e4e46dff8bc59ffde66b4b (diff) |
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Nerve_GIC@2424 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: 93dd0d2099fce5d2b5b05f2ddc22cfebecac14fc
Diffstat (limited to 'src/Nerve_GIC/include')
-rw-r--r-- | src/Nerve_GIC/include/gudhi/GIC.h | 201 |
1 files changed, 98 insertions, 103 deletions
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; - }*/ + } // ******************************************************************************************************************* // ******************************************************************************************************************* |