From 7f8b2a07a463fdcec12430abe2119d4f86a72517 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Fri, 29 May 2015 10:34:14 +0000 Subject: Alpha_shapes filtration value git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/alphashapes@599 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 303ed34da0718312e4bab3f50b7196acded94cec --- src/Alpha_shapes/example/CMakeLists.txt | 9 +- .../example/Delaunay_triangulation_off_rw.cpp | 2 +- .../Simplex_tree_from_delaunay_triangulation.cpp | 32 +----- src/Alpha_shapes/include/gudhi/Alpha_shapes.h | 89 +++++++++++----- .../Alpha_shapes/Delaunay_triangulation_off_io.h | 113 +++++++-------------- src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp | 2 +- src/Alpha_shapes/test/CMakeLists.txt | 6 +- 7 files changed, 112 insertions(+), 141 deletions(-) diff --git a/src/Alpha_shapes/example/CMakeLists.txt b/src/Alpha_shapes/example/CMakeLists.txt index dab3861c..582a9322 100644 --- a/src/Alpha_shapes/example/CMakeLists.txt +++ b/src/Alpha_shapes/example/CMakeLists.txt @@ -1,10 +1,10 @@ cmake_minimum_required(VERSION 2.6) project(GUDHIAlphaShapesExample) -# need CGAL 4.6 -# cmake -DCGAL_DIR=~/workspace/CGAL-4.6-beta1 ../../.. +# need CGAL 4.7 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.7-Ic-41 ../../.. if(CGAL_FOUND) - if (NOT CGAL_VERSION VERSION_LESS 4.6.0) + if (NOT CGAL_VERSION VERSION_LESS 4.7.0) message(STATUS "CGAL version: ${CGAL_VERSION}.") include( ${CGAL_USE_FILE} ) @@ -13,12 +13,13 @@ if(CGAL_FOUND) if (EIGEN3_FOUND) message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.") include( ${EIGEN3_USE_FILE} ) + + add_definitions(-DDEBUG_TRACES) add_executable ( dtoffrw Delaunay_triangulation_off_rw.cpp ) target_link_libraries(dtoffrw ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) add_test(dtoffrw_tore3D ${CMAKE_CURRENT_BINARY_DIR}/dtoffrw ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off 3) - #add_definitions(-DDEBUG_TRACES) add_executable ( stfromdt Simplex_tree_from_delaunay_triangulation.cpp ) target_link_libraries(stfromdt ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) else() diff --git a/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp b/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp index 03f6440d..a31c44ab 100644 --- a/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp +++ b/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp @@ -63,7 +63,7 @@ int main(int argc, char **argv) { T dt(dimension); std::string offFileName(argv[1]); - Gudhi::alphashapes::Delaunay_triangulation_off_reader off_reader(offFileName, dt, true, true); + Gudhi::alphashapes::Delaunay_triangulation_off_reader off_reader(offFileName, dt); if (!off_reader.is_valid()) { std::cerr << "Unable to read file " << offFileName << std::endl; exit(-1); // ----- >> diff --git a/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp b/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp index 3a17b519..11079a03 100644 --- a/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp +++ b/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp @@ -48,50 +48,24 @@ typedef CGAL::Delaunay_triangulation T; // TriangulationDataStructure template parameter void usage(char * const progName) { - std::cerr << "Usage: " << progName << " filename.off dimension" << std::endl; + std::cerr << "Usage: " << progName << " filename.off" << std::endl; exit(-1); // ----- >> } int main(int argc, char **argv) { - if (argc != 3) { + if (argc != 2) { std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl; usage(argv[0]); } - int dimension = 0; - int returnedScanValue = sscanf(argv[2], "%d", &dimension); - if ((returnedScanValue == EOF) || (dimension <= 0)) { - std::cerr << "Error: " << argv[2] << " is not correct" << std::endl; - usage(argv[0]); - } - - // ---------------------------------------------------------------------------- - // - // Init of an alpha-shape from a Delaunay triangulation - // - // ---------------------------------------------------------------------------- - T dt(dimension); std::string off_file_name(argv[1]); - Gudhi::alphashapes::Delaunay_triangulation_off_reader off_reader(off_file_name, dt, false, false); - if (!off_reader.is_valid()) { - std::cerr << "Unable to read file " << off_file_name << std::endl; - exit(-1); // ----- >> - } - - std::cout << "number of vertices=" << dt.number_of_vertices() << std::endl; - std::cout << "number of full cells=" << dt.number_of_full_cells() << std::endl; - std::cout << "number of finite full cells=" << dt.number_of_finite_full_cells() << std::endl; - - Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_dt(dt); - //std::cout << alpha_shapes_from_dt << std::endl; - // ---------------------------------------------------------------------------- // // Init of an alpha-shape from a OFF file // // ---------------------------------------------------------------------------- - Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_file(off_file_name, dimension); + Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_file(off_file_name); //std::cout << alpha_shapes_from_file << std::endl; std::cout << "alpha_shapes_from_file.dimension()=" << alpha_shapes_from_file.dimension() << std::endl; diff --git a/src/Alpha_shapes/include/gudhi/Alpha_shapes.h b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h index b8efdb4d..6f290938 100644 --- a/src/Alpha_shapes/include/gudhi/Alpha_shapes.h +++ b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h @@ -47,6 +47,8 @@ namespace Gudhi { namespace alphashapes { +#define Kinit(f) =k.f() + /** \defgroup alpha_shapes Alpha shapes in dimension N *
Implementations:
@@ -86,16 +88,22 @@ class Alpha_shapes { */ typedef CGAL::Delaunay_triangulation Delaunay_triangulation; + typedef typename Kernel::Compute_squared_radius_d Squared_Radius; + typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel; + + /** \brief Type required to insert into a simplex_tree (with or without subfaces).*/ + typedef std::vector typeVectorPoint; + private: /** \brief Upper bound on the simplex tree of the simplicial complex.*/ Gudhi::Simplex_tree<> _st; public: - Alpha_shapes(std::string off_file_name, int dimension) { - Delaunay_triangulation dt(dimension); - Gudhi::alphashapes::Delaunay_triangulation_off_reader - off_reader(off_file_name, dt, true, true); + Alpha_shapes(std::string off_file_name) { + // Construct a default Delaunay_triangulation (dim=0) - dim will be set in visitor reader init function + Delaunay_triangulation dt(3); + Gudhi::alphashapes::Delaunay_triangulation_off_reader off_reader(off_file_name, dt); if (!off_reader.is_valid()) { std::cerr << "Unable to read file " << off_file_name << std::endl; exit(-1); // ----- >> @@ -120,45 +128,70 @@ class Alpha_shapes { template void init(T triangulation) { _st.set_dimension(triangulation.maximal_dimension()); - _st.set_filtration(0.0); - // triangulation points list - for (auto vit = triangulation.finite_vertices_begin(); - vit != triangulation.finite_vertices_end(); ++vit) { - typeVectorVertex vertexVector; - Vertex_handle vertexHdl = std::distance(triangulation.finite_vertices_begin(), vit); - vertexVector.push_back(vertexHdl); + Filtration_value filtration_max = 0.0; - // Insert each point in the simplex tree - _st.insert_simplex(vertexVector, 0.0); + Kernel k; + Squared_Radius squared_radius Kinit(compute_squared_radius_d_object); + Is_Gabriel is_gabriel Kinit(side_of_bounded_sphere_d_object); + // triangulation full cells list + for (auto cit = triangulation.full_cells_begin(); cit != triangulation.full_cells_end(); ++cit) { + typeVectorVertex vertexVector; + typeVectorPoint pointVector; + for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { + if (!triangulation.is_infinite(*vit)) { + // Vector of vertex construction for simplex_tree structure + // Vertex handle is distance - 1 + Vertex_handle vertexHdl = std::distance(triangulation.vertices_begin(), *vit) - 1; + // infinite cell is -1 for infinite + vertexVector.push_back(vertexHdl); + // Vector of points for alpha_shapes filtration value computation + pointVector.push_back((*vit)->point()); #ifdef DEBUG_TRACES - std::cout << "P" << vertexHdl << ":"; - for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) { - std::cout << *Coord << " "; + std::cout << "Point "; + for (auto Coord = (*vit)->point().cartesian_begin(); Coord != (*vit)->point().cartesian_end(); ++Coord) { + std::cout << *Coord << " | "; + } + std::cout << std::endl; +#endif // DEBUG_TRACES + } } - std::cout << std::endl; + Filtration_value alpha_shapes_filtration = 0.0; + + if (!triangulation.is_infinite(cit)) { + alpha_shapes_filtration = squared_radius(pointVector.begin(), pointVector.end()); +#ifdef DEBUG_TRACES + std::cout << "Alpha_shape filtration value = " << alpha_shapes_filtration << std::endl; #endif // DEBUG_TRACES - } - // triangulation finite full cells list - for (auto cit = triangulation.finite_full_cells_begin(); - cit != triangulation.finite_full_cells_end(); ++cit) { - typeVectorVertex vertexVector; - for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { - // Vertex handle is distance - 1 - Vertex_handle vertexHdl = std::distance(triangulation.vertices_begin(), *vit) - 1; - vertexVector.push_back(vertexHdl); + } else { + Filtration_value tmp_filtration = 0.0; + bool is_gab = true; + for (auto vit = triangulation.finite_vertices_begin(); vit != triangulation.finite_vertices_end(); ++vit) { + if (CGAL::ON_UNBOUNDED_SIDE != is_gabriel(pointVector.begin(), pointVector.end(), vit->point())) { + is_gab = false; + // TODO(VR) : Compute minimum + + } + } + if (true == is_gab) { + alpha_shapes_filtration = squared_radius(pointVector.begin(), pointVector.end()); +#ifdef DEBUG_TRACES + std::cout << "Alpha_shape filtration value = " << alpha_shapes_filtration << std::endl; +#endif // DEBUG_TRACES + } } // Insert each point in the simplex tree - _st.insert_simplex_and_subfaces(vertexVector, 0.0); + _st.insert_simplex_and_subfaces(vertexVector, alpha_shapes_filtration); #ifdef DEBUG_TRACES - std::cout << "C" << std::distance(triangulation.finite_full_cells_begin(), cit) << ":"; + std::cout << "C" << std::distance(triangulation.full_cells_begin(), cit) << ":"; for (auto value : vertexVector) { std::cout << value << ' '; } std::cout << std::endl; #endif // DEBUG_TRACES } + _st.set_filtration(filtration_max); } public: diff --git a/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h b/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h index 693b393e..3215c8f6 100644 --- a/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h +++ b/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h @@ -36,27 +36,35 @@ namespace alphashapes { *@brief Off reader visitor with flag that can be passed to Off_reader to read a Delaunay_triangulation_complex. */ template -class Delaunay_triangulation_off_flag_visitor_reader { +class Delaunay_triangulation_off_visitor_reader { Complex& complex_; typedef typename Complex::Point Point; - const bool load_only_points_; - public: - explicit Delaunay_triangulation_off_flag_visitor_reader(Complex& complex, bool load_only_points = false) : - complex_(complex), - load_only_points_(load_only_points) { } + + explicit Delaunay_triangulation_off_visitor_reader(Complex& complex) : + complex_(complex) { } void init(int dim, int num_vertices, int num_faces, int num_edges) { #ifdef DEBUG_TRACES - std::cout << "init" << std::endl; + std::cout << "Delaunay_triangulation_off_visitor_reader::init - dim=" << dim << " - num_vertices=" << + num_vertices << " - num_faces=" << num_faces << " - num_edges=" << num_edges << std::endl; #endif // DEBUG_TRACES + if (num_faces > 0) { + std::cerr << "Delaunay_triangulation_off_visitor_reader::init faces are not taken into account from OFF " << + "file for Delaunay triangulation - faces are computed." << std::endl; + } + if (num_edges > 0) { + std::cerr << "Delaunay_triangulation_off_visitor_reader::init edges are not taken into account from OFF " << + "file for Delaunay triangulation - edges are computed." << std::endl; + } + //complex_.set_current_dimension(dim); } void point(const std::vector& point) { #ifdef DEBUG_TRACES - std::cout << "p "; - for (auto coordinate: point) { + std::cout << "Delaunay_triangulation_off_visitor_reader::point "; + for (auto coordinate : point) { std::cout << coordinate << " | "; } std::cout << std::endl; @@ -65,62 +73,19 @@ class Delaunay_triangulation_off_flag_visitor_reader { } void maximal_face(const std::vector& face) { - // For alpha shapes, only points are read - } - - void done() { + // For Delaunay Triangulation, only points are read #ifdef DEBUG_TRACES - std::cout << "done" << std::endl; -#endif // DEBUG_TRACES - } -}; - -/** - *@brief Off reader visitor that can be passed to Off_reader to read a Delaunay_triangulation_complex. - */ -template -class Delaunay_triangulation_off_visitor_reader { - Complex& complex_; - // typedef typename Complex::Vertex_handle Vertex_handle; - // typedef typename Complex::Simplex_handle Simplex_handle; - typedef typename Complex::Point Point; - - const bool load_only_points_; - std::vector points_; - // std::vector maximal_faces_; - - public: - explicit Delaunay_triangulation_off_visitor_reader(Complex& complex, bool load_only_points = false) : - complex_(complex), - load_only_points_(load_only_points) { } - - void init(int dim, int num_vertices, int num_faces, int num_edges) { -#ifdef DEBUG_TRACES - std::cout << "init - " << num_vertices << std::endl; -#endif // DEBUG_TRACES - // maximal_faces_.reserve(num_faces); - points_.reserve(num_vertices); - } - - void point(const std::vector& point) { -#ifdef DEBUG_TRACES - std::cout << "p "; - for (auto coordinate: point) { - std::cout << coordinate << " | "; + std::cout << "Delaunay_triangulation_off_visitor_reader::face "; + for (auto vertex : face) { + std::cout << vertex << " | "; } std::cout << std::endl; #endif // DEBUG_TRACES - points_.emplace_back(Point(point.size(), point.begin(), point.end())); - } - - void maximal_face(const std::vector& face) { - // For alpha shapes, only points are read } void done() { - complex_.insert(points_.begin(), points_.end()); #ifdef DEBUG_TRACES - std::cout << "done" << std::endl; + std::cout << "Delaunay_triangulation_off_visitor_reader::done" << std::endl; #endif // DEBUG_TRACES } }; @@ -131,27 +96,23 @@ class Delaunay_triangulation_off_visitor_reader { template class Delaunay_triangulation_off_reader { public: + /** * name_file : file to read * read_complex : complex that will receive the file content * read_only_points : specify true if only the points must be read */ - Delaunay_triangulation_off_reader(const std::string & name_file, Complex& read_complex, bool read_only_points = false, - bool is_flag = false) : valid_(false) { + Delaunay_triangulation_off_reader(const std::string & name_file, Complex& read_complex) : valid_(false) { std::ifstream stream(name_file); if (stream.is_open()) { - if (is_flag) { - // For alpha shapes, only points are read - Delaunay_triangulation_off_flag_visitor_reader off_visitor(read_complex, true); - Off_reader off_reader(stream); - valid_ = off_reader.read(off_visitor); - } else { - // For alpha shapes, only points are read - Delaunay_triangulation_off_visitor_reader off_visitor(read_complex, true); - Off_reader off_reader(stream); - valid_ = off_reader.read(off_visitor); - } + Delaunay_triangulation_off_visitor_reader off_visitor(read_complex); + Off_reader off_reader(stream); + valid_ = off_reader.read(off_visitor); + } else { + std::cerr << "Delaunay_triangulation_off_reader::Delaunay_triangulation_off_reader could not open file " << + name_file << std::endl; } + } /** @@ -160,7 +121,7 @@ class Delaunay_triangulation_off_reader { bool is_valid() const { return valid_; } - + private: bool valid_; }; @@ -168,6 +129,7 @@ class Delaunay_triangulation_off_reader { template class Delaunay_triangulation_off_writer { public: + /** * name_file : file where the off will be written * save_complex : complex that be outputted in the file @@ -191,7 +153,7 @@ class Delaunay_triangulation_off_writer { // Finite cells list for (auto cit = save_complex.finite_full_cells_begin(); cit != save_complex.finite_full_cells_end(); ++cit) { - stream << std::distance(cit->vertices_begin(), cit->vertices_end()) << " "; // Dimension + stream << std::distance(cit->vertices_begin(), cit->vertices_end()) << " "; // Dimension for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { auto vertexHdl = *vit; // auto vertexHdl = std::distance(save_complex.vertices_begin(), *vit) - 1; @@ -201,13 +163,14 @@ class Delaunay_triangulation_off_writer { } stream.close(); } else { - std::cerr << "could not open file " << name_file << std::endl; + std::cerr << "Delaunay_triangulation_off_writer::Delaunay_triangulation_off_writer could not open file " << + name_file << std::endl; } } }; -} // namespace alphashapes +} // namespace alphashapes -} // namespace Gudhi +} // namespace Gudhi #endif // SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_DELAUNAY_TRIANGULATION_OFF_IO_H_ diff --git a/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp b/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp index a90704b6..d5db3bfa 100644 --- a/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp +++ b/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp @@ -90,7 +90,7 @@ BOOST_AUTO_TEST_CASE( Delaunay_triangulation ) { std::string off_file_name("S8_10.off"); std::cout << "========== OFF FILE NAME = " << off_file_name << " ==========" << std::endl; - Gudhi::alphashapes::Delaunay_triangulation_off_reader off_reader(off_file_name, dt, true, true); + Gudhi::alphashapes::Delaunay_triangulation_off_reader off_reader(off_file_name, dt); std::cout << "off_reader.is_valid()=" << off_reader.is_valid() << std::endl; BOOST_CHECK(off_reader.is_valid()); diff --git a/src/Alpha_shapes/test/CMakeLists.txt b/src/Alpha_shapes/test/CMakeLists.txt index a48c1a8f..3cf97b71 100644 --- a/src/Alpha_shapes/test/CMakeLists.txt +++ b/src/Alpha_shapes/test/CMakeLists.txt @@ -1,10 +1,10 @@ cmake_minimum_required(VERSION 2.6) project(GUDHIAlphaShapesTest) -# need CGAL 4.6 -# cmake -DCGAL_DIR=~/workspace/CGAL-4.6-beta1 ../../.. +# need CGAL 4.7 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.7-Ic-41 ../../.. if(CGAL_FOUND) - if (NOT CGAL_VERSION VERSION_LESS 4.6.0) + if (NOT CGAL_VERSION VERSION_LESS 4.7.0) message(STATUS "CGAL version: ${CGAL_VERSION}.") include( ${CGAL_USE_FILE} ) -- cgit v1.2.3