summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/Cech_complex/benchmark/CMakeLists.txt18
-rw-r--r--src/Cech_complex/benchmark/cech_complex_benchmark.cpp27
-rw-r--r--src/Cech_complex/concept/SimplicialComplexForCech.h4
-rw-r--r--src/Cech_complex/example/CMakeLists.txt26
-rw-r--r--src/Cech_complex/example/cech_complex_example_from_points.cpp47
-rw-r--r--src/Cech_complex/example/cech_complex_step_by_step.cpp28
-rw-r--r--src/Cech_complex/include/gudhi/Cech_complex.h45
-rw-r--r--src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h107
-rw-r--r--src/Cech_complex/include/gudhi/Cech_complex_blocker.h98
-rw-r--r--src/Cech_complex/test/CMakeLists.txt19
-rw-r--r--src/Cech_complex/test/test_cech_complex.cpp73
-rw-r--r--src/Cech_complex/utilities/CMakeLists.txt24
-rw-r--r--src/Cech_complex/utilities/cech_persistence.cpp10
-rw-r--r--src/common/include/gudhi/distance_functions.h49
14 files changed, 382 insertions, 193 deletions
diff --git a/src/Cech_complex/benchmark/CMakeLists.txt b/src/Cech_complex/benchmark/CMakeLists.txt
index bc54c0f3..a6b3d70b 100644
--- a/src/Cech_complex/benchmark/CMakeLists.txt
+++ b/src/Cech_complex/benchmark/CMakeLists.txt
@@ -1,13 +1,15 @@
project(Cech_complex_benchmark)
-# Do not forget to copy test files in current binary dir
-file(COPY "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.0.1)
+ # Do not forget to copy test files in current binary dir
+ file(COPY "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
-if(TARGET Boost::filesystem)
- add_executable(cech_complex_benchmark cech_complex_benchmark.cpp)
- target_link_libraries(cech_complex_benchmark Boost::filesystem)
+ if(TARGET Boost::filesystem)
+ add_executable(cech_complex_benchmark cech_complex_benchmark.cpp)
+ target_link_libraries(cech_complex_benchmark Boost::filesystem)
- if (TBB_FOUND)
- target_link_libraries(cech_complex_benchmark ${TBB_LIBRARIES})
+ if (TBB_FOUND)
+ target_link_libraries(cech_complex_benchmark ${TBB_LIBRARIES})
+ endif()
endif()
-endif() \ No newline at end of file
+endif()
diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
index 2e4adce4..94c5fa4f 100644
--- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
+++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
@@ -17,6 +17,8 @@
#include <gudhi/Simplex_tree.h>
#include <gudhi/Miniball.hpp>
+#include <CGAL/Epick_d.h>
+
#include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations
#include <string>
@@ -30,7 +32,11 @@ using Point_cloud = std::vector<Point>;
using Points_off_reader = Gudhi::Points_off_reader<Point>;
using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>;
using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
-using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>;
+using Kernel = CGAL::Epick_d<CGAL::Dimension_tag<3>>;
+using Point_cgal = typename Kernel::Point_d;
+using Point_cloud_cgal = std::vector<Point_cgal>;
+using Points_off_reader_cgal = Gudhi::Points_off_reader<Point_cgal>;
+using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud_cgal, Kernel, Simplex_tree>;
class Minimal_enclosing_ball_radius {
public:
@@ -62,6 +68,7 @@ int main(int argc, char* argv[]) {
// Extract the points from the file filepoints
Points_off_reader off_reader(off_file_points);
+ Points_off_reader_cgal off_reader_cgal(off_file_points);
Gudhi::Clock euclidean_clock("Gudhi::Euclidean_distance");
// Compute the proximity graph of the points
@@ -76,16 +83,16 @@ int main(int argc, char* argv[]) {
off_reader.get_point_cloud(), threshold, Minimal_enclosing_ball_radius());
std::clog << miniball_clock << std::endl;
- Gudhi::Clock common_miniball_clock("Gudhi::Minimal_enclosing_ball_radius()");
+ Gudhi::Clock cgal_miniball_clock("Gudhi::Minimal_enclosing_ball_radius_cgal()");
// Compute the proximity graph of the points
- Proximity_graph common_miniball_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(
- off_reader.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius());
- std::clog << common_miniball_clock << std::endl;
+ Proximity_graph cgal_miniball_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(
+ off_reader_cgal.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius<Kernel>());
+ std::clog << cgal_miniball_clock << std::endl;
boost::filesystem::path full_path(boost::filesystem::current_path());
std::clog << "Current path is : " << full_path << std::endl;
- std::clog << "File name;Radius;Rips time;Cech time; Ratio Rips/Cech time;Rips nb simplices;Cech nb simplices;"
+ std::clog << "File name; Radius; Rips time; Cech time; Ratio Rips/Cech time; Rips nb simplices; Cech nb simplices;"
<< std::endl;
boost::filesystem::directory_iterator end_itr; // default construction yields past-the-end
for (boost::filesystem::directory_iterator itr(boost::filesystem::current_path()); itr != end_itr; ++itr) {
@@ -93,14 +100,16 @@ int main(int argc, char* argv[]) {
if (itr->path().extension() == ".off") // see below
{
Points_off_reader off_reader(itr->path().string());
+ Points_off_reader_cgal off_reader_cgal(itr->path().string());
+
Point p0 = off_reader.get_point_cloud()[0];
for (Filtration_value radius = 0.1; radius < 0.4; radius += 0.1) {
std::clog << itr->path().stem() << ";";
std::clog << radius << ";";
Gudhi::Clock rips_clock("Rips computation");
- Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), radius,
- Gudhi::Minimal_enclosing_ball_radius());
+ Rips_complex rips_complex_from_points(off_reader_cgal.get_point_cloud(), radius,
+ Gudhi::Minimal_enclosing_ball_radius<Kernel>());
Simplex_tree rips_stree;
rips_complex_from_points.create_complex(rips_stree, p0.size() - 1);
// ------------------------------------------
@@ -110,7 +119,7 @@ int main(int argc, char* argv[]) {
std::clog << rips_sec << ";";
Gudhi::Clock cech_clock("Cech computation");
- Cech_complex cech_complex_from_points(off_reader.get_point_cloud(), radius);
+ Cech_complex cech_complex_from_points(off_reader_cgal.get_point_cloud(), radius);
Simplex_tree cech_stree;
cech_complex_from_points.create_complex(cech_stree, p0.size() - 1);
// ------------------------------------------
diff --git a/src/Cech_complex/concept/SimplicialComplexForCech.h b/src/Cech_complex/concept/SimplicialComplexForCech.h
index 00c7df3a..6202fe92 100644
--- a/src/Cech_complex/concept/SimplicialComplexForCech.h
+++ b/src/Cech_complex/concept/SimplicialComplexForCech.h
@@ -47,8 +47,8 @@ struct SimplicialComplexForCech {
};
-} // namespace alpha_complex
+} // namespace cech_complex
} // namespace Gudhi
-#endif // CONCEPT_ALPHA_COMPLEX_SIMPLICIAL_COMPLEX_FOR_ALPHA_H_
+#endif // CONCEPT_CECH_COMPLEX_SIMPLICIAL_COMPLEX_FOR_CECH_H_
diff --git a/src/Cech_complex/example/CMakeLists.txt b/src/Cech_complex/example/CMakeLists.txt
index 1b08c7cb..4d11ace2 100644
--- a/src/Cech_complex/example/CMakeLists.txt
+++ b/src/Cech_complex/example/CMakeLists.txt
@@ -1,17 +1,19 @@
project(Cech_complex_examples)
-if (TARGET Boost::program_options)
- add_executable ( Cech_complex_example_step_by_step cech_complex_step_by_step.cpp )
- target_link_libraries(Cech_complex_example_step_by_step Boost::program_options)
- if (TBB_FOUND)
- target_link_libraries(Cech_complex_example_step_by_step ${TBB_LIBRARIES})
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.0.1)
+ if (TARGET Boost::program_options)
+ add_executable ( Cech_complex_example_step_by_step cech_complex_step_by_step.cpp )
+ target_link_libraries(Cech_complex_example_step_by_step Boost::program_options)
+ if (TBB_FOUND)
+ target_link_libraries(Cech_complex_example_step_by_step ${TBB_LIBRARIES})
+ endif()
+ add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:Cech_complex_example_step_by_step>
+ "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-d" "3")
endif()
- add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:Cech_complex_example_step_by_step>
- "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-d" "3")
-endif()
-add_executable ( Cech_complex_example_from_points cech_complex_example_from_points.cpp)
-if (TBB_FOUND)
- target_link_libraries(Cech_complex_example_from_points ${TBB_LIBRARIES})
+ add_executable ( Cech_complex_example_from_points cech_complex_example_from_points.cpp)
+ if (TBB_FOUND)
+ target_link_libraries(Cech_complex_example_from_points ${TBB_LIBRARIES})
+ endif()
+ add_test(NAME Cech_complex_example_from_points COMMAND $<TARGET_FILE:Cech_complex_example_from_points>)
endif()
-add_test(NAME Cech_complex_example_from_points COMMAND $<TARGET_FILE:Cech_complex_example_from_points>)
diff --git a/src/Cech_complex/example/cech_complex_example_from_points.cpp b/src/Cech_complex/example/cech_complex_example_from_points.cpp
index 1a1f708c..78861951 100644
--- a/src/Cech_complex/example/cech_complex_example_from_points.cpp
+++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp
@@ -1,6 +1,8 @@
#include <gudhi/Cech_complex.h>
#include <gudhi/Simplex_tree.h>
+#include <CGAL/Epeck_d.h> // For EXACT or SAFE version
+
#include <iostream>
#include <string>
#include <vector>
@@ -8,32 +10,47 @@
int main() {
// Type definitions
- using Point_cloud = std::vector<std::array<double, 2>>;
using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
using Filtration_value = Simplex_tree::Filtration_value;
- using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>;
+ using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+ using FT = typename Kernel::FT;
+ using Point = typename Kernel::Point_d;
+ using Point_cloud = std::vector<Point>;
+ using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel, Simplex_tree>;
Point_cloud points;
- points.push_back({1., 0.}); // 0
- points.push_back({0., 1.}); // 1
- points.push_back({2., 1.}); // 2
- points.push_back({3., 2.}); // 3
- points.push_back({0., 3.}); // 4
- points.push_back({3. + std::sqrt(3.), 3.}); // 5
- points.push_back({1., 4.}); // 6
- points.push_back({3., 4.}); // 7
- points.push_back({2., 4. + std::sqrt(3.)}); // 8
- points.push_back({0., 4.}); // 9
- points.push_back({-0.5, 2.}); // 10
+
+ std::vector<FT> point0({1., 0.});
+ points.emplace_back(point0.begin(), point0.end());
+ std::vector<FT> point1({0., 1.});
+ points.emplace_back(point1.begin(), point1.end());
+ std::vector<FT> point2({2., 1.});
+ points.emplace_back(point2.begin(), point2.end());
+ std::vector<FT> point3({3., 2.});
+ points.emplace_back(point3.begin(), point3.end());
+ std::vector<FT> point4({0., 3.});
+ points.emplace_back(point4.begin(), point4.end());
+ std::vector<FT> point5({3. + std::sqrt(3.), 3.});
+ points.emplace_back(point5.begin(), point5.end());
+ std::vector<FT> point6({1., 4.});
+ points.emplace_back(point6.begin(), point6.end());
+ std::vector<FT> point7({3., 4.});
+ points.emplace_back(point7.begin(), point7.end());
+ std::vector<FT> point8({2., 4. + std::sqrt(3.)});
+ points.emplace_back(point8.begin(), point8.end());
+ std::vector<FT> point9({0., 4.});
+ points.emplace_back(point9.begin(), point9.end());
+ std::vector<FT> point10({-0.5, 2.});
+ points.emplace_back(point10.begin(), point10.end());
// ----------------------------------------------------------------------------
// Init of a Cech complex from points
// ----------------------------------------------------------------------------
- Filtration_value max_radius = 1.;
+ Filtration_value max_radius = 100.; //100.;
Cech_complex cech_complex_from_points(points, max_radius);
Simplex_tree stree;
- cech_complex_from_points.create_complex(stree, 2);
+ cech_complex_from_points.create_complex(stree, 6); //6
// ----------------------------------------------------------------------------
// Display information about the one skeleton Cech complex
// ----------------------------------------------------------------------------
diff --git a/src/Cech_complex/example/cech_complex_step_by_step.cpp b/src/Cech_complex/example/cech_complex_step_by_step.cpp
index f59f0293..c8dd1585 100644
--- a/src/Cech_complex/example/cech_complex_step_by_step.cpp
+++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp
@@ -9,11 +9,11 @@
*/
#include <gudhi/graph_simplicial_complex.h>
-#include <gudhi/distance_functions.h>
+#include <gudhi/Cech_complex/Cech_kernel.h>
#include <gudhi/Simplex_tree.h>
#include <gudhi/Points_off_io.h>
-#include <gudhi/Miniball.hpp>
+#include <CGAL/Epeck_d.h>
#include <boost/program_options.hpp>
@@ -24,9 +24,9 @@
#include <map>
// ----------------------------------------------------------------------------
-// rips_persistence_step_by_step is an example of each step that is required to
-// build a Rips over a Simplex_tree. Please refer to rips_persistence to see
-// how to do the same thing with the Rips_complex wrapper for less detailed
+// cech_complex_step_by_step is an example of each step that is required to
+// build a Cech over a Simplex_tree. Please refer to cech_complex_example_from_points to see
+// how to do the same thing with the Cech complex wrapper for less detailed
// steps.
// ----------------------------------------------------------------------------
@@ -34,16 +34,14 @@
using Simplex_tree = Gudhi::Simplex_tree<>;
using Simplex_handle = Simplex_tree::Simplex_handle;
using Filtration_value = Simplex_tree::Filtration_value;
-using Point = std::vector<double>;
+using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+using Point = typename Kernel::Point_d;
using Points_off_reader = Gudhi::Points_off_reader<Point>;
using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>;
class Cech_blocker {
private:
using Point_cloud = std::vector<Point>;
- using Point_iterator = Point_cloud::const_iterator;
- using Coordinate_iterator = Point::const_iterator;
- using Min_sphere = Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
public:
bool operator()(Simplex_handle sh) {
@@ -54,7 +52,7 @@ class Cech_blocker {
std::clog << "#(" << vertex << ")#";
#endif // DEBUG_TRACES
}
- Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points);
+ Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius<Kernel>()(points);
#ifdef DEBUG_TRACES
std::clog << "radius = " << radius << " - " << (radius > max_radius_) << std::endl;
#endif // DEBUG_TRACES
@@ -63,14 +61,12 @@ class Cech_blocker {
}
Cech_blocker(Simplex_tree& simplex_tree, Filtration_value max_radius, const std::vector<Point>& point_cloud)
: simplex_tree_(simplex_tree), max_radius_(max_radius), point_cloud_(point_cloud) {
- dimension_ = point_cloud_[0].size();
}
private:
Simplex_tree simplex_tree_;
Filtration_value max_radius_;
std::vector<Point> point_cloud_;
- int dimension_;
};
void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& max_radius, int& dim_max);
@@ -87,9 +83,9 @@ int main(int argc, char* argv[]) {
// Compute the proximity graph of the points
Proximity_graph prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), max_radius,
- Gudhi::Minimal_enclosing_ball_radius());
+ Gudhi::Minimal_enclosing_ball_radius<Kernel>());
- // Construct the Rips complex in a Simplex Tree
+ // Construct the Cech complex in a Simplex Tree
Simplex_tree st;
// insert the proximity graph in the simplex tree
st.insert_graph(prox_graph);
@@ -129,9 +125,9 @@ void program_options(int argc, char* argv[], std::string& off_file_points, Filtr
visible.add_options()("help,h", "produce help message")(
"max-radius,r",
po::value<Filtration_value>(&max_radius)->default_value(std::numeric_limits<Filtration_value>::infinity()),
- "Maximal length of an edge for the Rips complex construction.")(
+ "Maximal length of an edge for the Cech complex construction.")(
"cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
- "Maximal dimension of the Rips complex we want to compute.");
+ "Maximal dimension of the Cech complex we want to compute.");
po::positional_options_description pos;
pos.add("input-file", 1);
diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h
index b0871e10..0031d861 100644
--- a/src/Cech_complex/include/gudhi/Cech_complex.h
+++ b/src/Cech_complex/include/gudhi/Cech_complex.h
@@ -11,14 +11,13 @@
#ifndef CECH_COMPLEX_H_
#define CECH_COMPLEX_H_
-#include <gudhi/distance_functions.h> // for Gudhi::Minimal_enclosing_ball_radius
+#include <gudhi/Cech_complex/Cech_kernel.h> // for Gudhi::Minimal_enclosing_ball_radius
#include <gudhi/graph_simplicial_complex.h> // for Gudhi::Proximity_graph
#include <gudhi/Debug_utils.h> // for GUDHI_CHECK
#include <gudhi/Cech_complex_blocker.h> // for Gudhi::cech_complex::Cech_blocker
#include <iostream>
#include <stdexcept> // for exception management
-#include <vector>
namespace Gudhi {
@@ -40,7 +39,7 @@ namespace cech_complex {
* \tparam ForwardPointRange must be a range for which `std::begin()` and `std::end()` methods return input
* iterators on a point. `std::begin()` and `std::end()` methods are also required for a point.
*/
-template <typename SimplicialComplexForProximityGraph, typename ForwardPointRange>
+template <typename SimplicialComplexForProximityGraph, typename ForwardPointRange, typename Kernel, typename SimplicialComplexForCechComplex>
class Cech_complex {
private:
// Required by compute_proximity_graph
@@ -48,16 +47,17 @@ class Cech_complex {
using Filtration_value = typename SimplicialComplexForProximityGraph::Filtration_value;
using Proximity_graph = Gudhi::Proximity_graph<SimplicialComplexForProximityGraph>;
- // Retrieve Coordinate type from ForwardPointRange
- using Point_from_range_iterator = typename boost::range_const_iterator<ForwardPointRange>::type;
- using Point_from_range = typename std::iterator_traits<Point_from_range_iterator>::value_type;
- using Coordinate_iterator = typename boost::range_const_iterator<Point_from_range>::type;
- using Coordinate = typename std::iterator_traits<Coordinate_iterator>::value_type;
-
public:
- // Point and Point_cloud type definition
- using Point = std::vector<Coordinate>;
- using Point_cloud = std::vector<Point>;
+
+ using cech_blocker = Cech_blocker<SimplicialComplexForCechComplex, Cech_complex, Kernel>;
+
+ using Point_d = typename cech_blocker::Point_d;
+ using Point_cloud = std::vector<Point_d>;
+
+ // Numeric type of coordinates in the kernel
+ using FT = typename cech_blocker::FT;
+ // Sphere is a pair of point and squared radius.
+ using Sphere = typename std::pair<Point_d, FT>;
public:
/** \brief Cech_complex constructor from a list of points.
@@ -70,11 +70,14 @@ class Cech_complex {
*/
Cech_complex(const ForwardPointRange& points, Filtration_value max_radius) : max_radius_(max_radius) {
// Point cloud deep copy
- point_cloud_.reserve(boost::size(points));
- for (auto&& point : points) point_cloud_.emplace_back(std::begin(point), std::end(point));
+
+// point_cloud_.reserve(boost::size(points));
+// for (auto&& point : points) point_cloud_.emplace_back(point.cartesian_begin(), point.cartesian_end());
+
+ point_cloud_.assign(points.begin(), points.end());
cech_skeleton_graph_ = Gudhi::compute_proximity_graph<SimplicialComplexForProximityGraph>(
- point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius());
+ point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius<Kernel>());
}
/** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal
@@ -85,7 +88,6 @@ class Cech_complex {
* @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0.
*
*/
- template <typename SimplicialComplexForCechComplex>
void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) {
GUDHI_CHECK(complex.num_vertices() == 0,
std::invalid_argument("Cech_complex::create_complex - simplicial complex is not empty"));
@@ -93,8 +95,7 @@ class Cech_complex {
// insert the proximity graph in the simplicial complex
complex.insert_graph(cech_skeleton_graph_);
// expand the graph until dimension dim_max
- complex.expansion_with_blockers(dim_max,
- Cech_blocker<SimplicialComplexForCechComplex, Cech_complex>(&complex, this));
+ complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this));
}
/** @return max_radius value given at construction. */
@@ -103,12 +104,18 @@ class Cech_complex {
/** @param[in] vertex Point position in the range.
* @return The point.
*/
- const Point& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; }
+ const Point_d& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; }
+
+ /**
+ * @return Vector of cached spheres.
+ */
+ std::vector<Sphere> & get_cache() { return cache_; }
private:
Proximity_graph cech_skeleton_graph_;
Filtration_value max_radius_;
Point_cloud point_cloud_;
+ std::vector<Sphere> cache_;
};
} // namespace cech_complex
diff --git a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h
new file mode 100644
index 00000000..89012206
--- /dev/null
+++ b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h
@@ -0,0 +1,107 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Hind Montassif
+ *
+ * Copyright (C) 2021 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef CECH_KERNEL_H_
+#define CECH_KERNEL_H_
+
+#include <CGAL/Epeck_d.h> // for #include <CGAL/NewKernel_d/KernelD_converter.h>
+
+#include <cmath> // for std::sqrt
+#include <vector>
+
+namespace Gudhi {
+
+// namespace cech_complex {
+
+/** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates.
+ * The points are assumed to have the same dimension. */
+template<typename Kernel>
+class Minimal_enclosing_ball_radius {
+ private:
+ Kernel kernel_;
+ public:
+ using Point = typename Kernel::Point_d;
+ using Point_cloud = typename std::vector<Point>;
+
+ /** \brief Enclosing ball radius from two points using CGAL.
+ *
+ * @param[in] point_1
+ * @param[in] point_2
+ * @return Enclosing ball radius for the two points.
+ * \tparam Point must be a Kernel::Point_d from CGAL.
+ *
+ */
+ double operator()(const Point& point_1, const Point& point_2) const {
+ return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.;
+ }
+
+
+ /** \brief Enclosing ball radius from a point cloud using CGAL.
+ *
+ * @param[in] point_cloud The points.
+ * @return Enclosing ball radius for the points.
+ * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL.
+ *
+ */
+ double operator()(const Point_cloud& point_cloud) const {
+ return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end())));
+ }
+
+};
+
+/**
+ * \class Cech_kernel
+ * \brief Cech complex kernel container.
+ *
+ * \details
+ * The Cech complex kernel container stores CGAL Kernel and dispatch basic computations.
+ */
+
+// template < typename Kernel >
+// class Cech_kernel<Kernel> {
+// private:
+// // Kernel for functions access.
+// Kernel kernel_;
+// public:
+// using Point_d = typename Kernel::Point_d;
+// // Numeric type of coordinates in the kernel
+// using FT = typename Kernel::FT;
+// // Sphere is a pair of point and squared radius.
+// using Sphere = typename std::pair<Point_d, FT>;
+//
+// int get_dimension(const Point_d& p0) const {
+// return kernel_.point_dimension_d_object()(p0);
+// }
+//
+// template<class PointIterator>
+// Sphere get_sphere(PointIterator begin, PointIterator end) const {
+// Point_d c = kernel_.construct_circumcenter_d_object()(begin, end);
+// FT r = kernel_.squared_distance_d_object()(c, *begin);
+// return std::make_pair(std::move(c), std::move(r));
+// }
+//
+// template<class PointIterator>
+// FT get_squared_radius(PointIterator begin, PointIterator end) const {
+// return kernel_.compute_squared_radius_d_object()(begin, end);
+// }
+//
+// FT get_squared_radius(const Sphere& sph) const {
+// return sph.second;
+// }
+// };
+
+
+//} // namespace cech_complex
+
+// namespace cechcomplex = cech_complex;
+
+} // namespace Gudhi
+
+#endif // CECH_KERNEL_H_
diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
index 31b9aab5..f7f86534 100644
--- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
+++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
@@ -11,10 +11,11 @@
#ifndef CECH_COMPLEX_BLOCKER_H_
#define CECH_COMPLEX_BLOCKER_H_
-#include <gudhi/distance_functions.h> // for Gudhi::Minimal_enclosing_ball_radius
+#include <CGAL/NT_converter.h> // for casting from FT to double
#include <iostream>
#include <vector>
+#include <set>
#include <cmath> // for std::sqrt
namespace Gudhi {
@@ -35,28 +36,104 @@ namespace cech_complex {
*
* \tparam Chech_complex is required by the blocker.
*/
-template <typename SimplicialComplexForCech, typename Cech_complex>
+template <typename SimplicialComplexForCech, typename Cech_complex, typename Kernel>
class Cech_blocker {
private:
- using Point_cloud = typename Cech_complex::Point_cloud;
using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle;
using Filtration_value = typename SimplicialComplexForCech::Filtration_value;
public:
+
+ using Point_d = typename Kernel::Point_d;
+ // Numeric type of coordinates in the kernel
+ using FT = typename Kernel::FT;
+ // Sphere is a pair of point and squared radius.
+ using Sphere = typename std::pair<Point_d, FT>;
+
+ template<class PointIterator>
+ FT get_squared_radius(PointIterator begin, PointIterator end) const {
+ return kernel_.compute_squared_radius_d_object()(begin, end);
+ }
+
+ template<class PointIterator>
+ Sphere get_sphere(PointIterator begin, PointIterator end) const {
+ Point_d c = kernel_.construct_circumcenter_d_object()(begin, end);
+ FT r = kernel_.squared_distance_d_object()(c, *begin);
+ return std::make_pair(std::move(c), std::move(r));
+ }
+
+
/** \internal \brief Čech complex blocker operator() - the oracle - assigns the filtration value from the simplex
* radius and returns if the simplex expansion must be blocked.
* \param[in] sh The Simplex_handle.
* \return true if the simplex radius is greater than the Cech_complex max_radius*/
bool operator()(Simplex_handle sh) {
- Point_cloud points;
- for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) {
- points.push_back(cc_ptr_->get_point(vertex));
-#ifdef DEBUG_TRACES
- std::clog << "#(" << vertex << ")#";
-#endif // DEBUG_TRACES
+ using Point_cloud = std::vector<Point_d>;
+ CGAL::NT_converter<FT, double> cast_to_double;
+ Filtration_value radius = 0.;
+
+ // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices
+ Sphere min_enclos_ball;
+ CGAL::NT_converter<double, FT> cast_to_FT;
+ min_enclos_ball.second = cast_to_FT(std::numeric_limits<double>::max());
+ Point_cloud face_points;
+ for (auto face : sc_ptr_->boundary_simplex_range(sh)) {
+ // Find which vertex of sh is missing in face. We rely on the fact that simplex_vertex_range is sorted.
+ auto longlist = sc_ptr_->simplex_vertex_range(sh);
+ auto shortlist = sc_ptr_->simplex_vertex_range(face);
+
+ auto longiter = std::begin(longlist);
+ auto shortiter = std::begin(shortlist);
+ auto enditer = std::end(shortlist);
+ while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; }
+ auto extra = *longiter; // Vertex_handle
+
+ for (auto vertex : sc_ptr_->simplex_vertex_range(face)) {
+ face_points.push_back(cc_ptr_->get_point(vertex));
+ #ifdef DEBUG_TRACES
+ std::clog << "#(" << vertex << ")#";
+ #endif // DEBUG_TRACES
+ }
+ Sphere sph;
+ auto k = sc_ptr_->key(face);
+ if(k != sc_ptr_->null_key()) {
+ sph = cc_ptr_->get_cache().at(k);
+ }
+ else {
+ sph = get_sphere(face_points.cbegin(), face_points.cend());
+ }
+ face_points.clear();
+
+ if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) {
+ radius = std::sqrt(cast_to_double(sph.second));
+ #ifdef DEBUG_TRACES
+ std::clog << "circumcenter: " << sph.first << ", radius: " << radius << std::endl;
+ #endif // DEBUG_TRACES
+ if (cast_to_double(sph.second) < cast_to_double(min_enclos_ball.second))
+ min_enclos_ball = sph;
+ }
+ }
+ // Get the minimal radius of all faces enclosing balls if exists
+ if(cast_to_double(min_enclos_ball.second) != std::numeric_limits<double>::max()) {
+ radius = std::sqrt(cast_to_double(min_enclos_ball.second));
+
+ sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size());
+ cc_ptr_->get_cache().push_back(min_enclos_ball);
}
- Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points);
+
+ if (radius == 0.) { // Spheres of each face don't contain the whole simplex
+ Point_cloud points;
+ for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) {
+ points.push_back(cc_ptr_->get_point(vertex));
+ }
+ Sphere sph = get_sphere(points.cbegin(), points.cend());
+ radius = std::sqrt(cast_to_double(sph.second));
+
+ sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size());
+ cc_ptr_->get_cache().push_back(sph);
+ }
+
#ifdef DEBUG_TRACES
if (radius > cc_ptr_->max_radius()) std::clog << "radius > max_radius => expansion is blocked\n";
#endif // DEBUG_TRACES
@@ -70,6 +147,7 @@ class Cech_blocker {
private:
SimplicialComplexForCech* sc_ptr_;
Cech_complex* cc_ptr_;
+ Kernel kernel_;
};
} // namespace cech_complex
diff --git a/src/Cech_complex/test/CMakeLists.txt b/src/Cech_complex/test/CMakeLists.txt
index e6a2a18f..2d736f27 100644
--- a/src/Cech_complex/test/CMakeLists.txt
+++ b/src/Cech_complex/test/CMakeLists.txt
@@ -1,11 +1,14 @@
-include(GUDHI_boost_test)
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.0.1)
+ include(GUDHI_boost_test)
-add_executable ( Cech_complex_test_unit test_cech_complex.cpp )
-if (TBB_FOUND)
- target_link_libraries(Cech_complex_test_unit ${TBB_LIBRARIES})
-endif()
+ add_executable ( Cech_complex_test_unit test_cech_complex.cpp )
+ if (TBB_FOUND)
+ target_link_libraries(Cech_complex_test_unit ${TBB_LIBRARIES})
+ endif()
+
+ # Do not forget to copy test files in current binary dir
+ file(COPY "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
-# Do not forget to copy test files in current binary dir
-file(COPY "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+ gudhi_add_boost_test(Cech_complex_test_unit)
-gudhi_add_boost_test(Cech_complex_test_unit)
+endif()
diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp
index 6e00d7b5..ca7a9778 100644
--- a/src/Cech_complex/test/test_cech_complex.cpp
+++ b/src/Cech_complex/test/test_cech_complex.cpp
@@ -22,21 +22,21 @@
// to construct Cech_complex from a OFF file of points
#include <gudhi/Points_off_io.h>
#include <gudhi/Simplex_tree.h>
-#include <gudhi/distance_functions.h>
+#include <gudhi/Cech_complex/Cech_kernel.h>
#include <gudhi/Unitary_tests_utils.h>
-#include <gudhi/Miniball.hpp>
+
+#include <CGAL/Epeck_d.h> // For EXACT or SAFE version
// Type definitions
using Simplex_tree = Gudhi::Simplex_tree<>;
using Filtration_value = Simplex_tree::Filtration_value;
-using Point = std::vector<Filtration_value>;
+using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+using FT = typename Kernel::FT;
+using Point = typename Kernel::Point_d;
+
using Point_cloud = std::vector<Point>;
using Points_off_reader = Gudhi::Points_off_reader<Point>;
-using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>;
-
-using Point_iterator = Point_cloud::const_iterator;
-using Coordinate_iterator = Point::const_iterator;
-using Min_sphere = Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
+using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel, Simplex_tree>;
BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) {
// ----------------------------------------------------------------------------
@@ -45,17 +45,29 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) {
//
// ----------------------------------------------------------------------------
Point_cloud points;
- points.push_back({1., 0.}); // 0
- points.push_back({0., 1.}); // 1
- points.push_back({2., 1.}); // 2
- points.push_back({3., 2.}); // 3
- points.push_back({0., 3.}); // 4
- points.push_back({3. + std::sqrt(3.), 3.}); // 5
- points.push_back({1., 4.}); // 6
- points.push_back({3., 4.}); // 7
- points.push_back({2., 4. + std::sqrt(3.)}); // 8
- points.push_back({0., 4.}); // 9
- points.push_back({-0.5, 2.}); // 10
+
+ std::vector<FT> point0({1., 0.});
+ points.emplace_back(point0.begin(), point0.end());
+ std::vector<FT> point1({0., 1.});
+ points.emplace_back(point1.begin(), point1.end());
+ std::vector<FT> point2({2., 1.});
+ points.emplace_back(point2.begin(), point2.end());
+ std::vector<FT> point3({3., 2.});
+ points.emplace_back(point3.begin(), point3.end());
+ std::vector<FT> point4({0., 3.});
+ points.emplace_back(point4.begin(), point4.end());
+ std::vector<FT> point5({3. + std::sqrt(3.), 3.});
+ points.emplace_back(point5.begin(), point5.end());
+ std::vector<FT> point6({1., 4.});
+ points.emplace_back(point6.begin(), point6.end());
+ std::vector<FT> point7({3., 4.});
+ points.emplace_back(point7.begin(), point7.end());
+ std::vector<FT> point8({2., 4. + std::sqrt(3.)});
+ points.emplace_back(point8.begin(), point8.end());
+ std::vector<FT> point9({0., 4.});
+ points.emplace_back(point9.begin(), point9.end());
+ std::vector<FT> point10({-0.5, 2.});
+ points.emplace_back(point10.begin(), point10.end());
Filtration_value max_radius = 1.0;
std::clog << "========== NUMBER OF POINTS = " << points.size() << " - Cech max_radius = " << max_radius
@@ -96,11 +108,11 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) {
std::clog << vertex << ",";
vp.push_back(points.at(vertex));
}
- std::clog << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1))
+ std::clog << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius<Kernel>()(vp.at(0), vp.at(1))
<< " - filtration =" << st.filtration(f_simplex) << std::endl;
BOOST_CHECK(vp.size() == 2);
GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex),
- Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1)));
+ Gudhi::Minimal_enclosing_ball_radius<Kernel>()(vp.at(0), vp.at(1)));
}
}
@@ -125,35 +137,34 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) {
for (std::size_t vertex = 0; vertex <= 2; vertex++) {
points012.push_back(cech_complex_for_doc.get_point(vertex));
}
- std::size_t dimension = points[0].end() - points[0].begin();
- Min_sphere ms012(dimension, points012.begin(), points012.end());
+ Kernel kern;
Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2}));
- std::clog << "f012= " << f012 << " | ms012_radius= " << std::sqrt(ms012.squared_radius()) << std::endl;
+ std::clog << "f012= " << f012 << std::endl;
- GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(ms012.squared_radius()));
+ CGAL::NT_converter<FT, double> cast_to_double;
+ GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points012.begin(), points012.end()))));
Point_cloud points1410;
points1410.push_back(cech_complex_for_doc.get_point(1));
points1410.push_back(cech_complex_for_doc.get_point(4));
points1410.push_back(cech_complex_for_doc.get_point(10));
- Min_sphere ms1410(dimension, points1410.begin(), points1410.end());
Simplex_tree::Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10}));
- std::clog << "f1410= " << f1410 << " | ms1410_radius= " << std::sqrt(ms1410.squared_radius()) << std::endl;
+ std::clog << "f1410= " << f1410 << std::endl;
- GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, std::sqrt(ms1410.squared_radius()));
+ // In this case, the computed sphere using CGAL kernel does not match the minimal enclosing ball; the filtration value check is therefore done against a hardcoded value
+ GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, 1.);
Point_cloud points469;
points469.push_back(cech_complex_for_doc.get_point(4));
points469.push_back(cech_complex_for_doc.get_point(6));
points469.push_back(cech_complex_for_doc.get_point(9));
- Min_sphere ms469(dimension, points469.begin(), points469.end());
Simplex_tree::Filtration_value f469 = st2.filtration(st2.find({4, 6, 9}));
- std::clog << "f469= " << f469 << " | ms469_radius= " << std::sqrt(ms469.squared_radius()) << std::endl;
+ std::clog << "f469= " << f469 << std::endl;
- GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(ms469.squared_radius()));
+ GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points469.begin(), points469.end()))));
BOOST_CHECK((st2.find({6, 7, 8}) == st2.null_simplex()));
BOOST_CHECK((st2.find({3, 5, 7}) == st2.null_simplex()));
diff --git a/src/Cech_complex/utilities/CMakeLists.txt b/src/Cech_complex/utilities/CMakeLists.txt
index b183c8d8..e80a698e 100644
--- a/src/Cech_complex/utilities/CMakeLists.txt
+++ b/src/Cech_complex/utilities/CMakeLists.txt
@@ -1,15 +1,17 @@
-project(Cech_complex_utilities)
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.0.1)
+ project(Cech_complex_utilities)
-if (TARGET Boost::program_options)
- add_executable(cech_persistence cech_persistence.cpp)
- target_link_libraries(cech_persistence Boost::program_options)
+ if (TARGET Boost::program_options)
+ add_executable(cech_persistence cech_persistence.cpp)
+ target_link_libraries(cech_persistence Boost::program_options)
- if (TBB_FOUND)
- target_link_libraries(cech_persistence ${TBB_LIBRARIES})
- endif()
+ if (TBB_FOUND)
+ target_link_libraries(cech_persistence ${TBB_LIBRARIES})
+ endif()
- add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:cech_persistence>
- "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3")
+ add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:cech_persistence>
+ "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3")
- install(TARGETS cech_persistence DESTINATION bin)
-endif() \ No newline at end of file
+ install(TARGETS cech_persistence DESTINATION bin)
+ endif()
+endif()
diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp
index daea08e2..ccf63e3e 100644
--- a/src/Cech_complex/utilities/cech_persistence.cpp
+++ b/src/Cech_complex/utilities/cech_persistence.cpp
@@ -9,13 +9,15 @@
*/
#include <gudhi/Cech_complex.h>
-#include <gudhi/distance_functions.h>
+#include <gudhi/Cech_complex/Cech_kernel.h>
#include <gudhi/Simplex_tree.h>
#include <gudhi/Persistent_cohomology.h>
#include <gudhi/Points_off_io.h>
#include <boost/program_options.hpp>
+#include <CGAL/Epeck_d.h> // For EXACT or SAFE version
+
#include <string>
#include <vector>
#include <limits> // infinity
@@ -23,10 +25,12 @@
// Types definition
using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
using Filtration_value = Simplex_tree::Filtration_value;
-using Point = std::vector<double>;
+
+using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+using Point = typename Kernel::Point_d;
using Point_cloud = std::vector<Point>;
using Points_off_reader = Gudhi::Points_off_reader<Point>;
-using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>;
+using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel, Simplex_tree>;
using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>;
diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h
index 9bbc62b7..5e5a1e31 100644
--- a/src/common/include/gudhi/distance_functions.h
+++ b/src/common/include/gudhi/distance_functions.h
@@ -13,8 +13,6 @@
#include <gudhi/Debug_utils.h>
-#include <gudhi/Miniball.hpp>
-
#include <boost/range/metafunctions.hpp>
#include <boost/range/size.hpp>
@@ -59,53 +57,6 @@ class Euclidean_distance {
}
};
-/** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates.
- * The points are assumed to have the same dimension. */
-class Minimal_enclosing_ball_radius {
- public:
- /** \brief Minimal_enclosing_ball_radius from two points.
- *
- * @param[in] point_1 First point.
- * @param[in] point_2 second point.
- * @return The minimal enclosing ball radius for the two points (aka. Euclidean distance / 2.).
- *
- * \tparam Point must be a range of Cartesian coordinates.
- *
- */
- template< typename Point >
- typename std::iterator_traits<typename boost::range_iterator<Point>::type>::value_type
- operator()(const Point& point_1, const Point& point_2) const {
- return Euclidean_distance()(point_1, point_2) / 2.;
- }
- /** \brief Minimal_enclosing_ball_radius from a point cloud.
- *
- * @param[in] point_cloud The points.
- * @return The minimal enclosing ball radius for the points.
- *
- * \tparam Point_cloud must be a range of points with Cartesian coordinates.
- * Point_cloud is a range over a range of Coordinate.
- *
- */
- template< typename Point_cloud,
- typename Point_iterator = typename boost::range_const_iterator<Point_cloud>::type,
- typename Point = typename std::iterator_traits<Point_iterator>::value_type,
- typename Coordinate_iterator = typename boost::range_const_iterator<Point>::type,
- typename Coordinate = typename std::iterator_traits<Coordinate_iterator>::value_type>
- Coordinate
- operator()(const Point_cloud& point_cloud) const {
- using Min_sphere = Miniball::Miniball<Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
-
- Min_sphere ms(boost::size(*point_cloud.begin()), point_cloud.begin(), point_cloud.end());
-#ifdef DEBUG_TRACES
- std::clog << "Minimal_enclosing_ball_radius = " << std::sqrt(ms.squared_radius()) << " | nb points = "
- << boost::size(point_cloud) << " | dimension = "
- << boost::size(*point_cloud.begin()) << std::endl;
-#endif // DEBUG_TRACES
-
- return std::sqrt(ms.squared_radius());
- }
-};
-
} // namespace Gudhi
#endif // DISTANCE_FUNCTIONS_H_