diff options
57 files changed, 3386 insertions, 805 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 005df6d7..f949a803 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -123,7 +123,7 @@ else() include_directories(src/common/include/) include_directories(src/Alpha_complex/include/) include_directories(src/Bitmap_cubical_complex/include/) - include_directories(src/Bottleneck/include/) + include_directories(src/Bottleneck_distance/include/) include_directories(src/Contraction/include/) include_directories(src/Hasse_complex/include/) include_directories(src/Persistent_cohomology/include/) @@ -156,6 +156,9 @@ else() add_subdirectory(src/Tangential_complex/example) add_subdirectory(src/Tangential_complex/test) add_subdirectory(src/Tangential_complex/benchmark) + add_subdirectory(src/Bottleneck_distance/example) + add_subdirectory(src/Bottleneck_distance/test) + add_subdirectory(src/Bottleneck_distance/benchmark) # data points generator add_subdirectory(data/points/generator) diff --git a/biblio/how_to_cite_gudhi.bib b/biblio/how_to_cite_gudhi.bib index 03c05728..d88ecb75 100644 --- a/biblio/how_to_cite_gudhi.bib +++ b/biblio/how_to_cite_gudhi.bib @@ -33,7 +33,7 @@ , year = 2015 } -@incollection{gudhi:SkeletonBlocker +@incollection{gudhi:Skeleton-Blocker , author = "David Salinas" , title = "Skeleton-Blocker" , publisher = "{GUDHI Editorial Board}" @@ -78,20 +78,11 @@ , year = 2016 } -@incollection{gudhi:SpatialSearching -, author = "Cl\'ement Jamin" -, title = "Spatial searching" +@incollection{gudhi:BottleneckDistance +, author = "Fran{{\c{c}}ois Godi" +, title = "Bottleneck distance" , publisher = "{GUDHI Editorial Board}" , booktitle = "{GUDHI} User and Reference Manual" -, url = "http://gudhi.gforge.inria.fr/doc/latest/group__spatial__searching.html" +, url = "http://gudhi.gforge.inria.fr/doc/latest/group__bottleneck__distance.html" , year = 2016 } - -@incollection{gudhi:TangentialComplex -, author = "Cl\'ement Jamin" -, title = "Tangential complex" -, publisher = "{GUDHI Editorial Board}" -, booktitle = "{GUDHI} User and Reference Manual" -, url = "http://gudhi.gforge.inria.fr/doc/latest/group__tangential__complex.html" -, year = 2016 -}
\ No newline at end of file diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 9d5a9bad..1ff95c3d 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -193,17 +193,17 @@ class Alpha_complex { triangulation_ = new Delaunay_triangulation(point_dimension(*first)); std::vector<Point_d> point_cloud(first, last); - + // Creates a vector {0, 1, ..., N-1} std::vector<std::ptrdiff_t> indices(boost::counting_iterator<std::ptrdiff_t>(0), boost::counting_iterator<std::ptrdiff_t>(point_cloud.size())); - + typedef boost::iterator_property_map<typename std::vector<Point_d>::iterator, CGAL::Identity_property_map<std::ptrdiff_t>> Point_property_map; typedef CGAL::Spatial_sort_traits_adapter_d<Kernel, Point_property_map> Search_traits_d; - + CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud))); - + typename Delaunay_triangulation::Full_cell_handle hint; for (auto index : indices) { typename Delaunay_triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint); diff --git a/src/Bottleneck/concept/Persistence_diagram.h b/src/Bottleneck/concept/Persistence_diagram.h deleted file mode 100644 index eaaf8bc5..00000000 --- a/src/Bottleneck/concept/Persistence_diagram.h +++ /dev/null @@ -1,7 +0,0 @@ -typedef typename std::pair<double,double> Diagram_point; - -struct Persistence_Diagram -{ - const_iterator<Diagram_point> cbegin() const; - const_iterator<Diagram_point> cend() const; -}; diff --git a/src/Bottleneck/example/CMakeLists.txt b/src/Bottleneck/example/CMakeLists.txt deleted file mode 100644 index 77797202..00000000 --- a/src/Bottleneck/example/CMakeLists.txt +++ /dev/null @@ -1,5 +0,0 @@ -cmake_minimum_required(VERSION 2.6) -project(Bottleneck_examples) - -add_executable ( RandomDiagrams random_diagrams.cpp ) -add_test(RandomDiagrams ${CMAKE_CURRENT_BINARY_DIR}/RandomDiagrams) diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h deleted file mode 100644 index ea47e1d5..00000000 --- a/src/Bottleneck/include/gudhi/Graph_matching.h +++ /dev/null @@ -1,197 +0,0 @@ -/* 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): Francois Godi - * - * Copyright (C) 2015 INRIA Saclay (France) - * - * 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 SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ - -#include <deque> -#include <list> -#include <vector> - -#include "gudhi/Layered_neighbors_finder.h" - -namespace Gudhi { - -namespace bottleneck { - -template<typename Persistence_diagram1, typename Persistence_diagram2> -double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.); - -class Graph_matching { - public: - Graph_matching(const Persistence_diagrams_graph& g); - Graph_matching& operator=(const Graph_matching& m); - bool perfect() const; - bool multi_augment(); - void set_r(double r); - - private: - const Persistence_diagrams_graph& g; - double r; - std::vector<int> v_to_u; - std::list<int> unmatched_in_u; - - Layered_neighbors_finder* layering() const; - bool augment(Layered_neighbors_finder* layered_nf, int u_start_index, int max_depth); - void update(std::deque<int>& path); -}; - -Graph_matching::Graph_matching(const Persistence_diagrams_graph& g) - : g(g), r(0), v_to_u(g.size()), unmatched_in_u() { - for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index) - unmatched_in_u.emplace_back(u_point_index); -} - -Graph_matching& Graph_matching::operator=(const Graph_matching& m) { - r = m.r; - v_to_u = m.v_to_u; - unmatched_in_u = m.unmatched_in_u; - return *this; -} - -inline bool Graph_matching::perfect() const { - return unmatched_in_u.empty(); -} - -inline bool Graph_matching::multi_augment() { - if (perfect()) - return false; - Layered_neighbors_finder* layered_nf = layering(); - double rn = sqrt(g.size()); - int nblmax = layered_nf->vlayers_number()*2 + 1; - // verification of a necessary criterion - if ((unmatched_in_u.size() > rn && nblmax > rn) || nblmax == 0) - return false; - bool successful = false; - std::list<int>* tries = new std::list<int>(unmatched_in_u); - for (auto it = tries->cbegin(); it != tries->cend(); it++) - successful = successful || augment(layered_nf, *it, nblmax); - delete tries; - delete layered_nf; - return successful; -} - -inline void Graph_matching::set_r(double r) { - this->r = r; -} - -Layered_neighbors_finder* Graph_matching::layering() const { - bool end = false; - int layer = 0; - std::list<int> u_vertices(unmatched_in_u); - std::list<int> v_vertices; - Neighbors_finder nf(g, r); - Layered_neighbors_finder* layered_nf = new Layered_neighbors_finder(g, r); - for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index) - nf.add(v_point_index); - while (!u_vertices.empty()) { - for (auto it = u_vertices.cbegin(); it != u_vertices.cend(); ++it) { - std::list<int>* u_succ = nf.pull_all_near(*it); - for (auto it = u_succ->cbegin(); it != u_succ->cend(); ++it) { - layered_nf->add(*it, layer); - v_vertices.emplace_back(*it); - } - delete u_succ; - } - u_vertices.clear(); - for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) { - if (v_to_u.at(*it) == null_point_index()) - end = true; - else - u_vertices.emplace_back(v_to_u.at(*it)); - } - if (end) - return layered_nf; - v_vertices.clear(); - layer++; - } - return layered_nf; -} - -bool Graph_matching::augment(Layered_neighbors_finder *layered_nf, int u_start_index, int max_depth) { - std::deque<int> path; - path.emplace_back(u_start_index); - // start is a point from U - do { - if (static_cast<int>(path.size()) > max_depth) { - path.pop_back(); - path.pop_back(); - } - if (path.empty()) - return false; - int w = path.back(); - path.emplace_back(layered_nf->pull_near(w, path.size() / 2)); - while (path.back() == null_point_index()) { - path.pop_back(); - path.pop_back(); - if (path.empty()) - return false; - path.pop_back(); - path.emplace_back(layered_nf->pull_near(path.back(), path.size() / 2)); - } - path.emplace_back(v_to_u.at(path.back())); - } while (path.back() != null_point_index()); - path.pop_back(); - update(path); - return true; -} - -void Graph_matching::update(std::deque<int>& path) { - unmatched_in_u.remove(path.front()); - for (auto it = path.cbegin(); it != path.cend(); ++it) { - int tmp = *it; - ++it; - v_to_u[*it] = tmp; - } -} - -template<typename Persistence_diagram1, typename Persistence_diagram2> -double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e) { - Persistence_diagrams_graph g(diag1, diag2, e); - std::vector<double>* sd = g.sorted_distances(); - int idmin = 0; - int idmax = sd->size() - 1; - double alpha = pow(sd->size(), 0.25); - Graph_matching m(g); - Graph_matching biggest_unperfect = m; - while (idmin != idmax) { - int pas = static_cast<int>((idmax - idmin) / alpha); - m.set_r(sd->at(idmin + pas)); - while (m.multi_augment()) {} - if (m.perfect()) { - idmax = idmin + pas; - m = biggest_unperfect; - } else { - biggest_unperfect = m; - idmin = idmin + pas + 1; - } - } - double b = sd->at(idmin); - delete sd; - return b; -} - -} // namespace bottleneck - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ diff --git a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h deleted file mode 100644 index de36e00b..00000000 --- a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h +++ /dev/null @@ -1,74 +0,0 @@ -/* 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): Francois Godi - * - * Copyright (C) 2015 INRIA Saclay (France) - * - * 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 SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ - -#include <vector> - -#include "Neighbors_finder.h" - -// Layered_neighbors_finder is a data structure used to find if a query point from U has neighbors in V in a given -// vlayer of the vlayered persistence diagrams graph. V's points have to be added manually using their index. -// A neighbor returned is automatically removed. - -namespace Gudhi { - -namespace bottleneck { - -class Layered_neighbors_finder { - public: - Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r); - void add(int v_point_index, int vlayer); - int pull_near(int u_point_index, int vlayer); - int vlayers_number() const; - - private: - const Persistence_diagrams_graph& g; - const double r; - std::vector<Neighbors_finder> neighbors_finder; -}; - -Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r) : - g(g), r(r), neighbors_finder() { } - -inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) { - for (int l = neighbors_finder.size(); l <= vlayer; l++) - neighbors_finder.emplace_back(Neighbors_finder(g, r)); - neighbors_finder.at(vlayer).add(v_point_index); -} - -inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) { - if (static_cast<int> (neighbors_finder.size()) <= vlayer) - return null_point_index(); - return neighbors_finder.at(vlayer).pull_near(u_point_index); -} - -inline int Layered_neighbors_finder::vlayers_number() const { - return neighbors_finder.size(); -} - -} // namespace bottleneck - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck/include/gudhi/Neighbors_finder.h b/src/Bottleneck/include/gudhi/Neighbors_finder.h deleted file mode 100644 index 98256571..00000000 --- a/src/Bottleneck/include/gudhi/Neighbors_finder.h +++ /dev/null @@ -1,96 +0,0 @@ -/* 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): Francois Godi - * - * Copyright (C) 2015 INRIA Saclay (France) - * - * 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 SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ - -#include <unordered_set> -#include <list> - -#include "gudhi/Planar_neighbors_finder.h" - -namespace Gudhi { - -namespace bottleneck { - -// Neighbors_finder is a data structure used to find if a query point from U has neighbors in V -// in the persistence diagrams graph. -// V's points have to be added manually using their index. A neighbor returned is automatically removed. - -class Neighbors_finder { - public: - Neighbors_finder(const Persistence_diagrams_graph& g, double r); - void add(int v_point_index); - int pull_near(int u_point_index); - std::list<int>* pull_all_near(int u_point_index); - - private: - const Persistence_diagrams_graph& g; - const double r; - Planar_neighbors_finder planar_neighbors_f; - std::unordered_set<int> projections_f; -}; - -Neighbors_finder::Neighbors_finder(const Persistence_diagrams_graph& g, double r) : - g(g), r(r), planar_neighbors_f(g, r), projections_f() { } - -inline void Neighbors_finder::add(int v_point_index) { - if (g.on_the_v_diagonal(v_point_index)) - projections_f.emplace(v_point_index); - else - planar_neighbors_f.add(v_point_index); -} - -inline int Neighbors_finder::pull_near(int u_point_index) { - int v_challenger = g.corresponding_point_in_v(u_point_index); - if (planar_neighbors_f.contains(v_challenger) && g.distance(u_point_index, v_challenger) < r) { - planar_neighbors_f.remove(v_challenger); - return v_challenger; - } - if (g.on_the_u_diagonal(u_point_index)) { - auto it = projections_f.cbegin(); - if (it != projections_f.cend()) { - int tmp = *it; - projections_f.erase(it); - return tmp; - } - } else { - return planar_neighbors_f.pull_near(u_point_index); - } - return null_point_index(); -} - -inline std::list<int>* Neighbors_finder::pull_all_near(int u_point_index) { - std::list<int>* all_pull = planar_neighbors_f.pull_all_near(u_point_index); - int last_pull = pull_near(u_point_index); - while (last_pull != null_point_index()) { - all_pull->emplace_back(last_pull); - last_pull = pull_near(u_point_index); - } - return all_pull; -} - -} // namespace bottleneck - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h deleted file mode 100644 index 73ad940b..00000000 --- a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h +++ /dev/null @@ -1,147 +0,0 @@ -/* 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): Francois Godi - * - * Copyright (C) 2015 INRIA Saclay (France) - * - * 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 SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ - -#include <vector> -#include <set> -#include <cmath> -#include <utility> // for pair<> -#include <algorithm> // for max - -namespace Gudhi { - -namespace bottleneck { - -// Diagram_point is the type of the persistence diagram's points -typedef std::pair<double, double> Diagram_point; - -// Return the used index for encoding none of the points -int null_point_index(); - -// Persistence_diagrams_graph is the interface beetwen any external representation of the two persistence diagrams and -// the bottleneck distance computation. An interface is necessary to ensure basic functions complexity. - -class Persistence_diagrams_graph { - public: - // Persistence_diagram1 and 2 are the types of any externals representations of persistence diagrams. - // They have to have an iterator over points, which have to have fields first (for birth) and second (for death). - template<typename Persistence_diagram1, typename Persistence_diagram2> - Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.); - Persistence_diagrams_graph(); - bool on_the_u_diagonal(int u_point_index) const; - bool on_the_v_diagonal(int v_point_index) const; - int corresponding_point_in_u(int v_point_index) const; - int corresponding_point_in_v(int u_point_index) const; - double distance(int u_point_index, int v_point_index) const; - int size() const; - std::vector<double>* sorted_distances(); - - private: - std::vector<Diagram_point> u; - std::vector<Diagram_point> v; - Diagram_point get_u_point(int u_point_index) const; - Diagram_point get_v_point(int v_point_index) const; -}; - -inline int null_point_index() { - return -1; -} - -template<typename Persistence_diagram1, typename Persistence_diagram2> -Persistence_diagrams_graph::Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e) - : u(), v() { - for (auto it = diag1.cbegin(); it != diag1.cend(); ++it) - if (it->second - it->first > e) - u.emplace_back(*it); - for (auto it = diag2.cbegin(); it != diag2.cend(); ++it) - if (it->second - it->first > e) - v.emplace_back(*it); - if (u.size() < v.size()) - swap(u, v); -} - -Persistence_diagrams_graph::Persistence_diagrams_graph() - : u(), v() { } - -inline bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) const { - return u_point_index >= static_cast<int> (u.size()); -} - -inline bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) const { - return v_point_index >= static_cast<int> (v.size()); -} - -inline int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) const { - return on_the_v_diagonal(v_point_index) ? - v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size()); -} - -inline int Persistence_diagrams_graph::corresponding_point_in_v(int u_point_index) const { - return on_the_u_diagonal(u_point_index) ? - u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size()); -} - -inline double Persistence_diagrams_graph::distance(int u_point_index, int v_point_index) const { - // could be optimized for the case where one point is the projection of the other - if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index)) - return 0; - Diagram_point p_u = get_u_point(u_point_index); - Diagram_point p_v = get_v_point(v_point_index); - return (std::max)(std::fabs(p_u.first - p_v.first), std::fabs(p_u.second - p_v.second)); -} - -inline int Persistence_diagrams_graph::size() const { - return static_cast<int> (u.size() + v.size()); -} - -inline std::vector<double>* Persistence_diagrams_graph::sorted_distances() { - // could be optimized - std::set<double> sorted_distances; - for (int u_point_index = 0; u_point_index < size(); ++u_point_index) - for (int v_point_index = 0; v_point_index < size(); ++v_point_index) - sorted_distances.emplace(distance(u_point_index, v_point_index)); - return new std::vector<double>(sorted_distances.cbegin(), sorted_distances.cend()); -} - -inline Diagram_point Persistence_diagrams_graph::get_u_point(int u_point_index) const { - if (!on_the_u_diagonal(u_point_index)) - return u.at(u_point_index); - Diagram_point projector = v.at(corresponding_point_in_v(u_point_index)); - double x = (projector.first + projector.second) / 2; - return Diagram_point(x, x); -} - -inline Diagram_point Persistence_diagrams_graph::get_v_point(int v_point_index) const { - if (!on_the_v_diagonal(v_point_index)) - return v.at(v_point_index); - Diagram_point projector = u.at(corresponding_point_in_u(v_point_index)); - double x = (projector.first + projector.second) / 2; - return Diagram_point(x, x); -} - -} // namespace bottleneck - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ diff --git a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h deleted file mode 100644 index 4af672e4..00000000 --- a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h +++ /dev/null @@ -1,119 +0,0 @@ -/* 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): Francois Godi - * - * Copyright (C) 2015 INRIA Saclay (France) - * - * 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 SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ -#define SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ - -#include <list> -#include <iostream> -#include <set> - -#include "Persistence_diagrams_graph.h" - -namespace Gudhi { - -namespace bottleneck { - -// Planar_neighbors_finder is a data structure used to find if a query point from U has planar neighbors in V with the -// planar distance. -// V's points have to be added manually using their index. A neighbor returned is automatically removed but we can also -// remove points manually using their index. - -class Abstract_planar_neighbors_finder { - public: - Abstract_planar_neighbors_finder(const Persistence_diagrams_graph& g, double r); - virtual ~Abstract_planar_neighbors_finder() = 0; - virtual void add(int v_point_index) = 0; - virtual void remove(int v_point_index) = 0; - virtual bool contains(int v_point_index) const = 0; - virtual int pull_near(int u_point_index) = 0; - virtual std::list<int>* pull_all_near(int u_point_index); - - protected: - const Persistence_diagrams_graph& g; - const double r; -}; - - -// Naive_pnf is a nave implementation of Abstract_planar_neighbors_finder - -class Naive_pnf : public Abstract_planar_neighbors_finder { - public: - Naive_pnf(const Persistence_diagrams_graph& g, double r); - void add(int v_point_index); - void remove(int v_point_index); - bool contains(int v_point_index) const; - int pull_near(int u_point_index); - - private: - std::set<int> candidates; -}; - - -// Planar_neighbors_finder is the used Abstract_planar_neighbors_finder's implementation -typedef Naive_pnf Planar_neighbors_finder; - -Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(const Persistence_diagrams_graph& g, double r) : - g(g), r(r) { } - -inline Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() { } - -inline std::list<int>* Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) { - std::list<int>* all_pull = new std::list<int>(); - int last_pull = pull_near(u_point_index); - while (last_pull != null_point_index()) { - all_pull->emplace_back(last_pull); - last_pull = pull_near(u_point_index); - } - return all_pull; -} - -Naive_pnf::Naive_pnf(const Persistence_diagrams_graph& g, double r) : - Abstract_planar_neighbors_finder(g, r), candidates() { } - -inline void Naive_pnf::add(int v_point_index) { - candidates.emplace(v_point_index); -} - -inline void Naive_pnf::remove(int v_point_index) { - candidates.erase(v_point_index); -} - -inline bool Naive_pnf::contains(int v_point_index) const { - return (candidates.count(v_point_index) > 0); -} - -inline int Naive_pnf::pull_near(int u_point_index) { - for (auto it = candidates.begin(); it != candidates.end(); ++it) - if (g.distance(u_point_index, *it) <= r) { - int tmp = *it; - candidates.erase(it); - return tmp; - } - return null_point_index(); -} - -} // namespace bottleneck - -} // namespace Gudhi - -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck/test/CMakeLists.txt b/src/Bottleneck/test/CMakeLists.txt deleted file mode 100644 index 9d88ab25..00000000 --- a/src/Bottleneck/test/CMakeLists.txt +++ /dev/null @@ -1,21 +0,0 @@ -cmake_minimum_required(VERSION 2.6) -project(Bottleneck_tests) - -if (GCOVR_PATH) - # for gcovr to make coverage reports - Corbera Jenkins plugin - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") -endif() -if (GPROF_PATH) - # for gprof to make coverage reports - Jenkins - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") -endif() - -add_executable ( BottleneckUT bottleneck_unit_test.cpp ) -target_link_libraries(BottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) - -# Unitary tests -add_test(NAME BottleneckUT - COMMAND ${CMAKE_CURRENT_BINARY_DIR}/BottleneckUT - # XML format for Jenkins xUnit plugin - --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/BottleneckUT.xml --log_level=test_suite --report_level=no) - diff --git a/src/Bottleneck/test/bottleneck_unit_test.cpp b/src/Bottleneck/test/bottleneck_unit_test.cpp deleted file mode 100644 index c60f5d8a..00000000 --- a/src/Bottleneck/test/bottleneck_unit_test.cpp +++ /dev/null @@ -1,26 +0,0 @@ -#define BOOST_TEST_DYN_LINK -#define BOOST_TEST_MODULE "bottleneck" -#include <boost/test/unit_test.hpp> - -#include "gudhi/Graph_matching.h" -#include <iostream> - -using namespace Gudhi::bottleneck; - -BOOST_AUTO_TEST_CASE(random_diagrams) { - int n = 100; - // Random construction - std::vector< std::pair<double, double> > v1, v2; - for (int i = 0; i < n; i++) { - int a = rand() % n; - v1.emplace_back(a, a + rand() % (n - a)); - int b = rand() % n; - v2.emplace_back(b, b + rand() % (n - b)); - } - // v1 and v2 are persistence diagrams containing each 100 randoms points. - double b = bottleneck_distance(v1, v2, 0); - // - std::cout << b << std::endl; - const double EXPECTED_DISTANCE = 98.5; - BOOST_CHECK(b == EXPECTED_DISTANCE); -} diff --git a/src/Bottleneck_distance/benchmark/CMakeLists.txt b/src/Bottleneck_distance/benchmark/CMakeLists.txt new file mode 100644 index 00000000..f70dd8ff --- /dev/null +++ b/src/Bottleneck_distance/benchmark/CMakeLists.txt @@ -0,0 +1,13 @@ +cmake_minimum_required(VERSION 2.6) +project(Bottleneck_distance_benchmark) + + +# requires CGAL 4.8 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.8 ../../.. +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.8.0) + if (EIGEN3_FOUND) + add_executable ( bottleneck_chrono bottleneck_chrono.cpp ) + endif() + endif () +endif() diff --git a/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp b/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp new file mode 100644 index 00000000..456c570b --- /dev/null +++ b/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp @@ -0,0 +1,62 @@ +/* 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: Francois Godi + * + * Copyright (C) 2015 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/>. + */ + +#include <gudhi/Bottleneck.h> +#include <chrono> +#include <fstream> +#include <random> + +using namespace Gudhi::persistence_diagram; + + +double upper_bound = 400.; // any real > 0 + +int main() { + std::ofstream result_file; + result_file.open("results.csv", std::ios::out); + + for (int n = 1000; n <= 10000; n += 1000) { + std::uniform_real_distribution<double> unif1(0., upper_bound); + std::uniform_real_distribution<double> unif2(upper_bound / 1000., upper_bound / 100.); + std::default_random_engine re; + std::vector< std::pair<double, double> > v1, v2; + for (int i = 0; i < n; i++) { + double a = unif1(re); + double b = unif1(re); + double x = unif2(re); + double y = unif2(re); + v1.emplace_back(std::min(a, b), std::max(a, b)); + v2.emplace_back(std::min(a, b) + std::min(x, y), std::max(a, b) + std::max(x, y)); + if (i % 5 == 0) + v1.emplace_back(std::min(a, b), std::min(a, b) + x); + if (i % 3 == 0) + v2.emplace_back(std::max(a, b), std::max(a, b) + y); + } + std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); + double b = bottleneck_distance(v1, v2); + std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); + typedef std::chrono::duration<int, std::milli> millisecs_t; + millisecs_t duration(std::chrono::duration_cast<millisecs_t>(end - start)); + result_file << n << ";" << duration.count() << ";" << b << std::endl; + } + result_file.close(); +} diff --git a/src/Bottleneck_distance/concept/Persistence_diagram.h b/src/Bottleneck_distance/concept/Persistence_diagram.h new file mode 100644 index 00000000..2706716b --- /dev/null +++ b/src/Bottleneck_distance/concept/Persistence_diagram.h @@ -0,0 +1,49 @@ +/* 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: François Godi + * + * Copyright (C) 2015 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 CONCEPT_BOTTLENECK_DISTANCE_PERSISTENCE_DIAGRAM_H_ +#define CONCEPT_BOTTLENECK_DISTANCE_PERSISTENCE_DIAGRAM_H_ + +namespace Gudhi { + +namespace bottleneck_distance { + +/** \brief Concept of Diagram_point. std::get<0>(point) must return the birth of the corresponding component and std::get<1>(point) its death. + * A valid implementation of this concept is std::pair<double,double>. + * Death should be larger than birth, death can be std::numeric_limits<double>::infinity() for components which stay alive. + * + * \ingroup bottleneck_distance + */ +typename Diagram_point; + +/** \brief Concept of persistence diagram. It's a range of Diagram_point. + * std::begin(diagram) and std::end(diagram) must return corresponding iterators. + * + * \ingroup bottleneck_distance + */ +typename Persistence_Diagram; + +} // namespace bottleneck_distance + +} // namespace Gudhi + +#endif // CONCEPT_BOTTLENECK_DISTANCE_PERSISTENCE_DIAGRAM_H_ diff --git a/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h new file mode 100644 index 00000000..21187f9c --- /dev/null +++ b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h @@ -0,0 +1,51 @@ +/* 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: François Godi + * + * Copyright (C) 2015 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 DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_ +#define DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_ + +// needs namespace for Doxygen to link on classes +namespace Gudhi { +// needs namespace for Doxygen to link on classes +namespace bottleneck_distance { + +/** \defgroup bottleneck_distance Bottleneck distance + * + * \author François Godi + * @{ + * + * \section bottleneckdefinition Definition + * + * The bottleneck distance measures the similarity between two persistence diagrams. It's the shortest distance b for which there exists a perfect matching between + * the points of the two diagrams (completed with all the points on the diagonal in order to ignore cardinality mismatchs) such that + * any couple of matched points are at distance at most b. + * + * \image html perturb_pd.png On this picture, the red edges represent the matching. The bottleneck distance is the length of the longest edge. + * + */ +/** @} */ // end defgroup bottleneck_distance + +} // namespace bottleneck_distance + +} // namespace Gudhi + +#endif // DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_ diff --git a/src/Bottleneck_distance/doc/perturb_pd.png b/src/Bottleneck_distance/doc/perturb_pd.png Binary files differnew file mode 100644 index 00000000..be638de0 --- /dev/null +++ b/src/Bottleneck_distance/doc/perturb_pd.png diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt new file mode 100644 index 00000000..c66623e9 --- /dev/null +++ b/src/Bottleneck_distance/example/CMakeLists.txt @@ -0,0 +1,15 @@ +cmake_minimum_required(VERSION 2.6) +project(Bottleneck_distance_examples) + +# requires CGAL 4.8 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.8 ../../.. +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.8.0) + if (EIGEN3_FOUND) + add_executable (bottleneck_read_file_example bottleneck_read_file_example.cpp) + add_executable (bottleneck_basic_example bottleneck_basic_example.cpp) + + add_test(bottleneck_basic_example ${CMAKE_CURRENT_BINARY_DIR}/bottleneck_basic_example) + endif() + endif () +endif() diff --git a/src/Bottleneck/example/random_diagrams.cpp b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp index 71f152a6..91a7302f 100644 --- a/src/Bottleneck/example/random_diagrams.cpp +++ b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Francois Godi + * Authors: Francois Godi, small modifications by Pawel Dlotko * - * Copyright (C) 2015 INRIA Saclay (France) + * Copyright (C) 2015 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 @@ -20,21 +20,29 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include "gudhi/Graph_matching.h" +#include <gudhi/Bottleneck.h> #include <iostream> -using namespace Gudhi::bottleneck; - int main() { - int n = 100; + std::vector< std::pair<double, double> > v1, v2; - for (int i = 0; i < n; i++) { - int a = rand() % n; - v1.emplace_back(a, a + rand() % (n - a)); - int b = rand() % n; - v2.emplace_back(b, b + rand() % (n - b)); - } - // v1 and v2 are persistence diagrams containing each 100 randoms points. - double b = bottleneck_distance(v1, v2, 0); - std::cout << b << std::endl; + + v1.emplace_back(2.7, 3.7); + v1.emplace_back(9.6, 14.); + v1.emplace_back(34.2, 34.974); + v1.emplace_back(3., std::numeric_limits<double>::infinity()); + + v2.emplace_back(2.8, 4.45); + v2.emplace_back(9.5, 14.1); + v2.emplace_back(3.2, std::numeric_limits<double>::infinity()); + + + double b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2); + + std::cout << "Bottleneck distance = " << b << std::endl; + + b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1); + + std::cout << "Approx bottleneck distance = " << b << std::endl; + } diff --git a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp new file mode 100644 index 00000000..4c74b66e --- /dev/null +++ b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp @@ -0,0 +1,70 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Authors: Francois Godi, small modifications by Pawel Dlotko + * + * Copyright (C) 2015 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/>. + */ + +#define CGAL_HAS_THREADS + +#include <gudhi/Bottleneck.h> +#include <iostream> +#include <fstream> +#include <sstream> +#include <string> + +std::vector< std::pair<double, double> > read_diagram_from_file(const char* filename) { + std::ifstream in; + in.open(filename); + std::vector< std::pair<double, double> > result; + if (!in.is_open()) { + std::cerr << "File : " << filename << " do not exist. The program will now terminate \n"; + throw "File do not exist \n"; + } + + std::string line; + while (!in.eof()) { + getline(in, line); + if (line.length() != 0) { + std::stringstream lineSS; + lineSS << line; + double beginn, endd; + lineSS >> beginn; + lineSS >> endd; + result.push_back(std::make_pair(beginn, endd)); + } + } + in.close(); + return result; +} // read_diagram_from_file + +int main(int argc, char** argv) { + if (argc < 3) { + std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file " << + "should contain a birth-death pair per line. Third, optional parameter is an error bound on a bottleneck" << + " distance (set by default to zero). The program will now terminate \n"; + } + std::vector< std::pair< double, double > > diag1 = read_diagram_from_file(argv[1]); + std::vector< std::pair< double, double > > diag2 = read_diagram_from_file(argv[2]); + double tolerance = 0.; + if (argc == 4) { + tolerance = atof(argv[3]); + } + double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance); + std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl; +} diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h new file mode 100644 index 00000000..2b7e4767 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -0,0 +1,99 @@ +/* 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: Francois Godi + * + * Copyright (C) 2015 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 BOTTLENECK_H_ +#define BOTTLENECK_H_ + +#include <gudhi/Graph_matching.h> +#include <cmath> + +namespace Gudhi { + +namespace persistence_diagram { + +double bottleneck_distance_approx(Persistence_graph& g, double e) { + double b_lower_bound = 0.; + double b_upper_bound = g.diameter_bound(); + const double alpha = std::pow(g.size(), 1. / 5.); + Graph_matching m(g); + Graph_matching biggest_unperfect(g); + while (b_upper_bound - b_lower_bound > 2 * e) { + double step = b_lower_bound + (b_upper_bound - b_lower_bound) / alpha; + if (step <= b_lower_bound || step >= b_upper_bound) // Avoid precision problem + break; + m.set_r(step); + while (m.multi_augment()); // compute a maximum matching (in the graph corresponding to the current r) + if (m.perfect()) { + m = biggest_unperfect; + b_upper_bound = step; + } else { + biggest_unperfect = m; + b_lower_bound = step; + } + } + return (b_lower_bound + b_upper_bound) / 2.; +} + +double bottleneck_distance_exact(Persistence_graph& g) { + std::vector<double> sd = g.sorted_distances(); + long lower_bound_i = 0; + long upper_bound_i = sd.size() - 1; + const double alpha = std::pow(g.size(), 1. / 5.); + Graph_matching m(g); + Graph_matching biggest_unperfect(g); + while (lower_bound_i != upper_bound_i) { + long step = lower_bound_i + static_cast<long> ((upper_bound_i - lower_bound_i - 1) / alpha); + m.set_r(sd.at(step)); + while (m.multi_augment()); // compute a maximum matching (in the graph corresponding to the current r) + if (m.perfect()) { + m = biggest_unperfect; + upper_bound_i = step; + } else { + biggest_unperfect = m; + lower_bound_i = step + 1; + } + } + return sd.at(lower_bound_i); +} + +/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams (see concepts). + * If the last parameter e is not 0, you get an additive e-approximation, which is a lot faster to compute whatever is + * e. + * Thus, by default, e is a very small positive double, actually the smallest double possible such that the + * floating-point inaccuracies don't lead to a failure of the algorithm. + * + * \ingroup bottleneck_distance + */ +template<typename Persistence_diagram1, typename Persistence_diagram2> +double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, + double e = std::numeric_limits<double>::min()) { + Persistence_graph g(diag1, diag2, e); + if (g.bottleneck_alive() == std::numeric_limits<double>::infinity()) + return std::numeric_limits<double>::infinity(); + return std::max(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e)); +} + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // BOTTLENECK_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h new file mode 100644 index 00000000..253c89b4 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -0,0 +1,179 @@ +/* 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: Francois Godi + * + * Copyright (C) 2015 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 GRAPH_MATCHING_H_ +#define GRAPH_MATCHING_H_ + +#include <gudhi/Neighbors_finder.h> + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Structure representing a graph matching. The graph is a Persistence_diagrams_graph. + * + * \ingroup bottleneck_distance + */ +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. */ + bool multi_augment(); + /** \internal \brief Sets the maximum length of the edges allowed to be added in the matching, 0 initially. */ + void set_r(double r); + + private: + Persistence_graph& g; + 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; + + /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ + Layered_neighbors_finder layering() const; + /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */ + bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth); + /** \internal \brief Update the matching with the simple augmenting path given as parameter. */ + void update(std::vector<int> & path); +}; + +inline Graph_matching::Graph_matching(Persistence_graph& g) + : g(g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u() { + 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; +} + +inline bool Graph_matching::perfect() const { + return unmatched_in_u.empty(); +} + +inline bool Graph_matching::multi_augment() { + if (perfect()) + return false; + Layered_neighbors_finder layered_nf(layering()); + int max_depth = layered_nf.vlayers_number()*2 - 1; + double rn = sqrt(g.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); + 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; + return successful; +} + +inline void Graph_matching::set_r(double r) { + this->r = r; +} + +inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth) { + // V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices. + std::vector<int> path; + path.emplace_back(u_start_index); + do { + if (static_cast<int> (path.size()) > max_depth) { + path.pop_back(); + path.pop_back(); + } + if (path.empty()) + return false; + path.emplace_back(layered_nf.pull_near(path.back(), static_cast<int> (path.size()) / 2)); + while (path.back() == null_point_index()) { + path.pop_back(); + path.pop_back(); + if (path.empty()) + return false; + path.pop_back(); + path.emplace_back(layered_nf.pull_near(path.back(), path.size() / 2)); + } + path.emplace_back(v_to_u.at(path.back())); + } while (path.back() != null_point_index()); + // if v_to_u.at(path.back()) has no successor, path.back() is an exposed vertex + path.pop_back(); + update(path); + return true; +} + +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) + nf.add(v_point_index); + Layered_neighbors_finder layered_nf(g, 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) { + std::vector<int> u_succ(nf.pull_all_near(*it1)); + for (auto it2 = u_succ.begin(); it2 != u_succ.end(); ++it2) { + layered_nf.add(*it2, layer); + v_vertices.emplace_back(*it2); + } + } + // When the above for finishes, we have progress of one half-step (from U to V) in the BFS + u_vertices.clear(); + bool end = false; + for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) + if (v_to_u.at(*it) == null_point_index()) + // we stop when a nearest exposed V vertex (from U exposed vertices) has been found + end = true; + else + u_vertices.emplace_back(v_to_u.at(*it)); + // When the above for finishes, we have progress of one half-step (from V to U) in the BFS + if (end) + return layered_nf; + v_vertices.clear(); + } + return layered_nf; +} + +inline void Graph_matching::update(std::vector<int>& path) { + unmatched_in_u.remove(path.front()); + for (auto it = path.cbegin(); it != path.cend(); ++it) { + // Be careful, the iterator is incremented twice each time + int tmp = *it; + v_to_u[*(++it)] = tmp; + } +} + + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // GRAPH_MATCHING_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h new file mode 100644 index 00000000..0b2d26fe --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h @@ -0,0 +1,91 @@ +/* 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: Francois Godi + * + * Copyright (C) 2015 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 INTERNAL_POINT_H_ +#define INTERNAL_POINT_H_ + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Returns the used index for encoding none of the points */ +int null_point_index(); + +/** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */ +struct Internal_point { + double vec[2]; + int point_index; + + Internal_point() { } + + Internal_point(double x, double y, int p_i) { + vec[0] = x; + vec[1] = y; + point_index = p_i; + } + + double x() const { + return vec[ 0 ]; + } + + double y() const { + return vec[ 1 ]; + } + + double& x() { + return vec[ 0 ]; + } + + double& y() { + return vec[ 1 ]; + } + + bool operator==(const Internal_point& p) const { + return point_index == p.point_index; + } + + bool operator!=(const Internal_point& p) const { + return !(*this == p); + } +}; + +inline int null_point_index() { + return -1; +} + +struct Construct_coord_iterator { + typedef const double* result_type; + + const double* operator()(const Internal_point& p) const { + return p.vec; + } + + const double* operator()(const Internal_point& p, int) const { + return p.vec + 2; + } +}; + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // INTERNAL_POINT_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h new file mode 100644 index 00000000..96ece360 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -0,0 +1,171 @@ +/* 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: Francois Godi + * + * Copyright (C) 2015 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 NEIGHBORS_FINDER_H_ +#define NEIGHBORS_FINDER_H_ + +// 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> +#include <gudhi/Internal_point.h> + +#include <unordered_set> + +namespace Gudhi { + +namespace persistence_diagram { + +/** \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). + * + * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically + * removed. + * + * \ingroup bottleneck_distance + */ +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; + + public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Neighbors_finder(const Persistence_graph& g, double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index); + /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if + * there isn't such a point. */ + int pull_near(int u_point_index); + /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */ + std::vector<int> pull_all_near(int u_point_index); + + private: + const Persistence_graph& g; + const double r; + Kd_tree kd_t; + std::unordered_set<int> projections_f; +}; + +/** \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) in a layered graph layer given as parmeter. + * + * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically + * removed. + * + * \ingroup bottleneck_distance + */ +class Layered_neighbors_finder { + public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Layered_neighbors_finder(const Persistence_graph& g, double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index, int vlayer); + /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if + * there isn't such a point. */ + int pull_near(int u_point_index, int vlayer); + /** \internal \brief Returns the number of layers. */ + int vlayers_number() const; + + private: + const Persistence_graph& g; + const double r; + std::vector<std::unique_ptr<Neighbors_finder>> neighbors_finder; +}; + +inline Neighbors_finder::Neighbors_finder(const Persistence_graph& g, double r) : + g(g), r(r), kd_t(), projections_f() { } + +inline void Neighbors_finder::add(int v_point_index) { + if (g.on_the_v_diagonal(v_point_index)) + projections_f.emplace(v_point_index); + else + kd_t.insert(g.get_v_point(v_point_index)); +} + +inline int Neighbors_finder::pull_near(int u_point_index) { + int tmp; + int c = g.corresponding_point_in_v(u_point_index); + if (g.on_the_u_diagonal(u_point_index) && !projections_f.empty()) { + // Any pair of projection is at distance 0 + tmp = *projections_f.cbegin(); + projections_f.erase(tmp); + } else if (projections_f.count(c) && (g.distance(u_point_index, c) <= r)) { + // Is the query point near to its projection ? + tmp = c; + projections_f.erase(tmp); + } 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) + return null_point_index(); + tmp = it->first.point_index; + kd_t.remove(g.get_v_point(tmp)); + } + return tmp; +} + +inline std::vector<int> Neighbors_finder::pull_all_near(int u_point_index) { + std::vector<int> all_pull; + int last_pull = pull_near(u_point_index); + while (last_pull != null_point_index()) { + all_pull.emplace_back(last_pull); + last_pull = pull_near(u_point_index); + } + return all_pull; +} + +inline Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_graph& g, double r) : + g(g), r(r), neighbors_finder() { } + +inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) { + for (int l = neighbors_finder.size(); l <= vlayer; l++) + neighbors_finder.emplace_back(std::unique_ptr<Neighbors_finder>(new Neighbors_finder(g, r))); + neighbors_finder.at(vlayer)->add(v_point_index); +} + +inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) { + if (static_cast<int> (neighbors_finder.size()) <= vlayer) + return null_point_index(); + return neighbors_finder.at(vlayer)->pull_near(u_point_index); +} + +inline int Layered_neighbors_finder::vlayers_number() const { + return static_cast<int> (neighbors_finder.size()); +} + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h new file mode 100644 index 00000000..3a4a5fec --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h @@ -0,0 +1,176 @@ +/* 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: Francois Godi + * + * Copyright (C) 2015 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 PERSISTENCE_GRAPH_H_ +#define PERSISTENCE_GRAPH_H_ + +#include <vector> +#include <algorithm> +#include <gudhi/Internal_point.h> + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Structure representing an euclidean bipartite graph containing + * the points from the two persistence diagrams (including the projections). + * + * \ingroup bottleneck_distance + */ +class Persistence_graph { + public: + /** \internal \brief Constructor taking 2 Persistence_Diagrams (concept) as parameters. */ + template<typename Persistence_diagram1, typename Persistence_diagram2> + Persistence_graph(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e); + /** \internal \brief Is the given point from U the projection of a point in V ? */ + bool on_the_u_diagonal(int u_point_index) const; + /** \internal \brief Is the given point from V the projection of a point in U ? */ + bool on_the_v_diagonal(int v_point_index) const; + /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */ + int corresponding_point_in_u(int v_point_index) const; + /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */ + int corresponding_point_in_v(int u_point_index) const; + /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */ + double distance(int u_point_index, int v_point_index) const; + /** \internal \brief Returns size = |U| = |V|. */ + int size() const; + /** \internal \brief Is there as many infinite points (alive components) in both diagrams ? */ + double bottleneck_alive() const; + /** \internal \brief Returns the O(n^2) sorted distances between the points. */ + std::vector<double> sorted_distances() const; + /** \internal \brief Returns an upper bound for the diameter of the convex hull of all non infinite points */ + double diameter_bound() const; + /** \internal \brief Returns the corresponding internal point */ + Internal_point get_u_point(int u_point_index) const; + /** \internal \brief Returns the corresponding internal point */ + Internal_point get_v_point(int v_point_index) const; + + private: + std::vector<Internal_point> u; + std::vector<Internal_point> v; + double b_alive; +}; + +template<typename Persistence_diagram1, typename Persistence_diagram2> +Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1, + const Persistence_diagram2 &diag2, double e) + : u(), v(), b_alive(0.) { + std::vector<double> u_alive; + std::vector<double> v_alive; + for (auto it = std::begin(diag1); it != std::end(diag1); ++it) { + if (std::get<1>(*it) == std::numeric_limits<double>::infinity()) + u_alive.push_back(std::get<0>(*it)); + else if (std::get<1>(*it) - std::get<0>(*it) > e) + u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size())); + } + for (auto it = std::begin(diag2); it != std::end(diag2); ++it) { + if (std::get<1>(*it) == std::numeric_limits<double>::infinity()) + v_alive.push_back(std::get<0>(*it)); + else if (std::get<1>(*it) - std::get<0>(*it) > e) + v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size())); + } + if (u.size() < v.size()) + swap(u, v); + std::sort(u_alive.begin(), u_alive.end()); + std::sort(v_alive.begin(), v_alive.end()); + if (u_alive.size() != v_alive.size()) + b_alive = std::numeric_limits<double>::infinity(); + else for (auto it_u = u_alive.cbegin(), it_v = v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v) + b_alive = std::max(b_alive, std::fabs(*it_u - *it_v)); +} + +inline bool Persistence_graph::on_the_u_diagonal(int u_point_index) const { + return u_point_index >= static_cast<int> (u.size()); +} + +inline bool Persistence_graph::on_the_v_diagonal(int v_point_index) const { + return v_point_index >= static_cast<int> (v.size()); +} + +inline int Persistence_graph::corresponding_point_in_u(int v_point_index) const { + return on_the_v_diagonal(v_point_index) ? + v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size()); +} + +inline int Persistence_graph::corresponding_point_in_v(int u_point_index) const { + return on_the_u_diagonal(u_point_index) ? + u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size()); +} + +inline double Persistence_graph::distance(int u_point_index, int v_point_index) const { + if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index)) + return 0.; + Internal_point p_u = get_u_point(u_point_index); + Internal_point p_v = get_v_point(v_point_index); + return std::max(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y())); +} + +inline int Persistence_graph::size() const { + return static_cast<int> (u.size() + v.size()); +} + +inline double Persistence_graph::bottleneck_alive() const { + return b_alive; +} + +inline std::vector<double> Persistence_graph::sorted_distances() const { + std::vector<double> distances; + distances.push_back(0.); // for empty diagrams + for (int u_point_index = 0; u_point_index < size(); ++u_point_index) { + distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index))); + for (int v_point_index = 0; v_point_index < size(); ++v_point_index) + distances.push_back(distance(u_point_index, v_point_index)); + } + std::sort(distances.begin(), distances.end()); + return distances; +} + +inline Internal_point Persistence_graph::get_u_point(int u_point_index) const { + if (!on_the_u_diagonal(u_point_index)) + return u.at(u_point_index); + Internal_point projector = v.at(corresponding_point_in_v(u_point_index)); + double m = (projector.x() + projector.y()) / 2.; + return Internal_point(m, m, u_point_index); +} + +inline Internal_point Persistence_graph::get_v_point(int v_point_index) const { + if (!on_the_v_diagonal(v_point_index)) + return v.at(v_point_index); + Internal_point projector = u.at(corresponding_point_in_u(v_point_index)); + double m = (projector.x() + projector.y()) / 2.; + return Internal_point(m, m, v_point_index); +} + +inline double Persistence_graph::diameter_bound() const { + double max = 0.; + for (auto it = u.cbegin(); it != u.cend(); it++) + max = std::max(max, it->y()); + for (auto it = v.cbegin(); it != v.cend(); it++) + max = std::max(max, it->y()); + return max; +} + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // PERSISTENCE_GRAPH_H_ diff --git a/src/Bottleneck_distance/test/CMakeLists.txt b/src/Bottleneck_distance/test/CMakeLists.txt new file mode 100644 index 00000000..a6979d3c --- /dev/null +++ b/src/Bottleneck_distance/test/CMakeLists.txt @@ -0,0 +1,28 @@ +cmake_minimum_required(VERSION 2.6) +project(Bottleneck_distance_tests) + + +if (GCOVR_PATH) + # for gcovr to make coverage reports - Corbera Jenkins plugin + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") +endif() +if (GPROF_PATH) + # for gprof to make coverage reports - Jenkins + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") +endif() + +# requires CGAL 4.8 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.8 ../../.. +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.8.0) + if (EIGEN3_FOUND) + add_executable ( bottleneckUT bottleneck_unit_test.cpp ) + target_link_libraries(bottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + + # Unitary tests + add_test(NAME bottleneckUT COMMAND ${CMAKE_CURRENT_BINARY_DIR}/bottleneckUT + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/bottleneckUT.xml --log_level=test_suite --report_level=no) + endif() + endif () +endif() diff --git a/src/Bottleneck/test/README b/src/Bottleneck_distance/test/README index 0e7b8673..0e7b8673 100644 --- a/src/Bottleneck/test/README +++ b/src/Bottleneck_distance/test/README diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp new file mode 100644 index 00000000..e39613b3 --- /dev/null +++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp @@ -0,0 +1,167 @@ +/* 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: Francois Godi + * + * Copyright (C) 2015 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/>. + */ + + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "bottleneck distance" +#include <boost/test/unit_test.hpp> + +#include <random> +#include <gudhi/Bottleneck.h> + +using namespace Gudhi::persistence_diagram; + +int n1 = 81; // a natural number >0 +int n2 = 180; // a natural number >0 +double upper_bound = 406.43; // any real >0 + + +std::uniform_real_distribution<double> unif(0., upper_bound); +std::default_random_engine re; +std::vector< std::pair<double, double> > v1, v2; + +BOOST_AUTO_TEST_CASE(persistence_graph) { + // Random construction + for (int i = 0; i < n1; i++) { + double a = unif(re); + double b = unif(re); + v1.emplace_back(std::min(a, b), std::max(a, b)); + } + for (int i = 0; i < n2; i++) { + double a = unif(re); + double b = unif(re); + v2.emplace_back(std::min(a, b), std::max(a, b)); + } + Persistence_graph g(v1, v2, 0.); + std::vector<double> d(g.sorted_distances()); + // + BOOST_CHECK(!g.on_the_u_diagonal(n1 - 1)); + BOOST_CHECK(!g.on_the_u_diagonal(n1)); + BOOST_CHECK(!g.on_the_u_diagonal(n2 - 1)); + BOOST_CHECK(g.on_the_u_diagonal(n2)); + BOOST_CHECK(!g.on_the_v_diagonal(n1 - 1)); + BOOST_CHECK(g.on_the_v_diagonal(n1)); + BOOST_CHECK(g.on_the_v_diagonal(n2 - 1)); + BOOST_CHECK(g.on_the_v_diagonal(n2)); + // + BOOST_CHECK(g.corresponding_point_in_u(0) == n2); + BOOST_CHECK(g.corresponding_point_in_u(n1) == 0); + BOOST_CHECK(g.corresponding_point_in_v(0) == n1); + BOOST_CHECK(g.corresponding_point_in_v(n2) == 0); + // + BOOST_CHECK(g.size() == (n1 + n2)); + // + BOOST_CHECK((int) d.size() == (n1 + n2)*(n1 + n2) + n1 + n2 + 1); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, 0)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n1 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n2 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n2)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, (n1 + n2) - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, 0)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n1 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n2 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n2)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, (n1 + n2) - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, 0)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n1 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n2 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n2)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, (n1 + n2) - 1)) > 0); +} + +BOOST_AUTO_TEST_CASE(neighbors_finder) { + Persistence_graph g(v1, v2, 0.); + Neighbors_finder nf(g, 1.); + for (int v_point_index = 1; v_point_index < ((n2 + n1)*9 / 10); v_point_index += 2) + nf.add(v_point_index); + // + int v_point_index_1 = nf.pull_near(n2 / 2); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.)); + std::vector<int> l = nf.pull_all_near(n2 / 2); + bool v = true; + for (auto it = l.cbegin(); it != l.cend(); ++it) + v = v && (g.distance(n2 / 2, *it) > 1.); + BOOST_CHECK(v); + int v_point_index_2 = nf.pull_near(n2 / 2); + BOOST_CHECK(v_point_index_2 == -1); +} + +BOOST_AUTO_TEST_CASE(layered_neighbors_finder) { + Persistence_graph g(v1, v2, 0.); + Layered_neighbors_finder lnf(g, 1.); + for (int v_point_index = 1; v_point_index < ((n2 + n1)*9 / 10); v_point_index += 2) + lnf.add(v_point_index, v_point_index % 7); + // + int v_point_index_1 = lnf.pull_near(n2 / 2, 6); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.)); + int v_point_index_2 = lnf.pull_near(n2 / 2, 6); + BOOST_CHECK(v_point_index_2 == -1); + v_point_index_1 = lnf.pull_near(n2 / 2, 0); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.)); + v_point_index_2 = lnf.pull_near(n2 / 2, 0); + BOOST_CHECK(v_point_index_2 == -1); +} + +BOOST_AUTO_TEST_CASE(graph_matching) { + Persistence_graph g(v1, v2, 0.); + Graph_matching m1(g); + m1.set_r(0.); + int e = 0; + while (m1.multi_augment()) + ++e; + BOOST_CHECK(e > 0); + BOOST_CHECK(e <= 2 * sqrt(2 * (n1 + n2))); + Graph_matching m2 = m1; + BOOST_CHECK(!m2.multi_augment()); + m2.set_r(upper_bound); + e = 0; + while (m2.multi_augment()) + ++e; + BOOST_CHECK(e <= 2 * sqrt(2 * (n1 + n2))); + BOOST_CHECK(m2.perfect()); + BOOST_CHECK(!m1.perfect()); +} + +BOOST_AUTO_TEST_CASE(global) { + std::uniform_real_distribution<double> unif1(0., upper_bound); + std::uniform_real_distribution<double> unif2(upper_bound / 10000., upper_bound / 100.); + std::default_random_engine re; + std::vector< std::pair<double, double> > v1, v2; + for (int i = 0; i < n1; i++) { + double a = unif1(re); + double b = unif1(re); + double x = unif2(re); + double y = unif2(re); + v1.emplace_back(std::min(a, b), std::max(a, b)); + v2.emplace_back(std::min(a, b) + std::min(x, y), std::max(a, b) + std::max(x, y)); + if (i % 5 == 0) + v1.emplace_back(std::min(a, b), std::min(a, b) + x); + if (i % 3 == 0) + v2.emplace_back(std::max(a, b), std::max(a, b) + y); + } + BOOST_CHECK(bottleneck_distance(v1, v2, 0.) <= upper_bound / 100.); + BOOST_CHECK(bottleneck_distance(v1, v2, upper_bound / 10000.) <= upper_bound / 100. + upper_bound / 10000.); + BOOST_CHECK(std::abs(bottleneck_distance(v1, v2, 0.) - bottleneck_distance(v1, v2, upper_bound / 10000.)) <= upper_bound / 10000.); +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e26b2d25..1fa62e35 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -115,6 +115,7 @@ else() add_subdirectory(example/Spatial_searching) add_subdirectory(example/Subsampling) add_subdirectory(example/Tangential_complex) + add_subdirectory(example/Bottleneck_distance) # data points generator add_subdirectory(data/points/generator) diff --git a/src/Contraction/example/Garland_heckbert.cpp b/src/Contraction/example/Garland_heckbert.cpp index 5347830c..4689519f 100644 --- a/src/Contraction/example/Garland_heckbert.cpp +++ b/src/Contraction/example/Garland_heckbert.cpp @@ -63,12 +63,10 @@ typedef Skeleton_blocker_contractor<Complex> Complex_contractor; * the point minimizing the cost of the quadric. */ class GH_placement : public Gudhi::contraction::Placement_policy<EdgeProfile> { - Complex& complex_; - public: typedef Gudhi::contraction::Placement_policy<EdgeProfile>::Placement_type Placement_type; - GH_placement(Complex& complex) : complex_(complex) { } + GH_placement(Complex& complex) { } Placement_type operator()(const EdgeProfile& profile) const override { auto sum_quad(profile.v0().quadric); @@ -87,12 +85,10 @@ class GH_placement : public Gudhi::contraction::Placement_policy<EdgeProfile> { * which expresses a squared distances with triangles planes. */ class GH_cost : public Gudhi::contraction::Cost_policy<EdgeProfile> { - Complex& complex_; - public: typedef Gudhi::contraction::Cost_policy<EdgeProfile>::Cost_type Cost_type; - GH_cost(Complex& complex) : complex_(complex) { } + GH_cost(Complex& complex) { } Cost_type operator()(EdgeProfile const& profile, boost::optional<Point> const& new_point) const override { Cost_type res; @@ -111,10 +107,8 @@ class GH_cost : public Gudhi::contraction::Cost_policy<EdgeProfile> { * and we update them when contracting an edge (the quadric become the sum of both quadrics). */ class GH_visitor : public Gudhi::contraction::Contraction_visitor<EdgeProfile> { - Complex& complex_; - public: - GH_visitor(Complex& complex) : complex_(complex) { } + GH_visitor(Complex& complex) { } // Compute quadrics for every vertex v // The quadric of v consists in the sum of quadric diff --git a/src/Doxyfile b/src/Doxyfile index 943869ad..71aa038d 100644 --- a/src/Doxyfile +++ b/src/Doxyfile @@ -848,7 +848,8 @@ IMAGE_PATH = doc/Skeleton_blocker/ \ doc/Bitmap_cubical_complex/ \ doc/Subsampling/ \ doc/Spatial_searching/ \ - doc/Tangential_complex/ + doc/Tangential_complex/ \ + doc/Bottleneck_distance/ # The INPUT_FILTER tag can be used to specify a program that doxygen should # invoke to filter for each input file. Doxygen will invoke the filter program diff --git a/src/GudhUI/gui/MainWindow.h b/src/GudhUI/gui/MainWindow.h index c8c3fcf6..3fe0d720 100644 --- a/src/GudhUI/gui/MainWindow.h +++ b/src/GudhUI/gui/MainWindow.h @@ -26,6 +26,9 @@ // Workaround for moc-qt4 not parsing boost headers #include <CGAL/config.h> +// Workaround https://svn.boost.org/trac/boost/ticket/12534 +#include <boost/container/flat_map.hpp> + #include <QMainWindow> #include "ui_main_window.h" #include "model/Model.h" diff --git a/src/GudhUI/utils/Critical_points.h b/src/GudhUI/utils/Critical_points.h index 3021a5fe..b88293e9 100644 --- a/src/GudhUI/utils/Critical_points.h +++ b/src/GudhUI/utils/Critical_points.h @@ -105,8 +105,6 @@ template<typename SkBlComplex> class Critical_points { if (link.empty()) return 0; - Edge_contractor<Complex> contractor(link, link.num_vertices() - 1); - if (link.num_connected_components() > 1) // one than more CC -> not contractible return 0; diff --git a/src/GudhUI/utils/Is_manifold.h b/src/GudhUI/utils/Is_manifold.h index 0640ea47..6dd7898e 100644 --- a/src/GudhUI/utils/Is_manifold.h +++ b/src/GudhUI/utils/Is_manifold.h @@ -76,7 +76,6 @@ template<typename SkBlComplex> class Is_manifold { bool is_k_sphere(Vertex_handle v, int k) { auto link = input_complex_.link(v); - Edge_contractor<Complex> contractor(link, link.num_vertices() - 1); return (is_sphere_simplex(link) == k); } diff --git a/src/GudhUI/utils/Vertex_collapsor.h b/src/GudhUI/utils/Vertex_collapsor.h index 2b36cb3a..3f0e7ffd 100644 --- a/src/GudhUI/utils/Vertex_collapsor.h +++ b/src/GudhUI/utils/Vertex_collapsor.h @@ -80,7 +80,6 @@ template<typename SkBlComplex> class Vertex_collapsor { if (link.empty()) return false; if (link.is_cone()) return true; if (link.num_connected_components() > 1) return false; - Edge_contractor<Complex> contractor(link, link.num_vertices() - 1); return (link.num_vertices() == 1); } }; diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index b31df6a4..b3339b7d 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -110,7 +110,7 @@ class Persistent_cohomology { cell_pool_() { if (cpx_->num_simplices() > std::numeric_limits<Simplex_key>::max()) { // num_simplices must be strictly lower than the limit, because a value is reserved for null_key. - throw std::out_of_range ("The number of simplices is more than Simplex_key type numeric limit."); + throw std::out_of_range("The number of simplices is more than Simplex_key type numeric limit."); } Simplex_key idx_fil = 0; for (auto sh : cpx_->filtration_simplex_range()) { @@ -300,8 +300,7 @@ class Persistent_cohomology { // with multiplicity. We used to sum the coefficients directly in // annotations_in_boundary by using a map, we now do it later. typedef std::pair<Column *, int> annotation_t; - // Danger: not thread-safe! - static std::vector<annotation_t> annotations_in_boundary; + thread_local std::vector<annotation_t> annotations_in_boundary; annotations_in_boundary.clear(); int sign = 1 - 2 * (dim_sigma % 2); // \in {-1,1} provides the sign in the // alternate sum in the boundary. @@ -690,6 +689,22 @@ class Persistent_cohomology { return persistent_pairs_; } + /** @brief Returns persistence intervals for a given dimension. + * @param[in] dimension Dimension to get the birth and death pairs from. + * @return A vector of persistence intervals (birth and death) on a fixed dimension. + */ + std::vector< std::pair< Filtration_value , Filtration_value > > + intervals_in_dimension(int dimension) { + std::vector< std::pair< Filtration_value , Filtration_value > > result; + // auto && pair, to avoid unnecessary copying + for (auto && pair : persistent_pairs_) { + if (cpx_->dimension(get<0>(pair)) == dimension) { + result.emplace_back(cpx_->filtration(get<0>(pair)), cpx_->filtration(get<1>(pair))); + } + } + return result; + } + private: /* * Structure representing a cocycle. diff --git a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp index 40221005..0ed3fddf 100644 --- a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp +++ b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp @@ -84,6 +84,8 @@ BOOST_AUTO_TEST_CASE( plain_homology_betti_numbers ) // 2 1 0 inf // means that in Z/2Z-homology, the Betti numbers are b0=2 and b1=1. + std::cout << "BETTI NUMBERS" << std::endl; + BOOST_CHECK(pcoh.betti_number(0) == 2); BOOST_CHECK(pcoh.betti_number(1) == 1); BOOST_CHECK(pcoh.betti_number(2) == 0); @@ -93,6 +95,8 @@ BOOST_AUTO_TEST_CASE( plain_homology_betti_numbers ) BOOST_CHECK(bns[0] == 2); BOOST_CHECK(bns[1] == 1); BOOST_CHECK(bns[2] == 0); + + std::cout << "GET PERSISTENT PAIRS" << std::endl; // Custom sort and output persistence cmp_intervals_by_dim_then_length<Mini_simplex_tree> cmp(&st); @@ -115,6 +119,33 @@ BOOST_AUTO_TEST_CASE( plain_homology_betti_numbers ) BOOST_CHECK(st.dimension(get<0>(persistent_pairs[2])) == 0); BOOST_CHECK(st.filtration(get<0>(persistent_pairs[2])) == 0); BOOST_CHECK(get<1>(persistent_pairs[2]) == st.null_simplex()); + + std::cout << "INTERVALS IN DIMENSION" << std::endl; + + auto intervals_in_dimension_0 = pcoh.intervals_in_dimension(0); + std::cout << "intervals_in_dimension_0.size() = " << intervals_in_dimension_0.size() << std::endl; + for (std::size_t i = 0; i < intervals_in_dimension_0.size(); i++) + std::cout << "intervals_in_dimension_0[" << i << "] = [" << intervals_in_dimension_0[i].first << "," << + intervals_in_dimension_0[i].second << "]" << std::endl; + BOOST_CHECK(intervals_in_dimension_0.size() == 2); + BOOST_CHECK(intervals_in_dimension_0[0].first == 0); + BOOST_CHECK(intervals_in_dimension_0[0].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + BOOST_CHECK(intervals_in_dimension_0[1].first == 0); + BOOST_CHECK(intervals_in_dimension_0[1].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + + + auto intervals_in_dimension_1 = pcoh.intervals_in_dimension(1); + std::cout << "intervals_in_dimension_1.size() = " << intervals_in_dimension_1.size() << std::endl; + for (std::size_t i = 0; i < intervals_in_dimension_1.size(); i++) + std::cout << "intervals_in_dimension_1[" << i << "] = [" << intervals_in_dimension_1[i].first << "," << + intervals_in_dimension_1[i].second << "]" << std::endl; + BOOST_CHECK(intervals_in_dimension_1.size() == 1); + BOOST_CHECK(intervals_in_dimension_1[0].first == 0); + BOOST_CHECK(intervals_in_dimension_1[0].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + + auto intervals_in_dimension_2 = pcoh.intervals_in_dimension(2); + std::cout << "intervals_in_dimension_2.size() = " << intervals_in_dimension_2.size() << std::endl; + BOOST_CHECK(intervals_in_dimension_2.size() == 0); } using Simplex_tree = Gudhi::Simplex_tree<>; @@ -231,4 +262,30 @@ BOOST_AUTO_TEST_CASE( betti_numbers ) BOOST_CHECK(st.dimension(get<0>(persistent_pairs[2])) == 0); BOOST_CHECK(st.filtration(get<0>(persistent_pairs[2])) == 1); BOOST_CHECK(get<1>(persistent_pairs[2]) == st.null_simplex()); + + std::cout << "INTERVALS IN DIMENSION" << std::endl; + + auto intervals_in_dimension_0 = pcoh.intervals_in_dimension(0); + std::cout << "intervals_in_dimension_0.size() = " << intervals_in_dimension_0.size() << std::endl; + for (std::size_t i = 0; i < intervals_in_dimension_0.size(); i++) + std::cout << "intervals_in_dimension_0[" << i << "] = [" << intervals_in_dimension_0[i].first << "," << + intervals_in_dimension_0[i].second << "]" << std::endl; + BOOST_CHECK(intervals_in_dimension_0.size() == 2); + BOOST_CHECK(intervals_in_dimension_0[0].first == 2); + BOOST_CHECK(intervals_in_dimension_0[0].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + BOOST_CHECK(intervals_in_dimension_0[1].first == 1); + BOOST_CHECK(intervals_in_dimension_0[1].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + + auto intervals_in_dimension_1 = pcoh.intervals_in_dimension(1); + std::cout << "intervals_in_dimension_1.size() = " << intervals_in_dimension_1.size() << std::endl; + for (std::size_t i = 0; i < intervals_in_dimension_1.size(); i++) + std::cout << "intervals_in_dimension_1[" << i << "] = [" << intervals_in_dimension_1[i].first << "," << + intervals_in_dimension_1[i].second << "]" << std::endl; + BOOST_CHECK(intervals_in_dimension_1.size() == 1); + BOOST_CHECK(intervals_in_dimension_1[0].first == 4); + BOOST_CHECK(intervals_in_dimension_1[0].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + + auto intervals_in_dimension_2 = pcoh.intervals_in_dimension(2); + std::cout << "intervals_in_dimension_2.size() = " << intervals_in_dimension_2.size() << std::endl; + BOOST_CHECK(intervals_in_dimension_2.size() == 0); } diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 63e3f0e5..317bce23 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1029,7 +1029,7 @@ class Simplex_tree { Dictionary_it next = siblings->members().begin(); ++next; - static std::vector<std::pair<Vertex_handle, Node> > inter; // static, not thread-safe. + thread_local std::vector<std::pair<Vertex_handle, Node> > inter; for (Dictionary_it s_h = siblings->members().begin(); s_h != siblings->members().end(); ++s_h, ++next) { Simplex_handle root_sh = find_vertex(s_h->first); diff --git a/src/Spatial_searching/example/CMakeLists.txt b/src/Spatial_searching/example/CMakeLists.txt index e73b201c..6238a0ec 100644 --- a/src/Spatial_searching/example/CMakeLists.txt +++ b/src/Spatial_searching/example/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 2.6) project(Spatial_searching_examples) if(CGAL_FOUND) - if (NOT CGAL_VERSION VERSION_LESS 4.9.0) + if (NOT CGAL_VERSION VERSION_LESS 4.8.1) if (EIGEN3_FOUND) add_executable( Spatial_searching_example_spatial_searching example_spatial_searching.cpp ) target_link_libraries(Spatial_searching_example_spatial_searching ${CGAL_LIBRARY}) diff --git a/src/Spatial_searching/test/CMakeLists.txt b/src/Spatial_searching/test/CMakeLists.txt index 7f443b79..2c685c72 100644 --- a/src/Spatial_searching/test/CMakeLists.txt +++ b/src/Spatial_searching/test/CMakeLists.txt @@ -11,7 +11,7 @@ if (GPROF_PATH) endif() if(CGAL_FOUND) - if (NOT CGAL_VERSION VERSION_LESS 4.9.0) + if (NOT CGAL_VERSION VERSION_LESS 4.8.1) if (EIGEN3_FOUND) add_executable( Spatial_searching_test_Kd_tree_search test_Kd_tree_search.cpp ) target_link_libraries(Spatial_searching_test_Kd_tree_search diff --git a/src/Subsampling/example/CMakeLists.txt b/src/Subsampling/example/CMakeLists.txt index 54349f0c..0fd3335c 100644 --- a/src/Subsampling/example/CMakeLists.txt +++ b/src/Subsampling/example/CMakeLists.txt @@ -6,6 +6,7 @@ if(CGAL_FOUND) if (EIGEN3_FOUND) add_executable(Subsampling_example_pick_n_random_points example_pick_n_random_points.cpp) add_executable(Subsampling_example_choose_n_farthest_points example_choose_n_farthest_points.cpp) + add_executable(Subsampling_example_custom_kernel example_custom_kernel.cpp) add_executable(Subsampling_example_sparsify_point_set example_sparsify_point_set.cpp) target_link_libraries(Subsampling_example_sparsify_point_set ${CGAL_LIBRARY}) diff --git a/src/Subsampling/example/example_custom_kernel.cpp b/src/Subsampling/example/example_custom_kernel.cpp new file mode 100644 index 00000000..25b5bf6c --- /dev/null +++ b/src/Subsampling/example/example_custom_kernel.cpp @@ -0,0 +1,63 @@ +#include <gudhi/choose_n_farthest_points.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Random.h> + +#include <vector> +#include <iterator> + + +/* The class Kernel contains a distance function defined on the set of points {0, 1, 2, 3} + * and computes a distance according to the matrix: + * 0 1 2 4 + * 1 0 4 2 + * 2 4 0 1 + * 4 2 1 0 + */ +class Kernel { + public: + typedef double FT; + typedef unsigned Point_d; + + // Class Squared_distance_d + class Squared_distance_d { + private: + std::vector<std::vector<FT>> matrix_; + + public: + Squared_distance_d() { + matrix_.push_back(std::vector<FT>({0, 1, 2, 4})); + matrix_.push_back(std::vector<FT>({1, 0, 4, 2})); + matrix_.push_back(std::vector<FT>({2, 4, 0, 1})); + matrix_.push_back(std::vector<FT>({4, 2, 1, 0})); + } + + FT operator()(Point_d p1, Point_d p2) { + return matrix_[p1][p2]; + } + }; + + // Constructor + Kernel() {} + + // Object of type Squared_distance_d + Squared_distance_d squared_distance_d_object() const { + return Squared_distance_d(); + } +}; + +int main(void) { + typedef Kernel K; + typedef typename K::Point_d Point_d; + + K k; + 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)); + 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"; + + return 0; +} diff --git a/src/Subsampling/include/gudhi/choose_n_farthest_points.h b/src/Subsampling/include/gudhi/choose_n_farthest_points.h index 9b45c640..5e908090 100644 --- a/src/Subsampling/include/gudhi/choose_n_farthest_points.h +++ b/src/Subsampling/include/gudhi/choose_n_farthest_points.h @@ -48,15 +48,28 @@ namespace 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`. + * \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[in] starting_point The seed in the farthest point algorithm. + * @param[out] output_it The output iterator. * */ template < typename Kernel, -typename Point_container, +typename Point_range, typename OutputIterator> void choose_n_farthest_points(Kernel const &k, - Point_container const &input_pts, + Point_range const &input_pts, std::size_t final_size, std::size_t starting_point, OutputIterator output_it) { @@ -101,15 +114,27 @@ void choose_n_farthest_points(Kernel const &k, * \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_container, +typename Point_range, typename OutputIterator> void choose_n_farthest_points(Kernel const& k, - Point_container const &input_pts, + Point_range const &input_pts, unsigned final_size, OutputIterator output_it) { // Tests to the limit diff --git a/src/Subsampling/include/gudhi/pick_n_random_points.h b/src/Subsampling/include/gudhi/pick_n_random_points.h index e89b2b2d..f0e3f1f1 100644 --- a/src/Subsampling/include/gudhi/pick_n_random_points.h +++ b/src/Subsampling/include/gudhi/pick_n_random_points.h @@ -57,7 +57,9 @@ void pick_n_random_points(Point_container const &points, #endif std::size_t nbP = boost::size(points); - assert(nbP >= final_size); + if (final_size > nbP) + final_size = nbP; + std::vector<int> landmarks(nbP); std::iota(landmarks.begin(), landmarks.end(), 0); diff --git a/src/Subsampling/include/gudhi/sparsify_point_set.h b/src/Subsampling/include/gudhi/sparsify_point_set.h index 7ff11b4c..507f8c79 100644 --- a/src/Subsampling/include/gudhi/sparsify_point_set.h +++ b/src/Subsampling/include/gudhi/sparsify_point_set.h @@ -64,8 +64,6 @@ sparsify_point_set( typedef typename Gudhi::spatial_searching::Kd_tree_search< Kernel, Point_range> Points_ds; - typename Kernel::Squared_distance_d sqdist = k.squared_distance_d_object(); - #ifdef GUDHI_SUBSAMPLING_PROFILING Gudhi::Clock t; #endif diff --git a/src/Tangential_complex/benchmark/CMakeLists.txt b/src/Tangential_complex/benchmark/CMakeLists.txt index 56dd8128..788c2b4d 100644 --- a/src/Tangential_complex/benchmark/CMakeLists.txt +++ b/src/Tangential_complex/benchmark/CMakeLists.txt @@ -1,15 +1,6 @@ cmake_minimum_required(VERSION 2.6) project(Tangential_complex_benchmark) -if (GCOVR_PATH) - # for gcovr to make coverage reports - Corbera Jenkins plugin - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") -endif() -if (GPROF_PATH) - # for gprof to make coverage reports - Jenkins - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") -endif() - # need CGAL 4.8 if(CGAL_FOUND) if (NOT CGAL_VERSION VERSION_LESS 4.8.0) diff --git a/src/Tangential_complex/benchmark/benchmark_tc.cpp b/src/Tangential_complex/benchmark/benchmark_tc.cpp index 943fcb54..6d6dd548 100644 --- a/src/Tangential_complex/benchmark/benchmark_tc.cpp +++ b/src/Tangential_complex/benchmark/benchmark_tc.cpp @@ -161,7 +161,7 @@ typename Kernel, typename OutputIteratorPoints> bool load_points_from_file( const std::string &filename, OutputIteratorPoints points, - std::size_t only_first_n_points = std::numeric_limits<std::size_t>::max()) { + std::size_t only_first_n_points = (std::numeric_limits<std::size_t>::max)()) { typedef typename Kernel::Point_d Point; std::ifstream in(filename); @@ -196,7 +196,7 @@ bool load_points_and_tangent_space_basis_from_file( OutputIteratorPoints points, OutputIteratorTS tangent_spaces, int intrinsic_dim, - std::size_t only_first_n_points = std::numeric_limits<std::size_t>::max()) { + std::size_t only_first_n_points = (std::numeric_limits<std::size_t>::max)()) { typedef typename Kernel::Point_d Point; typedef typename Kernel::Vector_d Vector; diff --git a/src/Tangential_complex/include/gudhi/Tangential_complex.h b/src/Tangential_complex/include/gudhi/Tangential_complex.h index e748d3b7..cfc82eb1 100644 --- a/src/Tangential_complex/include/gudhi/Tangential_complex.h +++ b/src/Tangential_complex/include/gudhi/Tangential_complex.h @@ -63,6 +63,7 @@ #include <iterator> #include <cmath> // for std::sqrt #include <string> +#include <cstddef> // for std::size_t #ifdef GUDHI_USE_TBB #include <tbb/parallel_for.h> @@ -82,7 +83,7 @@ using namespace internal; class Vertex_data { public: - Vertex_data(std::size_t data = std::numeric_limits<std::size_t>::max()) + Vertex_data(std::size_t data = (std::numeric_limits<std::size_t>::max)()) : m_data(data) { } operator std::size_t() { @@ -1048,7 +1049,7 @@ class Tangential_complex { #endif // GUDHI_USE_TBB bool is_infinite(Simplex const& s) const { - return *s.rbegin() == std::numeric_limits<std::size_t>::max(); + return *s.rbegin() == (std::numeric_limits<std::size_t>::max)(); } // Output: "triangulation" is a Regular Triangulation containing at least the @@ -1130,7 +1131,7 @@ class Tangential_complex { Tr_vertex_handle vh = triangulation.insert_if_in_star(proj_pt, center_vertex); // Tr_vertex_handle vh = triangulation.insert(proj_pt); - if (vh != Tr_vertex_handle()) { + if (vh != Tr_vertex_handle() && vh->data() == (std::numeric_limits<std::size_t>::max)()) { #ifdef GUDHI_TC_VERY_VERBOSE ++num_inserted_points; #endif @@ -1292,6 +1293,8 @@ class Tangential_complex { if (index != i) incident_simplex.insert(index); } + GUDHI_CHECK(incident_simplex.size() == cur_dim_plus_1 - 1, + std::logic_error("update_star: wrong size of incident simplex")); star.push_back(incident_simplex); } } @@ -1303,21 +1306,14 @@ class Tangential_complex { , bool normalize_basis = true , Orthogonal_space_basis *p_orth_space_basis = NULL ) { - unsigned int num_pts_for_pca = static_cast<unsigned int> (std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)); + unsigned int num_pts_for_pca = (std::min)(static_cast<unsigned int> (std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)), + static_cast<unsigned int> (m_points.size())); // Kernel functors typename K::Construct_vector_d constr_vec = m_k.construct_vector_d_object(); typename K::Compute_coordinate_d coord = m_k.compute_coordinate_d_object(); - typename K::Squared_length_d sqlen = - m_k.squared_length_d_object(); - typename K::Scaled_vector_d scaled_vec = - m_k.scaled_vector_d_object(); - typename K::Scalar_product_d scalar_pdct = - m_k.scalar_product_d_object(); - typename K::Difference_of_vectors_d diff_vec = - m_k.difference_of_vectors_d_object(); #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM KNS_range kns_range = m_points_ds_for_tse.query_k_nearest_neighbors( @@ -1392,7 +1388,8 @@ class Tangential_complex { // on it. Note that most points are duplicated. Tangent_space_basis compute_tangent_space(const Simplex &s, bool normalize_basis = true) { - unsigned int num_pts_for_pca = static_cast<unsigned int> (std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)); + unsigned int num_pts_for_pca = (std::min)(static_cast<unsigned int> (std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)), + static_cast<unsigned int> (m_points.size())); // Kernel functors typename K::Construct_vector_d constr_vec = @@ -1650,7 +1647,7 @@ class Tangential_complex { for (; it_point_idx != simplex.end(); ++it_point_idx) { std::size_t point_idx = *it_point_idx; // Don't check infinite simplices - if (point_idx == std::numeric_limits<std::size_t>::max()) + if (point_idx == (std::numeric_limits<std::size_t>::max)()) continue; Star const& star = m_stars[point_idx]; @@ -1689,7 +1686,7 @@ class Tangential_complex { for (; it_point_idx != s.end(); ++it_point_idx) { std::size_t point_idx = *it_point_idx; // Don't check infinite simplices - if (point_idx == std::numeric_limits<std::size_t>::max()) + if (point_idx == (std::numeric_limits<std::size_t>::max)()) continue; Star const& star = m_stars[point_idx]; @@ -1900,7 +1897,7 @@ class Tangential_complex { #ifdef GUDHI_TC_EXPORT_ALL_COORDS_IN_OFF int num_coords = m_ambient_dim; #else - int num_coords = std::min(m_ambient_dim, 3); + int num_coords = (std::min)(m_ambient_dim, 3); #endif #ifdef GUDHI_TC_EXPORT_NORMALS @@ -1955,7 +1952,7 @@ class Tangential_complex { Triangulation const& tr = it_tr->tr(); Tr_vertex_handle center_vh = it_tr->center_vertex(); - if (&tr == NULL || tr.current_dimension() < m_intrinsic_dim) + if (tr.current_dimension() < m_intrinsic_dim) continue; // Color for this star @@ -2154,7 +2151,7 @@ class Tangential_complex { typedef std::vector<Simplex> Triangles; Triangles triangles; - std::size_t num_vertices = c.size(); + int num_vertices = static_cast<int>(c.size()); // Do not export smaller dimension simplices if (num_vertices < m_intrinsic_dim + 1) continue; diff --git a/src/Tangential_complex/test/test_tangential_complex.cpp b/src/Tangential_complex/test/test_tangential_complex.cpp index f8b0d2fb..48156440 100644 --- a/src/Tangential_complex/test/test_tangential_complex.cpp +++ b/src/Tangential_complex/test/test_tangential_complex.cpp @@ -68,3 +68,61 @@ BOOST_AUTO_TEST_CASE(test_Spatial_tree_data_structure) { Gudhi::Simplex_tree<> stree; tc.create_complex(stree); } + +BOOST_AUTO_TEST_CASE(test_mini_tangential) { + typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> Kernel; + typedef Kernel::Point_d Point; + typedef tc::Tangential_complex<Kernel, CGAL::Dynamic_dimension_tag, CGAL::Parallel_tag> TC; + + + const int INTRINSIC_DIM = 1; + + // Generate points on a 2-sphere + std::vector<Point> points; + // [[0, 0], [1, 0], [0, 1], [1, 1]] + std::vector<double> point = {0.0, 0.0}; + points.push_back(Point(point.size(), point.begin(), point.end())); + point = {1.0, 0.0}; + points.push_back(Point(point.size(), point.begin(), point.end())); + point = {0.0, 1.0}; + points.push_back(Point(point.size(), point.begin(), point.end())); + point = {1.0, 1.0}; + points.push_back(Point(point.size(), point.begin(), point.end())); + std::cout << "points = " << points.size() << std::endl; + Kernel k; + + // Compute the TC + TC tc(points, INTRINSIC_DIM, k); + tc.compute_tangential_complex(); + TC::Num_inconsistencies num_inc = tc.number_of_inconsistent_simplices(); + std::cout << "TC vertices = " << tc.number_of_vertices() << " - simplices = " << num_inc.num_simplices << + " - inc simplices = " << num_inc.num_inconsistent_simplices << + " - inc stars = " << num_inc.num_inconsistent_stars << std::endl; + + BOOST_CHECK(tc.number_of_vertices() == 4); + BOOST_CHECK(num_inc.num_simplices == 4); + BOOST_CHECK(num_inc.num_inconsistent_simplices == 0); + BOOST_CHECK(num_inc.num_inconsistent_stars == 0); + + // Export the TC into a Simplex_tree + Gudhi::Simplex_tree<> stree; + tc.create_complex(stree); + std::cout << "ST vertices = " << stree.num_vertices() << " - simplices = " << stree.num_simplices() << std::endl; + + BOOST_CHECK(stree.num_vertices() == 4); + BOOST_CHECK(stree.num_simplices() == 6); + + tc.fix_inconsistencies_using_perturbation(0.01, 30.0); + + BOOST_CHECK(tc.number_of_vertices() == 4); + BOOST_CHECK(num_inc.num_simplices == 4); + BOOST_CHECK(num_inc.num_inconsistent_simplices == 0); + BOOST_CHECK(num_inc.num_inconsistent_stars == 0); + + // Export the TC into a Simplex_tree + tc.create_complex(stree); + std::cout << "ST vertices = " << stree.num_vertices() << " - simplices = " << stree.num_simplices() << std::endl; + + BOOST_CHECK(stree.num_vertices() == 4); + BOOST_CHECK(stree.num_simplices() == 6); +} diff --git a/src/cmake/modules/GUDHI_user_version_target.txt b/src/cmake/modules/GUDHI_user_version_target.txt index 51553e7e..c91d5659 100644 --- a/src/cmake/modules/GUDHI_user_version_target.txt +++ b/src/cmake/modules/GUDHI_user_version_target.txt @@ -48,7 +48,7 @@ if (NOT CMAKE_VERSION VERSION_LESS 2.8.11) add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_SOURCE_DIR}/src/GudhUI ${GUDHI_USER_VERSION_DIR}/GudhUI) - set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Contraction;Hasse_complex;Persistent_cohomology;Simplex_tree;Skeleton_blocker;Spatial_searching;Subsampling;Tangential_complex;Witness_complex") + set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Bottleneck_distance;Contraction;Hasse_complex;Persistent_cohomology;Simplex_tree;Skeleton_blocker;Spatial_searching;Subsampling;Tangential_complex;Witness_complex") foreach(GUDHI_MODULE ${GUDHI_MODULES}) # doc files @@ -81,21 +81,21 @@ if (NOT CMAKE_VERSION VERSION_LESS 2.8.11) endif() endforeach() - # include files - file(GLOB GUDHI_INCLUDE_FILES ${CMAKE_SOURCE_DIR}/src/${GUDHI_MODULE}/include/gudhi/*) + # include files and patches + file(GLOB GUDHI_INCLUDE_FILES ${CMAKE_SOURCE_DIR}/src/${GUDHI_MODULE}/include/*) foreach(GUDHI_INCLUDE_FILE ${GUDHI_INCLUDE_FILES}) get_filename_component(GUDHI_INCLUDE_FILE_NAME ${GUDHI_INCLUDE_FILE} NAME) # GUDHI_INCLUDE_FILE can be a file or a directory if(IS_DIRECTORY ${GUDHI_INCLUDE_FILE}) add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E - copy_directory ${GUDHI_INCLUDE_FILE} ${GUDHI_USER_VERSION_DIR}/include/gudhi/${GUDHI_INCLUDE_FILE_NAME}) + copy_directory ${GUDHI_INCLUDE_FILE} ${GUDHI_USER_VERSION_DIR}/include/${GUDHI_INCLUDE_FILE_NAME}) else() add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E - copy ${GUDHI_INCLUDE_FILE} ${GUDHI_USER_VERSION_DIR}/include/gudhi/${GUDHI_INCLUDE_FILE_NAME}) + copy ${GUDHI_INCLUDE_FILE} ${GUDHI_USER_VERSION_DIR}/include/${GUDHI_INCLUDE_FILE_NAME}) endif() endforeach() - + # concept files file(GLOB GUDHI_CONCEPT_FILES ${CMAKE_SOURCE_DIR}/src/${GUDHI_MODULE}/concept/*.h) diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h index 46c43a90..9fed4106 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.h @@ -28,6 +28,7 @@ <b>Author:</b> Vincent Rouvreau<br> <b>Introduced in:</b> GUDHI 1.3.0<br> <b>Copyright:</b> GPL v3<br> + <b>Requires:</b> \ref cgal ≥ 4.7.0 and \ref eigen3 </td> <td width="75%"> Alpha_complex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation.<br> @@ -130,6 +131,26 @@ </table> \section Toolbox Toolbox + \subsection BottleneckDistanceToolbox Bottleneck distance + \image html "perturb_pd.png" "Bottleneck distance is the length of the longest edge" +<table border="0"> + <tr> + <td width="25%"> + <b>Author:</b> François Godi<br> + <b>Introduced in:</b> GUDHI 1.4.0<br> + <b>Copyright:</b> GPL v3<br> + <b>Requires:</b> \ref cgal ≥ 4.8.0 and \ref eigen3 + </td> + <td width="75%"> + Bottleneck distance measures the similarity between two persistence diagrams. + It's the shortest distance b for which there exists a perfect matching between + the points of the two diagrams (+ all the diagonal points) such that + any couple of matched points are at distance at most b. + <br> + <b>User manual:</b> \ref bottleneck_distance + </td> + </tr> +</table> \subsection ContractionToolbox Contraction \image html "sphere_contraction_representation.png" "Sphere contraction example" <table border="0"> @@ -196,7 +217,7 @@ make \endverbatim * \verbatim make test \endverbatim * * \section optionallibrary Optional third-party library - * \subsection gmp GMP: + * \subsection gmp GMP * The multi-field persistent homology algorithm requires GMP which is a free library for arbitrary-precision * arithmetic, operating on signed integers, rational numbers, and floating point numbers. * @@ -209,11 +230,11 @@ make \endverbatim * * Having GMP version 4.2 or higher installed is recommended. * - * \subsection cgal CGAL: - * The \ref alpha_complex data structure and few examples requires CGAL, which is a C++ library which provides easy - * access to efficient and reliable geometric algorithms. + * \subsection cgal CGAL + * The \ref alpha_complex data structure, \ref bottleneck_distance, and few examples requires CGAL, which is a C++ + * library which provides easy access to efficient and reliable geometric algorithms. * - * Having CGAL version 4.4 or higher installed is recommended. The procedure to install this library according to + * Having CGAL version 4.4.0 or higher installed is recommended. The procedure to install this library according to * your operating system is detailed here http://doc.cgal.org/latest/Manual/installation.html * * The following examples require the <a target="_blank" href="http://www.cgal.org/">Computational Geometry Algorithms @@ -223,11 +244,11 @@ make \endverbatim * \li <a href="_simplex_tree_2example_alpha_shapes_3_simplex_tree_from_off_file_8cpp-example.html"> * Simplex_tree/example_alpha_shapes_3_simplex_tree_from_off_file.cpp</a> * - * The following example requires CGAL version ≥ 4.6: + * The following example requires CGAL version ≥ 4.6.0: * \li <a href="_witness_complex_2witness_complex_sphere_8cpp-example.html"> * Witness_complex/witness_complex_sphere.cpp</a> * - * The following example requires CGAL version ≥ 4.7: + * The following example requires CGAL version ≥ 4.7.0: * \li <a href="_alpha_complex_2_alpha_complex_from_off_8cpp-example.html"> * Alpha_complex/Alpha_complex_from_off.cpp</a> * \li <a href="_alpha_complex_2_alpha_complex_from_points_8cpp-example.html"> @@ -239,7 +260,11 @@ make \endverbatim * \li <a href="_persistent_cohomology_2custom_persistence_sort_8cpp-example.html"> * Persistent_cohomology/custom_persistence_sort.cpp</a> * - * \subsection eigen3 Eigen3: + * The following example requires CGAL version ≥ 4.8.0: + * \li <a href="_bottleneck_distance_2_bottleneck_example_8cpp-example.html"> + * + * Bottleneck_distance/bottleneck_example.cpp</a> + * \subsection eigen3 Eigen3 * The \ref alpha_complex data structure and few examples requires * <a target="_blank" href="http://eigen.tuxfamily.org/">Eigen3</a> is a C++ template library for linear algebra: * matrices, vectors, numerical solvers, and related algorithms. @@ -247,9 +272,11 @@ make \endverbatim * The following example requires the <a target="_blank" href="http://eigen.tuxfamily.org/">Eigen3</a> and will not be * built if Eigen3 is not installed: * \li <a href="_alpha_complex_2_alpha_complex_from_off_8cpp-example.html"> - * Alpha_complex/Alpha_complex_from_off.cpp</a> (requires also Eigen3) + * Alpha_complex/Alpha_complex_from_off.cpp</a> * \li <a href="_alpha_complex_2_alpha_complex_from_points_8cpp-example.html"> - * Alpha_complex/Alpha_complex_from_points.cpp</a> (requires also Eigen3) + * Alpha_complex/Alpha_complex_from_points.cpp</a> + * \li <a href="_bottleneck_distance_2_bottleneck_example_8cpp-example.html"> + * Bottleneck_distance/bottleneck_example.cpp</a> * \li <a href="_persistent_cohomology_2alpha_complex_persistence_8cpp-example.html"> * Persistent_cohomology/alpha_complex_persistence.cpp</a> * \li <a href="_persistent_cohomology_2periodic_alpha_complex_3d_persistence_8cpp-example.html"> @@ -257,7 +284,7 @@ make \endverbatim * \li <a href="_persistent_cohomology_2custom_persistence_sort_8cpp-example.html"> * Persistent_cohomology/custom_persistence_sort.cpp</a> * - * \subsection tbb Threading Building Blocks: + * \subsection tbb Threading Building Blocks * <a target="_blank" href="https://www.threadingbuildingblocks.org/">Intel® TBB</a> lets you easily write parallel * C++ programs that take full advantage of multicore performance, that are portable and composable, and that have * future-proof scalability. @@ -331,6 +358,7 @@ make \endverbatim /*! @file Examples * @example Alpha_complex/Alpha_complex_from_off.cpp * @example Alpha_complex/Alpha_complex_from_points.cpp + * @example Bottleneck_distance/bottleneck_example.cpp * @example Bitmap_cubical_complex/Bitmap_cubical_complex.cpp * @example Bitmap_cubical_complex/Bitmap_cubical_complex_periodic_boundary_conditions.cpp * @example Bitmap_cubical_complex/Random_bitmap_cubical_complex.cpp diff --git a/src/common/include/gudhi/random_point_generators.h b/src/common/include/gudhi/random_point_generators.h index 3050b7ea..2ec465ef 100644 --- a/src/common/include/gudhi/random_point_generators.h +++ b/src/common/include/gudhi/random_point_generators.h @@ -192,10 +192,9 @@ static void generate_uniform_points_on_torus_d(const Kernel &k, int dim, std::si double radius_noise_percentage = 0., std::vector<typename Kernel::FT> current_point = std::vector<typename Kernel::FT>()) { CGAL::Random rng; - if (current_point.size() == 2 * dim) { - *out++ = k.construct_point_d_object()( - static_cast<int> (current_point.size()), - current_point.begin(), current_point.end()); + int point_size = static_cast<int>(current_point.size()); + if (point_size == 2 * dim) { + *out++ = k.construct_point_d_object()(point_size, current_point.begin(), current_point.end()); } else { for (std::size_t slice_idx = 0; slice_idx < num_slices; ++slice_idx) { double radius_noise_ratio = 1.; @@ -338,8 +337,6 @@ std::vector<typename Kernel::Point_d> generate_points_on_3sphere_and_circle(std: std::vector<Point> points; points.reserve(num_points); - typename Kernel::Translated_point_d k_transl = - k.translated_point_d_object(); typename Kernel::Compute_coordinate_d k_coord = k.compute_coordinate_d_object(); for (std::size_t i = 0; i < num_points;) { diff --git a/src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt b/src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt new file mode 100644 index 00000000..a588d113 --- /dev/null +++ b/src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt @@ -0,0 +1,3 @@ +CGAL/Kd_tree.h +CGAL/Kd_tree_node.h +CGAL/Orthogonal_incremental_neighbor_search.h diff --git a/src/common/include/gudhi_patches/CGAL/Kd_tree.h b/src/common/include/gudhi_patches/CGAL/Kd_tree.h new file mode 100644 index 00000000..f085b0da --- /dev/null +++ b/src/common/include/gudhi_patches/CGAL/Kd_tree.h @@ -0,0 +1,582 @@ +// Copyright (c) 2002,2011,2014 Utrecht University (The Netherlands), Max-Planck-Institute Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// 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. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Hans Tangelder (<hanst@cs.uu.nl>), +// : Waqar Khan <wkhan@mpi-inf.mpg.de> + +#ifndef CGAL_KD_TREE_H +#define CGAL_KD_TREE_H + +#include "Kd_tree_node.h" + +#include <CGAL/basic.h> +#include <CGAL/assertions.h> +#include <vector> + +#include <CGAL/algorithm.h> +#include <CGAL/internal/Get_dimension_tag.h> +#include <CGAL/Search_traits.h> + + +#include <deque> +#include <boost/container/deque.hpp> +#include <boost/optional.hpp> + +#ifdef CGAL_HAS_THREADS +#include <CGAL/mutex.h> +#endif + +namespace CGAL { + +//template <class SearchTraits, class Splitter_=Median_of_rectangle<SearchTraits>, class UseExtendedNode = Tag_true > +template <class SearchTraits, class Splitter_=Sliding_midpoint<SearchTraits>, class UseExtendedNode = Tag_true > +class Kd_tree { + +public: + typedef SearchTraits Traits; + typedef Splitter_ Splitter; + typedef typename SearchTraits::Point_d Point_d; + typedef typename Splitter::Container Point_container; + + typedef typename SearchTraits::FT FT; + typedef Kd_tree_node<SearchTraits, Splitter, UseExtendedNode > Node; + typedef Kd_tree_leaf_node<SearchTraits, Splitter, UseExtendedNode > Leaf_node; + typedef Kd_tree_internal_node<SearchTraits, Splitter, UseExtendedNode > Internal_node; + typedef Kd_tree<SearchTraits, Splitter> Tree; + typedef Kd_tree<SearchTraits, Splitter,UseExtendedNode> Self; + + typedef Node* Node_handle; + typedef const Node* Node_const_handle; + typedef Leaf_node* Leaf_node_handle; + typedef const Leaf_node* Leaf_node_const_handle; + typedef Internal_node* Internal_node_handle; + typedef const Internal_node* Internal_node_const_handle; + typedef typename std::vector<const Point_d*>::const_iterator Point_d_iterator; + typedef typename std::vector<const Point_d*>::const_iterator Point_d_const_iterator; + typedef typename Splitter::Separator Separator; + typedef typename std::vector<Point_d>::const_iterator iterator; + typedef typename std::vector<Point_d>::const_iterator const_iterator; + + typedef typename std::vector<Point_d>::size_type size_type; + + typedef typename internal::Get_dimension_tag<SearchTraits>::Dimension D; + +private: + SearchTraits traits_; + Splitter split; + + + // wokaround for https://svn.boost.org/trac/boost/ticket/9332 +#if (_MSC_VER == 1800) && (BOOST_VERSION == 105500) + std::deque<Internal_node> internal_nodes; + std::deque<Leaf_node> leaf_nodes; +#else + boost::container::deque<Internal_node> internal_nodes; + boost::container::deque<Leaf_node> leaf_nodes; +#endif + + Node_handle tree_root; + + Kd_tree_rectangle<FT,D>* bbox; + std::vector<Point_d> pts; + + // Instead of storing the points in arrays in the Kd_tree_node + // we put all the data in a vector in the Kd_tree. + // and we only store an iterator range in the Kd_tree_node. + // + std::vector<const Point_d*> data; + + + #ifdef CGAL_HAS_THREADS + mutable CGAL_MUTEX building_mutex;//mutex used to protect const calls inducing build() + #endif + bool built_; + bool removed_; + + // protected copy constructor + Kd_tree(const Tree& tree) + : traits_(tree.traits_),built_(tree.built_) + {}; + + + // Instead of the recursive construction of the tree in the class Kd_tree_node + // we do this in the tree class. The advantage is that we then can optimize + // the allocation of the nodes. + + // The leaf node + Node_handle + create_leaf_node(Point_container& c) + { + Leaf_node node(true , static_cast<unsigned int>(c.size())); + std::ptrdiff_t tmp = c.begin() - data.begin(); + node.data = pts.begin() + tmp; + + leaf_nodes.push_back(node); + Leaf_node_handle nh = &leaf_nodes.back(); + + + return nh; + } + + + // The internal node + + Node_handle + create_internal_node(Point_container& c, const Tag_true&) + { + return create_internal_node_use_extension(c); + } + + Node_handle + create_internal_node(Point_container& c, const Tag_false&) + { + return create_internal_node(c); + } + + + + // TODO: Similiar to the leaf_init function above, a part of the code should be + // moved to a the class Kd_tree_node. + // It is not proper yet, but the goal was to see if there is + // a potential performance gain through the Compact_container + Node_handle + create_internal_node_use_extension(Point_container& c) + { + Internal_node node(false); + internal_nodes.push_back(node); + Internal_node_handle nh = &internal_nodes.back(); + + Separator sep; + Point_container c_low(c.dimension(),traits_); + split(sep, c, c_low); + nh->set_separator(sep); + + int cd = nh->cutting_dimension(); + if(!c_low.empty()){ + nh->lower_low_val = c_low.tight_bounding_box().min_coord(cd); + nh->lower_high_val = c_low.tight_bounding_box().max_coord(cd); + } + else{ + nh->lower_low_val = nh->cutting_value(); + nh->lower_high_val = nh->cutting_value(); + } + if(!c.empty()){ + nh->upper_low_val = c.tight_bounding_box().min_coord(cd); + nh->upper_high_val = c.tight_bounding_box().max_coord(cd); + } + else{ + nh->upper_low_val = nh->cutting_value(); + nh->upper_high_val = nh->cutting_value(); + } + + CGAL_assertion(nh->cutting_value() >= nh->lower_low_val); + CGAL_assertion(nh->cutting_value() <= nh->upper_high_val); + + if (c_low.size() > split.bucket_size()){ + nh->lower_ch = create_internal_node_use_extension(c_low); + }else{ + nh->lower_ch = create_leaf_node(c_low); + } + if (c.size() > split.bucket_size()){ + nh->upper_ch = create_internal_node_use_extension(c); + }else{ + nh->upper_ch = create_leaf_node(c); + } + + + + + return nh; + } + + + // Note also that I duplicated the code to get rid if the if's for + // the boolean use_extension which was constant over the construction + Node_handle + create_internal_node(Point_container& c) + { + Internal_node node(false); + internal_nodes.push_back(node); + Internal_node_handle nh = &internal_nodes.back(); + Separator sep; + + Point_container c_low(c.dimension(),traits_); + split(sep, c, c_low); + nh->set_separator(sep); + + if (c_low.size() > split.bucket_size()){ + nh->lower_ch = create_internal_node(c_low); + }else{ + nh->lower_ch = create_leaf_node(c_low); + } + if (c.size() > split.bucket_size()){ + nh->upper_ch = create_internal_node(c); + }else{ + nh->upper_ch = create_leaf_node(c); + } + + + + return nh; + } + + + +public: + + Kd_tree(Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) + : traits_(traits),split(s), built_(false), removed_(false) + {} + + template <class InputIterator> + Kd_tree(InputIterator first, InputIterator beyond, + Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) + : traits_(traits),split(s), built_(false), removed_(false) + { + pts.insert(pts.end(), first, beyond); + } + + bool empty() const { + return pts.empty(); + } + + void + build() + { + // This function is not ready to be called when a tree already exists, one + // must call invalidate_built() first. + CGAL_assertion(!is_built()); + CGAL_assertion(!removed_); + const Point_d& p = *pts.begin(); + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits_.construct_cartesian_const_iterator_d_object(); + int dim = static_cast<int>(std::distance(ccci(p), ccci(p,0))); + + data.reserve(pts.size()); + for(unsigned int i = 0; i < pts.size(); i++){ + data.push_back(&pts[i]); + } + Point_container c(dim, data.begin(), data.end(),traits_); + bbox = new Kd_tree_rectangle<FT,D>(c.bounding_box()); + if (c.size() <= split.bucket_size()){ + tree_root = create_leaf_node(c); + }else { + tree_root = create_internal_node(c, UseExtendedNode()); + } + + //Reorder vector for spatial locality + std::vector<Point_d> ptstmp; + ptstmp.resize(pts.size()); + for (std::size_t i = 0; i < pts.size(); ++i){ + ptstmp[i] = *data[i]; + } + for(std::size_t i = 0; i < leaf_nodes.size(); ++i){ + std::ptrdiff_t tmp = leaf_nodes[i].begin() - pts.begin(); + leaf_nodes[i].data = ptstmp.begin() + tmp; + } + pts.swap(ptstmp); + + data.clear(); + + built_ = true; + } + +private: + //any call to this function is for the moment not threadsafe + void const_build() const { + #ifdef CGAL_HAS_THREADS + //this ensure that build() will be called once + CGAL_SCOPED_LOCK(building_mutex); + if(!is_built()) + #endif + const_cast<Self*>(this)->build(); //THIS IS NOT THREADSAFE + } +public: + + bool is_built() const + { + return built_; + } + + void invalidate_built() + { + if(removed_){ + // Walk the tree to collect the remaining points. + // Writing directly to pts would likely work, but better be safe. + std::vector<Point_d> ptstmp; + //ptstmp.resize(root()->num_items()); + root()->tree_items(std::back_inserter(ptstmp)); + pts.swap(ptstmp); + removed_=false; + CGAL_assertion(is_built()); // the rest of the cleanup must happen + } + if(is_built()){ + internal_nodes.clear(); + leaf_nodes.clear(); + data.clear(); + delete bbox; + built_ = false; + } + } + + void clear() + { + invalidate_built(); + pts.clear(); + removed_ = false; + } + + void + insert(const Point_d& p) + { + invalidate_built(); + pts.push_back(p); + } + + template <class InputIterator> + void + insert(InputIterator first, InputIterator beyond) + { + invalidate_built(); + pts.insert(pts.end(),first, beyond); + } + +private: + struct Equal_by_coordinates { + SearchTraits const* traits; + Point_d const* pp; + bool operator()(Point_d const&q) const { + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits->construct_cartesian_const_iterator_d_object(); + return std::equal(ccci(*pp), ccci(*pp,0), ccci(q)); + } + }; + Equal_by_coordinates equal_by_coordinates(Point_d const&p){ + Equal_by_coordinates ret = { &traits(), &p }; + return ret; + } + +public: + void + remove(const Point_d& p) + { + remove(p, equal_by_coordinates(p)); + } + + template<class Equal> + void + remove(const Point_d& p, Equal const& equal_to_p) + { +#if 0 + // This code could have quadratic runtime. + if (!is_built()) { + std::vector<Point_d>::iterator pi = std::find(pts.begin(), pts.end(), p); + // Precondition: the point must be there. + CGAL_assertion (pi != pts.end()); + pts.erase(pi); + return; + } +#endif + bool success = remove_(p, 0, false, 0, false, root(), equal_to_p); + CGAL_assertion(success); + + // Do not set the flag is the tree has been cleared. + if(is_built()) + removed_ |= success; + } +private: + template<class Equal> + bool remove_(const Point_d& p, + Internal_node_handle grandparent, bool parent_islower, + Internal_node_handle parent, bool islower, + Node_handle node, Equal const& equal_to_p) { + // Recurse to locate the point + if (!node->is_leaf()) { + Internal_node_handle newparent = static_cast<Internal_node_handle>(node); + // FIXME: This should be if(x<y) remove low; else remove up; + if (traits().construct_cartesian_const_iterator_d_object()(p)[newparent->cutting_dimension()] <= newparent->cutting_value()) { + if (remove_(p, parent, islower, newparent, true, newparent->lower(), equal_to_p)) + return true; + } + //if (traits().construct_cartesian_const_iterator_d_object()(p)[newparent->cutting_dimension()] >= newparent->cutting_value()) + return remove_(p, parent, islower, newparent, false, newparent->upper(), equal_to_p); + + CGAL_assertion(false); // Point was not found + } + + // Actual removal + Leaf_node_handle lnode = static_cast<Leaf_node_handle>(node); + if (lnode->size() > 1) { + iterator pi = std::find_if(lnode->begin(), lnode->end(), equal_to_p); + // FIXME: we should ensure this never happens + if (pi == lnode->end()) return false; + iterator lasti = lnode->end() - 1; + if (pi != lasti) { + // Hack to get a non-const iterator + std::iter_swap(pts.begin()+(pi-pts.begin()), pts.begin()+(lasti-pts.begin())); + } + lnode->drop_last_point(); + } else if (!equal_to_p(*lnode->begin())) { + // FIXME: we should ensure this never happens + return false; + } else if (grandparent) { + Node_handle brother = islower ? parent->upper() : parent->lower(); + if (parent_islower) + grandparent->set_lower(brother); + else + grandparent->set_upper(brother); + } else if (parent) { + tree_root = islower ? parent->upper() : parent->lower(); + } else { + clear(); + } + return true; + } + +public: + //For efficiency; reserve the size of the points vectors in advance (if the number of points is already known). + void reserve(size_t size) + { + pts.reserve(size); + } + + //Get the capacity of the underlying points vector. + size_t capacity() + { + return pts.capacity(); + } + + + template <class OutputIterator, class FuzzyQueryItem> + OutputIterator + search(OutputIterator it, const FuzzyQueryItem& q) const + { + if(! pts.empty()){ + + if(! is_built()){ + const_build(); + } + Kd_tree_rectangle<FT,D> b(*bbox); + return tree_root->search(it,q,b); + } + return it; + } + + + template <class FuzzyQueryItem> + boost::optional<Point_d> + search_any_point(const FuzzyQueryItem& q) const + { + if(! pts.empty()){ + + if(! is_built()){ + const_build(); + } + Kd_tree_rectangle<FT,D> b(*bbox); + return tree_root->search_any_point(q,b); + } + return boost::none; + } + + + ~Kd_tree() { + if(is_built()){ + delete bbox; + } + } + + + const SearchTraits& + traits() const + { + return traits_; + } + + Node_const_handle + root() const + { + if(! is_built()){ + const_build(); + } + return tree_root; + } + + Node_handle + root() + { + if(! is_built()){ + build(); + } + return tree_root; + } + + void + print() const + { + if(! is_built()){ + const_build(); + } + root()->print(); + } + + const Kd_tree_rectangle<FT,D>& + bounding_box() const + { + if(! is_built()){ + const_build(); + } + return *bbox; + } + + const_iterator + begin() const + { + return pts.begin(); + } + + const_iterator + end() const + { + return pts.end(); + } + + size_type + size() const + { + return pts.size(); + } + + // Print statistics of the tree. + std::ostream& + statistics(std::ostream& s) const + { + if(! is_built()){ + const_build(); + } + s << "Tree statistics:" << std::endl; + s << "Number of items stored: " + << root()->num_items() << std::endl; + s << "Number of nodes: " + << root()->num_nodes() << std::endl; + s << " Tree depth: " << root()->depth() << std::endl; + return s; + } + + +}; + +} // namespace CGAL + +#endif // CGAL_KD_TREE_H diff --git a/src/common/include/gudhi_patches/CGAL/Kd_tree_node.h b/src/common/include/gudhi_patches/CGAL/Kd_tree_node.h new file mode 100644 index 00000000..909ee260 --- /dev/null +++ b/src/common/include/gudhi_patches/CGAL/Kd_tree_node.h @@ -0,0 +1,586 @@ +// Copyright (c) 2002,2011 Utrecht University (The Netherlands). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// 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. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Authors : Hans Tangelder (<hanst@cs.uu.nl>) + +#ifndef CGAL_KD_TREE_NODE_H +#define CGAL_KD_TREE_NODE_H + +#include "CGAL/Splitters.h" + +#include <CGAL/Compact_container.h> +#include <boost/cstdint.hpp> + +namespace CGAL { + + template <class SearchTraits, class Splitter, class UseExtendedNode> + class Kd_tree; + + template < class TreeTraits, class Splitter, class UseExtendedNode > + class Kd_tree_node { + + friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>; + + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_handle Node_handle; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_const_handle Node_const_handle; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Internal_node_handle Internal_node_handle; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Internal_node_const_handle Internal_node_const_handle; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Leaf_node_handle Leaf_node_handle; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Leaf_node_const_handle Leaf_node_const_handle; + typedef typename TreeTraits::Point_d Point_d; + + typedef typename TreeTraits::FT FT; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Separator Separator; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Point_d_iterator Point_d_iterator; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::iterator iterator; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::D D; + + bool leaf; + + public : + Kd_tree_node(bool leaf_) + :leaf(leaf_){} + + bool is_leaf() const{ + return leaf; + } + + std::size_t + num_items() const + { + if (is_leaf()){ + Leaf_node_const_handle node = + static_cast<Leaf_node_const_handle>(this); + return node->size(); + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + return node->lower()->num_items() + node->upper()->num_items(); + } + } + + std::size_t + num_nodes() const + { + if (is_leaf()) return 1; + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + return node->lower()->num_nodes() + node->upper()->num_nodes(); + } + } + + int + depth(const int current_max_depth) const + { + if (is_leaf()){ + return current_max_depth; + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + return + (std::max)( node->lower()->depth(current_max_depth + 1), + node->upper()->depth(current_max_depth + 1)); + } + } + + int + depth() const + { + return depth(1); + } + + template <class OutputIterator> + OutputIterator + tree_items(OutputIterator it) const { + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast<Leaf_node_const_handle>(this); + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + {*it=*i; ++it;} + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + it=node->lower()->tree_items(it); + it=node->upper()->tree_items(it); + } + return it; + } + + + boost::optional<Point_d> + any_tree_item() const { + boost::optional<Point_d> result = boost::none; + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast<Leaf_node_const_handle>(this); + if (node->size()>0){ + return boost::make_optional(*(node->begin())); + } + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + result = node->lower()->any_tree_item(); + if(! result){ + result = node->upper()->any_tree_item(); + } + } + return result; + } + + + void + indent(int d) const + { + for(int i = 0; i < d; i++){ + std::cout << " "; + } + } + + + void + print(int d = 0) const + { + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast<Leaf_node_const_handle>(this); + indent(d); + std::cout << "leaf" << std::endl; + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + {indent(d);std::cout << *i << std::endl;} + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + indent(d); + std::cout << "lower tree" << std::endl; + node->lower()->print(d+1); + indent(d); + std::cout << "separator: dim = " << node->cutting_dimension() << " val = " << node->cutting_value() << std::endl; + indent(d); + std::cout << "upper tree" << std::endl; + node->upper()->print(d+1); + } + } + + + template <class OutputIterator, class FuzzyQueryItem> + OutputIterator + search(OutputIterator it, const FuzzyQueryItem& q, + Kd_tree_rectangle<FT,D>& b) const + { + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast<Leaf_node_const_handle>(this); + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + if (q.contains(*i)) + {*it++=*i;} + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + // after splitting b denotes the lower part of b + Kd_tree_rectangle<FT,D> b_upper(b); + b.split(b_upper, node->cutting_dimension(), + node->cutting_value()); + + if (q.outer_range_contains(b)) + it=node->lower()->tree_items(it); + else + if (q.inner_range_intersects(b)) + it=node->lower()->search(it,q,b); + if (q.outer_range_contains(b_upper)) + it=node->upper()->tree_items(it); + else + if (q.inner_range_intersects(b_upper)) + it=node->upper()->search(it,q,b_upper); + }; + return it; + } + + + template <class FuzzyQueryItem> + boost::optional<Point_d> + search_any_point(const FuzzyQueryItem& q, + Kd_tree_rectangle<FT,D>& b) const + { + boost::optional<Point_d> result = boost::none; + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast<Leaf_node_const_handle>(this); + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + if (q.contains(*i)) + { result = *i; break; } + } + else { + Internal_node_const_handle node = + static_cast<Internal_node_const_handle>(this); + // after splitting b denotes the lower part of b + Kd_tree_rectangle<FT,D> b_upper(b); + b.split(b_upper, node->cutting_dimension(), + node->cutting_value()); + + if (q.outer_range_contains(b)){ + result = node->lower()->any_tree_item(); + }else{ + if (q.inner_range_intersects(b)){ + result = node->lower()->search_any_point(q,b); + } + } + if(result){ + return result; + } + if (q.outer_range_contains(b_upper)){ + result = node->upper()->any_tree_item(); + }else{ + if (q.inner_range_intersects(b_upper)) + result = node->upper()->search_any_point(q,b_upper); + } + } + return result; + } + + }; + + + template < class TreeTraits, class Splitter, class UseExtendedNode > + class Kd_tree_leaf_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ + + friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>; + + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::iterator iterator; + typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; + typedef typename TreeTraits::Point_d Point_d; + + private: + + // private variables for leaf nodes + boost::int32_t n; // denotes number of items in a leaf node + iterator data; // iterator to data in leaf node + + + public: + + // default constructor + Kd_tree_leaf_node() + {} + + Kd_tree_leaf_node(bool leaf_ ) + : Base(leaf_) + {} + + Kd_tree_leaf_node(bool leaf_,unsigned int n_ ) + : Base(leaf_), n(n_) + {} + + // members for all nodes + + // members for leaf nodes only + inline + unsigned int + size() const + { + return n; + } + + inline + iterator + begin() const + { + return data; + } + + inline + iterator + end() const + { + return data + n; + } + + inline + void + drop_last_point() + { + --n; + } + + }; //leaf node + + + + template < class TreeTraits, class Splitter, class UseExtendedNode> + class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ + + friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>; + + typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_handle Node_handle; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_const_handle Node_const_handle; + + typedef typename TreeTraits::FT FT; + typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Separator Separator; + + private: + + // private variables for internal nodes + boost::int32_t cut_dim; + FT cut_val; + Node_handle lower_ch, upper_ch; + + + // private variables for extended internal nodes + FT upper_low_val; + FT upper_high_val; + FT lower_low_val; + FT lower_high_val; + + + public: + + // default constructor + Kd_tree_internal_node() + {} + + Kd_tree_internal_node(bool leaf_) + : Base(leaf_) + {} + + + // members for internal node and extended internal node + + inline + Node_const_handle + lower() const + { + return lower_ch; + } + + inline + Node_const_handle + upper() const + { + return upper_ch; + } + + inline + Node_handle + lower() + { + return lower_ch; + } + + inline + Node_handle + upper() + { + return upper_ch; + } + + inline + void + set_lower(Node_handle nh) + { + lower_ch = nh; + } + + inline + void + set_upper(Node_handle nh) + { + upper_ch = nh; + } + + // inline Separator& separator() {return sep; } + // use instead + inline + void set_separator(Separator& sep){ + cut_dim = sep.cutting_dimension(); + cut_val = sep.cutting_value(); + } + + inline + FT + cutting_value() const + { + return cut_val; + } + + inline + int + cutting_dimension() const + { + return cut_dim; + } + + // members for extended internal node only + inline + FT + upper_low_value() const + { + return upper_low_val; + } + + inline + FT + upper_high_value() const + { + return upper_high_val; + } + + inline + FT + lower_low_value() const + { + return lower_low_val; + } + + inline + FT + lower_high_value() const + { + return lower_high_val; + } + + /*Separator& + separator() + { + return Separator(cutting_dimension,cutting_value); + }*/ + + + };//internal node + + template < class TreeTraits, class Splitter> + class Kd_tree_internal_node<TreeTraits,Splitter,Tag_false> : public Kd_tree_node< TreeTraits, Splitter, Tag_false >{ + + friend class Kd_tree<TreeTraits,Splitter,Tag_false>; + + typedef Kd_tree_node< TreeTraits, Splitter, Tag_false> Base; + typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Node_handle Node_handle; + typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Node_const_handle Node_const_handle; + + typedef typename TreeTraits::FT FT; + typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Separator Separator; + + private: + + // private variables for internal nodes + boost::uint8_t cut_dim; + FT cut_val; + + Node_handle lower_ch, upper_ch; + + public: + + // default constructor + Kd_tree_internal_node() + {} + + Kd_tree_internal_node(bool leaf_) + : Base(leaf_) + {} + + + // members for internal node and extended internal node + + inline + Node_const_handle + lower() const + { + return lower_ch; + } + + inline + Node_const_handle + upper() const + { + return upper_ch; + } + + inline + Node_handle + lower() + { + return lower_ch; + } + + inline + Node_handle + upper() + { + return upper_ch; + } + + inline + void + set_lower(Node_handle nh) + { + lower_ch = nh; + } + + inline + void + set_upper(Node_handle nh) + { + upper_ch = nh; + } + + // inline Separator& separator() {return sep; } + // use instead + + inline + void set_separator(Separator& sep){ + cut_dim = sep.cutting_dimension(); + cut_val = sep.cutting_value(); + } + + inline + FT + cutting_value() const + { + return cut_val; + } + + inline + int + cutting_dimension() const + { + return cut_dim; + } + + /* Separator& + separator() + { + return Separator(cutting_dimension,cutting_value); + }*/ + + + };//internal node + + + +} // namespace CGAL +#endif // CGAL_KDTREE_NODE_H diff --git a/src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h b/src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h new file mode 100644 index 00000000..e29ce14f --- /dev/null +++ b/src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h @@ -0,0 +1,620 @@ +// Copyright (c) 2002,2011 Utrecht University (The Netherlands). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// 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. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Hans Tangelder (<hanst@cs.uu.nl>) + +#ifndef CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH +#define CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH + +#include <CGAL/Kd_tree.h> +#include <cstring> +#include <list> +#include <queue> +#include <memory> +#include <CGAL/Euclidean_distance.h> +#include <CGAL/tuple.h> + +namespace CGAL { + + template <class SearchTraits, + class Distance_= typename internal::Spatial_searching_default_distance<SearchTraits>::type, + class Splitter_ = Sliding_midpoint<SearchTraits>, + class Tree_= Kd_tree<SearchTraits, Splitter_, Tag_true> > + class Orthogonal_incremental_neighbor_search { + + public: + typedef Splitter_ Splitter; + typedef Tree_ Tree; + typedef Distance_ Distance; + typedef typename SearchTraits::Point_d Point_d; + typedef typename Distance::Query_item Query_item; + typedef typename SearchTraits::FT FT; + typedef typename Tree::Point_d_iterator Point_d_iterator; + typedef typename Tree::Node_const_handle Node_const_handle; + + typedef std::pair<Point_d,FT> Point_with_transformed_distance; + typedef CGAL::cpp11::tuple<Node_const_handle,FT,std::vector<FT> > Node_with_distance; + typedef std::vector<Node_with_distance*> Node_with_distance_vector; + typedef std::vector<Point_with_transformed_distance*> Point_with_transformed_distance_vector; + + template<class T> + struct Object_wrapper + { + T object; + Object_wrapper(const T& t):object(t){} + const T& operator* () const { return object; } + const T* operator-> () const { return &object; } + }; + + class Iterator_implementation { + SearchTraits traits; + public: + + int number_of_neighbours_computed; + int number_of_internal_nodes_visited; + int number_of_leaf_nodes_visited; + int number_of_items_visited; + + private: + + typedef std::vector<FT> Distance_vector; + + Distance_vector dists; + + Distance Orthogonal_distance_instance; + + FT multiplication_factor; + + Query_item query_point; + + FT distance_to_root; + + bool search_nearest_neighbour; + + FT rd; + + + class Priority_higher { + public: + + bool search_nearest; + + Priority_higher(bool search_the_nearest_neighbour) + : search_nearest(search_the_nearest_neighbour) + {} + + //highest priority is smallest distance + bool + operator() (Node_with_distance* n1, Node_with_distance* n2) const + { + return (search_nearest) ? (CGAL::cpp11::get<1>(*n1) > CGAL::cpp11::get<1>(*n2)) : (CGAL::cpp11::get<1>(*n2) > CGAL::cpp11::get<1>(*n1)); + } + }; + + class Distance_smaller { + + public: + + bool search_nearest; + + Distance_smaller(bool search_the_nearest_neighbour) + : search_nearest(search_the_nearest_neighbour) + {} + + //highest priority is smallest distance + bool operator() (Point_with_transformed_distance* p1, Point_with_transformed_distance* p2) const + { + return (search_nearest) ? (p1->second > p2->second) : (p2->second > p1->second); + } + }; + + + std::priority_queue<Node_with_distance*, Node_with_distance_vector, + Priority_higher> PriorityQueue; + + public: + std::priority_queue<Point_with_transformed_distance*, Point_with_transformed_distance_vector, + Distance_smaller> Item_PriorityQueue; + + + public: + + int reference_count; + + + + // constructor + Iterator_implementation(const Tree& tree,const Query_item& q, const Distance& tr, + FT Eps=FT(0.0), bool search_nearest=true) + : traits(tree.traits()),number_of_neighbours_computed(0), number_of_internal_nodes_visited(0), + number_of_leaf_nodes_visited(0), number_of_items_visited(0), + Orthogonal_distance_instance(tr), multiplication_factor(Orthogonal_distance_instance.transformed_distance(FT(1.0)+Eps)), + query_point(q), search_nearest_neighbour(search_nearest), + PriorityQueue(Priority_higher(search_nearest)), Item_PriorityQueue(Distance_smaller(search_nearest)), + reference_count(1) + + + { + if (tree.empty()) return; + + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits.construct_cartesian_const_iterator_d_object(); + int dim = static_cast<int>(std::distance(ccci(q), ccci(q,0))); + + dists.resize(dim); + for(int i=0 ; i<dim ; ++i){ + dists[i] = 0; + } + + if (search_nearest){ + distance_to_root= + Orthogonal_distance_instance.min_distance_to_rectangle(q, tree.bounding_box(),dists); + Node_with_distance *The_Root = new Node_with_distance(tree.root(), + distance_to_root, dists); + PriorityQueue.push(The_Root); + + // rd is the distance of the top of the priority queue to q + rd=CGAL::cpp11::get<1>(*The_Root); + Compute_the_next_nearest_neighbour(); + } + else{ + distance_to_root= + Orthogonal_distance_instance.max_distance_to_rectangle(q, + tree.bounding_box(), dists); + Node_with_distance *The_Root = new Node_with_distance(tree.root(), + distance_to_root, dists); + PriorityQueue.push(The_Root); + + // rd is the distance of the top of the priority queue to q + rd=CGAL::cpp11::get<1>(*The_Root); + Compute_the_next_furthest_neighbour(); + } + + + } + + // * operator + const Point_with_transformed_distance& + operator* () const + { + return *(Item_PriorityQueue.top()); + } + + // prefix operator + Iterator_implementation& + operator++() + { + Delete_the_current_item_top(); + if(search_nearest_neighbour) + Compute_the_next_nearest_neighbour(); + else + Compute_the_next_furthest_neighbour(); + return *this; + } + + // postfix operator + Object_wrapper<Point_with_transformed_distance> + operator++(int) + { + Object_wrapper<Point_with_transformed_distance> result( *(Item_PriorityQueue.top()) ); + ++*this; + return result; + } + + // Print statistics of the general priority search process. + std::ostream& + statistics (std::ostream& s) const { + s << "Orthogonal priority search statistics:" + << std::endl; + s << "Number of internal nodes visited:" + << number_of_internal_nodes_visited << std::endl; + s << "Number of leaf nodes visited:" + << number_of_leaf_nodes_visited << std::endl; + s << "Number of items visited:" + << number_of_items_visited << std::endl; + s << "Number of neighbours computed:" + << number_of_neighbours_computed << std::endl; + return s; + } + + + //destructor + ~Iterator_implementation() + { + while (!PriorityQueue.empty()) { + Node_with_distance* The_top=PriorityQueue.top(); + PriorityQueue.pop(); + delete The_top; + } + while (!Item_PriorityQueue.empty()) { + Point_with_transformed_distance* The_top=Item_PriorityQueue.top(); + Item_PriorityQueue.pop(); + delete The_top; + } + } + + private: + + void + Delete_the_current_item_top() + { + Point_with_transformed_distance* The_item_top=Item_PriorityQueue.top(); + Item_PriorityQueue.pop(); + delete The_item_top; + } + + void + Compute_the_next_nearest_neighbour() + { + // compute the next item + bool next_neighbour_found=false; + if (!(Item_PriorityQueue.empty())) { + next_neighbour_found= + (multiplication_factor*rd > Item_PriorityQueue.top()->second); + } + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point); + // otherwise browse the tree further + while ((!next_neighbour_found) && (!PriorityQueue.empty())) { + Node_with_distance* The_node_top=PriorityQueue.top(); + Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top); + dists = CGAL::cpp11::get<2>(*The_node_top); + PriorityQueue.pop(); + delete The_node_top; + FT copy_rd=rd; + while (!(N->is_leaf())) { // compute new distance + typename Tree::Internal_node_const_handle node = + static_cast<typename Tree::Internal_node_const_handle>(N); + number_of_internal_nodes_visited++; + int new_cut_dim=node->cutting_dimension(); + FT new_rd,dst = dists[new_cut_dim]; + FT val = *(query_point_it + new_cut_dim); + FT diff1 = val - node->upper_low_value(); + FT diff2 = val - node->lower_high_value(); + if (diff1 + diff2 < FT(0.0)) { + new_rd= + Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); + + CGAL_assertion(new_rd >= copy_rd); + dists[new_cut_dim] = diff1; + Node_with_distance *Upper_Child = + new Node_with_distance(node->upper(), new_rd, dists); + PriorityQueue.push(Upper_Child); + dists[new_cut_dim] = dst; + N=node->lower(); + + } + else { // compute new distance + new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); + CGAL_assertion(new_rd >= copy_rd); + dists[new_cut_dim] = diff2; + Node_with_distance *Lower_Child = + new Node_with_distance(node->lower(), new_rd, dists); + PriorityQueue.push(Lower_Child); + dists[new_cut_dim] = dst; + N=node->upper(); + } + } + // n is a leaf + typename Tree::Leaf_node_const_handle node = + static_cast<typename Tree::Leaf_node_const_handle>(N); + number_of_leaf_nodes_visited++; + if (node->size() > 0) { + for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { + number_of_items_visited++; + FT distance_to_query_point= + Orthogonal_distance_instance.transformed_distance(query_point,*it); + Point_with_transformed_distance *NN_Candidate= + new Point_with_transformed_distance(*it,distance_to_query_point); + Item_PriorityQueue.push(NN_Candidate); + } + // old top of PriorityQueue has been processed, + // hence update rd + + if (!(PriorityQueue.empty())) { + rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); + next_neighbour_found = + (multiplication_factor*rd > + Item_PriorityQueue.top()->second); + } + else // priority queue empty => last neighbour found + { + next_neighbour_found=true; + } + + number_of_neighbours_computed++; + } + } // next_neighbour_found or priority queue is empty + // in the latter case also the item priority quee is empty + } + + + void + Compute_the_next_furthest_neighbour() + { + // compute the next item + bool next_neighbour_found=false; + if (!(Item_PriorityQueue.empty())) { + next_neighbour_found= + (rd < multiplication_factor*Item_PriorityQueue.top()->second); + } + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point); + // otherwise browse the tree further + while ((!next_neighbour_found) && (!PriorityQueue.empty())) { + Node_with_distance* The_node_top=PriorityQueue.top(); + Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top); + dists = CGAL::cpp11::get<2>(*The_node_top); + PriorityQueue.pop(); + delete The_node_top; + FT copy_rd=rd; + while (!(N->is_leaf())) { // compute new distance + typename Tree::Internal_node_const_handle node = + static_cast<typename Tree::Internal_node_const_handle>(N); + number_of_internal_nodes_visited++; + int new_cut_dim=node->cutting_dimension(); + FT new_rd,dst = dists[new_cut_dim]; + FT val = *(query_point_it + new_cut_dim); + FT diff1 = val - node->upper_low_value(); + FT diff2 = val - node->lower_high_value(); + if (diff1 + diff2 < FT(0.0)) { + diff1 = val - node->upper_high_value(); + new_rd= + Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); + Node_with_distance *Lower_Child = + new Node_with_distance(node->lower(), copy_rd, dists); + PriorityQueue.push(Lower_Child); + N=node->upper(); + dists[new_cut_dim] = diff1; + copy_rd=new_rd; + + } + else { // compute new distance + diff2 = val - node->lower_low_value(); + new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); + Node_with_distance *Upper_Child = + new Node_with_distance(node->upper(), copy_rd, dists); + PriorityQueue.push(Upper_Child); + N=node->lower(); + dists[new_cut_dim] = diff2; + copy_rd=new_rd; + } + } + // n is a leaf + typename Tree::Leaf_node_const_handle node = + static_cast<typename Tree::Leaf_node_const_handle>(N); + number_of_leaf_nodes_visited++; + if (node->size() > 0) { + for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { + number_of_items_visited++; + FT distance_to_query_point= + Orthogonal_distance_instance.transformed_distance(query_point,*it); + Point_with_transformed_distance *NN_Candidate= + new Point_with_transformed_distance(*it,distance_to_query_point); + Item_PriorityQueue.push(NN_Candidate); + } + // old top of PriorityQueue has been processed, + // hence update rd + + if (!(PriorityQueue.empty())) { + rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); + next_neighbour_found = + (multiplication_factor*rd < + Item_PriorityQueue.top()->second); + } + else // priority queue empty => last neighbour found + { + next_neighbour_found=true; + } + + number_of_neighbours_computed++; + } + } // next_neighbour_found or priority queue is empty + // in the latter case also the item priority quee is empty + } + }; // class Iterator_implementaion + + + + + + + + + + public: + class iterator; + typedef iterator const_iterator; + + // constructor + Orthogonal_incremental_neighbor_search(const Tree& tree, + const Query_item& q, FT Eps = FT(0.0), + bool search_nearest=true, const Distance& tr=Distance()) + : m_tree(tree),m_query(q),m_dist(tr),m_Eps(Eps),m_search_nearest(search_nearest) + {} + + iterator + begin() const + { + return iterator(m_tree,m_query,m_dist,m_Eps,m_search_nearest); + } + + iterator + end() const + { + return iterator(); + } + + std::ostream& + statistics(std::ostream& s) + { + begin()->statistics(s); + return s; + } + + + + + class iterator { + + public: + + typedef std::input_iterator_tag iterator_category; + typedef Point_with_transformed_distance value_type; + typedef Point_with_transformed_distance* pointer; + typedef const Point_with_transformed_distance& reference; + typedef std::size_t size_type; + typedef std::ptrdiff_t difference_type; + typedef int distance_type; + + //class Iterator_implementation; + Iterator_implementation *Ptr_implementation; + + + public: + + // default constructor + iterator() + : Ptr_implementation(0) + {} + + int + the_number_of_items_visited() + { + return Ptr_implementation->number_of_items_visited; + } + + // constructor + iterator(const Tree& tree,const Query_item& q, const Distance& tr=Distance(), FT eps=FT(0.0), + bool search_nearest=true) + : Ptr_implementation(new Iterator_implementation(tree, q, tr, eps, search_nearest)) + {} + + // copy constructor + iterator(const iterator& Iter) + { + Ptr_implementation = Iter.Ptr_implementation; + if (Ptr_implementation != 0) Ptr_implementation->reference_count++; + } + + iterator& operator=(const iterator& Iter) + { + if (Ptr_implementation != Iter.Ptr_implementation){ + if (Ptr_implementation != 0 && --(Ptr_implementation->reference_count)==0) { + delete Ptr_implementation; + } + Ptr_implementation = Iter.Ptr_implementation; + if (Ptr_implementation != 0) Ptr_implementation->reference_count++; + } + return *this; + } + + + const Point_with_transformed_distance& + operator* () const + { + return *(*Ptr_implementation); + } + + // -> operator + const Point_with_transformed_distance* + operator-> () const + { + return &*(*Ptr_implementation); + } + + // prefix operator + iterator& + operator++() + { + ++(*Ptr_implementation); + return *this; + } + + // postfix operator + Object_wrapper<Point_with_transformed_distance> + operator++(int) + { + return (*Ptr_implementation)++; + } + + + bool + operator==(const iterator& It) const + { + if ( + ((Ptr_implementation == 0) || + Ptr_implementation->Item_PriorityQueue.empty()) && + ((It.Ptr_implementation == 0) || + It.Ptr_implementation->Item_PriorityQueue.empty()) + ) + return true; + // else + return (Ptr_implementation == It.Ptr_implementation); + } + + bool + operator!=(const iterator& It) const + { + return !(*this == It); + } + + std::ostream& + statistics (std::ostream& s) + { + Ptr_implementation->statistics(s); + return s; + } + + ~iterator() + { + if (Ptr_implementation != 0) { + Ptr_implementation->reference_count--; + if (Ptr_implementation->reference_count==0) { + delete Ptr_implementation; + Ptr_implementation = 0; + } + } + } + + + }; // class iterator + + //data members + const Tree& m_tree; + Query_item m_query; + Distance m_dist; + FT m_Eps; + bool m_search_nearest; + }; // class + + template <class Traits, class Query_item, class Distance> + void swap (typename Orthogonal_incremental_neighbor_search<Traits, + Query_item, Distance>::iterator& x, + typename Orthogonal_incremental_neighbor_search<Traits, + Query_item, Distance>::iterator& y) + { + typename Orthogonal_incremental_neighbor_search<Traits, + Query_item, Distance>::iterator::Iterator_implementation + *tmp = x.Ptr_implementation; + x.Ptr_implementation = y.Ptr_implementation; + y.Ptr_implementation = tmp; + } + +} // namespace CGAL + +#endif // CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH_H diff --git a/src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt b/src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt new file mode 100644 index 00000000..5b9581a0 --- /dev/null +++ b/src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt @@ -0,0 +1,82 @@ +CGAL/Regular_triangulation_traits_adapter.h +CGAL/Triangulation_ds_vertex.h +CGAL/Triangulation_data_structure.h +CGAL/transforming_pair_iterator.h +CGAL/NewKernel_d/static_int.h +CGAL/NewKernel_d/Cartesian_LA_functors.h +CGAL/NewKernel_d/Cartesian_change_FT.h +CGAL/NewKernel_d/Wrapper/Vector_d.h +CGAL/NewKernel_d/Wrapper/Hyperplane_d.h +CGAL/NewKernel_d/Wrapper/Ref_count_obj.h +CGAL/NewKernel_d/Wrapper/Cartesian_wrap.h +CGAL/NewKernel_d/Wrapper/Point_d.h +CGAL/NewKernel_d/Wrapper/Segment_d.h +CGAL/NewKernel_d/Wrapper/Weighted_point_d.h +CGAL/NewKernel_d/Wrapper/Sphere_d.h +CGAL/NewKernel_d/Cartesian_per_dimension.h +CGAL/NewKernel_d/Kernel_object_converter.h +CGAL/NewKernel_d/KernelD_converter.h +CGAL/NewKernel_d/Vector/sse2.h +CGAL/NewKernel_d/Vector/avx4.h +CGAL/NewKernel_d/Vector/determinant_of_vectors_small_dim_internal.h +CGAL/NewKernel_d/Vector/determinant_of_iterator_to_points_from_points.h +CGAL/NewKernel_d/Vector/determinant_of_points_from_vectors.h +CGAL/NewKernel_d/Vector/array.h +CGAL/NewKernel_d/Vector/determinant_of_iterator_to_points_from_iterator_to_vectors.h +CGAL/NewKernel_d/Vector/determinant_of_iterator_to_vectors_from_vectors.h +CGAL/NewKernel_d/Vector/determinant_of_vectors_small_dim.h +CGAL/NewKernel_d/Vector/vector.h +CGAL/NewKernel_d/Vector/v2int.h +CGAL/NewKernel_d/Vector/mix.h +CGAL/NewKernel_d/Cartesian_static_filters.h +CGAL/NewKernel_d/Cartesian_LA_base.h +CGAL/NewKernel_d/Lazy_cartesian.h +CGAL/NewKernel_d/Coaffine.h +CGAL/NewKernel_d/store_kernel.h +CGAL/NewKernel_d/Dimension_base.h +CGAL/NewKernel_d/Kernel_3_interface.h +CGAL/NewKernel_d/Cartesian_complete.h +CGAL/NewKernel_d/Cartesian_base.h +CGAL/NewKernel_d/Cartesian_filter_K.h +CGAL/NewKernel_d/functor_tags.h +CGAL/NewKernel_d/Filtered_predicate2.h +CGAL/NewKernel_d/functor_properties.h +CGAL/NewKernel_d/Define_kernel_types.h +CGAL/NewKernel_d/LA_eigen/LA.h +CGAL/NewKernel_d/LA_eigen/constructors.h +CGAL/NewKernel_d/Types/Aff_transformation.h +CGAL/NewKernel_d/Types/Sphere.h +CGAL/NewKernel_d/Types/Hyperplane.h +CGAL/NewKernel_d/Types/Line.h +CGAL/NewKernel_d/Types/Ray.h +CGAL/NewKernel_d/Types/Iso_box.h +CGAL/NewKernel_d/Types/Weighted_point.h +CGAL/NewKernel_d/Types/Segment.h +CGAL/NewKernel_d/Kernel_d_interface.h +CGAL/NewKernel_d/utils.h +CGAL/NewKernel_d/Kernel_2_interface.h +CGAL/NewKernel_d/Cartesian_filter_NT.h +CGAL/NewKernel_d/function_objects_cartesian.h +CGAL/Convex_hull.h +CGAL/Triangulation_ds_full_cell.h +CGAL/Regular_triangulation.h +CGAL/Epick_d.h +CGAL/transforming_iterator.h +CGAL/iterator_from_indices.h +CGAL/Delaunay_triangulation.h +CGAL/IO/Triangulation_off_ostream.h +CGAL/typeset.h +CGAL/Triangulation_full_cell.h +CGAL/Triangulation.h +CGAL/internal/Static_or_dynamic_array.h +CGAL/internal/Combination_enumerator.h +CGAL/internal/Triangulation/utilities.h +CGAL/internal/Triangulation/Triangulation_ds_iterators.h +CGAL/internal/Triangulation/Dummy_TDS.h +CGAL/argument_swaps.h +CGAL/Epeck_d.h +CGAL/determinant_of_vectors.h +CGAL/TDS_full_cell_default_storage_policy.h +CGAL/TDS_full_cell_mirror_storage_policy.h +CGAL/Triangulation_face.h +CGAL/Triangulation_vertex.h |
