From e9dd788438bff7289ddff1e0ade2de0f031a2f9b Mon Sep 17 00:00:00 2001 From: fgodi Date: Mon, 28 Nov 2016 09:57:43 +0000 Subject: Bug fix git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1789 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: f2c3ef0dab1a0d80cfcf0018da83dd2b6b9a9ef1 --- .../include/gudhi/CGAL/Kd_tree.h | 14 +- .../include/gudhi/CGAL/Kd_tree_node.h | 2 +- .../include/gudhi/CGAL/Splitters.h | 313 --------------------- .../include/gudhi/Neighbors_finder.h | 63 +++-- .../include/gudhi/Persistence_diagrams_graph.h | 2 +- .../include/gudhi/Planar_neighbors_finder.h | 156 ---------- 6 files changed, 48 insertions(+), 502 deletions(-) delete mode 100644 src/Bottleneck_distance/include/gudhi/CGAL/Splitters.h delete mode 100644 src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h (limited to 'src/Bottleneck_distance/include') diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h index dbdf5259..f085b0da 100644 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h +++ b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h @@ -392,24 +392,26 @@ public: bool success = remove_(p, 0, false, 0, false, root(), equal_to_p); CGAL_assertion(success); - removed_ |= 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 islower, - Internal_node_handle parent, bool islower2, + 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, islower2, newparent, true, newparent->lower(), equal_to_p)) + 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, islower2, newparent, false, newparent->upper(), equal_to_p); + return remove_(p, parent, islower, newparent, false, newparent->upper(), equal_to_p); CGAL_assertion(false); // Point was not found } @@ -431,7 +433,7 @@ private: return false; } else if (grandparent) { Node_handle brother = islower ? parent->upper() : parent->lower(); - if (islower2) + if (parent_islower) grandparent->set_lower(brother); else grandparent->set_upper(brother); diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h index 49b0c022..909ee260 100644 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h +++ b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h @@ -21,7 +21,7 @@ #ifndef CGAL_KD_TREE_NODE_H #define CGAL_KD_TREE_NODE_H -#include "Splitters.h" +#include "CGAL/Splitters.h" #include #include diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Splitters.h b/src/Bottleneck_distance/include/gudhi/CGAL/Splitters.h deleted file mode 100644 index e58a593d..00000000 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Splitters.h +++ /dev/null @@ -1,313 +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 () - -// Defines rules used for constructing a split node. That is, it implements, -// in several ways, the concept -// Boxtree_splitter. - -#ifndef CGAL_SPLITTERS_H -#define CGAL_SPLITTERS_H - -#include -#include - -namespace CGAL { - - template - class Splitter_base { - private: - unsigned int the_bucket_size; - FT the_aspect_ratio; - - public: - //default bucket_size should be 10 - Splitter_base(unsigned int bucket_size = 10, - FT aspect_ratio = FT(3)) - : the_bucket_size(bucket_size), - the_aspect_ratio(aspect_ratio) - {} - - FT - aspect_ratio() const - { - return the_aspect_ratio; - } - - unsigned int - bucket_size() const - { - return the_bucket_size; - } - - }; - - - template > - class Median_of_max_spread - : public Splitter_base - { - - typedef Splitter_base Base; - public: - typedef typename SearchTraits::FT FT; - typedef Point_container Container; - typedef Separator_ Separator; - - Median_of_max_spread() - : Base() - {} - - Median_of_max_spread(unsigned int bucket_size) - : Base(bucket_size) - {} - - void - operator() (Separator& sep, - Container& c0, - Container& c1) const - { - sep=Separator(c0.max_tight_span_coord(),FT(0)); - sep.set_cutting_value(c0.median(sep.cutting_dimension())); - c0.split(c1,sep,true); - } - }; - - template > - class Fair - : public Splitter_base - { - - typedef Splitter_base Base; - public: - typedef typename SearchTraits::FT FT; - typedef Point_container Container; - typedef Separator_ Separator; - - Fair() - : Base() - {} - - Fair(unsigned int bucket_size, - FT aspect_ratio=FT(3)) - : Base(bucket_size, aspect_ratio) - {} - - void - operator()(Separator& sep, Container& c0, Container& c1) const - { - // find legal cut with max spread - sep=Separator(c0.max_tight_span_coord_balanced(this->aspect_ratio()), - FT(0)); - sep.set_cutting_value(c0.balanced_fair(sep.cutting_dimension(), - this->aspect_ratio())); - c0.split(c1,sep); - } - }; - - template > - class Sliding_fair - : public Splitter_base - { - - typedef Splitter_base Base; - - public: - typedef typename SearchTraits::FT FT; - typedef Point_container Container; - typedef Separator_ Separator; - - Sliding_fair() - : Base() - {} - - Sliding_fair(unsigned int bucket_size, - FT aspect_ratio=FT(3)) - : Base(bucket_size, aspect_ratio) - {} - - void - operator() (Separator& sep, Container& c0, Container& c1) const - { - // find legal cut with max spread - - sep = Separator(c0.max_tight_span_coord_balanced(this->aspect_ratio()), - FT(0)); - - sep.set_cutting_value(c0.balanced_sliding_fair(sep.cutting_dimension(), - this->aspect_ratio())); - c0.split(c1,sep,true); - } - }; - - - template > - class Sliding_midpoint - : public Splitter_base - { - - typedef Splitter_base Base; - - public: - typedef typename SearchTraits::FT FT; - typedef Point_container Container; - typedef Separator_ Separator; - - Sliding_midpoint() - : Base() - {} - - Sliding_midpoint(unsigned int bucket_size) - : Base(bucket_size) - {} - - void - operator()(Separator& sep, Container& c0, Container& c1) const - { - CGAL_assertion(c0.is_valid()); - CGAL_assertion(c1.is_valid()); - int cutdim = c0.max_span_coord(); - - //Bugfix: avoid linear tree in degenerated cases - if(c0.tight_bounding_box().min_coord(cutdim) != c0.tight_bounding_box().max_coord(cutdim)){ - sep = Separator(cutdim, - (c0.max_span_upper() + c0.max_span_lower())/FT(2)); - } - else{ - cutdim = c0.max_tight_span_coord(); - sep = Separator(cutdim, - (c0.max_tight_span_upper() + c0.max_tight_span_lower())/FT(2)); - } - - FT max_span_lower = - c0.tight_bounding_box().min_coord(cutdim); - - CGAL_assertion(max_span_lower >= c0.bounding_box().min_coord(cutdim)); - FT max_span_upper = - c0.tight_bounding_box().max_coord(cutdim); - - CGAL_assertion(max_span_upper <= c0.bounding_box().max_coord(cutdim)); - if (max_span_upper <= sep.cutting_value()) { - sep.set_cutting_value(max_span_upper); - }; - if (max_span_lower >= sep.cutting_value()) { - sep.set_cutting_value(max_span_lower); - }; - c0.split(c1,sep,true); - } - }; - - template > - class Median_of_rectangle - : public Splitter_base - { - - typedef Splitter_base Base; - - public: - typedef typename SearchTraits::FT FT; - typedef Point_container Container; - typedef Separator_ Separator; - - - Median_of_rectangle() - : Base() - {} - - Median_of_rectangle(unsigned int bucket_size) - : Base(bucket_size) - {} - - void - operator() (Separator& sep, Container& c0, Container& c1) const - { - sep = Separator(c0.max_span_coord(),FT(0)); - sep.set_cutting_value(c0.median(sep.cutting_dimension())); - c0.split(c1,sep,true); - } - }; - - template > - class Midpoint_of_max_spread - : public Splitter_base - { - typedef Splitter_base Base; - - public: - typedef typename SearchTraits::FT FT; - typedef Point_container Container; - typedef Separator_ Separator; - - - Midpoint_of_max_spread() - : Base() - {} - - Midpoint_of_max_spread(unsigned int bucket_size) - : Base(bucket_size) - {} - - void - operator()(Separator& sep, Container& c0, Container& c1) const - { - sep = Separator(c0.max_tight_span_coord(), - (c0.max_tight_span_upper() + c0.max_tight_span_lower())/FT(2)); - c0.split(c1,sep); - } - }; - - template > - class Midpoint_of_rectangle - : public Splitter_base - { - - typedef Splitter_base Base; - public: - typedef typename SearchTraits::FT FT; - typedef Point_container Container; - typedef Separator_ Separator; - - - Midpoint_of_rectangle() - : Base() - {} - - Midpoint_of_rectangle(unsigned int bucket_size) - : Base(bucket_size) - {} - - void - operator()(Separator& sep, Container& c0, Container& c1) const - { - sep = Separator(c0.max_span_coord(), - (c0.max_span_upper() + c0.max_span_lower())/FT(2)); - c0.split(c1,sep); - } - - }; - -} // namespace CGAL -#endif // CGAL_SPLITTERS diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h index 20b47d0b..f2d9abb2 100644 --- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -23,8 +23,18 @@ #ifndef NEIGHBORS_FINDER_H_ #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_incremental_neighbor_search.h" + +#include +#include + +#include +#include + #include -#include namespace Gudhi { @@ -38,6 +48,13 @@ namespace bottleneck_distance { * \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_incremental_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(double r); @@ -50,10 +67,8 @@ public: private: const double r; - Planar_neighbors_finder planar_neighbors_f; + Kd_tree kd_t; std::unordered_set projections_f; - void remove(int v_point_index); - bool contains(int v_point_index); }; /** \internal \brief data structure used to find any point (including projections) in V near to a query point from U @@ -80,46 +95,44 @@ private: }; inline Neighbors_finder::Neighbors_finder(double r) : - r(r), planar_neighbors_f(r), projections_f() { } + 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 - planar_neighbors_f.add(v_point_index); -} - -inline void Neighbors_finder::remove(int v_point_index) { - if(v_point_index == null_point_index()) - return; - if (G::on_the_v_diagonal(v_point_index)) - projections_f.erase(v_point_index); - else - planar_neighbors_f.remove(v_point_index); -} - -inline bool Neighbors_finder::contains(int v_point_index) { - return planar_neighbors_f.contains(v_point_index) || (projections_f.count(v_point_index)>0); + 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()) + if (G::on_the_u_diagonal(u_point_index) && !projections_f.empty()){ //Any pair of projection is at distance 0 tmp = *projections_f.cbegin(); - else if (contains(c) && (G::distance(u_point_index, c) <= r)) + 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; - else + projections_f.erase(tmp); + } + else{ //Is the query point near to a V point in the plane ? - tmp = planar_neighbors_f.pull_near(u_point_index); - remove(tmp); + Internal_point u_point = G::get_u_point(u_point_index); + std::vector w = {1., 1.}; + K_neighbor_search search(kd_t, u_point, 0., true, Distance(0, 2, w)); + 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(planar_neighbors_f.pull_all_near(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); diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h index 8242ce2b..fd17bf4a 100644 --- a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h +++ b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h @@ -65,7 +65,7 @@ private: static Internal_point get_u_point(int u_point_index); static Internal_point get_v_point(int v_point_index); - friend class Planar_neighbors_finder; + friend class Neighbors_finder; }; /** \internal \typedef \brief Shorter alias */ diff --git a/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h deleted file mode 100644 index c1f4d908..00000000 --- a/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h +++ /dev/null @@ -1,156 +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: Francois Godi - * - * Copyright (C) 2015 INRIA (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 . - */ - -#ifndef PLANAR_NEIGHBORS_FINDER_H_ -#define PLANAR_NEIGHBORS_FINDER_H_ - -// Inclusion order is important for CGAL patch -#include "CGAL/Kd_tree_node.h" -#include "CGAL/Kd_tree.h" -#include "CGAL/Orthogonal_incremental_neighbor_search.h" - -#include -#include - -#include -#include - - -namespace Gudhi { - -namespace bottleneck_distance { - -/** \internal \brief Structure used to find any point in V near (according to the planar distance) to a query point from U. - * - * V points have to be added manually using their index and before the first remove/pull. A neighbor pulled is automatically removed. but we can also - * remove points manually using their index. - * - * \ingroup bottleneck_distance - */ -class Naive_pnf { -public: - /** \internal \brief Constructor taking the near distance definition as parameter. */ - Naive_pnf(double r_); - /** \internal \brief A point added will be possibly pulled. */ - void add(int v_point_index); - /** \internal \brief A point manually removed will no longer be possibly pulled. */ - void remove(int v_point_index); - /** \internal \brief Can the point given as parameter be returned ? */ - bool contains(int v_point_index) const; - /** \internal \brief Provide 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 Provide and remove all the V points near to the U point given as parameter. */ - std::vector pull_all_near(int u_point_index); - -private: - double r; - std::pair get_v_key(int v_point_index) const; - std::multimap,int> grid; -}; - -class Planar_neighbors_finder { - - typedef CGAL::Dimension_tag<2> D; - typedef CGAL::Search_traits Traits; - typedef CGAL::Weighted_Minkowski_distance Distance; - typedef CGAL::Orthogonal_incremental_neighbor_search K_neighbor_search; - typedef K_neighbor_search::Tree Kd_tree; - - -public: - /** \internal \brief Constructor taking the near distance definition as parameter. */ - Planar_neighbors_finder(double r_); - /** \internal \brief A point added will be possibly pulled. */ - void add(int v_point_index); - /** \internal \brief A point manually removed will no longer be possibly pulled. */ - void remove(int v_point_index); - /** \internal \brief Can the point given as parameter be returned ? */ - bool contains(int v_point_index) const; - /** \internal \brief Provide 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 Provide and remove all the V points near to the U point given as parameter. */ - virtual std::vector pull_all_near(int u_point_index); - -private: - double r; - std::set contents; - Kd_tree kd_t; -}; - - -/** \internal \brief Constructor taking the near distance definition as parameter. */ -inline Planar_neighbors_finder::Planar_neighbors_finder(double r_) - : r(r_), contents(), kd_t() {} - - -/** \internal \brief A point added will be possibly pulled. */ -inline void Planar_neighbors_finder::add(int v_point_index){ - if(v_point_index == null_point_index()) - return; - contents.insert(v_point_index); - kd_t.insert(G::get_v_point(v_point_index)); -} - -/** \internal \brief A point manually removed will no longer be possibly pulled. */ -inline void Planar_neighbors_finder::remove(int v_point_index){ - contents.erase(v_point_index); - kd_t.remove(G::get_v_point(v_point_index)); -} - -/** \internal \brief Can the point given as parameter be returned ? */ -inline bool Planar_neighbors_finder::contains(int v_point_index) const{ - if(v_point_index == null_point_index()) - return false; - return contents.count(v_point_index)>0; -} - -/** \internal \brief Provide and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ -inline int Planar_neighbors_finder::pull_near(int u_point_index){ - Internal_point u_point = G::get_u_point(u_point_index); - std::vector w = {1., 1.}; - K_neighbor_search search(kd_t, u_point, 0., true, Distance(0, 2, w)); - auto it = search.begin(); - if(it==search.end() || G::distance(u_point_index, it->first.point_index) > r) - return null_point_index(); - int tmp = it->first.point_index; - if(!contains(tmp)) - std::cout << "!! A kd_tree returns a point (Point_index:" << tmp << ") previously removed !!" << std::endl; - remove(tmp); - return tmp; -} - -inline std::vector Planar_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; -} - - -} // namespace bottleneck_distance - -} // namespace Gudhi - -#endif // PLANAR_NEIGHBORS_FINDER_H_ -- cgit v1.2.3