From 3477d83ef23505e920c459157e9110acbb276b63 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Tue, 26 Jun 2018 13:22:45 +0000 Subject: weighted_alpha_complex_3d_persistence.cpp using Alpha_complex_3d new module git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/alpha_complex_3d_module_vincent@3639 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: b7e5f4674a56bef537ba5f4dda62e6605bf8c249 --- .../weighted_alpha_complex_3d_persistence.cpp | 181 ++------------------- 1 file changed, 13 insertions(+), 168 deletions(-) diff --git a/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp index cdeeabfc..56fd6a7c 100644 --- a/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp @@ -23,23 +23,12 @@ #include #include +#include +#include #include #include #include -#include -#include -#include -#include -#include -#include -#include - -// For CGAL < 4.11 -#if CGAL_VERSION_NR < 1041100000 -#include -#endif // CGAL_VERSION_NR < 1041100000 - #include #include #include @@ -51,58 +40,12 @@ #include "alpha_complex_3d_helper.h" -// Alpha_shape_3 templates type definitions -using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; - -// For CGAL < 4.11 -#if CGAL_VERSION_NR < 1041100000 -using Gt = CGAL::Regular_triangulation_euclidean_traits_3; -using Vb = CGAL::Alpha_shape_vertex_base_3; -using Fb = CGAL::Alpha_shape_cell_base_3; -using Tds = CGAL::Triangulation_data_structure_3; -using Triangulation_3 = CGAL::Regular_triangulation_3; - -// From file type definition -using Point_3 = Gt::Bare_point; -using Weighted_point_3 = Gt::Weighted_point; - -// For CGAL >= 4.11 -#else // CGAL_VERSION_NR < 1041100000 -using Rvb = CGAL::Regular_triangulation_vertex_base_3; -using Vb = CGAL::Alpha_shape_vertex_base_3; -using Rcb = CGAL::Regular_triangulation_cell_base_3; -using Cb = CGAL::Alpha_shape_cell_base_3; -using Tds = CGAL::Triangulation_data_structure_3; -using Triangulation_3 = CGAL::Regular_triangulation_3; - -// From file type definition -using Point_3 = Triangulation_3::Bare_point; -using Weighted_point_3 = Triangulation_3::Weighted_point; -#endif // CGAL_VERSION_NR < 1041100000 - -using Alpha_shape_3 = CGAL::Alpha_shape_3; - -// filtration with alpha values needed type definition -using Alpha_value_type = Alpha_shape_3::FT; -using Object = CGAL::Object; -using Dispatch = - CGAL::Dispatch_output_iterator, - CGAL::cpp11::tuple >, - std::back_insert_iterator > > >; -using Cell_handle = Alpha_shape_3::Cell_handle; -using Facet = Alpha_shape_3::Facet; -using Edge_3 = Alpha_shape_3::Edge; -using Vertex_handle = Alpha_shape_3::Vertex_handle; -using Vertex_list = std::vector; - // gudhi type definition -using ST = Gudhi::Simplex_tree; -using Filtration_value = ST::Filtration_value; -using Simplex_tree_vertex = ST::Vertex_handle; -using Alpha_shape_simplex_tree_map = std::map; -using Simplex_tree_vector_vertex = std::vector; +using Alpha_complex_3d = Gudhi::alpha_complex::Alpha_complex_3d; +using Simplex_tree = Gudhi::Simplex_tree; +using Filtration_value = Simplex_tree::Filtration_value; using Persistent_cohomology = - Gudhi::persistent_cohomology::Persistent_cohomology; + Gudhi::persistent_cohomology::Persistent_cohomology; void program_options(int argc, char *argv[], std::string &off_file_points, std::string &weight_file, std::string &output_file_diag, int &coeff_field_characteristic, Filtration_value &min_persistence); @@ -118,131 +61,33 @@ int main(int argc, char **argv) { min_persistence); // Read the OFF file (input file name given as parameter) and triangulate points - Gudhi::Points_3D_off_reader off_reader(off_file_points); + Gudhi::Points_3D_off_reader 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; exit(-1); } - // Retrieve the points - std::vector lp = off_reader.get_point_cloud(); - // Read weights information from file std::ifstream weights_ifstr(weight_file); - std::vector wp; + std::vector weights; if (weights_ifstr.good()) { double weight = 0.0; - std::size_t index = 0; - wp.reserve(lp.size()); // Attempt read the weight in a double format, return false if it fails - while ((weights_ifstr >> weight) && (index < lp.size())) { - wp.push_back(Weighted_point_3(lp[index], weight)); - index++; - } - if (index != lp.size()) { - std::cerr << "Bad number of weights in file " << weight_file << std::endl; - exit(-1); + while (weights_ifstr >> weight) { + weights.push_back(weight); } } else { std::cerr << "Unable to read weights file " << weight_file << std::endl; exit(-1); } - // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. - Alpha_shape_3 as(wp.begin(), wp.end(), 0, Alpha_shape_3::GENERAL); -#ifdef DEBUG_TRACES - std::cout << "Alpha shape computed in GENERAL mode" << std::endl; -#endif // DEBUG_TRACES - - // 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 - - 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; + Alpha_complex_3d alpha_complex(off_reader.get_point_cloud(), weights); - // Loop on objects vector - Vertex_list vertex_list; - ST simplex_tree; - Alpha_shape_simplex_tree_map map_cgal_simplex_tree; - 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); - count_cells++; - } else if (const Facet *facet = CGAL::object_cast(&object_iterator)) { - vertex_list = from_facet(*facet); - count_facets++; - } else if (const Edge_3 *edge = CGAL::object_cast(&object_iterator)) { - vertex_list = from_edge(*edge); - count_edges++; - } else if (const Vertex_handle *vertex = CGAL::object_cast(&object_iterator)) { - count_vertices++; - vertex_list = from_vertex(*vertex); - } - // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex - Simplex_tree_vector_vertex the_simplex; - for (auto the_alpha_shape_vertex : vertex_list) { - Alpha_shape_simplex_tree_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 - Simplex_tree_vertex 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 - Simplex_tree_vertex 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); - if (the_alpha_value_iterator != the_alpha_values.end()) - ++the_alpha_value_iterator; - else - std::cout << "This shall not happen" << std::endl; - } + Simplex_tree simplex_tree; -#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; + alpha_complex.create_complex(simplex_tree); - 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 // Sort the simplices in the order of the filtration simplex_tree.initialize_filtration(); -- cgit v1.2.3