From da5114d23e2a42e45805773f60520059c7f8cf22 Mon Sep 17 00:00:00 2001 From: skachano Date: Tue, 26 May 2015 08:51:43 +0000 Subject: witness complex in a torus git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@594 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: af4ef3a1233aec08d21f635bb3606cfc5d19c4d7 --- src/Witness_complex/example/CMakeLists.txt | 6 +- .../example/witness_complex_flat_torus.cpp | 303 ++++++++++++++++++--- .../example/witness_complex_perturbations.cpp | 4 +- 3 files changed, 270 insertions(+), 43 deletions(-) (limited to 'src/Witness_complex') diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt index 49c8c0bb..e3768138 100644 --- a/src/Witness_complex/example/CMakeLists.txt +++ b/src/Witness_complex/example/CMakeLists.txt @@ -80,9 +80,9 @@ if(CGAL_FOUND) add_executable ( witness_complex_perturbations witness_complex_perturbations.cpp ) target_link_libraries(witness_complex_perturbations ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) add_test(witness_complex_perturbations ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_perturbations ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) - #add_executable ( range_queries range_queries.cpp ) - #target_link_libraries(range_queries ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) - #add_test(range_queries ${CMAKE_CURRENT_BINARY_DIR}/range_queries ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + add_executable ( witness_complex_flat_torus witness_complex_flat_torus.cpp ) + target_link_libraries(witness_complex_flat_torus ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test(witness_complex_flat_torus ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_flat_torus ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) else() message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") endif() diff --git a/src/Witness_complex/example/witness_complex_flat_torus.cpp b/src/Witness_complex/example/witness_complex_flat_torus.cpp index afb42807..c1fb5082 100644 --- a/src/Witness_complex/example/witness_complex_flat_torus.cpp +++ b/src/Witness_complex/example/witness_complex_flat_torus.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -50,7 +51,8 @@ #include #include #include -#include +#include + #include #include @@ -60,11 +62,6 @@ using namespace Gudhi; //using namespace boost::filesystem; -typedef std::vector< Vertex_handle > typeVectorVertex; - -//typedef std::pair typeSimplex; -//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; - typedef CGAL::Epick_d K; typedef K::FT FT; typedef K::Point_d Point_d; @@ -72,10 +69,197 @@ typedef CGAL::Search_traits< FT, Point_d, typename K::Cartesian_const_iterator_d, typename K::Construct_cartesian_const_iterator_d> Traits_base; +typedef CGAL::Euclidean_distance Euclidean_distance; + +/** + * \brief Class of distance in a flat torus in dimaension D + * + */ +class Torus_distance : public Euclidean_distance { + +public: + typedef FT FT; + typedef Point_d Point_d; + typedef Point_d Query_item; + typedef typename CGAL::Dynamic_dimension_tag D; + + FT transformed_distance(Query_item q, Point_d p) const + { + FT distance = FT(0); + FT coord = FT(0); + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1), pit = construct_it(p); + for(; qit != qe; qit++, pit++) + { + coord = sqrt(((*qit)-(*pit))*((*qit)-(*pit))); + if (coord*coord <= (1-coord)*(1-coord)) + distance += coord*coord; + else + distance += (1-coord)*(1-coord); + } + return distance; + } + + FT min_distance_to_rectangle(const Query_item& q, + const CGAL::Kd_tree_rectangle& r) const { + FT distance = FT(0); + FT dist1, dist2; + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + for(unsigned int i = 0;qit != qe; i++, qit++) + { + if((*qit) < r.min_coord(i)) + { + dist1 = (r.min_coord(i)-(*qit)); + dist2 = (1 - r.max_coord(i)+(*qit)); + if (dist1 < dist2) + distance += dist1*dist1; + else + distance += dist2*dist2; + } + else if ((*qit) > r.max_coord(i)) + { + dist1 = (1 - (*qit)+r.min_coord(i)); + dist2 = ((*qit) - r.max_coord(i)); + if (dist1 < dist2) + distance += dist1*dist1; + else + distance += dist2*dist2; + } + } + return distance; + } + + FT min_distance_to_rectangle(const Query_item& q, + const CGAL::Kd_tree_rectangle& r, + std::vector& dists) const { + FT distance = FT(0); + FT dist1, dist2; + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + for(unsigned int i = 0;qit != qe; i++, qit++) + { + if((*qit) < r.min_coord(i)) + { + dist1 = (r.min_coord(i)-(*qit)); + dist2 = (1 - r.max_coord(i)+(*qit)); + if (dist1 < dist2) + { + dists[i] = dist1; + distance += dist1*dist1; + } + else + { + dists[i] = dist2; + distance += dist2*dist2; + } + } + else if ((*qit) > r.max_coord(i)) + { + dist1 = (1 - (*qit)+r.min_coord(i)); + dist2 = ((*qit) - r.max_coord(i)); + if (dist1 < dist2) + { + dists[i] = dist1; + distance += dist1*dist1; + } + else + { + dists[i] = dist2; + distance += dist2*dist2; + } + } + }; + return distance; + } + + FT max_distance_to_rectangle(const Query_item& q, + const CGAL::Kd_tree_rectangle& r) const { + FT distance=FT(0); + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + for(unsigned int i = 0;qit != qe; i++, qit++) + { + if (FT(1) <= (r.min_coord(i)+r.max_coord(i))) + if ((r.max_coord(i)+r.min_coord(i)-FT(1))/FT(2.0) <= (*qit) && + (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) + distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); + else + distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); + else + if ((FT(1)-r.max_coord(i)-r.min_coord(i))/FT(2.0) <= (*qit) || + (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) + distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); + else + distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); + } + return distance; + } + + FT max_distance_to_rectangle(const Query_item& q, + const CGAL::Kd_tree_rectangle& r, + std::vector& dists) const { + FT distance=FT(0); + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + for(unsigned int i = 0;qit != qe; i++, qit++) + { + if (FT(1) <= (r.min_coord(i)+r.max_coord(i))) + if ((r.max_coord(i)+r.min_coord(i)-FT(1))/FT(2.0) <= (*qit) && + (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) + { + dists[i] = r.max_coord(i)-(*qit); + distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); + } + else + { + dists[i] = sqrt(((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i))); + distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); + } + else + if ((FT(1)-r.max_coord(i)-r.min_coord(i))/FT(2.0) <= (*qit) || + (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) + { + dists[i] = sqrt((r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit))); + distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); + + } + else + { + dists[i] = (*qit)-r.min_coord(i); + distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); + } + } + return distance; + } + + inline FT transformed_distance(FT d) const { + return d*d; + } + + inline FT inverse_of_transformed_distance(FT d) const { + return sqrt(d); + } + +}; + + +typedef std::vector< Vertex_handle > typeVectorVertex; + +//typedef std::pair typeSimplex; +//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; + typedef CGAL::Search_traits_adapter< std::ptrdiff_t, Point_d*, Traits_base> STraits; //typedef K TreeTraits; -typedef CGAL::Orthogonal_k_neighbor_search K_neighbor_search; +//typedef CGAL::Distance_adapter Euclidean_adapter; +//typedef CGAL::Kd_tree Kd_tree; +typedef CGAL::Orthogonal_k_neighbor_search> K_neighbor_search; typedef K_neighbor_search::Tree Tree; typedef K_neighbor_search::Distance Distance; typedef K_neighbor_search::iterator KNS_iterator; @@ -87,9 +271,12 @@ typedef CGAL::Fuzzy_sphere Fuzzy_sphere; typedef std::vector Point_Vector; -typedef CGAL::Euclidean_distance Euclidean_distance; //typedef K::Equal_d Equal_d; +typedef CGAL::Random_points_in_cube_d Random_cube_iterator; typedef CGAL::Random_points_in_ball_d Random_point_iterator; + + + /** * \brief Customized version of read_points * which takes into account a possible nbP first line @@ -118,6 +305,20 @@ read_points_cust ( std::string file_name , Point_Vector & points) in_file.close(); } +void generate_points_grid(Point_Vector& W, int width, int D) +{ + +} + +void generate_points_random_box(Point_Vector& W, int nbP, int dim) +{ + Random_cube_iterator rp(dim, 1); + for (int i = 0; i < nbP; i++) + { + W.push_back(*rp++); + } +} + /* void read_points_to_tree (std::string file_name, Tree& tree) { @@ -169,12 +370,14 @@ void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, int chosen_landmark; //std::pair res = std::make_pair(L_i.begin(),false); Point_d* p; - srand(24660); + CGAL::Random rand; for (int i = 0; i < nbL; i++) { // while (!res.second) // { - chosen_landmark = rand()%nbP; + chosen_landmark = rand.get_int(0,nbP); + //rand++; + //std::cout << "Chose " << chosen_landmark << std::endl; p = &W[chosen_landmark]; //L_i.emplace(chosen_landmark,i); // } @@ -191,9 +394,11 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< int nbP = W.size(); int nbL = landmarks.size(); //Point_Vector landmarks_ = landmarks; - Euclidean_distance ed; + Torus_distance ed; //Equal_d ed; + //Point_d p1(std::vector({0.8,0.8})), p2(std::vector({0.1,0.1})); FT lambda = ed.transformed_distance(landmarks[0],landmarks[1]); + //std::cout << "Lambda=" << lambda << std::endl; //FT lambda = 0.1;//Euclidean_distance(); std::vector< std::vector > WL(nbP); Tree L(boost::counting_iterator(0), @@ -215,7 +420,8 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< //std::cout << "Safely constructed a point\n"; ////Search D+1 nearest neighbours from the tree of landmarks L K_neighbor_search search(L, w, D+1, FT(0), true, - CGAL::Distance_adapter>(&(landmarks[0])) ); + //CGAL::Distance_adapter(&(landmarks[0])) ); + CGAL::Distance_adapter(&(landmarks[0])) ); //std::cout << "Safely found nearest landmarks\n"; for(K_neighbor_search::iterator it = search.begin(); it != search.end(); ++it) { @@ -236,10 +442,10 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< } } //std::cout << "\n"; - /* + std::string out_file = "wl_result"; write_wl(out_file,WL); - */ + //******************** Constructng a witness complex std::cout << "Entered witness complex construction\n"; Witness_complex<> witnessComplex; @@ -249,16 +455,16 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< std::cout << "Entered bad links\n"; std::set< int > perturbL; int count_badlinks = 0; - std::cout << "Bad links around "; + //std::cout << "Bad links around "; for (auto u: witnessComplex.complex_vertex_range()) if (!witnessComplex.has_good_link(u)) { //std::cout << "Landmark " << u << " start!" << std::endl; //perturbL.insert(u); count_badlinks++; - std::cout << u << " "; + //std::cout << u << " "; Point_d& l = landmarks[u]; - Fuzzy_sphere fs(l, sqrt(lambda)*2, 0, STraits(&(landmarks[0]))); + Fuzzy_sphere fs(l, sqrt(lambda)*3, 0, STraits(&(landmarks[0]))); L.search(std::insert_iterator>(perturbL,perturbL.begin()),fs); //L.search(std::inserter(perturbL,perturbL.begin()),fs); //L.search(std::ostream_iterator(std::cout,"\n"),fs); @@ -270,13 +476,20 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< for (auto u: perturbL) { - Random_point_iterator rp(D,sqrt(lambda)/2); + Random_point_iterator rp(D,sqrt(lambda)/16); //std::cout << landmarks[u] << std::endl; std::vector point; for (int i = 0; i < D; i++) { - point.push_back(W[landmarks_ind[u]][i] + (*rp)[i]); + //FT coord = W[landmarks_ind[u]][i] + (*rp)[i]; + FT coord = landmarks[u][i] + (*rp)[i]; + if (coord > 1) + point.push_back(coord-1); + else if (coord < -1) + point.push_back(coord+1); + else + point.push_back(coord); } landmarks[u] = Point_d(point); //std::cout << landmarks[u] << std::endl; @@ -284,14 +497,13 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< //std::cout << "landmark[0][0] after" << landmarks[0][0] << std::endl; std::cout << "lambda=" << lambda << std::endl; - // Write the WL matrix in a file + + //std::cout << "WL size" << WL.size() << std::endl; /* - mkdir("output", S_IRWXU); - const size_t last_slash_idx = file_name.find_last_of("/"); - if (std::string::npos != last_slash_idx) - { - file_name.erase(0, last_slash_idx + 1); - } + std::cout << "L:" << std::endl; + for (int i = 0; i < landmarks.size(); i++) + std::cout << landmarks[i] << std::endl; + */ char buffer[100]; int i = sprintf(buffer,"stree_result.txt"); @@ -303,39 +515,47 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< witnessComplex.st_to_file(ofs); ofs.close(); } - */ + return count_badlinks; } int main (int argc, char * const argv[]) { - if (argc != 3) + + if (argc != 4) { std::cerr << "Usage: " << argv[0] - << " path_to_point_file nbL \n"; + << " nbP nbL dim\n"; return 0; } /* boost::filesystem::path p; - for (; argc > 2; --argc, ++argv) p /= argv[1]; */ - std::string file_name = argv[1]; - int nbL = atoi(argv[2]); + int nbP = atoi(argv[1]); + int nbL = atoi(argv[2]); + int dim = atoi(argv[3]); //clock_t start, end; //Construct the Simplex Tree //Witness_complex<> witnessComplex; std::cout << "Let the carnage begin!\n"; Point_Vector point_vector; - read_points_cust(file_name, point_vector); + //read_points_cust(file_name, point_vector); + generate_points_random_box(point_vector, nbP, dim); + /* + for (auto &p: point_vector) + { + assert(std::count(point_vector.begin(),point_vector.end(),p) == 1); + } + */ //std::cout << "Successfully read the points\n"; //witnessComplex.setNbL(nbL); // witnessComplex.witness_complex_from_points(point_vector); - int nbP = point_vector.size(); + //int nbP = point_vector.size(); //std::vector > WL(nbP); //std::set L; Point_Vector L; @@ -344,11 +564,18 @@ int main (int argc, char * const argv[]) //start = clock(); //witnessComplex.landmark_choice_by_furthest_points(point_vector, point_vector.size(), WL); landmark_choice(point_vector, nbP, nbL, L, chosen_landmarks); - int bl = 1; - for (int i = 0; bl != 0; i++) + for (auto i: chosen_landmarks) + { + assert(std::count(chosen_landmarks.begin(),chosen_landmarks.end(),i) == 1); + } + int bl = nbL, curr_min = bl; + + for (int i = 0; bl > 0; i++) { - std::cout << "========== Start iteration " << i << " ========\n"; - bl = landmark_perturbation(point_vector, L, chosen_landmarks); + std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; + bl=landmark_perturbation(point_vector, L, chosen_landmarks); + if (bl < curr_min) + curr_min=bl; } //end = clock(); diff --git a/src/Witness_complex/example/witness_complex_perturbations.cpp b/src/Witness_complex/example/witness_complex_perturbations.cpp index afb42807..514b115b 100644 --- a/src/Witness_complex/example/witness_complex_perturbations.cpp +++ b/src/Witness_complex/example/witness_complex_perturbations.cpp @@ -236,10 +236,10 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< } } //std::cout << "\n"; - /* + std::string out_file = "wl_result"; write_wl(out_file,WL); - */ + //******************** Constructng a witness complex std::cout << "Entered witness complex construction\n"; Witness_complex<> witnessComplex; -- cgit v1.2.3