From 82f15e71f43a5547627699b625830126dad1c4ab Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Tue, 19 Jun 2018 21:26:17 +0000 Subject: First create_complex version weighted periodic option git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/alpha_complex_3d_module_vincent@3624 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 7d22e307b12aac1390593541f77cabcccffbc6ad --- src/Alpha_complex/include/gudhi/Alpha_complex_3d.h | 243 ++++++++++++++++++++- .../include/gudhi/Alpha_complex_3d_options.h | 26 +++ 2 files changed, 268 insertions(+), 1 deletion(-) (limited to 'src/Alpha_complex/include') diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h index 21707e5e..79f8593d 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h @@ -23,6 +23,12 @@ #ifndef ALPHA_COMPLEX_3D_H_ #define ALPHA_COMPLEX_3D_H_ +#include +#include + +#if BOOST_VERSION >= 105400 +#include +#endif #include #include @@ -49,7 +55,13 @@ private: using Cell_handle = typename Alpha_shape_3::Cell_handle; using Facet = typename Alpha_shape_3::Facet; using Edge_3 = typename Alpha_shape_3::Edge; - using Vertex_handle = typename Alpha_shape_3::Vertex_handle; + using Alpha_vertex_handle = typename Alpha_shape_3::Vertex_handle; + +#if BOOST_VERSION >= 105400 + using Vertex_list = boost::container::static_vector; +#else + using Vertex_list = std::vector; +#endif public: /** \brief Alpha_complex constructor from a list of points. @@ -166,6 +178,235 @@ public: } + /** \brief Alpha_complex constructor from a list of points. +* +* Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous. +* +* @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d +* +* The type InputPointRange must be a range for which std::begin and +* std::end return input iterators on a Kernel::Point_d. +*/ + template + Alpha_complex_3d(const InputPointRange& points, WeightRange weights, + Alpha_value_type x_min, Alpha_value_type y_min, Alpha_value_type z_min, + Alpha_value_type x_max, Alpha_value_type y_max, Alpha_value_type z_max) { + static_assert(AlphaComplex3dOptions::weighted, + "This constructor is not available for non-weighted versions of Alpha_complex_3d"); + static_assert(AlphaComplex3dOptions::periodic, + "This constructor is not available for non-periodic versions of Alpha_complex_3d"); + GUDHI_CHECK((weights.size() == points.size()), + std::invalid_argument("Alpha_complex_3d constructor with weights requires points number to be equal with points number")); + // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it. + GUDHI_CHECK((x_max - x_min != y_max - y_min) || (x_max - x_min != z_max - z_min) || (z_max - z_min != y_max - y_min), + std::invalid_argument("The size of the cuboid in every directions is not the same.")); + + double maximal_possible_weight = 0.015625 * (x_max - x_min) * (x_max - x_min); + + using Weighted_point_3 = typename AlphaComplex3dOptions::Weighted_point_3; + std::vector weighted_points_3; + + std::size_t index = 0; + weighted_points_3.reserve(points.size()); + while ((index < weights.size()) && (index < points.size())) { + GUDHI_CHECK((weights[index] >= maximal_possible_weight) || (weights[index] < 0), + std::invalid_argument("Invalid weight at line" + std::to_string(index + 1) + + ". Must be positive and less than maximal possible weight = 1/64*cuboid length " + "squared, which is not an acceptable input.")); + weighted_points_3.push_back(Weighted_point_3(points[index], weights[index])); + index++; + } + + using Periodic_delaunay_triangulation_3 = typename AlphaComplex3dOptions::Periodic_delaunay_triangulation_3; + using Iso_cuboid_3 = typename AlphaComplex3dOptions::Iso_cuboid_3; + // Define the periodic cube + Periodic_delaunay_triangulation_3 pdt(Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max)); + // Heuristic for inserting large point sets (if pts is reasonably large) + pdt.insert(std::begin(weighted_points_3), std::end(weighted_points_3), true); + // As pdt won't be modified anymore switch to 1-sheeted cover if possible + GUDHI_CHECK(pdt.is_triangulation_in_1_sheet(), + std::invalid_argument("Uable to construct a triangulation within a single periodic domain.")); + + // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode + // Maybe need to set it to GENERAL mode + Alpha_shape_3 as(pdt, 0, Alpha_shape_3::GENERAL); + + // filtration with alpha values from alpha shape + std::vector the_objects; + std::vector the_alpha_values; + + Dispatch disp = CGAL::dispatch_output(std::back_inserter(the_objects), + std::back_inserter(the_alpha_values)); + + as.filtration_with_alpha_values(disp); +#ifdef DEBUG_TRACES + std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl; +#endif // DEBUG_TRACES + } + + template + void create_complex(SimplicialComplexForAlpha3d& complex) { + using Filtration_value = typename SimplicialComplexForAlpha3d::Filtration_value; + create_complex(complex, std::numeric_limits::infinity()); + } + + /** \brief Inserts all Delaunay triangulation into the simplicial complex. + * It also computes the filtration values accordingly to the \ref createcomplexalgorithm + * + * \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$. + * + * @return true if creation succeeds, false otherwise. + * + * @pre Delaunay triangulation must be already constructed with dimension strictly greater than 0. + * @pre The simplicial complex must be empty (no vertices) + * + * Initialization can be launched once. + */ + template + void create_complex(SimplicialComplexForAlpha3d& complex, + typename SimplicialComplexForAlpha3d::Filtration_value max_alpha_square) { + using Filtration_value = typename SimplicialComplexForAlpha3d::Filtration_value; + using Vertex_handle = typename SimplicialComplexForAlpha3d::Vertex_handle; +#ifdef DEBUG_TRACES + Alpha_shape_3::size_type count_vertices = 0; + Alpha_shape_3::size_type count_edges = 0; + Alpha_shape_3::size_type count_facets = 0; + Alpha_shape_3::size_type count_cells = 0; +#endif // DEBUG_TRACES + + // Loop on objects vector + Vertex_list vertex_list; + SimplicialComplexForAlpha3d simplex_tree; + std::map map_cgal_simplex_tree; + typename std::vector::iterator the_alpha_value_iterator = the_alpha_values.begin(); + for (auto object_iterator : the_objects) { + // Retrieve Alpha shape vertex list from object + if (const Cell_handle *cell = CGAL::object_cast(&object_iterator)) { + vertex_list = from_cell(*cell); +#ifdef DEBUG_TRACES + count_cells++; +#endif // DEBUG_TRACES + } else if (const Facet *facet = CGAL::object_cast(&object_iterator)) { + vertex_list = from_facet(*facet); +#ifdef DEBUG_TRACES + count_facets++; +#endif // DEBUG_TRACES + } else if (const Edge_3 *edge = CGAL::object_cast(&object_iterator)) { + vertex_list = from_edge(*edge); +#ifdef DEBUG_TRACES + count_edges++; +#endif // DEBUG_TRACES + } else if (const Vertex_handle *vertex = CGAL::object_cast(&object_iterator)) { + vertex_list = from_vertex(*vertex); +#ifdef DEBUG_TRACES + count_vertices++; +#endif // DEBUG_TRACES + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + std::vector the_simplex; + for (auto the_alpha_shape_vertex : vertex_list) { + typename std::map::iterator the_map_iterator = + map_cgal_simplex_tree.find(the_alpha_shape_vertex); + if (the_map_iterator == map_cgal_simplex_tree.end()) { + // alpha shape not found + 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; +#endif // DEBUG_TRACES + the_simplex.push_back(vertex); + map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex); + } else { + // alpha shape found + Vertex_handle vertex = the_map_iterator->second; +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex.push_back(vertex); + } + } + // Construction of the simplex_tree + Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator); +#ifdef DEBUG_TRACES + std::cout << "filtration = " << filtr << std::endl; +#endif // DEBUG_TRACES + simplex_tree.insert_simplex(the_simplex, filtr); + GUDHI_CHECK(the_alpha_value_iterator != the_alpha_values.end(), "CGAL provided more simplices than values"); + ++the_alpha_value_iterator; + } + +#ifdef DEBUG_TRACES + std::cout << "vertices \t\t" << count_vertices << std::endl; + std::cout << "edges \t\t" << count_edges << std::endl; + std::cout << "facets \t\t" << count_facets << std::endl; + std::cout << "cells \t\t" << count_cells << std::endl; + + std::cout << "Information of the Simplex Tree: " << std::endl; + std::cout << " Number of vertices = " << simplex_tree.num_vertices() << " "; + std::cout << " Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl; + std::cout << " Dimension = " << simplex_tree.dimension() << " "; +#endif // DEBUG_TRACES + +#ifdef DEBUG_TRACES + std::cout << "Iterator on vertices: " << std::endl; + for (auto vertex : simplex_tree.complex_vertex_range()) { + std::cout << vertex << " "; + } +#endif // DEBUG_TRACES + + } + +private: + template + Vertex_list from_cell(const Cell_handle& ch) { + Vertex_list the_list; + for (auto i = 0; i < 4; i++) { +#ifdef DEBUG_TRACES + std::cout << "from cell[" << i << "]=" << ch->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + the_list.push_back(ch->vertex(i)); + } + return the_list; + } + + template + Vertex_list from_facet(const Facet& fct) { + Vertex_list the_list; + for (auto i = 0; i < 4; i++) { + if (fct.second != i) { +#ifdef DEBUG_TRACES + std::cout << "from facet=[" << i << "]" << fct.first->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + the_list.push_back(fct.first->vertex(i)); + } + } + return the_list; + } + + template + Vertex_list from_edge(const Edge_3& edg) { + Vertex_list the_list; + for (auto i : {edg.second, edg.third}) { +#ifdef DEBUG_TRACES + std::cout << "from edge[" << i << "]=" << edg.first->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + the_list.push_back(edg.first->vertex(i)); + } + return the_list; + } + + template + Vertex_list from_vertex(const Vertex_handle& vh) { + Vertex_list the_list; +#ifdef DEBUG_TRACES + std::cout << "from vertex=" << vh->point() << std::endl; +#endif // DEBUG_TRACES + the_list.push_back(vh); + return the_list; + } + private: std::vector the_objects; std::vector the_alpha_values; diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d_options.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d_options.h index e05e5a7f..9433afdc 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d_options.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d_options.h @@ -28,6 +28,8 @@ #include #include #include +#include +#include #include #include #include @@ -119,6 +121,30 @@ public: static const bool periodic = true; }; +class Weighted_periodic_alpha_shapes_3d { +private: + using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; + using Periodic_kernel = CGAL::Periodic_3_regular_triangulation_traits_3; + using DsVb = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>; + using Vb = CGAL::Regular_triangulation_vertex_base_3; + using AsVb = CGAL::Alpha_shape_vertex_base_3; + using DsCb = CGAL::Periodic_3_triangulation_ds_cell_base_3<>; + using Cb = CGAL::Regular_triangulation_cell_base_3; + using AsCb = CGAL::Alpha_shape_cell_base_3; + using Tds = CGAL::Triangulation_data_structure_3; + +public: + using Periodic_delaunay_triangulation_3 = CGAL::Periodic_3_regular_triangulation_3; + using Alpha_shape_3 = CGAL::Alpha_shape_3; + using Point_3 = Periodic_delaunay_triangulation_3::Bare_point; + using Weighted_point_3 = Periodic_delaunay_triangulation_3::Weighted_point; + using Alpha_value_type = Alpha_shape_3::FT; + using Iso_cuboid_3 = Periodic_kernel::Iso_cuboid_3; + + static const bool weighted = true; + static const bool periodic = true; +}; + } // namespace alpha_complex } // namespace Gudhi -- cgit v1.2.3