summaryrefslogtreecommitdiff
path: root/src/Alpha_complex
diff options
context:
space:
mode:
Diffstat (limited to 'src/Alpha_complex')
-rw-r--r--src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp4
-rw-r--r--src/Alpha_complex/doc/Intro_alpha_complex.h31
-rw-r--r--src/Alpha_complex/doc/alpha_complex_representation.ipe6
-rw-r--r--src/Alpha_complex/doc/alpha_complex_representation.pngbin14606 -> 19568 bytes
-rw-r--r--src/Alpha_complex/example/Alpha_complex_from_off.cpp13
-rw-r--r--src/Alpha_complex/example/Alpha_complex_from_points.cpp5
-rw-r--r--src/Alpha_complex/example/CMakeLists.txt16
-rw-r--r--src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp65
-rw-r--r--src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp12
-rw-r--r--src/Alpha_complex/include/gudhi/Alpha_complex.h75
-rw-r--r--src/Alpha_complex/include/gudhi/Alpha_complex_3d.h81
-rw-r--r--src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp53
-rw-r--r--src/Alpha_complex/test/Alpha_complex_unit_test.cpp73
-rw-r--r--src/Alpha_complex/test/CMakeLists.txt22
-rw-r--r--src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp37
-rw-r--r--src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp50
-rw-r--r--src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp24
-rw-r--r--src/Alpha_complex/utilities/CMakeLists.txt30
-rw-r--r--src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp4
-rw-r--r--src/Alpha_complex/utilities/alpha_complex_persistence.cpp107
-rw-r--r--src/Alpha_complex/utilities/alphacomplex.md2
21 files changed, 498 insertions, 212 deletions
diff --git a/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp b/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp
index 005a712a..99ad94b9 100644
--- a/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp
+++ b/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp
@@ -115,7 +115,7 @@ void benchmark_weighted_points_on_torus_3D(const std::string& msg) {
std::cout << " Alpha complex 3d on torus with " << nb_points << " points." << std::endl;
std::vector<K::Point_d> points_on_torus = Gudhi::generate_points_on_torus_3D<K>(nb_points, 1.0, 0.5);
- using Point = typename Weighted_alpha_complex_3d::Point_3;
+ using Point = typename Weighted_alpha_complex_3d::Bare_point_3;
using Weighted_point = typename Weighted_alpha_complex_3d::Weighted_point_3;
std::vector<Weighted_point> points;
@@ -206,7 +206,7 @@ void benchmark_weighted_periodic_points(const std::string& msg) {
std::cout << " Weighted periodic alpha complex 3d with " << nb_points * nb_points * nb_points << " points."
<< std::endl;
- using Point = typename Weighted_periodic_alpha_complex_3d::Point_3;
+ using Point = typename Weighted_periodic_alpha_complex_3d::Bare_point_3;
using Weighted_point = typename Weighted_periodic_alpha_complex_3d::Weighted_point_3;
std::vector<Weighted_point> points;
diff --git a/src/Alpha_complex/doc/Intro_alpha_complex.h b/src/Alpha_complex/doc/Intro_alpha_complex.h
index b075d1fc..a8b1a106 100644
--- a/src/Alpha_complex/doc/Intro_alpha_complex.h
+++ b/src/Alpha_complex/doc/Intro_alpha_complex.h
@@ -31,26 +31,37 @@ namespace alpha_complex {
* circumsphere is empty (the simplex is then said to be Gabriel), and as the minimum of the filtration
* values of the codimension 1 cofaces that make it not Gabriel otherwise.
*
- * All simplices that have a filtration value strictly greater than a given alpha squared value are not inserted into
- * the complex.
+ * All simplices that have a filtration value \f$ > \alpha^2 \f$ are removed from the Delaunay complex
+ * when creating the simplicial complex if it is specified.
*
* \image html "alpha_complex_representation.png" "Alpha-complex representation"
*
* Alpha_complex is constructing a <a target="_blank"
* href="http://doc.cgal.org/latest/Triangulation/index.html#Chapter_Triangulations">Delaunay Triangulation</a>
- * \cite cgal:hdj-t-15b from <a target="_blank" href="http://www.cgal.org/">CGAL</a> (the Computational Geometry
- * Algorithms Library \cite cgal:eb-15b) and is able to create a `SimplicialComplexForAlpha`.
+ * \cite cgal:hdj-t-19b from <a target="_blank" href="http://www.cgal.org/">CGAL</a> (the Computational Geometry
+ * Algorithms Library \cite cgal:eb-19b) and is able to create a `SimplicialComplexForAlpha`.
*
* The complex is a template class requiring an Epick_d <a target="_blank"
* href="http://doc.cgal.org/latest/Kernel_d/index.html#Chapter_dD_Geometry_Kernel">dD Geometry Kernel</a>
- * \cite cgal:s-gkd-15b from CGAL as template parameter.
+ * \cite cgal:s-gkd-19b from CGAL as template parameter.
*
* \remark
- * - When the simplicial complex is constructed with an infinite value of alpha, the complex is a Delaunay
- * complex.
+ * - 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).
* - 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">
+ * 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
+ * exact, but the computation of filtration values can exceptionally be arbitrarily bad. In all cases, we still
+ * guarantee that the output is a valid filtration (faces have a filtration value no larger than their cofaces).
+ * - For performances reasons, it is advised to use `Alpha_complex` with \ref cgal &ge; 5.0.0.
*
* \section pointsexample Example from points
*
@@ -124,13 +135,13 @@ namespace alpha_complex {
*
* \subsubsection nondecreasing Non decreasing filtration values
*
- * As the squared radii computed by CGAL are an approximation, it might happen that these alpha squared values do not
- * quite define a proper filtration (i.e. non-decreasing with respect to inclusion).
+ * As the squared radii computed by CGAL are an approximation, it might happen that these \f$ \alpha^2 \f$ values do
+ * not quite define a proper filtration (i.e. non-decreasing with respect to inclusion).
* We fix that up by calling `SimplicialComplexForAlpha::make_filtration_non_decreasing()`.
*
* \subsubsection pruneabove Prune above given filtration value
*
- * The simplex tree is pruned from the given maximum alpha squared value (cf.
+ * The simplex tree is pruned from the given maximum \f$ \alpha^2 \f$ value (cf.
* `SimplicialComplexForAlpha::prune_above_filtration()`).
* In the following example, the value is given by the user as argument of the program.
*
diff --git a/src/Alpha_complex/doc/alpha_complex_representation.ipe b/src/Alpha_complex/doc/alpha_complex_representation.ipe
index e8096b93..40ff1d0f 100644
--- a/src/Alpha_complex/doc/alpha_complex_representation.ipe
+++ b/src/Alpha_complex/doc/alpha_complex_representation.ipe
@@ -1,7 +1,7 @@
<?xml version="1.0"?>
<!DOCTYPE ipe SYSTEM "ipe.dtd">
-<ipe version="70107" creator="Ipe 7.1.10">
-<info created="D:20150603143945" modified="D:20160404172133"/>
+<ipe version="70206" creator="Ipe 7.2.7">
+<info created="D:20150603143945" modified="D:20200110100102"/>
<ipestyle name="basic">
<symbol name="arrow/arc(spx)">
<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen">
@@ -305,7 +305,7 @@ h
108.275 743.531 m
166.45 743.531 l
</path>
-<text matrix="1 0 0 1 142.618 -109.867" transformations="translations" pos="127.397 746.763" stroke="darkgray" type="label" width="6.41" height="4.289" depth="0" valign="baseline">$\alpha$</text>
+<text matrix="1 0 0 1 126.618 -109.867" transformations="translations" pos="127.397 746.763" stroke="darkgray" type="label" width="45.707" height="9.041" depth="1.32" valign="baseline" style="math">\alpha = \sqrt{32.0}</text>
<use matrix="1 0 0 1 -209.478 12.0238" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="black" fill="white"/>
<use matrix="1 0 0 1 -210.178 22.1775" name="mark/fdisk(sfx)" pos="280 660" size="normal" stroke="black" fill="white"/>
<use matrix="1 0 0 1 -210.178 22.1775" name="mark/fdisk(sfx)" pos="370 690" size="normal" stroke="black" fill="white"/>
diff --git a/src/Alpha_complex/doc/alpha_complex_representation.png b/src/Alpha_complex/doc/alpha_complex_representation.png
index 7b81cd69..5ebb1e75 100644
--- a/src/Alpha_complex/doc/alpha_complex_representation.png
+++ b/src/Alpha_complex/doc/alpha_complex_representation.png
Binary files differ
diff --git a/src/Alpha_complex/example/Alpha_complex_from_off.cpp b/src/Alpha_complex/example/Alpha_complex_from_off.cpp
index d411e90a..220a66de 100644
--- a/src/Alpha_complex/example/Alpha_complex_from_off.cpp
+++ b/src/Alpha_complex/example/Alpha_complex_from_off.cpp
@@ -2,8 +2,6 @@
// to construct a simplex_tree from alpha complex
#include <gudhi/Simplex_tree.h>
-#include <CGAL/Epick_d.h>
-
#include <iostream>
#include <string>
@@ -23,22 +21,21 @@ int main(int argc, char **argv) {
// ----------------------------------------------------------------------------
// Init of an alpha complex from an OFF file
// ----------------------------------------------------------------------------
- using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >;
- Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_name);
+ Gudhi::alpha_complex::Alpha_complex<> alpha_complex_from_file(off_file_name);
- std::streambuf* streambufffer;
+ std::streambuf* streambuffer;
std::ofstream ouput_file_stream;
if (argc == 4) {
ouput_file_stream.open(std::string(argv[3]));
- streambufffer = ouput_file_stream.rdbuf();
+ streambuffer = ouput_file_stream.rdbuf();
} else {
- streambufffer = std::cout.rdbuf();
+ streambuffer = std::cout.rdbuf();
}
Gudhi::Simplex_tree<> simplex;
if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) {
- std::ostream output_stream(streambufffer);
+ std::ostream output_stream(streambuffer);
// ----------------------------------------------------------------------------
// Display information about the alpha complex
diff --git a/src/Alpha_complex/example/Alpha_complex_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_from_points.cpp
index 981aa470..6526ca3a 100644
--- a/src/Alpha_complex/example/Alpha_complex_from_points.cpp
+++ b/src/Alpha_complex/example/Alpha_complex_from_points.cpp
@@ -2,12 +2,13 @@
// to construct a simplex_tree from alpha complex
#include <gudhi/Simplex_tree.h>
-#include <CGAL/Epick_d.h>
+#include <CGAL/Epeck_d.h>
#include <iostream>
#include <vector>
-using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<2> >;
+// Explicit dimension 2 Epeck_d kernel
+using Kernel = CGAL::Epeck_d< CGAL::Dimension_tag<2> >;
using Point = Kernel::Point_d;
using Vector_of_points = std::vector<Point>;
diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt
index b069b443..b0337934 100644
--- a/src/Alpha_complex/example/CMakeLists.txt
+++ b/src/Alpha_complex/example/CMakeLists.txt
@@ -5,9 +5,12 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
target_link_libraries(Alpha_complex_example_from_points ${CGAL_LIBRARY})
add_executable ( Alpha_complex_example_from_off Alpha_complex_from_off.cpp )
target_link_libraries(Alpha_complex_example_from_off ${CGAL_LIBRARY})
+ add_executable ( Alpha_complex_example_fast_from_off Fast_alpha_complex_from_off.cpp )
+ target_link_libraries(Alpha_complex_example_fast_from_off ${CGAL_LIBRARY})
if (TBB_FOUND)
target_link_libraries(Alpha_complex_example_from_points ${TBB_LIBRARIES})
target_link_libraries(Alpha_complex_example_from_off ${TBB_LIBRARIES})
+ target_link_libraries(Alpha_complex_example_fast_from_off ${TBB_LIBRARIES})
endif()
add_test(NAME Alpha_complex_example_from_points COMMAND $<TARGET_FILE:Alpha_complex_example_from_points>)
@@ -16,7 +19,13 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
"${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "60.0" "${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_60.txt")
add_test(NAME Alpha_complex_example_from_off_32 COMMAND $<TARGET_FILE:Alpha_complex_example_from_off>
"${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "32.0" "${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_32.txt")
- if (DIFF_PATH)
+
+ add_test(NAME Alpha_complex_example_fast_from_off_60 COMMAND $<TARGET_FILE:Alpha_complex_example_fast_from_off>
+ "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "60.0" "${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_60.txt")
+ add_test(NAME Alpha_complex_example_fast_from_off_32 COMMAND $<TARGET_FILE:Alpha_complex_example_fast_from_off>
+ "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "32.0" "${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_32.txt")
+
+if (DIFF_PATH)
# Do not forget to copy test results files in current binary dir
file(COPY "alphaoffreader_for_doc_32.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
file(COPY "alphaoffreader_for_doc_60.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
@@ -25,6 +34,11 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_60.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_60.txt)
add_test(Alpha_complex_example_from_off_32_diff_files ${DIFF_PATH}
${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt)
+
+ add_test(Alpha_complex_example_fast_from_off_60_diff_files ${DIFF_PATH}
+ ${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_60.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_60.txt)
+ add_test(Alpha_complex_example_fast_from_off_32_diff_files ${DIFF_PATH}
+ ${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt)
endif()
add_executable ( Alpha_complex_example_weighted_3d_from_points Weighted_alpha_complex_3d_from_points.cpp )
diff --git a/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp b/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp
new file mode 100644
index 00000000..f181005a
--- /dev/null
+++ b/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp
@@ -0,0 +1,65 @@
+#include <gudhi/Alpha_complex.h>
+// to construct a simplex_tree from alpha complex
+#include <gudhi/Simplex_tree.h>
+
+#include <CGAL/Epick_d.h>
+
+#include <iostream>
+#include <string>
+
+void usage(int nbArgs, char * const progName) {
+ std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n";
+ std::cerr << "Usage: " << progName << " filename.off alpha_square_max_value [ouput_file.txt]\n";
+ std::cerr << " i.e.: " << progName << " ../../data/points/alphacomplexdoc.off 60.0\n";
+ exit(-1); // ----- >>
+}
+
+int main(int argc, char **argv) {
+ if ((argc != 3) && (argc != 4)) usage(argc, (argv[0] - 1));
+
+ std::string off_file_name {argv[1]};
+ double alpha_square_max_value {atof(argv[2])};
+
+ // WARNING : CGAL::Epick_d is fast but not safe (unlike CGAL::Epeck_d)
+ // (i.e. when the points are on a grid)
+ using Fast_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
+ // ----------------------------------------------------------------------------
+ // Init of an alpha complex from an OFF file
+ // ----------------------------------------------------------------------------
+ Gudhi::alpha_complex::Alpha_complex<Fast_kernel> alpha_complex_from_file(off_file_name);
+
+ std::streambuf* streambuffer;
+ std::ofstream ouput_file_stream;
+
+ if (argc == 4) {
+ ouput_file_stream.open(std::string(argv[3]));
+ streambuffer = ouput_file_stream.rdbuf();
+ } else {
+ streambuffer = std::cout.rdbuf();
+ }
+
+ Gudhi::Simplex_tree<> simplex;
+ if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) {
+ std::ostream output_stream(streambuffer);
+
+ // ----------------------------------------------------------------------------
+ // Display information about the alpha complex
+ // ----------------------------------------------------------------------------
+ output_stream << "Alpha complex is of dimension " << simplex.dimension() <<
+ " - " << simplex.num_simplices() << " simplices - " <<
+ simplex.num_vertices() << " vertices." << std::endl;
+
+ output_stream << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" <<
+ std::endl;
+ for (auto f_simplex : simplex.filtration_simplex_range()) {
+ output_stream << " ( ";
+ for (auto vertex : simplex.simplex_vertex_range(f_simplex)) {
+ output_stream << vertex << " ";
+ }
+ output_stream << ") -> " << "[" << simplex.filtration(f_simplex) << "] ";
+ output_stream << std::endl;
+ }
+ }
+ ouput_file_stream.close();
+ return 0;
+}
diff --git a/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp b/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp
index ac11b68c..fcf80802 100644
--- a/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp
+++ b/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp
@@ -10,7 +10,7 @@
// Complexity = FAST, weighted = true, periodic = false
using Weighted_alpha_complex_3d =
Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, false>;
-using Point = Weighted_alpha_complex_3d::Point_3;
+using Bare_point = Weighted_alpha_complex_3d::Bare_point_3;
using Weighted_point = Weighted_alpha_complex_3d::Weighted_point_3;
int main(int argc, char **argv) {
@@ -18,11 +18,11 @@ int main(int argc, char **argv) {
// Init of a list of points and weights from a small molecule
// ----------------------------------------------------------------------------
std::vector<Weighted_point> weighted_points;
- weighted_points.push_back(Weighted_point(Point(1, -1, -1), 4.));
- weighted_points.push_back(Weighted_point(Point(-1, 1, -1), 4.));
- weighted_points.push_back(Weighted_point(Point(-1, -1, 1), 4.));
- weighted_points.push_back(Weighted_point(Point(1, 1, 1), 4.));
- weighted_points.push_back(Weighted_point(Point(2, 2, 2), 1.));
+ weighted_points.push_back(Weighted_point(Bare_point(1, -1, -1), 4.));
+ weighted_points.push_back(Weighted_point(Bare_point(-1, 1, -1), 4.));
+ weighted_points.push_back(Weighted_point(Bare_point(-1, -1, 1), 4.));
+ weighted_points.push_back(Weighted_point(Bare_point(1, 1, 1), 4.));
+ weighted_points.push_back(Weighted_point(Bare_point(2, 2, 2), 1.));
// ----------------------------------------------------------------------------
// Init of an alpha complex from the list of points
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h
index 8919cdb9..f2a05e95 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h
@@ -20,11 +20,12 @@
#include <math.h> // isnan, fmax
#include <CGAL/Delaunay_triangulation.h>
-#include <CGAL/Epick_d.h>
+#include <CGAL/Epeck_d.h> // For EXACT or SAFE version
+#include <CGAL/Epick_d.h> // For FAST version
#include <CGAL/Spatial_sort_traits_adapter_d.h>
#include <CGAL/property_map.h> // for CGAL::Identity_property_map
-#include <CGAL/NT_converter.h>
#include <CGAL/version.h> // for CGAL_VERSION_NR
+#include <CGAL/NT_converter.h>
#include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST
@@ -39,17 +40,20 @@
// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
#if CGAL_VERSION_NR < 1041101000
-# error Alpha_complex_3d is only available for CGAL >= 4.11
+# error Alpha_complex is only available for CGAL >= 4.11
#endif
#if !EIGEN_VERSION_AT_LEAST(3,1,0)
-# error Alpha_complex_3d is only available for Eigen3 >= 3.1.0 installed with CGAL
+# error Alpha_complex is only available for Eigen3 >= 3.1.0 installed with CGAL
#endif
namespace Gudhi {
namespace alpha_complex {
+template<typename D> struct Is_Epeck_D { static const bool value = false; };
+template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool value = true; };
+
/**
* \class Alpha_complex Alpha_complex.h gudhi/Alpha_complex.h
* \brief Alpha complex data structure.
@@ -63,17 +67,31 @@ namespace alpha_complex {
*
* Please refer to \ref alpha_complex for examples.
*
- * The complex is a template class requiring an Epick_d <a target="_blank"
+ * The complex is a template class requiring an <a target="_blank"
+ * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>,
+ * or an <a target="_blank"
+ * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epick__d.html">CGAL::Epick_d</a> <a target="_blank"
* href="http://doc.cgal.org/latest/Kernel_d/index.html#Chapter_dD_Geometry_Kernel">dD Geometry Kernel</a>
- * \cite cgal:s-gkd-15b from CGAL as template, default value is <a target="_blank"
- * href="http://doc.cgal.org/latest/Kernel_d/classCGAL_1_1Epick__d.html">CGAL::Epick_d</a>
+ * \cite cgal:s-gkd-19b from CGAL as template, default value is <a target="_blank"
+ * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>
* < <a target="_blank" href="http://doc.cgal.org/latest/Kernel_23/classCGAL_1_1Dynamic__dimension__tag.html">
* CGAL::Dynamic_dimension_tag </a> >
*
- * \remark When Alpha_complex is constructed with an infinite value of alpha, the complex is a Delaunay complex.
- *
+ * \remark
+ * - When Alpha_complex is constructed with an infinite value of alpha, the complex is a Delaunay complex.
+ * - 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">
+ * 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
+ * exact, but the computation of filtration values can exceptionally be arbitrarily bad. In all cases, we still
+ * guarantee that the output is a valid filtration (faces have a filtration value no larger than their cofaces).
+ * - For performances reasons, it is advised to use `Alpha_complex` with \ref cgal &ge; 5.0.0.
*/
-template<class Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>>
+template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>
class Alpha_complex {
public:
// Add an int in TDS to save point index in the structure
@@ -103,8 +121,8 @@ class Alpha_complex {
// size_type type from CGAL.
typedef typename Delaunay_triangulation::size_type size_type;
- // Map type to switch from simplex tree vertex handle to CGAL vertex iterator.
- typedef typename std::map< std::size_t, CGAL_vertex_iterator > Vector_vertex_iterator;
+ // Structure to switch from simplex tree vertex handle to CGAL vertex iterator.
+ typedef typename std::vector< CGAL_vertex_iterator > Vector_vertex_iterator;
private:
/** \brief Vertex iterator vector to switch from simplex tree vertex handle to CGAL vertex iterator.
@@ -173,17 +191,15 @@ class Alpha_complex {
return vertex_handle_to_iterator_.at(vertex)->point();
}
- /** \brief number_of_vertices returns the number of vertices (same as the number of points).
- *
- * @return The number of vertices.
- */
- std::size_t number_of_vertices() const {
- return vertex_handle_to_iterator_.size();
- }
-
private:
template<typename InputPointRange >
void init_from_range(const InputPointRange& points) {
+ #if CGAL_VERSION_NR < 1050000000
+ if (Is_Epeck_D<Kernel>::value)
+ std::cerr << "It is strongly advised to use a CGAL version >= 5.0 with Epeck_d Kernel for performance reasons."
+ << std::endl;
+ #endif
+
auto first = std::begin(points);
auto last = std::end(points);
@@ -214,14 +230,16 @@ class Alpha_complex {
hint = pos->full_cell();
}
// --------------------------------------------------------------------------------------------
- // double map to retrieve simplex tree vertex handles from CGAL vertex iterator and vice versa
+ // structure to retrieve CGAL points from vertex handle - one vertex handle per point.
+ // Needs to be constructed before as vertex handles arrives in no particular order.
+ vertex_handle_to_iterator_.resize(point_cloud.size());
// Loop on triangulation vertices list
for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) {
if (!triangulation_->is_infinite(*vit)) {
#ifdef DEBUG_TRACES
std::cout << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl;
#endif // DEBUG_TRACES
- vertex_handle_to_iterator_.emplace(vit->data(), vit);
+ vertex_handle_to_iterator_[vit->data()] = vit;
}
}
// --------------------------------------------------------------------------------------------
@@ -237,6 +255,8 @@ class Alpha_complex {
* @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.
+ * @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>.
*
* @return true if creation succeeds, false otherwise.
*
@@ -248,7 +268,8 @@ class Alpha_complex {
template <typename SimplicialComplexForAlpha,
typename Filtration_value = typename SimplicialComplexForAlpha::Filtration_value>
bool create_complex(SimplicialComplexForAlpha& complex,
- Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity()) {
+ Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(),
+ bool exact = 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;
@@ -324,9 +345,13 @@ class Alpha_complex {
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;
- alpha_complex_filtration = cv(squared_radius(pointVector.begin(), pointVector.end()));
+ 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
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
index 13ebb9c1..7f96c94c 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
@@ -43,7 +43,7 @@
#include <vector>
#include <unordered_map>
#include <stdexcept>
-#include <cstddef>
+#include <cstddef> // for std::size_t
#include <memory> // for std::unique_ptr
#include <type_traits> // for std::conditional and std::enable_if
#include <limits> // for numeric_limits<>
@@ -97,7 +97,7 @@ struct Value_from_iterator<complexity::EXACT> {
* \details
* The data structure is constructing a <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html">CGAL 3D Alpha
* Shapes</a> from a range of points (can be read from an OFF file, cf. Points_off_reader).
- * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous.
+ * Duplicate points are inserted once in the Alpha_complex.
*
* \tparam Complexity shall be `Gudhi::alpha_complex::complexity` type. Default value is
* `Gudhi::alpha_complex::complexity::SAFE`.
@@ -225,23 +225,23 @@ class Alpha_complex_3d {
* Must be compatible with double. */
using FT = typename Alpha_shape_3::FT;
- /** \brief Gives public access to the Point_3 type. Here is a Point_3 constructor example:
+ /** \brief Gives public access to the Bare_point_3 (bare aka. unweighed) type.
+ * Here is a Bare_point_3 constructor example:
\code{.cpp}
using Alpha_complex_3d = Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, false, false>;
// x0 = 1., y0 = -1.1, z0 = -1..
-Alpha_complex_3d::Point_3 p0(1., -1.1, -1.);
+Alpha_complex_3d::Bare_point_3 p0(1., -1.1, -1.);
\endcode
* */
- using Point_3 = typename Kernel::Point_3;
+ using Bare_point_3 = typename Kernel::Point_3;
/** \brief Gives public access to the Weighted_point_3 type. A Weighted point can be constructed as follows:
\code{.cpp}
-using Weighted_alpha_complex_3d =
- Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, false>;
+using Weighted_alpha_complex_3d = Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, false>;
// x0 = 1., y0 = -1.1, z0 = -1., weight = 4.
-Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point_3(1., -1.1, -1.), 4.);
+Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_point_3(1., -1.1, -1.), 4.);
\endcode
*
* Note: This type is defined to void if Alpha complex is not weighted.
@@ -249,6 +249,11 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
* */
using Weighted_point_3 = typename Triangulation_3<Kernel, Tds, Weighted, Periodic>::Weighted_point_3;
+ /** \brief `Alpha_complex_3d::Point_3` type is either a `Alpha_complex_3d::Bare_point_3` (Weighted = false) or a
+ * `Alpha_complex_3d::Weighted_point_3` (Weighted = true).
+ */
+ using Point_3 = typename Alpha_shape_3::Point;
+
private:
using Dispatch =
CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<CGAL::Object, FT>,
@@ -264,13 +269,12 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
public:
/** \brief Alpha_complex constructor from a list of points.
*
- * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3` or
- * `Alpha_complex_3d::Weighted_point_3`.
+ * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3`.
*
* @pre Available if Alpha_complex_3d is not Periodic.
*
* The type InputPointRange must be a range for which std::begin and std::end return input iterators on a
- * `Alpha_complex_3d::Point_3` or a `Alpha_complex_3d::Weighted_point_3`.
+ * `Alpha_complex_3d::Point_3`.
*/
template <typename InputPointRange>
Alpha_complex_3d(const InputPointRange& points) {
@@ -284,13 +288,13 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
*
* @exception std::invalid_argument In debug mode, if points and weights do not have the same size.
*
- * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3`.
+ * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Bare_point_3`.
* @param[in] weights Range of weights on points. Weights shall be in double.
*
* @pre Available if Alpha_complex_3d is Weighted and not Periodic.
*
* The type InputPointRange must be a range for which std::begin and
- * std::end return input iterators on a `Alpha_complex_3d::Point_3`.
+ * std::end return input iterators on a `Alpha_complex_3d::Bare_point_3`.
* The type WeightRange must be a range for which std::begin and
* std::end return an input iterator on a double.
*/
@@ -318,8 +322,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
*
* @exception std::invalid_argument In debug mode, if the size of the cuboid in every directions is not the same.
*
- * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3` or
- * `Alpha_complex_3d::Weighted_point_3`.
+ * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3`.
* @param[in] x_min Iso-oriented cuboid x_min.
* @param[in] y_min Iso-oriented cuboid y_min.
* @param[in] z_min Iso-oriented cuboid z_min.
@@ -330,7 +333,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
* @pre Available if Alpha_complex_3d is Periodic.
*
* The type InputPointRange must be a range for which std::begin and std::end return input iterators on a
- * `Alpha_complex_3d::Point_3` or a `Alpha_complex_3d::Weighted_point_3`.
+ * `Alpha_complex_3d::Point_3`.
*
* @note In weighted version, please check weights are greater than zero, and lower than 1/64*cuboid length
* squared.
@@ -366,7 +369,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
* @exception std::invalid_argument In debug mode, if a weight is negative, zero, or greater than 1/64*cuboid length
* squared.
*
- * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3`.
+ * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Bare_point_3`.
* @param[in] weights Range of weights on points. Weights shall be in double.
* @param[in] x_min Iso-oriented cuboid x_min.
* @param[in] y_min Iso-oriented cuboid y_min.
@@ -378,7 +381,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
* @pre Available if Alpha_complex_3d is Weighted and Periodic.
*
* The type InputPointRange must be a range for which std::begin and
- * std::end return input iterators on a `Alpha_complex_3d::Point_3`.
+ * std::end return input iterators on a `Alpha_complex_3d::Bare_point_3`.
* The type WeightRange must be a range for which std::begin and
* std::end return an input iterator on a double.
* The type of x_min, y_min, z_min, x_max, y_max and z_max must be a double.
@@ -452,9 +455,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
return false; // ----- >>
}
- // using Filtration_value = typename SimplicialComplexForAlpha3d::Filtration_value;
using Complex_vertex_handle = typename SimplicialComplexForAlpha3d::Vertex_handle;
- using Alpha_shape_simplex_tree_map = std::unordered_map<Alpha_vertex_handle, Complex_vertex_handle>;
using Simplex_tree_vector_vertex = std::vector<Complex_vertex_handle>;
#ifdef DEBUG_TRACES
@@ -474,7 +475,6 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
std::cout << "filtration_with_alpha_values returns : " << objects.size() << " objects" << std::endl;
#endif // DEBUG_TRACES
- Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
using Alpha_value_iterator = typename std::vector<FT>::const_iterator;
Alpha_value_iterator alpha_value_iterator = alpha_values.begin();
for (auto object_iterator : objects) {
@@ -484,7 +484,8 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
for (auto i = 0; i < 4; i++) {
#ifdef DEBUG_TRACES
- std::cout << "from cell[" << i << "]=" << (*cell)->vertex(i)->point() << std::endl;
+ std::cout << "from cell[" << i << "] - Point coordinates (" << (*cell)->vertex(i)->point() << ")"
+ << std::endl;
#endif // DEBUG_TRACES
vertex_list.push_back((*cell)->vertex(i));
}
@@ -495,7 +496,8 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
for (auto i = 0; i < 4; i++) {
if ((*facet).second != i) {
#ifdef DEBUG_TRACES
- std::cout << "from facet=[" << i << "]" << (*facet).first->vertex(i)->point() << std::endl;
+ std::cout << "from facet=[" << i << "] - Point coordinates (" << (*facet).first->vertex(i)->point() << ")"
+ << std::endl;
#endif // DEBUG_TRACES
vertex_list.push_back((*facet).first->vertex(i));
}
@@ -506,7 +508,8 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
} else if (const Edge* edge = CGAL::object_cast<Edge>(&object_iterator)) {
for (auto i : {(*edge).second, (*edge).third}) {
#ifdef DEBUG_TRACES
- std::cout << "from edge[" << i << "]=" << (*edge).first->vertex(i)->point() << std::endl;
+ std::cout << "from edge[" << i << "] - Point coordinates (" << (*edge).first->vertex(i)->point() << ")"
+ << std::endl;
#endif // DEBUG_TRACES
vertex_list.push_back((*edge).first->vertex(i));
}
@@ -516,7 +519,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
} else if (const Alpha_vertex_handle* vertex = CGAL::object_cast<Alpha_vertex_handle>(&object_iterator)) {
#ifdef DEBUG_TRACES
count_vertices++;
- std::cout << "from vertex=" << (*vertex)->point() << std::endl;
+ std::cout << "from vertex - Point coordinates (" << (*vertex)->point() << ")" << std::endl;
#endif // DEBUG_TRACES
vertex_list.push_back((*vertex));
}
@@ -528,7 +531,8 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
// alpha shape not found
Complex_vertex_handle vertex = map_cgal_simplex_tree.size();
#ifdef DEBUG_TRACES
- std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
+ std::cout << "Point (" << the_alpha_shape_vertex->point() << ") not found - insert new vertex id " << vertex
+ << std::endl;
#endif // DEBUG_TRACES
the_simplex.push_back(vertex);
map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex);
@@ -536,7 +540,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
// alpha shape found
Complex_vertex_handle vertex = the_map_iterator->second;
#ifdef DEBUG_TRACES
- std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
+ std::cout << "Point (" << the_alpha_shape_vertex->point() << ") found as vertex id " << vertex << std::endl;
#endif // DEBUG_TRACES
the_simplex.push_back(vertex);
}
@@ -567,9 +571,32 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point
return true;
}
+ /** \brief get_point returns the point corresponding to the vertex given as parameter.
+ *
+ * @param[in] vertex Vertex handle of the point to retrieve.
+ * @return The point found.
+ * @exception std::out_of_range In case vertex is not found (cf. std::vector::at).
+ */
+ const Point_3& get_point(std::size_t vertex) {
+ if (map_cgal_simplex_tree.size() != cgal_vertex_iterator_vector.size()) {
+ cgal_vertex_iterator_vector.resize(map_cgal_simplex_tree.size());
+ for (auto map_iterator : map_cgal_simplex_tree) {
+ cgal_vertex_iterator_vector[map_iterator.second] = map_iterator.first;
+ }
+
+ }
+ auto cgal_vertex_iterator = cgal_vertex_iterator_vector.at(vertex);
+ return cgal_vertex_iterator->point();
+ }
+
private:
// use of a unique_ptr on cgal Alpha_shape_3, as copy and default constructor is not available - no need to be freed
std::unique_ptr<Alpha_shape_3> alpha_shape_3_ptr_;
+
+ // Map type to switch from CGAL vertex iterator to simplex tree vertex handle.
+ std::unordered_map<Alpha_vertex_handle, std::size_t> map_cgal_simplex_tree;
+ // Vector type to switch from simplex tree vertex handle to CGAL vertex iterator.
+ std::vector<Alpha_vertex_handle> cgal_vertex_iterator_vector;
};
} // namespace alpha_complex
diff --git a/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp
index 1102838a..cd698a27 100644
--- a/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp
+++ b/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp
@@ -56,21 +56,52 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) {
// -----------------
std::cout << "Fast alpha complex 3d" << std::endl;
- Fast_alpha_complex_3d alpha_complex(get_points<Fast_alpha_complex_3d::Point_3>());
+ std::vector<Fast_alpha_complex_3d::Bare_point_3> points = get_points<Fast_alpha_complex_3d::Bare_point_3>();
+ Fast_alpha_complex_3d alpha_complex(points);
Gudhi::Simplex_tree<> stree;
alpha_complex.create_complex(stree);
+ for (std::size_t index = 0; index < points.size(); index++) {
+ bool found = false;
+ for (auto point : points) {
+ if (point == alpha_complex.get_point(index)) {
+ found = true;
+ break;
+ }
+ }
+ // Check all points from alpha complex are found in the input point cloud
+ BOOST_CHECK(found);
+ }
+ // Exception if we go out of range
+ BOOST_CHECK_THROW(alpha_complex.get_point(points.size()), std::out_of_range);
+
// -----------------
// Exact version
// -----------------
std::cout << "Exact alpha complex 3d" << std::endl;
- Exact_alpha_complex_3d exact_alpha_complex(get_points<Exact_alpha_complex_3d::Point_3>());
+ std::vector<Exact_alpha_complex_3d::Bare_point_3> exact_points = get_points<Exact_alpha_complex_3d::Bare_point_3>();
+ Exact_alpha_complex_3d exact_alpha_complex(exact_points);
Gudhi::Simplex_tree<> exact_stree;
exact_alpha_complex.create_complex(exact_stree);
+ for (std::size_t index = 0; index < exact_points.size(); index++) {
+ bool found = false;
+ Exact_alpha_complex_3d::Bare_point_3 ap = exact_alpha_complex.get_point(index);
+ for (auto point : points) {
+ if ((point.x() == ap.x()) && (point.y() == ap.y()) && (point.z() == ap.z())) {
+ found = true;
+ break;
+ }
+ }
+ // Check all points from alpha complex are found in the input point cloud
+ BOOST_CHECK(found);
+ }
+ // Exception if we go out of range
+ BOOST_CHECK_THROW(exact_alpha_complex.get_point(exact_points.size()), std::out_of_range);
+
// ---------------------
// Compare both versions
// ---------------------
@@ -110,11 +141,27 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) {
// -----------------
std::cout << "Safe alpha complex 3d" << std::endl;
- Safe_alpha_complex_3d safe_alpha_complex(get_points<Safe_alpha_complex_3d::Point_3>());
+ std::vector<Safe_alpha_complex_3d::Bare_point_3> safe_points = get_points<Safe_alpha_complex_3d::Bare_point_3>();
+ Safe_alpha_complex_3d safe_alpha_complex(safe_points);
Gudhi::Simplex_tree<> safe_stree;
safe_alpha_complex.create_complex(safe_stree);
+ for (std::size_t index = 0; index < safe_points.size(); index++) {
+ bool found = false;
+ Safe_alpha_complex_3d::Bare_point_3 ap = safe_alpha_complex.get_point(index);
+ for (auto point : points) {
+ if ((point.x() == ap.x()) && (point.y() == ap.y()) && (point.z() == ap.z())) {
+ found = true;
+ break;
+ }
+ }
+ // Check all points from alpha complex are found in the input point cloud
+ BOOST_CHECK(found);
+ }
+ // Exception if we go out of range
+ BOOST_CHECK_THROW(safe_alpha_complex.get_point(safe_points.size()), std::out_of_range);
+
// ---------------------
// Compare both versions
// ---------------------
diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
index 01e4cee3..27b671dd 100644
--- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
+++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
@@ -15,6 +15,7 @@
#include <CGAL/Delaunay_triangulation.h>
#include <CGAL/Epick_d.h>
+#include <CGAL/Epeck_d.h>
#include <cmath> // float comparison
#include <limits>
@@ -28,12 +29,16 @@
#include <gudhi/Unitary_tests_utils.h>
// Use dynamic_dimension_tag for the user to be able to set dimension
-typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel_d;
+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::Epick_d< CGAL::Dimension_tag<3> > Kernel_s;
+typedef CGAL::Epeck_d< CGAL::Dimension_tag<3> > 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<3> > Inexact_kernel_s;
// The triangulation uses the default instantiation of the TriangulationDataStructure template parameter
-typedef boost::mpl::list<Kernel_d, Kernel_s> list_of_kernel_variants;
+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) {
// ----------------------------------------------------------------------------
@@ -48,20 +53,12 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of
Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_file(off_file_name);
- std::cout << "alpha_complex_from_points.number_of_vertices()=" << alpha_complex_from_file.number_of_vertices()
- << std::endl;
- BOOST_CHECK(alpha_complex_from_file.number_of_vertices() == 7);
-
Gudhi::Simplex_tree<> simplex_tree_60;
BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_60, max_alpha_square_value));
std::cout << "simplex_tree_60.dimension()=" << simplex_tree_60.dimension() << std::endl;
BOOST_CHECK(simplex_tree_60.dimension() == 2);
- std::cout << "alpha_complex_from_points.number_of_vertices()=" << alpha_complex_from_file.number_of_vertices()
- << std::endl;
- BOOST_CHECK(alpha_complex_from_file.number_of_vertices() == 7);
-
std::cout << "simplex_tree_60.num_vertices()=" << simplex_tree_60.num_vertices() << std::endl;
BOOST_CHECK(simplex_tree_60.num_vertices() == 7);
@@ -86,7 +83,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of
}
// Use static dimension_tag for the user not to be able to set dimension
-typedef CGAL::Epick_d< CGAL::Dimension_tag<4> > Kernel_4;
+typedef CGAL::Epeck_d< CGAL::Dimension_tag<4> > Kernel_4;
typedef Kernel_4::Point_d Point_4;
typedef std::vector<Point_4> Vector_4_Points;
@@ -123,10 +120,6 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
Gudhi::Simplex_tree<> simplex_tree;
BOOST_CHECK(alpha_complex_from_points.create_complex(simplex_tree));
- std::cout << "alpha_complex_from_points.number_of_vertices()=" << alpha_complex_from_points.number_of_vertices()
- << std::endl;
- BOOST_CHECK(alpha_complex_from_points.number_of_vertices() == points.size());
-
// Another way to check num_simplices
std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
int num_simplices = 0;
@@ -146,7 +139,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
std::cout << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl;
BOOST_CHECK(simplex_tree.dimension() == 3);
std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl;
- BOOST_CHECK(simplex_tree.num_vertices() == 4);
+ BOOST_CHECK(simplex_tree.num_vertices() == points.size());
for (auto f_simplex : simplex_tree.filtration_simplex_range()) {
switch (simplex_tree.dimension(f_simplex)) {
@@ -256,10 +249,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_empty_points, TestedKernel, lis
Gudhi::Simplex_tree<> simplex_tree;
BOOST_CHECK(!alpha_complex_from_points.create_complex(simplex_tree));
- std::cout << "alpha_complex_from_points.number_of_vertices()=" << alpha_complex_from_points.number_of_vertices()
- << std::endl;
- BOOST_CHECK(alpha_complex_from_points.number_of_vertices() == points.size());
-
std::cout << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl;
BOOST_CHECK(simplex_tree.num_simplices() == 0);
@@ -267,5 +256,45 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_empty_points, TestedKernel, lis
BOOST_CHECK(simplex_tree.dimension() == -1);
std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl;
- BOOST_CHECK(simplex_tree.num_vertices() == 0);
+ BOOST_CHECK(simplex_tree.num_vertices() == points.size());
+}
+
+using Inexact_kernel_2 = CGAL::Epick_d< CGAL::Dimension_tag<2> >;
+using Exact_kernel_2 = CGAL::Epeck_d< CGAL::Dimension_tag<2> >;
+using list_of_kernel_2_variants = boost::mpl::list<Inexact_kernel_2, Exact_kernel_2>;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_with_duplicated_points, TestedKernel, list_of_kernel_2_variants) {
+ std::cout << "========== Alpha_complex_with_duplicated_points ==========" << std::endl;
+
+ using Point = typename TestedKernel::Point_d;
+ using Vector_of_points = std::vector<Point>;
+
+ // ----------------------------------------------------------------------------
+ // Init of a list of points
+ // ----------------------------------------------------------------------------
+ Vector_of_points points;
+ points.push_back(Point(1.0, 1.0));
+ points.push_back(Point(7.0, 0.0));
+ points.push_back(Point(4.0, 6.0));
+ points.push_back(Point(9.0, 6.0));
+ points.push_back(Point(0.0, 14.0));
+ points.push_back(Point(2.0, 19.0));
+ points.push_back(Point(9.0, 17.0));
+ // duplicated points
+ points.push_back(Point(1.0, 1.0));
+ points.push_back(Point(7.0, 0.0));
+
+ // ----------------------------------------------------------------------------
+ // Init of an alpha complex from the list of points
+ // ----------------------------------------------------------------------------
+ std::cout << "Init" << std::endl;
+ Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_points(points);
+
+ Gudhi::Simplex_tree<> simplex_tree;
+ std::cout << "create_complex" << std::endl;
+ BOOST_CHECK(alpha_complex_from_points.create_complex(simplex_tree));
+
+ std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices()
+ << std::endl;
+ BOOST_CHECK(simplex_tree.num_vertices() < points.size());
}
diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt
index ad5b6314..0476c6d4 100644
--- a/src/Alpha_complex/test/CMakeLists.txt
+++ b/src/Alpha_complex/test/CMakeLists.txt
@@ -1,27 +1,27 @@
project(Alpha_complex_tests)
-include(GUDHI_test_coverage)
+include(GUDHI_boost_test)
if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
# 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}/)
add_executable ( Alpha_complex_test_unit Alpha_complex_unit_test.cpp )
- target_link_libraries(Alpha_complex_test_unit ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+ target_link_libraries(Alpha_complex_test_unit ${CGAL_LIBRARY})
if (TBB_FOUND)
target_link_libraries(Alpha_complex_test_unit ${TBB_LIBRARIES})
endif()
- gudhi_add_coverage_test(Alpha_complex_test_unit)
+ gudhi_add_boost_test(Alpha_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} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+ target_link_libraries(Alpha_complex_3d_test_unit ${CGAL_LIBRARY})
add_executable ( Weighted_alpha_complex_3d_test_unit Weighted_alpha_complex_3d_unit_test.cpp )
- target_link_libraries(Weighted_alpha_complex_3d_test_unit ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+ target_link_libraries(Weighted_alpha_complex_3d_test_unit ${CGAL_LIBRARY})
add_executable ( Periodic_alpha_complex_3d_test_unit Periodic_alpha_complex_3d_unit_test.cpp )
- target_link_libraries(Periodic_alpha_complex_3d_test_unit ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+ target_link_libraries(Periodic_alpha_complex_3d_test_unit ${CGAL_LIBRARY})
add_executable ( Weighted_periodic_alpha_complex_3d_test_unit Weighted_periodic_alpha_complex_3d_unit_test.cpp )
- target_link_libraries(Weighted_periodic_alpha_complex_3d_test_unit ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+ target_link_libraries(Weighted_periodic_alpha_complex_3d_test_unit ${CGAL_LIBRARY})
if (TBB_FOUND)
target_link_libraries(Alpha_complex_3d_test_unit ${TBB_LIBRARIES})
target_link_libraries(Weighted_alpha_complex_3d_test_unit ${TBB_LIBRARIES})
@@ -29,9 +29,9 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
target_link_libraries(Weighted_periodic_alpha_complex_3d_test_unit ${TBB_LIBRARIES})
endif()
- gudhi_add_coverage_test(Alpha_complex_3d_test_unit)
- gudhi_add_coverage_test(Weighted_alpha_complex_3d_test_unit)
- gudhi_add_coverage_test(Periodic_alpha_complex_3d_test_unit)
- gudhi_add_coverage_test(Weighted_periodic_alpha_complex_3d_test_unit)
+ gudhi_add_boost_test(Alpha_complex_3d_test_unit)
+ gudhi_add_boost_test(Weighted_alpha_complex_3d_test_unit)
+ gudhi_add_boost_test(Periodic_alpha_complex_3d_test_unit)
+ gudhi_add_boost_test(Weighted_periodic_alpha_complex_3d_test_unit)
endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
diff --git a/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp
index ac3791a4..731763fa 100644
--- a/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp
+++ b/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp
@@ -44,11 +44,11 @@ typedef boost::mpl::list<Fast_periodic_alpha_complex_3d, Safe_periodic_alpha_com
BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_periodic_throw, Periodic_alpha_complex_3d, periodic_variants_type_list) {
std::cout << "Periodic alpha complex 3d exception throw" << std::endl;
- using Point_3 = typename Periodic_alpha_complex_3d::Point_3;
- std::vector<Point_3> p_points;
+ using Bare_point_3 = typename Periodic_alpha_complex_3d::Bare_point_3;
+ std::vector<Bare_point_3> p_points;
// Not important, this is not what we want to check
- p_points.push_back(Point_3(0.0, 0.0, 0.0));
+ p_points.push_back(Bare_point_3(0.0, 0.0, 0.0));
std::cout << "Check exception throw in debug mode" << std::endl;
// Check it throws an exception when the cuboid is not iso
@@ -73,13 +73,13 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_periodic) {
// ---------------------
std::cout << "Fast periodic alpha complex 3d" << std::endl;
- using Creator = CGAL::Creator_uniform_3<double, Fast_periodic_alpha_complex_3d::Point_3>;
+ using Creator = CGAL::Creator_uniform_3<double, Fast_periodic_alpha_complex_3d::Bare_point_3>;
CGAL::Random random(7);
- CGAL::Random_points_in_cube_3<Fast_periodic_alpha_complex_3d::Point_3, Creator> in_cube(1, random);
- std::vector<Fast_periodic_alpha_complex_3d::Point_3> p_points;
+ CGAL::Random_points_in_cube_3<Fast_periodic_alpha_complex_3d::Bare_point_3, Creator> in_cube(1, random);
+ std::vector<Fast_periodic_alpha_complex_3d::Bare_point_3> p_points;
for (int i = 0; i < 50; i++) {
- Fast_periodic_alpha_complex_3d::Point_3 p = *in_cube++;
+ Fast_periodic_alpha_complex_3d::Bare_point_3 p = *in_cube++;
p_points.push_back(p);
}
@@ -88,15 +88,30 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_periodic) {
Gudhi::Simplex_tree<> stree;
periodic_alpha_complex.create_complex(stree);
+ for (std::size_t index = 0; index < p_points.size(); index++) {
+ bool found = false;
+ Fast_periodic_alpha_complex_3d::Bare_point_3 ap = periodic_alpha_complex.get_point(index);
+ for (auto point : p_points) {
+ if ((point.x() == ap.x()) && (point.y() == ap.y()) && (point.z() == ap.z())) {
+ found = true;
+ break;
+ }
+ }
+ // Check all points from alpha complex are found in the input point cloud
+ BOOST_CHECK(found);
+ }
+ // Exception if we go out of range
+ BOOST_CHECK_THROW(periodic_alpha_complex.get_point(p_points.size()), std::out_of_range);
+
// ----------------------
// Exact periodic version
// ----------------------
std::cout << "Exact periodic alpha complex 3d" << std::endl;
- std::vector<Exact_periodic_alpha_complex_3d::Point_3> e_p_points;
+ std::vector<Exact_periodic_alpha_complex_3d::Bare_point_3> e_p_points;
for (auto p : p_points) {
- e_p_points.push_back(Exact_periodic_alpha_complex_3d::Point_3(p[0], p[1], p[2]));
+ e_p_points.push_back(Exact_periodic_alpha_complex_3d::Bare_point_3(p[0], p[1], p[2]));
}
Exact_periodic_alpha_complex_3d exact_alpha_complex(e_p_points, -1., -1., -1., 1., 1., 1.);
@@ -142,10 +157,10 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_periodic) {
// ----------------------
std::cout << "Safe periodic alpha complex 3d" << std::endl;
- std::vector<Safe_periodic_alpha_complex_3d::Point_3> s_p_points;
+ std::vector<Safe_periodic_alpha_complex_3d::Bare_point_3> s_p_points;
for (auto p : p_points) {
- s_p_points.push_back(Safe_periodic_alpha_complex_3d::Point_3(p[0], p[1], p[2]));
+ s_p_points.push_back(Safe_periodic_alpha_complex_3d::Bare_point_3(p[0], p[1], p[2]));
}
Safe_periodic_alpha_complex_3d safe_alpha_complex(s_p_points, -1., -1., -1., 1., 1., 1.);
diff --git a/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp
index 44deb930..8035f6e8 100644
--- a/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp
+++ b/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp
@@ -43,14 +43,14 @@ typedef boost::mpl::list<Fast_weighted_alpha_complex_3d, Safe_weighted_alpha_com
#ifdef GUDHI_DEBUG
BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_throw, Weighted_alpha_complex_3d, weighted_variants_type_list) {
- using Point_3 = typename Weighted_alpha_complex_3d::Point_3;
- std::vector<Point_3> w_points;
- w_points.push_back(Point_3(0.0, 0.0, 0.0));
- w_points.push_back(Point_3(0.0, 0.0, 0.2));
- w_points.push_back(Point_3(0.2, 0.0, 0.2));
- // w_points.push_back(Point_3(0.6, 0.6, 0.0));
- // w_points.push_back(Point_3(0.8, 0.8, 0.2));
- // w_points.push_back(Point_3(0.2, 0.8, 0.6));
+ using Bare_point_3 = typename Weighted_alpha_complex_3d::Bare_point_3;
+ std::vector<Bare_point_3> w_points;
+ w_points.push_back(Bare_point_3(0.0, 0.0, 0.0));
+ w_points.push_back(Bare_point_3(0.0, 0.0, 0.2));
+ w_points.push_back(Bare_point_3(0.2, 0.0, 0.2));
+ // w_points.push_back(Bare_point_3(0.6, 0.6, 0.0));
+ // w_points.push_back(Bare_point_3(0.8, 0.8, 0.2));
+ // w_points.push_back(Bare_point_3(0.2, 0.8, 0.6));
// weights size is different from w_points size to make weighted Alpha_complex_3d throw in debug mode
std::vector<double> weights = {0.01, 0.005, 0.006, 0.01, 0.009, 0.001};
@@ -62,14 +62,14 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_throw, Weighted_alpha_compl
BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted, Weighted_alpha_complex_3d, weighted_variants_type_list) {
std::cout << "Weighted alpha complex 3d from points and weights" << std::endl;
- using Point_3 = typename Weighted_alpha_complex_3d::Point_3;
- std::vector<Point_3> w_points;
- w_points.push_back(Point_3(0.0, 0.0, 0.0));
- w_points.push_back(Point_3(0.0, 0.0, 0.2));
- w_points.push_back(Point_3(0.2, 0.0, 0.2));
- w_points.push_back(Point_3(0.6, 0.6, 0.0));
- w_points.push_back(Point_3(0.8, 0.8, 0.2));
- w_points.push_back(Point_3(0.2, 0.8, 0.6));
+ using Bare_point_3 = typename Weighted_alpha_complex_3d::Bare_point_3;
+ std::vector<Bare_point_3> w_points;
+ w_points.push_back(Bare_point_3(0.0, 0.0, 0.0));
+ w_points.push_back(Bare_point_3(0.0, 0.0, 0.2));
+ w_points.push_back(Bare_point_3(0.2, 0.0, 0.2));
+ w_points.push_back(Bare_point_3(0.6, 0.6, 0.0));
+ w_points.push_back(Bare_point_3(0.8, 0.8, 0.2));
+ w_points.push_back(Bare_point_3(0.2, 0.8, 0.6));
// weights size is different from w_points size to make weighted Alpha_complex_3d throw in debug mode
std::vector<double> weights = {0.01, 0.005, 0.006, 0.01, 0.009, 0.001};
@@ -91,6 +91,24 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted, Weighted_alpha_complex_3d,
Gudhi::Simplex_tree<> stree_bis;
alpha_complex_w_p.create_complex(stree_bis);
+ for (std::size_t index = 0; index < weighted_points.size(); index++) {
+ bool found = false;
+ Weighted_point_3 awp = alpha_complex_w_p.get_point(index);
+ for (auto weighted_point : weighted_points) {
+ if ((weighted_point.weight() == awp.weight()) &&
+ (weighted_point.x() == awp.x()) &&
+ (weighted_point.y() == awp.y()) &&
+ (weighted_point.z() == awp.z())) {
+ found = true;
+ break;
+ }
+ }
+ // Check all points from alpha complex are found in the input point cloud
+ BOOST_CHECK(found);
+ }
+ // Exception if we go out of range
+ BOOST_CHECK_THROW(alpha_complex_w_p.get_point(weighted_points.size()), std::out_of_range);
+
// ---------------------
// Compare both versions
// ---------------------
diff --git a/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp
index 670c7799..b09e92d5 100644
--- a/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp
+++ b/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp
@@ -47,13 +47,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_periodic_throw, Weighted_pe
wp_variants_type_list) {
std::cout << "Weighted periodic alpha complex 3d exception throw" << std::endl;
- using Creator = CGAL::Creator_uniform_3<double, typename Weighted_periodic_alpha_complex_3d::Point_3>;
+ using Creator = CGAL::Creator_uniform_3<double, typename Weighted_periodic_alpha_complex_3d::Bare_point_3>;
CGAL::Random random(7);
- CGAL::Random_points_in_cube_3<typename Weighted_periodic_alpha_complex_3d::Point_3, Creator> in_cube(1, random);
- std::vector<typename Weighted_periodic_alpha_complex_3d::Point_3> wp_points;
+ CGAL::Random_points_in_cube_3<typename Weighted_periodic_alpha_complex_3d::Bare_point_3, Creator> in_cube(1, random);
+ std::vector<typename Weighted_periodic_alpha_complex_3d::Bare_point_3> wp_points;
for (int i = 0; i < 50; i++) {
- typename Weighted_periodic_alpha_complex_3d::Point_3 p = *in_cube++;
+ typename Weighted_periodic_alpha_complex_3d::Bare_point_3 p = *in_cube++;
wp_points.push_back(p);
}
std::vector<double> p_weights;
@@ -117,13 +117,13 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_weighted_periodic) {
// ---------------------
std::cout << "Fast weighted periodic alpha complex 3d" << std::endl;
- using Creator = CGAL::Creator_uniform_3<double, Fast_weighted_periodic_alpha_complex_3d::Point_3>;
+ using Creator = CGAL::Creator_uniform_3<double, Fast_weighted_periodic_alpha_complex_3d::Bare_point_3>;
CGAL::Random random(7);
- CGAL::Random_points_in_cube_3<Fast_weighted_periodic_alpha_complex_3d::Point_3, Creator> in_cube(1, random);
- std::vector<Fast_weighted_periodic_alpha_complex_3d::Point_3> p_points;
+ CGAL::Random_points_in_cube_3<Fast_weighted_periodic_alpha_complex_3d::Bare_point_3, Creator> in_cube(1, random);
+ std::vector<Fast_weighted_periodic_alpha_complex_3d::Bare_point_3> p_points;
for (int i = 0; i < 50; i++) {
- Fast_weighted_periodic_alpha_complex_3d::Point_3 p = *in_cube++;
+ Fast_weighted_periodic_alpha_complex_3d::Bare_point_3 p = *in_cube++;
p_points.push_back(p);
}
std::vector<double> p_weights;
@@ -142,10 +142,10 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_weighted_periodic) {
// ----------------------
std::cout << "Exact weighted periodic alpha complex 3d" << std::endl;
- std::vector<Exact_weighted_periodic_alpha_complex_3d::Point_3> e_p_points;
+ std::vector<Exact_weighted_periodic_alpha_complex_3d::Bare_point_3> e_p_points;
for (auto p : p_points) {
- e_p_points.push_back(Exact_weighted_periodic_alpha_complex_3d::Point_3(p[0], p[1], p[2]));
+ e_p_points.push_back(Exact_weighted_periodic_alpha_complex_3d::Bare_point_3(p[0], p[1], p[2]));
}
Exact_weighted_periodic_alpha_complex_3d exact_alpha_complex(e_p_points, p_weights, -1., -1., -1., 1., 1., 1.);
@@ -191,10 +191,10 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_weighted_periodic) {
// ----------------------
std::cout << "Safe weighted periodic alpha complex 3d" << std::endl;
- std::vector<Safe_weighted_periodic_alpha_complex_3d::Point_3> s_p_points;
+ std::vector<Safe_weighted_periodic_alpha_complex_3d::Bare_point_3> s_p_points;
for (auto p : p_points) {
- s_p_points.push_back(Safe_weighted_periodic_alpha_complex_3d::Point_3(p[0], p[1], p[2]));
+ s_p_points.push_back(Safe_weighted_periodic_alpha_complex_3d::Bare_point_3(p[0], p[1], p[2]));
}
Safe_weighted_periodic_alpha_complex_3d safe_alpha_complex(s_p_points, p_weights, -1., -1., -1., 1., 1., 1.);
diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt
index 5295f3cd..57b92942 100644
--- a/src/Alpha_complex/utilities/CMakeLists.txt
+++ b/src/Alpha_complex/utilities/CMakeLists.txt
@@ -7,9 +7,19 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
if (TBB_FOUND)
target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES})
endif(TBB_FOUND)
- add_test(NAME Alpha_complex_utilities_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence>
- "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45")
-
+ add_test(NAME Alpha_complex_utilities_safe_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence>
+ "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "safe.pers")
+ add_test(NAME Alpha_complex_utilities_fast_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence>
+ "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "fast.pers" "-f")
+ add_test(NAME Alpha_complex_utilities_exact_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence>
+ "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "exact.pers" "-e")
+ if (DIFF_PATH)
+ add_test(Alpha_complex_utilities_diff_exact_alpha_complex ${DIFF_PATH}
+ "exact.pers" "safe.pers")
+ add_test(Alpha_complex_utilities_diff_fast_alpha_complex ${DIFF_PATH}
+ "fast.pers" "safe.pers")
+ endif()
+
install(TARGETS alpha_complex_persistence DESTINATION bin)
add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp)
@@ -20,21 +30,21 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
add_test(NAME Alpha_complex_utilities_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
"${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off"
- "-p" "2" "-m" "0.45" "-o" "safe.pers")
+ "-p" "2" "-m" "0.45" "-o" "safe_3d.pers")
add_test(NAME Alpha_complex_utilities_exact_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
"${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off"
- "-p" "2" "-m" "0.45" "-o" "exact.pers" "-e")
+ "-p" "2" "-m" "0.45" "-o" "exact_3d.pers" "-e")
add_test(NAME Alpha_complex_utilities_safe_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
"${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off"
- "-p" "2" "-m" "0.45" "-o" "fast.pers" "-f")
+ "-p" "2" "-m" "0.45" "-o" "fast_3d.pers" "-f")
if (DIFF_PATH)
- add_test(Alpha_complex_utilities_diff_alpha_complex_3d ${DIFF_PATH}
- "exact.pers" "safe.pers")
- add_test(Alpha_complex_utilities_diff_alpha_complex_3d ${DIFF_PATH}
- "fast.pers" "safe.pers")
+ add_test(Alpha_complex_utilities_diff_exact_alpha_complex_3d ${DIFF_PATH}
+ "exact_3d.pers" "safe_3d.pers")
+ add_test(Alpha_complex_utilities_diff_fast_alpha_complex_3d ${DIFF_PATH}
+ "fast_3d.pers" "safe_3d.pers")
endif()
add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
index 2272576e..929fc2e8 100644
--- a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
+++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
@@ -62,9 +62,9 @@ bool read_cuboid_file(const std::string &cuboid_file, double &x_min, double &y_m
}
template <typename AlphaComplex3d>
-std::vector<typename AlphaComplex3d::Point_3> read_off(const std::string &off_file_points) {
+std::vector<typename AlphaComplex3d::Bare_point_3> read_off(const std::string &off_file_points) {
// Read the OFF file (input file name given as parameter) and triangulate points
- Gudhi::Points_3D_off_reader<typename AlphaComplex3d::Point_3> off_reader(off_file_points);
+ Gudhi::Points_3D_off_reader<typename AlphaComplex3d::Bare_point_3> off_reader(off_file_points);
// Check the read operation was correct
if (!off_reader.is_valid()) {
std::cerr << "Unable to read OFF file " << off_file_points << std::endl;
diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
index fab7bd30..486347cc 100644
--- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
+++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
@@ -24,63 +24,84 @@
using Simplex_tree = Gudhi::Simplex_tree<>;
using Filtration_value = Simplex_tree::Filtration_value;
-void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag,
- Filtration_value &alpha_square_max_value, int &coeff_field_characteristic,
- Filtration_value &min_persistence);
+void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast,
+ std::string &output_file_diag, Filtration_value &alpha_square_max_value,
+ int &coeff_field_characteristic, Filtration_value &min_persistence);
int main(int argc, char **argv) {
std::string off_file_points;
std::string output_file_diag;
+ bool exact_version = false;
+ bool fast_version = false;
Filtration_value alpha_square_max_value;
int coeff_field_characteristic;
Filtration_value min_persistence;
- program_options(argc, argv, off_file_points, output_file_diag, alpha_square_max_value, coeff_field_characteristic,
- min_persistence);
+ program_options(argc, argv, off_file_points, exact_version, fast_version, output_file_diag, alpha_square_max_value,
+ coeff_field_characteristic, min_persistence);
- // ----------------------------------------------------------------------------
- // Init of an alpha complex from an OFF file
- // ----------------------------------------------------------------------------
- using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
- Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_points);
+ if ((exact_version) && (fast_version)) {
+ std::cerr << "You cannot set the exact and the fast version." << std::endl;
+ exit(-1);
+ }
Simplex_tree simplex;
- if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) {
- // ----------------------------------------------------------------------------
- // Display information about the alpha complex
- // ----------------------------------------------------------------------------
- std::cout << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
- << " simplices - " << simplex.num_vertices() << " vertices." << std::endl;
-
- // Sort the simplices in the order of the filtration
- simplex.initialize_filtration();
-
- std::cout << "Simplex_tree dim: " << simplex.dimension() << std::endl;
- // Compute the persistence diagram of the complex
- Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(
- simplex);
- // initializes the coefficient field for homology
- pcoh.init_coefficients(coeff_field_characteristic);
-
- pcoh.compute_persistent_cohomology(min_persistence);
-
- // Output the diagram in filediag
- if (output_file_diag.empty()) {
- pcoh.output_diagram();
- } else {
- std::cout << "Result in file: " << output_file_diag << std::endl;
- std::ofstream out(output_file_diag);
- pcoh.output_diagram(out);
- out.close();
+ if (fast_version) {
+ // WARNING : CGAL::Epick_d is fast but not safe (unlike CGAL::Epeck_d)
+ // (i.e. when the points are on a grid)
+ using Fast_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
+
+ // Init of an alpha complex from an OFF file
+ Gudhi::alpha_complex::Alpha_complex<Fast_kernel> alpha_complex_from_file(off_file_points);
+
+ if (!alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) {
+ std::cerr << "Fast Alpha complex simplicial complex creation failed." << std::endl;
+ exit(-1);
}
- }
+ } else {
+ using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+ // Init of an alpha complex from an OFF file
+ Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_points);
+
+ if (!alpha_complex_from_file.create_complex(simplex, alpha_square_max_value, exact_version)) {
+ std::cerr << "Alpha complex simplicial complex creation failed." << std::endl;
+ exit(-1);
+ }
+ }
+ // ----------------------------------------------------------------------------
+ // Display information about the alpha complex
+ // ----------------------------------------------------------------------------
+ std::cout << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
+ << " simplices - " << simplex.num_vertices() << " vertices." << std::endl;
+
+ // Sort the simplices in the order of the filtration
+ simplex.initialize_filtration();
+
+ std::cout << "Simplex_tree dim: " << simplex.dimension() << std::endl;
+ // Compute the persistence diagram of the complex
+ Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(
+ simplex);
+ // initializes the coefficient field for homology
+ pcoh.init_coefficients(coeff_field_characteristic);
+
+ pcoh.compute_persistent_cohomology(min_persistence);
+
+ // Output the diagram in filediag
+ if (output_file_diag.empty()) {
+ pcoh.output_diagram();
+ } else {
+ std::cout << "Result in file: " << output_file_diag << std::endl;
+ std::ofstream out(output_file_diag);
+ pcoh.output_diagram(out);
+ out.close();
+ }
return 0;
}
-void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag,
- Filtration_value &alpha_square_max_value, int &coeff_field_characteristic,
- Filtration_value &min_persistence) {
+void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast,
+ std::string &output_file_diag, Filtration_value &alpha_square_max_value,
+ int &coeff_field_characteristic, Filtration_value &min_persistence) {
namespace po = boost::program_options;
po::options_description hidden("Hidden options");
hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
@@ -88,6 +109,10 @@ void program_options(int argc, char *argv[], std::string &off_file_points, std::
po::options_description visible("Allowed options", 100);
visible.add_options()("help,h", "produce help message")(
+ "exact,e", po::bool_switch(&exact),
+ "To activate exact version of Alpha complex (default is false, not available if fast is set)")(
+ "fast,f", po::bool_switch(&fast),
+ "To activate fast version of Alpha complex (default is false, not available if exact is set)")(
"output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
"Name of file in which the persistence diagram is written. Default print in std::cout")(
"max-alpha-square-value,r", po::value<Filtration_value>(&alpha_square_max_value)
diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md
index fcd16a3b..527598a9 100644
--- a/src/Alpha_complex/utilities/alphacomplex.md
+++ b/src/Alpha_complex/utilities/alphacomplex.md
@@ -46,6 +46,8 @@ for the Alpha complex construction.
coefficient field Z/pZ for computing homology.
* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature
to be recorded. Enter a negative value to see zero length intervals.
+* `-e [ --exact ]` for the exact computation version.
+* `-f [ --fast ]` for the fast computation version.
**Example**
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-26 06:22:04 +0000