summaryrefslogtreecommitdiff
path: root/src/Alpha_complex
diff options
context:
space:
mode:
Diffstat (limited to 'src/Alpha_complex')
-rw-r--r--src/Alpha_complex/doc/Intro_alpha_complex.h15
-rw-r--r--src/Alpha_complex/include/gudhi/Alpha_complex.h117
-rw-r--r--src/Alpha_complex/test/CMakeLists.txt4
-rw-r--r--src/Alpha_complex/test/Delaunay_complex_unit_test.cpp72
4 files changed, 146 insertions, 62 deletions
diff --git a/src/Alpha_complex/doc/Intro_alpha_complex.h b/src/Alpha_complex/doc/Intro_alpha_complex.h
index a8b1a106..60da7169 100644
--- a/src/Alpha_complex/doc/Intro_alpha_complex.h
+++ b/src/Alpha_complex/doc/Intro_alpha_complex.h
@@ -46,16 +46,17 @@ namespace alpha_complex {
* \cite cgal:s-gkd-19b from CGAL as template parameter.
*
* \remark
- * - When an \f$\alpha\f$-complex is constructed with an infinite value of \f$ \alpha^2 \f$, the complex is a Delaunay
- * complex (with special filtration values).
+ * - When the simplicial complex is constructed with an infinite value of \f$ \alpha^2 \f$, the complex is a Delaunay
+ * complex with special filtration values. The Delaunay complex without filtration values is also available by passing
+ * `default_filtration_value=true` to `Alpha_complex::create_complex`.
* - For people only interested in the topology of the \ref alpha_complex (for instance persistence),
* \ref alpha_complex is equivalent to the \ref cech_complex and much smaller if you do not bound the radii.
* \ref cech_complex can still make sense in higher dimension precisely because you can bound the radii.
- * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass exact=true to create_complex, the
- * filtration values are the exact ones converted to the filtration value type of the simplicial complex. This can be
- * very slow. If you pass exact=false (the default), the filtration values are only guaranteed to have a small
- * multiplicative error compared to the exact value, see <code><a class="el" target="_blank"
- * href="https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html">
+ * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass `exact=true` to
+ * `Alpha_complex::create_complex`, the filtration values are the exact ones converted to the filtration value type of
+ * the simplicial complex. This can be very slow. If you pass `exact=false` (the default), the filtration values are
+ * only guaranteed to have a small multiplicative error compared to the exact value, see <code>
+ * <a class="el" target="_blank" href="https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html">
* CGAL::Lazy_exact_nt<NT>::set_relative_precision_of_to_double</a></code> for details. A drawback, when computing
* persistence, is that an empty exact interval [10^12,10^12] may become a non-empty approximate interval
* [10^12,10^12+10^6]. Using `CGAL::Epick_d` makes the computations slightly faster, and the combinatorics are still
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h
index 0839ae6c..65c969d2 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h
@@ -248,16 +248,20 @@ class Alpha_complex {
public:
/** \brief Inserts all Delaunay triangulation into the simplicial complex.
- * It also computes the filtration values accordingly to the \ref createcomplexalgorithm
+ * It also computes the filtration values accordingly to the \ref createcomplexalgorithm if default_filtration_value
+ * is not set.
*
* \tparam SimplicialComplexForAlpha must meet `SimplicialComplexForAlpha` concept.
*
* @param[in] complex SimplicialComplexForAlpha to be created.
* @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$, and there is very
- * little point using anything else since it does not save time.
+ * little point using anything else since it does not save time. Useless if `default_filtration_value` is set to
+ * `true`.
* @param[in] exact Exact filtration values computation. Not exact if `Kernel` is not <a target="_blank"
* href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>.
- *
+ * @param[in] default_filtration_value Set this value to `true` if filtration values are not needed to be computed.
+ * Default value is `false` (which means compute the filtration values).
+ *
* @return true if creation succeeds, false otherwise.
*
* @pre Delaunay triangulation must be already constructed with dimension strictly greater than 0.
@@ -269,7 +273,8 @@ class Alpha_complex {
typename Filtration_value = typename SimplicialComplexForAlpha::Filtration_value>
bool create_complex(SimplicialComplexForAlpha& complex,
Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(),
- bool exact = false) {
+ bool exact = false,
+ bool default_filtration_value = false) {
// From SimplicialComplexForAlpha type required to insert into a simplicial complex (with or without subfaces).
typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle;
typedef typename SimplicialComplexForAlpha::Simplex_handle Simplex_handle;
@@ -316,62 +321,64 @@ class Alpha_complex {
}
// --------------------------------------------------------------------------------------------
- // --------------------------------------------------------------------------------------------
- // Will be re-used many times
- Vector_of_CGAL_points pointVector;
- // ### For i : d -> 0
- for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) {
- // ### Foreach Sigma of dim i
- for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) {
- int f_simplex_dim = complex.dimension(f_simplex);
- if (decr_dim == f_simplex_dim) {
- pointVector.clear();
-#ifdef DEBUG_TRACES
- std::clog << "Sigma of dim " << decr_dim << " is";
-#endif // DEBUG_TRACES
- for (auto vertex : complex.simplex_vertex_range(f_simplex)) {
- pointVector.push_back(get_point(vertex));
-#ifdef DEBUG_TRACES
- std::clog << " " << vertex;
-#endif // DEBUG_TRACES
- }
-#ifdef DEBUG_TRACES
- std::clog << std::endl;
-#endif // DEBUG_TRACES
- // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma)
- if (std::isnan(complex.filtration(f_simplex))) {
- Filtration_value alpha_complex_filtration = 0.0;
- // No need to compute squared_radius on a single point - alpha is 0.0
- if (f_simplex_dim > 0) {
- // squared_radius function initialization
- Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object();
-
- CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv;
- auto sqrad = squared_radius(pointVector.begin(), pointVector.end());
-#if CGAL_VERSION_NR >= 1050000000
- if(exact) CGAL::exact(sqrad);
-#endif
- alpha_complex_filtration = cv(sqrad);
+ if (!default_filtration_value) {
+ // --------------------------------------------------------------------------------------------
+ // Will be re-used many times
+ Vector_of_CGAL_points pointVector;
+ // ### For i : d -> 0
+ for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) {
+ // ### Foreach Sigma of dim i
+ for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) {
+ int f_simplex_dim = complex.dimension(f_simplex);
+ if (decr_dim == f_simplex_dim) {
+ pointVector.clear();
+ #ifdef DEBUG_TRACES
+ std::clog << "Sigma of dim " << decr_dim << " is";
+ #endif // DEBUG_TRACES
+ for (auto vertex : complex.simplex_vertex_range(f_simplex)) {
+ pointVector.push_back(get_point(vertex));
+ #ifdef DEBUG_TRACES
+ std::clog << " " << vertex;
+ #endif // DEBUG_TRACES
}
- complex.assign_filtration(f_simplex, alpha_complex_filtration);
-#ifdef DEBUG_TRACES
- std::clog << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl;
-#endif // DEBUG_TRACES
+ #ifdef DEBUG_TRACES
+ std::clog << std::endl;
+ #endif // DEBUG_TRACES
+ // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma)
+ if (std::isnan(complex.filtration(f_simplex))) {
+ Filtration_value alpha_complex_filtration = 0.0;
+ // No need to compute squared_radius on a single point - alpha is 0.0
+ if (f_simplex_dim > 0) {
+ // squared_radius function initialization
+ Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object();
+
+ CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv;
+ auto sqrad = squared_radius(pointVector.begin(), pointVector.end());
+ #if CGAL_VERSION_NR >= 1050000000
+ if(exact) CGAL::exact(sqrad);
+ #endif
+ alpha_complex_filtration = cv(sqrad);
+ }
+ complex.assign_filtration(f_simplex, alpha_complex_filtration);
+ #ifdef DEBUG_TRACES
+ std::clog << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl;
+ #endif // DEBUG_TRACES
+ }
+ // No need to propagate further, unweighted points all have value 0
+ if (decr_dim > 1)
+ propagate_alpha_filtration(complex, f_simplex);
}
- // No need to propagate further, unweighted points all have value 0
- if (decr_dim > 1)
- propagate_alpha_filtration(complex, f_simplex);
}
}
+ // --------------------------------------------------------------------------------------------
+
+ // --------------------------------------------------------------------------------------------
+ // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension
+ complex.make_filtration_non_decreasing();
+ // Remove all simplices that have a filtration value greater than max_alpha_square
+ complex.prune_above_filtration(max_alpha_square);
+ // --------------------------------------------------------------------------------------------
}
- // --------------------------------------------------------------------------------------------
-
- // --------------------------------------------------------------------------------------------
- // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension
- complex.make_filtration_non_decreasing();
- // Remove all simplices that have a filtration value greater than max_alpha_square
- complex.prune_above_filtration(max_alpha_square);
- // --------------------------------------------------------------------------------------------
return true;
}
diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt
index 0476c6d4..fe4b23e4 100644
--- a/src/Alpha_complex/test/CMakeLists.txt
+++ b/src/Alpha_complex/test/CMakeLists.txt
@@ -8,11 +8,15 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
add_executable ( Alpha_complex_test_unit Alpha_complex_unit_test.cpp )
target_link_libraries(Alpha_complex_test_unit ${CGAL_LIBRARY})
+ add_executable ( Delaunay_complex_test_unit Delaunay_complex_unit_test.cpp )
+ target_link_libraries(Delaunay_complex_test_unit ${CGAL_LIBRARY})
if (TBB_FOUND)
target_link_libraries(Alpha_complex_test_unit ${TBB_LIBRARIES})
+ target_link_libraries(Delaunay_complex_test_unit ${TBB_LIBRARIES})
endif()
gudhi_add_boost_test(Alpha_complex_test_unit)
+ gudhi_add_boost_test(Delaunay_complex_test_unit)
add_executable ( Alpha_complex_3d_test_unit Alpha_complex_3d_unit_test.cpp )
target_link_libraries(Alpha_complex_3d_test_unit ${CGAL_LIBRARY})
diff --git a/src/Alpha_complex/test/Delaunay_complex_unit_test.cpp b/src/Alpha_complex/test/Delaunay_complex_unit_test.cpp
new file mode 100644
index 00000000..71164705
--- /dev/null
+++ b/src/Alpha_complex/test/Delaunay_complex_unit_test.cpp
@@ -0,0 +1,72 @@
+/* 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): Vincent Rouvreau
+ *
+ * Copyright (C) 2020 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "delaunay_complex"
+#include <boost/test/unit_test.hpp>
+#include <boost/mpl/list.hpp>
+
+#include <CGAL/Epick_d.h>
+#include <CGAL/Epeck_d.h>
+
+#include <vector>
+#include <limits> // NaN
+#include <cmath>
+
+#include <gudhi/Alpha_complex.h>
+// to construct a simplex_tree from Delaunay_triangulation
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Unitary_tests_utils.h>
+#include <gudhi/random_point_generators.h>
+
+// Use dynamic_dimension_tag for the user to be able to set dimension
+typedef CGAL::Epeck_d< CGAL::Dynamic_dimension_tag > Exact_kernel_d;
+// Use static dimension_tag for the user not to be able to set dimension
+typedef CGAL::Epeck_d< CGAL::Dimension_tag<5> > Exact_kernel_s;
+// Use dynamic_dimension_tag for the user to be able to set dimension
+typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Inexact_kernel_d;
+// Use static dimension_tag for the user not to be able to set dimension
+typedef CGAL::Epick_d< CGAL::Dimension_tag<5> > Inexact_kernel_s;
+// The triangulation uses the default instantiation of the TriangulationDataStructure template parameter
+
+typedef boost::mpl::list<Exact_kernel_d, Exact_kernel_s, Inexact_kernel_d, Inexact_kernel_s> list_of_kernel_variants;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of_kernel_variants) {
+ std::cout << "*****************************************************************************************************";
+ using Point = typename TestedKernel::Point_d;
+ std::vector<Point> points;
+ // 50 points on a 4-sphere
+ points = Gudhi::generate_points_on_sphere_d<TestedKernel>(10, 5, 1.);
+
+ Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex(points);
+
+ // Alpha complex
+ Gudhi::Simplex_tree<> stree_from_alpha_complex;
+ BOOST_CHECK(alpha_complex.create_complex(stree_from_alpha_complex));
+ stree_from_alpha_complex.initialize_filtration();
+
+ // Delaunay complex
+ Gudhi::Simplex_tree<> stree_from_delaunay_complex;
+ BOOST_CHECK(alpha_complex.create_complex(stree_from_delaunay_complex, 0., false, true));
+
+ // Check all the simplices from alpha complex are in the Delaunay complex
+ for (auto f_simplex : stree_from_alpha_complex.filtration_simplex_range()) {
+ std::vector<Gudhi::Simplex_tree<>::Vertex_handle> simplex;
+ for (Gudhi::Simplex_tree<>::Vertex_handle vertex : stree_from_alpha_complex.simplex_vertex_range(f_simplex)) {
+ std::cout << "(" << vertex << ")";
+ simplex.push_back(vertex);
+ }
+ std::cout << std::endl;
+ Gudhi::Simplex_tree<>::Simplex_handle sh = stree_from_delaunay_complex.find(simplex);
+ BOOST_CHECK(std::isnan(stree_from_delaunay_complex.filtration(sh)));
+ BOOST_CHECK(sh != stree_from_delaunay_complex.null_simplex());
+ }
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-28 16:08:23 +0000