From 64dd2a2b0eec1374ed23ca079f86b312125d03f7 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Wed, 4 Jan 2017 08:39:18 +0000 Subject: Move bottleneck_chrono in benchmark Add test in cmake for basic example Move CGAL gudhi patches for bottleneck in dedicated directory Fix cpplint syntax issue git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1920 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 24671facf791de93dc6fd94bb39ca7362bb22959 --- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 93 +-- .../include/gudhi/CGAL/Kd_tree.h | 582 ------------------- .../include/gudhi/CGAL/Kd_tree_node.h | 586 ------------------- .../CGAL/Orthogonal_incremental_neighbor_search.h | 621 --------------------- .../include/gudhi/Graph_matching.h | 220 ++++---- .../include/gudhi/Internal_point.h | 69 ++- .../include/gudhi/Neighbors_finder.h | 172 +++--- .../include/gudhi/Persistence_graph.h | 191 ++++--- 8 files changed, 384 insertions(+), 2150 deletions(-) delete mode 100644 src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h delete mode 100644 src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h delete mode 100644 src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h (limited to 'src/Bottleneck_distance/include') diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index 42a0d444..2b7e4767 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (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 @@ -31,62 +31,65 @@ 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; - } + 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.; + } + return (b_lower_bound + b_upper_bound) / 2.; } double bottleneck_distance_exact(Persistence_graph& g) { - std::vector 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((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; - } + std::vector 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 ((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); + } + 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. + * 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 -double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=std::numeric_limits::min()) { - Persistence_graph g(diag1, diag2, e); - if(g.bottleneck_alive() == std::numeric_limits::infinity()) - return std::numeric_limits::infinity(); - return std::max(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e)); +double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, + double e = std::numeric_limits::min()) { + Persistence_graph g(diag1, diag2, e); + if (g.bottleneck_alive() == std::numeric_limits::infinity()) + return std::numeric_limits::infinity(); + return std::max(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e)); } } // namespace persistence_diagram diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h deleted file mode 100644 index f085b0da..00000000 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h +++ /dev/null @@ -1,582 +0,0 @@ -// 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 (), -// : Waqar Khan - -#ifndef CGAL_KD_TREE_H -#define CGAL_KD_TREE_H - -#include "Kd_tree_node.h" - -#include -#include -#include - -#include -#include -#include - - -#include -#include -#include - -#ifdef CGAL_HAS_THREADS -#include -#endif - -namespace CGAL { - -//template , class UseExtendedNode = Tag_true > -template , 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 Node; - typedef Kd_tree_leaf_node Leaf_node; - typedef Kd_tree_internal_node Internal_node; - typedef Kd_tree Tree; - typedef Kd_tree 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_iterator Point_d_iterator; - typedef typename std::vector::const_iterator Point_d_const_iterator; - typedef typename Splitter::Separator Separator; - typedef typename std::vector::const_iterator iterator; - typedef typename std::vector::const_iterator const_iterator; - - typedef typename std::vector::size_type size_type; - - typedef typename internal::Get_dimension_tag::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_nodes; - std::deque leaf_nodes; -#else - boost::container::deque internal_nodes; - boost::container::deque leaf_nodes; -#endif - - Node_handle tree_root; - - Kd_tree_rectangle* bbox; - std::vector 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 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(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 - 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(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(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 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(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 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 - 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 - void - remove(const Point_d& p, Equal const& equal_to_p) - { -#if 0 - // This code could have quadratic runtime. - if (!is_built()) { - std::vector::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 - 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(node); - // FIXME: This should be if(xcutting_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(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 - OutputIterator - search(OutputIterator it, const FuzzyQueryItem& q) const - { - if(! pts.empty()){ - - if(! is_built()){ - const_build(); - } - Kd_tree_rectangle b(*bbox); - return tree_root->search(it,q,b); - } - return it; - } - - - template - boost::optional - search_any_point(const FuzzyQueryItem& q) const - { - if(! pts.empty()){ - - if(! is_built()){ - const_build(); - } - Kd_tree_rectangle 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& - 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/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h deleted file mode 100644 index 909ee260..00000000 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h +++ /dev/null @@ -1,586 +0,0 @@ -// 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 () - -#ifndef CGAL_KD_TREE_NODE_H -#define CGAL_KD_TREE_NODE_H - -#include "CGAL/Splitters.h" - -#include -#include - -namespace CGAL { - - template - class Kd_tree; - - template < class TreeTraits, class Splitter, class UseExtendedNode > - class Kd_tree_node { - - friend class Kd_tree; - - typedef typename Kd_tree::Node_handle Node_handle; - typedef typename Kd_tree::Node_const_handle Node_const_handle; - typedef typename Kd_tree::Internal_node_handle Internal_node_handle; - typedef typename Kd_tree::Internal_node_const_handle Internal_node_const_handle; - typedef typename Kd_tree::Leaf_node_handle Leaf_node_handle; - typedef typename Kd_tree::Leaf_node_const_handle Leaf_node_const_handle; - typedef typename TreeTraits::Point_d Point_d; - - typedef typename TreeTraits::FT FT; - typedef typename Kd_tree::Separator Separator; - typedef typename Kd_tree::Point_d_iterator Point_d_iterator; - typedef typename Kd_tree::iterator iterator; - typedef typename Kd_tree::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(this); - return node->size(); - } - else { - Internal_node_const_handle node = - static_cast(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(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(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 - OutputIterator - tree_items(OutputIterator it) const { - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(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(this); - it=node->lower()->tree_items(it); - it=node->upper()->tree_items(it); - } - return it; - } - - - boost::optional - any_tree_item() const { - boost::optional result = boost::none; - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(this); - if (node->size()>0){ - return boost::make_optional(*(node->begin())); - } - } - else { - Internal_node_const_handle node = - static_cast(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(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(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 - OutputIterator - search(OutputIterator it, const FuzzyQueryItem& q, - Kd_tree_rectangle& b) const - { - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(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(this); - // after splitting b denotes the lower part of b - Kd_tree_rectangle 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 - boost::optional - search_any_point(const FuzzyQueryItem& q, - Kd_tree_rectangle& b) const - { - boost::optional result = boost::none; - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(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(this); - // after splitting b denotes the lower part of b - Kd_tree_rectangle 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; - - typedef typename Kd_tree::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; - - typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; - typedef typename Kd_tree::Node_handle Node_handle; - typedef typename Kd_tree::Node_const_handle Node_const_handle; - - typedef typename TreeTraits::FT FT; - typedef typename Kd_tree::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 : public Kd_tree_node< TreeTraits, Splitter, Tag_false >{ - - friend class Kd_tree; - - typedef Kd_tree_node< TreeTraits, Splitter, Tag_false> Base; - typedef typename Kd_tree::Node_handle Node_handle; - typedef typename Kd_tree::Node_const_handle Node_const_handle; - - typedef typename TreeTraits::FT FT; - typedef typename Kd_tree::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/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h b/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h deleted file mode 100644 index dbe707ed..00000000 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h +++ /dev/null @@ -1,621 +0,0 @@ -// 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 () - -#ifndef CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH -#define CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH - -#include "CGAL/Kd_tree.h" - -#include -#include -#include -#include -#include -#include - -namespace CGAL { - - template ::type, - class Splitter_ = Sliding_midpoint, - class Tree_= Kd_tree > - 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_with_transformed_distance; - typedef CGAL::cpp11::tuple > Node_with_distance; - typedef std::vector Node_with_distance_vector; - typedef std::vector Point_with_transformed_distance_vector; - - template - 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 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 PriorityQueue; - - public: - std::priority_queue 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(std::distance(ccci(q), ccci(q,0))); - - dists.resize(dim); - for(int i=0 ; i(*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 - operator++(int) - { - Object_wrapper 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(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(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(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(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 - 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 - void swap (typename Orthogonal_incremental_neighbor_search::iterator& x, - typename Orthogonal_incremental_neighbor_search::iterator& y) - { - typename Orthogonal_incremental_neighbor_search::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/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h index e9f455d7..253c89b4 100644 --- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (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 @@ -34,141 +34,141 @@ namespace persistence_diagram { * \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 v_to_u; - /** \internal \brief All the unmatched points in U. */ - std::list 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 & path); + 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 v_to_u; + /** \internal \brief All the unmatched points in U. */ + std::list 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 & 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); + 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; + 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(); + 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 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; + 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 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; + 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 path; - path.emplace_back(u_start_index); - do { - if (static_cast(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(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; + // V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices. + std::vector path; + path.emplace_back(u_start_index); + do { + if (static_cast (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 (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 u_vertices(unmatched_in_u); - std::list 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 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(); + std::list u_vertices(unmatched_in_u); + std::list 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 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); + } } - return layered_nf; + // 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& 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; - } + 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; + } } diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h index 70342d0c..0b2d26fe 100644 --- a/src/Bottleneck_distance/include/gudhi/Internal_point.h +++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (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 @@ -32,39 +32,60 @@ 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); } + 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; + 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; } + 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 index 792925b7..96ece360 100644 --- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (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 @@ -24,10 +24,9 @@ #define NEIGHBORS_FINDER_H_ // Inclusion order is important for CGAL patch -#include "CGAL/Kd_tree_node.h" -#include "CGAL/Kd_tree.h" -#include "CGAL/Orthogonal_k_neighbor_search.h" - +#include +#include +#include #include #include @@ -43,123 +42,126 @@ 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. + * 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 Traits; - typedef CGAL::Weighted_Minkowski_distance Distance; - typedef CGAL::Orthogonal_k_neighbor_search 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 pull_all_near(int u_point_index); - -private: - const Persistence_graph& g; - const double r; - Kd_tree kd_t; - std::unordered_set projections_f; + typedef CGAL::Dimension_tag<2> D; + typedef CGAL::Search_traits Traits; + typedef CGAL::Weighted_Minkowski_distance Distance; + typedef CGAL::Orthogonal_k_neighbor_search 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 pull_all_near(int u_point_index); + + private: + const Persistence_graph& g; + const double r; + Kd_tree kd_t; + std::unordered_set 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. + * 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> 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> 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)); + 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 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; + 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 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 Neighbors_finder::pull_all_near(int u_point_index) { - std::vector 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; + std::vector 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(new Neighbors_finder(g, r))); - neighbors_finder.at(vlayer)->add(v_point_index); + for (int l = neighbors_finder.size(); l <= vlayer; l++) + neighbors_finder.emplace_back(std::unique_ptr(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 (neighbors_finder.size()) <= vlayer) - return null_point_index(); - return neighbors_finder.at(vlayer)->pull_near(u_point_index); + if (static_cast (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(neighbors_finder.size()); + return static_cast (neighbors_finder.size()); } } // namespace persistence_diagram diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h index 45a4d586..3a4a5fec 100644 --- a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h +++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (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 @@ -31,147 +31,144 @@ 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 - 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 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 u; - std::vector v; - double b_alive; + public: + /** \internal \brief Constructor taking 2 Persistence_Diagrams (concept) as parameters. */ + template + 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 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 u; + std::vector v; + double b_alive; }; template Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) - : u(), v(), b_alive(0.) -{ - std::vector u_alive; - std::vector v_alive; - for (auto it = std::begin(diag1); it != std::end(diag1); ++it){ - if(std::get<1>(*it) == std::numeric_limits::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::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::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)); + : u(), v(), b_alive(0.) { + std::vector u_alive; + std::vector v_alive; + for (auto it = std::begin(diag1); it != std::end(diag1); ++it) { + if (std::get<1>(*it) == std::numeric_limits::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::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::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 (u.size()); + return u_point_index >= static_cast (u.size()); } inline bool Persistence_graph::on_the_v_diagonal(int v_point_index) const { - return v_point_index >= static_cast (v.size()); + return v_point_index >= static_cast (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 (v.size()) : v_point_index + static_cast (u.size()); + return on_the_v_diagonal(v_point_index) ? + v_point_index - static_cast (v.size()) : v_point_index + static_cast (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 (u.size()) : u_point_index + static_cast (v.size()); + return on_the_u_diagonal(u_point_index) ? + u_point_index - static_cast (u.size()) : u_point_index + static_cast (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())); + 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 (u.size() + v.size()); + return static_cast (u.size() + v.size()); } -inline double Persistence_graph::bottleneck_alive() const{ - return b_alive; +inline double Persistence_graph::bottleneck_alive() const { + return b_alive; } inline std::vector Persistence_graph::sorted_distances() const { - std::vector 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; + std::vector 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); + 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); + 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; + 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 -- cgit v1.2.3