diff options
author | mcarrier <mcarrier@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2017-01-31 16:34:26 +0000 |
---|---|---|
committer | mcarrier <mcarrier@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2017-01-31 16:34:26 +0000 |
commit | 3625c2fc5c5b3a0f0d3c00abe41d856be147a214 (patch) | |
tree | f2ec6ac52cfebddb501b0be4dbee2b1e8ab46f42 /src/Nerve_GIC | |
parent | ab2e89665e16be2dcd36e55df8543d4ff5974675 (diff) |
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Nerve_GIC@2036 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: f4998927b72b677f51727fb4efa5ea9fe0c0916f
Diffstat (limited to 'src/Nerve_GIC')
-rw-r--r-- | src/Nerve_GIC/example/simple_GIC.cpp | 14 | ||||
-rw-r--r-- | src/Nerve_GIC/include/gudhi/GIC.h | 173 |
2 files changed, 169 insertions, 18 deletions
diff --git a/src/Nerve_GIC/example/simple_GIC.cpp b/src/Nerve_GIC/example/simple_GIC.cpp index d9297d04..a2b021f6 100644 --- a/src/Nerve_GIC/example/simple_GIC.cpp +++ b/src/Nerve_GIC/example/simple_GIC.cpp @@ -2,17 +2,19 @@ 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 cover [ouput_file.txt]\n"; + std::cerr << "Usage: " << progName << " filename.off threshold function resolution gain [ouput_file.txt]\n"; std::cerr << " i.e.: " << progName << " ../../data/points/test.off 1.5 test_cov \n"; exit(-1); // ----- >> } int main(int argc, char **argv) { - if ((argc != 4) && (argc != 5)) usage(argc, (argv[0] - 1)); + if ((argc != 6) && (argc != 7)) usage(argc, (argv[0] - 1)); std::string off_file_name(argv[1]); double threshold = atof(argv[2]); - std::string cover_file_name(argv[3]); + std::string function_file_name(argv[3]); + double resolution = atof(argv[4]); + double gain = atof(argv[5]); // Type definitions using Graph_t = boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS,\ @@ -24,13 +26,15 @@ int main(int argc, char **argv) { // ---------------------------------------------------------------------------- Gudhi::graph_induced_complex::Graph_induced_complex GIC; - GIC.set_graph_simplex_tree(threshold, off_file_name); GIC.set_cover(cover_file_name); GIC.find_simplices(); + GIC.set_graph_simplex_tree(threshold, off_file_name); + GIC.set_cover_from_function(function_file_name,resolution,gain,0); + GIC.find_simplices(); Simplex_tree stree; GIC.create_complex(stree); std::streambuf* streambufffer; std::ofstream ouput_file_stream; - if (argc == 5) { + if (argc == 7) { ouput_file_stream.open(std::string(argv[4])); streambufffer = ouput_file_stream.rdbuf(); } else { diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index 91f0327e..1c88f973 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -50,12 +50,7 @@ using Filtration_value = Simplex_tree::Filtration_value; using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; using Point = std::vector<float>; -bool simplex_comp(const std::vector<int>& s1, const std::vector<int>& s2){ - if(s1.size() == s2.size()){ - return s1[0] < s2[0]; - } - else return s1.size() < s2.size(); -} +std::map<int, double> func; namespace Gudhi { @@ -82,7 +77,7 @@ class Graph_induced_complex { std::vector<std::vector<Cover_t> > cliques; private: - std::map<int, std::vector<Cover_t> > Cov; + std::map<int, std::vector<Cover_t> > cover; private: int maximal_dim; @@ -90,6 +85,33 @@ class Graph_induced_complex { private: Simplex_tree<> st; + + // Simplex comparator + private: + bool simplex_comp(const std::vector<Cover_t>& s1, const std::vector<Cover_t>& s2){ + if(s1.size() == s2.size()){ + return s1[0] < s2[0]; + } + else return s1.size() < s2.size(); + } + + // Point comparator + private: + static bool functional_comp(const int& a, const int& b){ + if(func[a] == func[b]) return a < b; + else return func[a] < func[b]; + } + + private: + void dfs(std::map<int,std::vector<int> >& G, const int& p, std::vector<int>& cc, std::map<int,bool>& visit){ + cc.push_back(p); + visit[p] = true; int neighb = G[p].size(); + for (int i = 0; i < neighb; i++) + if ( visit.find(G[p][i]) != visit.end() ) + if( !(visit[G[p][i]]) ) + dfs(G,G[p][i],cc,visit); + } + public: void set_cover(const std::string& cover_file_name){ int vertex_id = 0; Cover_t cov; std::vector<Cover_t> cov_elts, cov_number; @@ -97,7 +119,7 @@ class Graph_induced_complex { while(std::getline(input,line)){ cov_elts.clear(); std::stringstream stream(line); while(stream >> cov){cov_elts.push_back(cov); cov_number.push_back(cov);} - Cov.insert(std::pair<int, std::vector<Cover_t> >(vertex_id, cov_elts)); vertex_id++; + cover.insert(std::pair<int, std::vector<Cover_t> >(vertex_id, cov_elts)); vertex_id++; } std::vector<Cover_t>::iterator it; std::sort(cov_number.begin(),cov_number.end()); it = std::unique(cov_number.begin(),cov_number.end()); @@ -106,6 +128,129 @@ class Graph_induced_complex { } public: + void set_cover_from_function(const std::string& func_file_name, const double& resolution, const double& gain, const bool& token){ + + // Read function values and compute min and max + int vertex_id = 0; double f; std::vector<double> range; double maxf, minf; minf = std::numeric_limits<float>::max(); maxf = std::numeric_limits<float>::min(); + std::ifstream input(func_file_name); std::string line; + while(std::getline(input,line)){ + std::stringstream stream(line); + stream >> f; range.push_back(f); minf = std::min(minf, f); maxf = std::max(maxf, f); + func.insert(std::pair<int,double>(vertex_id, f)); vertex_id++; + } + int num_pts = func.size(); + + // Compute cover of im(f) + std::vector<std::pair<double,double> > intervals; int res; + if(!token){ + double incr = (maxf-minf)/resolution; double x = minf; double alpha = (incr*gain)/(2-2*gain); + double y = minf + incr + alpha; std::pair<double, double> interm(x,y); intervals.push_back(interm); + for(int i = 1; i < resolution-1; i++){ + x = minf + i*incr - alpha; + y = minf + (i+1)*incr + alpha; + std::pair<double, double> inter(x,y); intervals.push_back(inter); + } + x = minf + (resolution-1)*incr - alpha; y = maxf; + std::pair<double, double> interM(x,y); intervals.push_back(interM); res = intervals.size(); + } + else{ + double x = minf; double y = x + resolution; + while(y <= maxf && maxf - (y-gain*resolution) >= resolution){ + std::pair<double, double> inter(x,y); intervals.push_back(inter); + x = y - gain*resolution; + y = x + resolution; + } + std::pair<double, double> interM(x,maxf); intervals.push_back(interM); res = intervals.size(); + } + /*for (int i = 0; i < res; i++) + std::cout << " " << intervals[i].first << " " << intervals[i].second << " " << std::endl;*/ + //std::cout << maxf << " " << minf << std::endl; + + // Sort points according to function values + std::vector<int> points(num_pts); for(int i = 0; i < num_pts; i++) points[i] = i; + std::sort(points.begin(),points.end(),functional_comp); + + /*for(int i = 0; i < num_pts; i++){ + std::cout << points[i] << "--" << func[points[i]] << " "; + } + std::cout << std::endl;*/ + + // Build adjacency matrix + std::map<int,std::vector<int> > G; std::vector<int> empty; + for(int i = 0; i < num_pts; i++) G.insert(std::pair<int,std::vector<int> >(points[i],empty)); + for (auto simplex : st.complex_simplex_range()) { + if(st.dimension(simplex) == 1){ + std::vector<int> vertices; + for(auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); + G[vertices[0]].push_back(vertices[1]); G[vertices[1]].push_back(vertices[0]); + } + } + + int id = 0; int pos = 0; + for(int i = 0; i < res; i++){ + + // Find points in the preimage + std::map<int,std::vector<int> > prop; prop.clear(); + std::pair<double, double> inter1 = intervals[i]; + //std::map<int,std::vector<int> >::iterator pos = G.begin(); + //std::map<int,std::vector<int> >::iterator tmp = pos; + int tmp = pos; + + if(i != res-1){ + if(i != 0){ + std::pair<double, double> inter3 = intervals[i-1]; + while(func[points[tmp]] < inter3.second && tmp != num_pts){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],G[points[tmp]])); + tmp++; + } + } + std::pair<double, double> inter2 = intervals[i+1]; + while(func[points[tmp]] < inter2.first && tmp != num_pts){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],G[points[tmp]])); //std::cout << func[points[tmp]] << " "; + tmp++; + } + pos = tmp; + while(func[points[tmp]] < inter1.second && tmp != num_pts){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],G[points[tmp]])); //std::cout << func[points[tmp]] << " "; + tmp++; + } + + } + else{ + std::pair<double, double> inter3 = intervals[i-1]; + while(func[points[tmp]] < inter3.second && tmp != num_pts){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],G[points[tmp]])); + tmp++; + } + while(tmp != num_pts){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],G[points[tmp]])); //std::cout << func[points[tmp]] << " "; + tmp++; + } + + } + + //std::cout << i << std::endl; + //std::cout << prop.empty() << std::endl; + // Compute the connected components with DFS + std::map<int,bool> visit; + for(std::map<int, std::vector<int> >::iterator it = prop.begin(); it != prop.end(); it++) + visit.insert(std::pair<int,bool>(it->first, false)); + if (!(prop.empty())){ + for(std::map<int, std::vector<int> >::iterator it = prop.begin(); it != prop.end(); it++){ + if ( !(visit[it->first]) ){ + std::vector<int> cc; cc.clear(); + dfs(prop,it->first,cc,visit); int cci = cc.size(); + for(int i = 0; i < cci; i++) cover[cc[i]].push_back(id); + //std::cout << " " << id << std::endl; id++; + } + } + } + + } + + } + + public: void set_graph_simplex_tree(const double& threshold, const std::string& off_file_name){ Points_off_reader<Point> off_reader(off_file_name); Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), threshold, Euclidean_distance()); @@ -134,13 +279,13 @@ class Graph_induced_complex { std::vector<Cover_t>::iterator it; std::sort(simplex.begin(),simplex.end()); it = std::unique(simplex.begin(),simplex.end()); simplex.resize(std::distance(simplex.begin(),it)); - cliques.push_back(simplex); + cliques.push_back(simplex);std::cout << "test2" << std::endl; } } else{ int num_clus = cover_elts[token].size(); for(int i = 0; i < num_clus; i++){ - std::vector<Cover_t> simplex = simplex_tmp; simplex.push_back(cover_elts[token][i]); + std::vector<Cover_t> simplex = simplex_tmp; simplex.push_back(cover_elts[token][i]);std::cout << "test1" << std::endl; find_all_simplices(cliques, cover_elts, ++token, simplex); } } @@ -149,14 +294,14 @@ class Graph_induced_complex { public: void find_simplices() { + std::cout << "Removing edges" << std::endl; // Find IDs of edges to remove std::vector<int> simplex_to_remove; int simplex_id = 0; for (auto simplex : st.complex_simplex_range()) { if(st.dimension(simplex) == 1){ std::vector<std::vector<Cover_t> > comp; - for(auto vertex : st.simplex_vertex_range(simplex)) comp.push_back(Cov[vertex]); + for(auto vertex : st.simplex_vertex_range(simplex)) comp.push_back(cover[vertex]); if(comp[0].size() == 1 && comp[1].size() == 1 && comp[0][0] == comp[1][0]) simplex_to_remove.push_back(simplex_id); - } simplex_id++; } @@ -171,14 +316,16 @@ class Graph_induced_complex { simplex++; } + std::cout << "Find cliques" << std::endl; // Build the Simplex Tree corresponding to the graph st.expansion(maximal_dim); + std::cout << "Find simplices" << std::endl; // Find simplices of GIC cliques.clear(); for (auto simplex : st.complex_simplex_range()) { std::vector<std::vector<Cover_t> > cover_elts; int token = 0; std::vector<Cover_t> sim; - for (auto vertex : st.simplex_vertex_range(simplex)) cover_elts.push_back(Cov[vertex]); + for (auto vertex : st.simplex_vertex_range(simplex)) cover_elts.push_back(cover[vertex]); find_all_simplices(cliques,cover_elts,token,sim); } std::vector<std::vector<Cover_t> >::iterator it; |