diff --git a/.gitignore b/.gitignore index 4bfdf2b..66de543 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,9 @@ # ignore all files * +# tags file +tags + # unignore directories !*/ diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..7356599 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,9 @@ +[submodule "KDTree"] + path = KDTree + url = git@github.com:crvs/KDTree.git +[submodule "tinyply"] + path = tinyply + url = https://github.com/ddiakopoulos/tinyply +[submodule "tinyobjloader"] + path = tinyobjloader + url = https://github.com/syoyo/tinyobjloader diff --git a/CMakeLists.txt b/CMakeLists.txt index 7e4c837..bbbbcaf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,22 +1,18 @@ -cmake_minimum_required(VERSION 2.8) +cmake_minimum_required(VERSION 3.0) project(coefficient-flow CXX) # requirements: # * Gudhi under /usr/include/gudhi # * Eigen3 under /usr/include/eigen3 # * boost under /usr/include/boost -# * CGAL under /usr/include/CGAL # * libyaml under /usr/local/include/libyaml -# * pcl 1.8 under /usr/include/pcl-1.8 # -# in case Eigen, and/or CGAL, and/or GUDHI are not found by find_package +# in case Eigen, and/or GUDHI are not found by find_package # try setting the paths to them manually # # set(EIGEN3_INCLUDE_DIR "/usr/local/include/eigen3") # include_directories(${EIGEN3_INCLUDE_DIR}) -# set(CGAL_INCLUDE_DIR "/usr/local/include/eigen3") -# include_directories(${CGAL_INCLUDE_DIR}) # # GUDHI version should be 1.3 or higher, provided they don't break the @@ -25,33 +21,32 @@ project(coefficient-flow CXX) # libyaml-cpp this should actually be found in "/usr/include" so there is no # real need for the find_package, but it's always nice to make sure # -# pcl-1.8 version 1.7 doesn't let one read ply files into meshes which is -# rather strange, but it is a thing. -# -find_package(Boost REQUIRED QUIET) -message(STATUS "Boost version: ${BOOST_VERSION}") -message(STATUS "Boost directory: ${BOOST_INCLUDE_DIRS}") +# TODO remove +add_subdirectory(KDTree) + + +add_subdirectory(tinyply) +include_directories("./tinyply/source/") + +add_subdirectory(tinyobjloader) +include_directories("./tinyobjloader") + +find_package(Boost REQUIRED) +message(STATUS "Boost version: ${Boost_VERSION}") +message(STATUS "Boost directory: ${Boost_INCLUDE_DIRS}") +message(STATUS "boost libraries: ${Boost_LIBRARIES}") find_package(Eigen3 3.1.2 REQUIRED QUIET) message(STATUS "Eigen3 version: ${EIGEN3_VERSION}") message(STATUS "Eigen3 directory: ${EIGEN3_INCLUDE_DIR}") include_directories(${EIGEN3_INCLUDE_DIR}) -find_package(CGAL REQUIRED QUIET) -message(STATUS "CGAL version: ${CGAL_VERSION}") -message(STATUS "CGAL directory: ${CGAL_INCLUDE_DIRS}") - find_package(yaml-cpp REQUIRED QUIET) message(STATUS "yaml-cpp version: ${YAML_CPP_VERSION}") message(STATUS "yaml-cpp directory: ${YAML_CPP_INCLUDE_DIR}") include_directories(${YAML_CPP_INCLUDE_DIR}) -find_package(PCL 1.8 REQUIRED QUIET COMPONENTS io) -message(STATUS "PCL version: ${PCL_VERSION}") -message(STATUS "PCL directory: ${PCL_INCLUDE_DIRS}") -include_directories(${PCL_INCLUDE_DIRS}) - find_package(GUDHI 1.3 REQUIRED QUIET) message(STATUS "Gudhi version: ${GUDHI_VERSION}") string(REGEX MATCH "^/home/" LOCAL_INSTALLED "${GUDHI_DIR}") @@ -79,16 +74,22 @@ set(CMAKE_OBJDUMP "/usr/bin/llvm-objdump") set(CMAKE_RANLIB "/usr/bin/llvm-ranlib") set(CMAKE_BUILD_TYPE Release) +# set(CMAKE_BUILD_TYPE Debug) + +include_directories( + "./lib" # where all the library files are + "KDTree" # for the KDTree library + ) -include_directories("./lib") # where all the library files are add_library(scomplex SHARED "./lib/scomplex/simplicial_complex.cpp") add_library(pathsnap SHARED "./lib/scomplex/path_snapper.cpp") +target_link_libraries(pathsnap KDTree) add_executable(qhull2ply "src/make_mesh.cpp") add_executable(yamltest "src/testing_facility.cpp") -target_link_libraries(yamltest yaml-cpp scomplex pathsnap ${PCL_LIBRARIES}) +target_link_libraries(yamltest yaml-cpp scomplex pathsnap tinyply tinyobjloader) message(INFO ${CMAKE_CURRENT_SOURCE_DIR}) message(INFO ${CMAKE_CURRENT_BINARY_DIR}) @@ -98,3 +99,4 @@ add_custom_target(timing_test COMMAND "cp" "${CMAKE_CURRENT_SOURCE_DIR}/test/dumbexample.yaml" "${CMAKE_CURRENT_BINARY_DIR}" COMMAND "${CMAKE_CURRENT_BINARY_DIR}/run_dumb_tests.sh" "1" DEPENDS yamltest qhull2ply "test/run_dumb_tests.sh" "test/dumbexample.yaml") + diff --git a/KDTree b/KDTree new file mode 160000 index 0000000..a881293 --- /dev/null +++ b/KDTree @@ -0,0 +1 @@ +Subproject commit a88129389956023daa04f9e650c7e738b539a3ed diff --git a/README.md b/README.md index cefadda..4d0e808 100644 --- a/README.md +++ b/README.md @@ -14,14 +14,11 @@ In order to compile this repository you will need the following packages: - `clang (v == 3.8)` - `Eigen (v == 3.3.2)` - `Boost (v >= 1.54.0)` -- `CGAL (v >= 4.9.1)` - `Gudhi (v >= 1.3.1)` -- `pcl (v >= 1.8.0)` - `libyaml` **Note:** - In case you are using a different version of `clang`, everything should still compile without problems, but you will need to change the appropriate lines in `CMakeLists.txt`. -- pcl needs to be version 1.8 or higher because 1.7 does not support reading ply files into `pcl::PolygonMesh` elements. ## Running diff --git a/lib/scomplex/chain_calc.hpp b/lib/scomplex/chain_calc.hpp index a67c853..32ceb69 100644 --- a/lib/scomplex/chain_calc.hpp +++ b/lib/scomplex/chain_calc.hpp @@ -41,13 +41,13 @@ bounding_chain::bounding_chain(std::shared_ptr sc) { } bounding_chain::bounding_chain(simplicial_complex& sc) { - s_comp.reset(new simplicial_complex(sc)); + s_comp = std::make_shared(sc); populate_matrices(); } bounding_chain::bounding_chain(std::vector& points, std::vector& tris) { - s_comp.reset(new simplicial_complex(points, tris)); + s_comp = std::make_shared(points,tris); populate_matrices(); } diff --git a/lib/scomplex/graph_utils.hpp b/lib/scomplex/graph_utils.hpp index aa961a2..9ef3513 100644 --- a/lib/scomplex/graph_utils.hpp +++ b/lib/scomplex/graph_utils.hpp @@ -110,7 +110,7 @@ std::vector shortest_path(graph_t& g, vertex_t s, vertex_t t) { std::vector complete_path(graph_t& g, std::vector vec) { - std::vector::iterator s{vec.begin()}; + std::vector::iterator s = vec.begin(); std::vector::iterator t = std::next(s); // start calculatinng the full path diff --git a/lib/scomplex/nn_utils.hpp b/lib/scomplex/nn_utils.hpp deleted file mode 100644 index 86294c4..0000000 --- a/lib/scomplex/nn_utils.hpp +++ /dev/null @@ -1,98 +0,0 @@ -#pragma once - -#include - -#include -#include -#include - -#include -#include -#include -#include - -#include - -namespace gsimp { - -/* -types used outside: - tree_t -functions provided: - make_tree(tree_t& tree, std::vector& point_list); - nearest_neighbor(tree_t& tree, point_t point); - nearest_neighbor_index(tree_t& tree, point_t point); - snap_points(tree_t& tree, std::vector points); - snap_points_to_indexes( tree_t& tree, std::vector points); -*/ - -// kd-tree point type, d-dimensional point -typedef CGAL::Cartesian_d::Point_d point_d; - -typedef CGAL::Search_traits_d> cartesian_traits; - -typedef CGAL::Search_traits_adapter< // define the search tratis - boost::tuple, // pairing d-point with its index - CGAL::Nth_of_tuple_property_map< // custom kernel - 0, // position of the points to compare - boost::tuple>, // type tha will be used - cartesian_traits> // traits to use - traits; // name - -// type of search tree -typedef CGAL::Orthogonal_k_neighbor_search neighbor_search_t; - -typedef neighbor_search_t::Tree tree_t; - -void make_tree(tree_t& tree, std::vector& point_list) { - // tree.reserve(point_list.size()); - size_t ind = 0; - for (auto pt : point_list){ - point_d pt_d(pt.size(), pt.begin(), pt.end()); - tree.insert(boost::make_tuple(pt_d, ind)); - ++ind; - } -} - -// given a point snap it to the tree -point_t nearest_neighbor(tree_t& tree, point_t point) { - point_d pt(point.size(), point.begin(), point.end()); - // "1" refers to the number of nearest neighbors to find - neighbor_search_t search(tree, pt, 1); - point_d result_pt; - size_t ind; - boost::tie(result_pt, ind) = search.begin()->first; - point_t point_vec(result_pt.cartesian_begin(), result_pt.cartesian_end()); - return point_vec; -} - -// given a vector of points snap them to the tree -std::vector snap_points(tree_t& tree, std::vector points) { - std::vector snapped_points; - for (auto point : points) { - snapped_points.push_back(nearest_neighbor(tree, point)); - } - return snapped_points; -} - -size_t nearest_neighbor_index(tree_t& tree, point_t point) { - point_d pt(point.size(), point.begin(), point.end()); - // "1" refers to the number of nearest neighbors to find - neighbor_search_t search(tree, pt, 1); - point_d result_pt; - size_t ind; - boost::tie(result_pt, ind) = search.begin()->first; - return ind; -} - -std::vector snap_points_to_indexes( // - tree_t& tree, // - std::vector points) { - std::vector snapped_points; - for (auto point : points) { - snapped_points.push_back( // - nearest_neighbor_index(tree, point)); - } - return snapped_points; -} -}; diff --git a/lib/scomplex/path_snapper.cpp b/lib/scomplex/path_snapper.cpp index 5653bb9..8303578 100644 --- a/lib/scomplex/path_snapper.cpp +++ b/lib/scomplex/path_snapper.cpp @@ -1,46 +1,52 @@ -#include #include #include -#include -#include -#include -#include + +#include "scomplex/graph_utils.hpp" +#include "scomplex/path_snapper.hpp" +#include "scomplex/simplicial_complex.hpp" +#include "scomplex/types.hpp" + +#include "KDTree.hpp" namespace gsimp { struct path_snapper::impl { friend class clusterer; - tree_t point_tree; // defined in nn_utils.hpp + KDTree point_tree; graph_t vertex_graph; // defined in graph_utils.hpp - std::shared_ptr s_comp; + std::shared_ptr< simplicial_complex > s_comp; - impl(std::shared_ptr p_sc) { + impl(std::shared_ptr< simplicial_complex > p_sc) { s_comp = p_sc; auto points = s_comp->get_points(); - make_tree(point_tree, points); + point_tree = KDTree(points); vertex_graph = calculate_one_skelleton_graph(*s_comp); }; impl(simplicial_complex& sc) { - s_comp.reset(new simplicial_complex(sc)); + s_comp = std::make_shared< simplicial_complex >(sc); auto points = s_comp->get_points(); - make_tree(point_tree, points); + point_tree = KDTree(points); vertex_graph = calculate_one_skelleton_graph(*s_comp); } - impl(std::vector& pts, std::vector& cells) { - s_comp = std::shared_ptr( - new simplicial_complex(pts, cells)); - make_tree(point_tree, pts); + impl(std::vector< point_t >& pts, std::vector< cell_t >& cells) { + s_comp = std::make_shared< simplicial_complex >(pts, cells); + point_tree = KDTree(pts); vertex_graph = calculate_one_skelleton_graph(*s_comp); } ~impl(){}; - std::vector snap_path(std::vector path) { - auto way_points = snap_points_to_indexes(point_tree, path); + std::vector< size_t > snap_path(std::vector< point_t > path) { + std::vector< size_t > way_points; + for (point_t pt : path) { + size_t pti = point_tree.nearest_index(pt); + way_points.push_back(pti); + } + // std::cout << "processed path: \n"; // for (size_t i = 0 ; i < way_points.size() ; ++i ) { // for (auto c : path[i] ) std::cout << c << " "; @@ -53,39 +59,40 @@ struct path_snapper::impl { return snapped_path; } - std::vector> index_pairs(std::vector path) { + std::vector< std::pair< cell_t, int > > index_pairs( + std::vector< point_t > path) { auto vertex_path = snap_path(path); - std::vector> pair_seq; + std::vector< std::pair< cell_t, int > > pair_seq; for (auto it = vertex_path.begin(); std::next(it) != vertex_path.end(); ++it) { // // - size_t src(*it), trg(*(std::next(it))); // - if (src < trg) { // - pair_seq.push_back( // - std::make_pair, int>( // - std::vector{src, trg}, 1)); // - } else { // - pair_seq.push_back( // - std::make_pair, int>( // - std::vector{trg, src}, -1)); // - } // + size_t src(*it), trg(*(std::next(it))); // + if (src < trg) { // + pair_seq.push_back( // + std::make_pair< std::vector< size_t >, int >( // + std::vector< size_t >{src, trg}, 1)); // + } else { // + pair_seq.push_back( // + std::make_pair< std::vector< size_t >, int >( // + std::vector< size_t >{trg, src}, -1)); // + } // } return pair_seq; } }; -path_snapper::path_snapper(std::shared_ptr sc) { - p_impl = std::shared_ptr(new impl(sc)); +path_snapper::path_snapper(std::shared_ptr< simplicial_complex > sc) { + p_impl = std::make_shared< impl >(sc); }; path_snapper::path_snapper(simplicial_complex& sc) { - p_impl = std::shared_ptr(new impl(sc)); + p_impl = std::make_shared< impl >(sc); } -path_snapper::path_snapper(std::vector& pts, - std::vector& cells) { - p_impl = std::shared_ptr(new impl(pts, cells)); +path_snapper::path_snapper(std::vector< point_t >& pts, + std::vector< cell_t >& cells) { + p_impl = std::make_shared< impl >(pts, cells); } path_snapper::~path_snapper() {} @@ -97,61 +104,64 @@ path_snapper& path_snapper::operator=(const path_snapper& other) { return *this; } -std::vector path_snapper::snap_path_to_points( - std::vector path) { +std::vector< point_t > path_snapper::snap_path_to_points( + std::vector< point_t > path) { auto index_path = p_impl->snap_path(path); - std::vector point_path; + std::vector< point_t > point_path; for (size_t p : index_path) point_path.push_back(p_impl->s_comp->get_points().at(p)); return point_path; } -std::vector path_snapper::snap_path_to_indices( - std::vector path) { +std::vector< size_t > path_snapper::snap_path_to_indices( + std::vector< point_t > path) { return p_impl->snap_path(path); } -chain_v path_snapper::snap_path_to_v_chain(std::vector path) { +chain_v path_snapper::snap_path_to_v_chain(std::vector< point_t > path) { auto pair_seq = p_impl->index_pairs(path); chain_v rep = p_impl->s_comp->new_v_chain(1); for (auto pair : pair_seq) { - auto cell = std::get<0>(pair); + auto cell = std::get< 0 >(pair); if (cell[0] != cell[1]) { - size_t index = p_impl->s_comp->cell_to_index(std::get<0>(pair)); - chain_val(rep, index) += std::get<1>(pair); + size_t index = p_impl->s_comp->cell_to_index(std::get< 0 >(pair)); + chain_val(rep, index) += std::get< 1 >(pair); } } return rep; } -chain_t path_snapper::snap_path_to_chain(std::vector path) { +chain_t path_snapper::snap_path_to_chain(std::vector< point_t > path) { auto pair_seq = p_impl->index_pairs(path); chain_t rep = p_impl->s_comp->new_chain(1); for (auto pair : pair_seq) { - auto cell = std::get<0>(pair); + auto cell = std::get< 0 >(pair); if (cell[0] != cell[1]) { - size_t index = p_impl->s_comp->cell_to_index(std::get<0>(pair)); - chain_val(rep, index) += std::get<1>(pair); + size_t index = p_impl->s_comp->cell_to_index(std::get< 0 >(pair)); + chain_val(rep, index) += std::get< 1 >(pair); } } return rep; } -std::vector path_snapper::index_sequence_to_point( - std::vector ind_path) { - std::vector pt_path; +std::vector< point_t > path_snapper::index_sequence_to_point( + std::vector< size_t > ind_path) { + std::vector< point_t > pt_path; for (size_t ind : ind_path) pt_path.push_back(p_impl->s_comp->get_points().at(ind)); return pt_path; } -std::vector path_snapper::point_sequence_to_index( - std::vector pt_path) { - std::vector ind_path = - snap_points_to_indexes(p_impl->point_tree, pt_path); +std::vector< size_t > path_snapper::point_sequence_to_index( + std::vector< point_t > pt_path) { + std::vector< size_t > ind_path; + for (point_t pt : pt_path) { + ind_path.push_back(p_impl->point_tree.nearest_index(pt)); + } return ind_path; } -chain_v path_snapper::index_sequence_to_v_chain(std::vector ind_path) { +chain_v path_snapper::index_sequence_to_v_chain( + std::vector< size_t > ind_path) { auto it = ind_path.begin(); chain_v chain = p_impl->s_comp->new_v_chain(1); for (; std::next(it) != ind_path.end(); ++it) { @@ -161,7 +171,7 @@ chain_v path_snapper::index_sequence_to_v_chain(std::vector ind_path) { return chain; } -chain_t path_snapper::index_sequence_to_chain(std::vector ind_path) { +chain_t path_snapper::index_sequence_to_chain(std::vector< size_t > ind_path) { auto it = ind_path.begin(); chain_t chain = p_impl->s_comp->new_chain(1); for (; std::next(it) != ind_path.end(); ++it) { @@ -171,11 +181,11 @@ chain_t path_snapper::index_sequence_to_chain(std::vector ind_path) { return chain; } -chain_t path_snapper::point_sequence_to_chain(std::vector pt_path) { +chain_t path_snapper::point_sequence_to_chain(std::vector< point_t > pt_path) { return index_sequence_to_chain(point_sequence_to_index(pt_path)); } -std::shared_ptr path_snapper::get_underlying_complex() { +std::shared_ptr< simplicial_complex > path_snapper::get_underlying_complex() { return p_impl->s_comp; } -}; +}; // namespace gsimp diff --git a/lib/scomplex/path_snapper.hpp b/lib/scomplex/path_snapper.hpp index da03524..4d19970 100644 --- a/lib/scomplex/path_snapper.hpp +++ b/lib/scomplex/path_snapper.hpp @@ -1,34 +1,35 @@ #pragma once -#include -#include +#include "scomplex/simplicial_complex.hpp" +#include "scomplex/types.hpp" namespace gsimp { class path_snapper { private: struct impl; - std::shared_ptr p_impl; + std::shared_ptr< impl > p_impl; public: // constructors and such - path_snapper(std::vector&, std::vector&); - path_snapper(std::shared_ptr); - path_snapper(simplicial_complex&); + path_snapper(std::vector< point_t >&, std::vector< cell_t >&); + explicit path_snapper(std::shared_ptr< simplicial_complex >); + explicit path_snapper(simplicial_complex&); ~path_snapper(); path_snapper(path_snapper&); path_snapper& operator=(const path_snapper&); // the things we want to do - std::vector snap_path_to_points(std::vector); - std::vector snap_path_to_indices(std::vector); - chain_t snap_path_to_chain(std::vector); - chain_v snap_path_to_v_chain(std::vector); + std::vector< point_t > snap_path_to_points(std::vector< point_t >); + std::vector< size_t > snap_path_to_indices(std::vector< point_t >); + chain_t snap_path_to_chain(std::vector< point_t >); + chain_v snap_path_to_v_chain(std::vector< point_t >); // interconversion - std::vector index_sequence_to_point(std::vector); - std::vector point_sequence_to_index(std::vector); - chain_t index_sequence_to_chain(std::vector); - chain_v index_sequence_to_v_chain(std::vector); - chain_t point_sequence_to_chain(std::vector); - std::shared_ptr get_underlying_complex(); -}; -}; + std::vector< point_t > index_sequence_to_point(std::vector< size_t >); + std::vector< size_t > point_sequence_to_index(std::vector< point_t >); + chain_t index_sequence_to_chain(std::vector< size_t >); + chain_v index_sequence_to_v_chain(std::vector< size_t >); + chain_t point_sequence_to_chain(std::vector< point_t >); + std::shared_ptr< simplicial_complex > get_underlying_complex(); +}; // class path_snapper + +}; // namespace gsimp diff --git a/lib/scomplex/simplicial_complex.cpp b/lib/scomplex/simplicial_complex.cpp index b9d580b..ef18154 100644 --- a/lib/scomplex/simplicial_complex.cpp +++ b/lib/scomplex/simplicial_complex.cpp @@ -3,34 +3,34 @@ #include // for debuging purposes -#include #include -#include +#include #include +#include -#include -#include // smart pointers -#include #include +#include +#include // smart pointers +#include #include namespace gsimp { // hasse diagram structures struct hasse_node { - std::pair handle; - std::vector> cofaces; + std::pair< int, size_t > handle; + std::vector< std::shared_ptr< hasse_node > > cofaces; - hasse_node(std::pair _handle) { handle = _handle; } + hasse_node(std::pair< int, size_t > _handle) { handle = _handle; } - void add_coface(std::shared_ptr& node) { + void add_coface(std::shared_ptr< hasse_node >& node) { cofaces.push_back(node); } }; struct hasse_diag { - std::vector>> cells; - std::vector> cells_c; + std::vector< std::vector< std::shared_ptr< hasse_node > > > cells; + std::vector< std::vector< bool > > cells_c; hasse_diag(){}; @@ -38,24 +38,24 @@ struct hasse_diag { // vector of dimension d for (int d = 0; d <= s_comp.dimension(); d++) { size_t lv_s = s_comp.get_level_size(d); - std::vector> lv(lv_s); + std::vector< std::shared_ptr< hasse_node > > lv(lv_s); cells.push_back(lv); - std::vector c(lv_s, false); + std::vector< bool > c(lv_s, false); cells_c.push_back(c); } for (int d = s_comp.dimension(); d > 0; d--) { for (size_t face_i = 0; face_i < s_comp.get_level_size(d); face_i++) { - std::shared_ptr face_ptr = + std::shared_ptr< hasse_node > face_ptr = get_face(d, face_i, true); - std::vector> faces = + std::vector< std::pair< int, size_t > > faces = s_comp.get_bdry_and_ind_index(d, face_i); for (auto b_face_i : faces) { - size_t b_face = std::get<1>(b_face_i); - std::shared_ptr b_face_ptr = + size_t b_face = std::get< 1 >(b_face_i); + std::shared_ptr< hasse_node > b_face_ptr = get_face(d - 1, b_face, true); b_face_ptr->add_coface(face_ptr); } @@ -65,25 +65,26 @@ struct hasse_diag { cells_c.clear(); } - std::shared_ptr get_face(int d, size_t face_i, - bool building = false) { - std::shared_ptr face; + std::shared_ptr< hasse_node > get_face(int d, size_t face_i, + bool building = false) { + std::shared_ptr< hasse_node > face; if (building && cells_c[d][face_i]) face = cells[d][face_i]; else { if (building) cells_c[d][face_i] = true; - cells[d][face_i] = std::shared_ptr( - new hasse_node(std::pair(d, face_i))); + cells[d][face_i] = std::make_shared< hasse_node >( + std::pair< int, size_t >(d, face_i)); face = cells[d][face_i]; } return face; } - std::vector get_coface_i(int d, size_t face_i) { + std::vector< size_t > get_coface_i(int d, size_t face_i) { auto node = cells[d][face_i]; - std::vector cofaces; - for (auto v : node->cofaces) cofaces.push_back(std::get<1>(v->handle)); + std::vector< size_t > cofaces; + for (auto v : node->cofaces) + cofaces.push_back(std::get< 1 >(v->handle)); return cofaces; } }; @@ -94,22 +95,21 @@ struct simplicial_complex::impl { typedef size_t Vertex_handle; }; - typedef Gudhi::Simplex_tree simp_tree; - typedef Gudhi::Simplex_tree::Simplex_handle simp_handle; - typedef std::vector level_t; - typedef std::vector> levels_t; + typedef Gudhi::Simplex_tree< SimpleOptions > simp_tree; + typedef Gudhi::Simplex_tree< SimpleOptions >::Simplex_handle simp_handle; + typedef std::vector< simp_handle* > level_t; + typedef std::vector< std::unique_ptr< level_t > > levels_t; // member variables - std::vector points; + std::vector< point_t > points; simp_tree simplices; - std::vector boundary_matrices; + std::vector< matrix_t > boundary_matrices; levels_t levels; bool has_hasse; hasse_diag incidence; - - impl(std::vector& arg_points, std::vector& arg_tris) + impl(std::vector< point_t >& arg_points, std::vector< cell_t >& arg_tris) : points(arg_points) { // create the simplex tree for (auto tri : arg_tris) { @@ -121,10 +121,10 @@ struct simplicial_complex::impl { } // assign a key to each simplex in each level - std::vector count(simplices.dimension() + 1, 0); + std::vector< size_t > count(simplices.dimension() + 1, 0); for (int i = 0; i < simplices.dimension() + 1; ++i) { level_t* level = new level_t(); - levels.push_back(std::unique_ptr(level)); + levels.push_back(std::unique_ptr< level_t >(level)); } auto simplex_range = simplices.complex_simplex_range(); @@ -164,9 +164,9 @@ struct simplicial_complex::impl { return pow(-1, orient); } - std::vector dedupe_vec(std::vector& vec) { - std::set no_reps; - std::vector no_reps_list; + std::vector< size_t > dedupe_vec(std::vector< size_t >& vec) { + std::set< size_t > no_reps; + std::vector< size_t > no_reps_list; for (size_t el : vec) { no_reps.insert(el); } @@ -177,7 +177,7 @@ struct simplicial_complex::impl { } void calculate_matrices() { - boundary_matrices = std::vector(); + boundary_matrices = std::vector< matrix_t >(); for (int k = 0; k < simplices.dimension(); k++) { boundary_matrices.push_back( matrix_t(get_level_size(k), get_level_size(k + 1))); @@ -192,8 +192,8 @@ struct simplicial_complex::impl { } } - std::vector get_level(int level) { - std::vector level_cells; + std::vector< cell_t > get_level(int level) { + std::vector< cell_t > level_cells; for (auto simp : *levels[level]) { cell_t v_simp; for (auto v : simplices.simplex_vertex_range(*simp)) { @@ -223,45 +223,45 @@ struct simplicial_complex::impl { }; // struct impl -std::vector> simplicial_complex::get_bdry_and_ind( +std::vector< std::pair< int, cell_t > > simplicial_complex::get_bdry_and_ind( cell_t cell) { impl::simp_handle simp = p_impl->cell_to_handle(cell); auto c_boundary = p_impl->simplices.boundary_simplex_range(simp); - std::vector> boundary_and_indices; + std::vector< std::pair< int, cell_t > > boundary_and_indices; for (auto face : c_boundary) boundary_and_indices.push_back( // - std::make_pair( // + std::make_pair< int, cell_t >( // p_impl->boundary_index(face, simp), // p_impl->handle_to_cell(face))); // return boundary_and_indices; }; -std::vector> simplicial_complex::get_bdry_and_ind_index( - int d, size_t cell) { +std::vector< std::pair< int, size_t > > +simplicial_complex::get_bdry_and_ind_index(int d, size_t cell) { impl::simp_handle simp = p_impl->index_to_handle(d, cell); auto c_boundary = p_impl->simplices.boundary_simplex_range(simp); - std::vector> boundary_and_indices; + std::vector< std::pair< int, size_t > > boundary_and_indices; for (auto face : c_boundary) boundary_and_indices.push_back( // - std::make_pair( // + std::make_pair< int, size_t >( // p_impl->boundary_index(face, simp), // p_impl->handle_to_index(face))); // return boundary_and_indices; }; -std::vector simplicial_complex::cell_boundary_index(int d, - size_t cell) { +std::vector< size_t > simplicial_complex::cell_boundary_index(int d, + size_t cell) { impl::simp_handle simp = p_impl->index_to_handle(d, cell); auto c_boundary = p_impl->simplices.boundary_simplex_range(simp); - std::vector s_boundary; + std::vector< size_t > s_boundary; for (auto c : c_boundary) s_boundary.push_back(p_impl->handle_to_index(c)); return s_boundary; } -std::vector simplicial_complex::cell_boundary(cell_t cell) { +std::vector< cell_t > simplicial_complex::cell_boundary(cell_t cell) { impl::simp_handle simp = p_impl->cell_to_handle(cell); auto c_boundary = p_impl->simplices.boundary_simplex_range(simp); - std::vector s_boundary; + std::vector< cell_t > s_boundary; for (auto c : c_boundary) s_boundary.push_back(p_impl->handle_to_cell(c)); return s_boundary; } @@ -278,9 +278,9 @@ int simplicial_complex::boundary_inclusion_index(int d1, size_t s1, // p_impl->cell_to_handle(c2)); }; -std::vector> simplicial_complex::get_cof_and_ind( +std::vector< std::pair< int, cell_t > > simplicial_complex::get_cof_and_ind( cell_t cell) { - std::vector> c_cofaces; + std::vector< std::pair< int, cell_t > > c_cofaces; for (auto face : get_cofaces(cell)) c_cofaces.push_back( // std::make_pair( // @@ -288,9 +288,9 @@ std::vector> simplicial_complex::get_cof_and_ind( return c_cofaces; } -std::vector> simplicial_complex::get_cof_and_ind_index( - int d, size_t c) { - std::vector> c_cofaces; +std::vector< std::pair< int, size_t > > +simplicial_complex::get_cof_and_ind_index(int d, size_t c) { + std::vector< std::pair< int, size_t > > c_cofaces; for (auto face : get_cofaces_index(d, c)) c_cofaces.push_back( // std::make_pair( // @@ -302,9 +302,14 @@ int simplicial_complex::get_level_size(int level) { return p_impl->get_level_size(level); } -simplicial_complex::simplicial_complex(std::vector& arg_points, - std::vector& arg_tris) { - p_impl = std::shared_ptr(new impl(arg_points, arg_tris)); +simplicial_complex::simplicial_complex(std::vector< cell_t >& arg_tris) { + std::vector< point_t > points = {}; + p_impl = std::make_shared< impl >(points, arg_tris); +} + +simplicial_complex::simplicial_complex(std::vector< point_t >& arg_points, + std::vector< cell_t >& arg_tris) { + p_impl = std::make_shared< impl >(arg_points, arg_tris); } simplicial_complex::simplicial_complex(const simplicial_complex& other) { @@ -319,13 +324,15 @@ simplicial_complex& simplicial_complex::operator=( simplicial_complex::~simplicial_complex() {} -std::vector simplicial_complex::get_points() { return p_impl->points; } +std::vector< point_t > simplicial_complex::get_points() { + return p_impl->points; +} point_t simplicial_complex::get_point(size_t index) { return p_impl->points[index]; } -std::vector simplicial_complex::get_level(int level) { +std::vector< cell_t > simplicial_complex::get_level(int level) { return p_impl->get_level(level); } @@ -352,7 +359,8 @@ size_t simplicial_complex::cell_to_index(cell_t simp) { return p_impl->simplices.key(sh); } -std::vector simplicial_complex::get_cofaces_index(int d, size_t face) { +std::vector< size_t > simplicial_complex::get_cofaces_index(int d, + size_t face) { // codimension 1 faces if (!p_impl->has_hasse) { calculate_hasse(); @@ -362,12 +370,12 @@ std::vector simplicial_complex::get_cofaces_index(int d, size_t face) { return s_cofaces; } -std::vector simplicial_complex::get_cofaces(cell_t face) { +std::vector< cell_t > simplicial_complex::get_cofaces(cell_t face) { if (!p_impl->has_hasse) { calculate_hasse(); p_impl->has_hasse = true; } - std::vector s_cofaces; + std::vector< cell_t > s_cofaces; auto face_i = cell_to_index(face); int d = face.size() - 1; // codimension 1 faces @@ -380,9 +388,8 @@ std::vector simplicial_complex::get_cofaces(cell_t face) { } chain_v simplicial_complex::new_v_chain(int d) { - std::vector v(get_level_size(d),0); - return chain_v(d,v); - + std::vector< double > v(get_level_size(d), 0); + return chain_v(d, v); } chain_t simplicial_complex::new_chain(int d) { vector_t v(get_level_size(d)); @@ -393,4 +400,4 @@ void simplicial_complex::calculate_hasse() { p_impl->has_hasse = true; p_impl->incidence = hasse_diag(*this); } -}; +}; // namespace gsimp diff --git a/lib/scomplex/simplicial_complex.hpp b/lib/scomplex/simplicial_complex.hpp index b21180d..3bdb28e 100644 --- a/lib/scomplex/simplicial_complex.hpp +++ b/lib/scomplex/simplicial_complex.hpp @@ -20,6 +20,7 @@ class simplicial_complex { void calculate_hasse(); // constructor (no default) + simplicial_complex(std::vector&); simplicial_complex(std::vector&, std::vector&); simplicial_complex(const simplicial_complex&); simplicial_complex& operator=(const simplicial_complex&); diff --git a/src/testing_facility.cpp b/src/testing_facility.cpp index 17e3366..c8e0905 100644 --- a/src/testing_facility.cpp +++ b/src/testing_facility.cpp @@ -1,101 +1,124 @@ -#include -#include - -#include #include +#include #include #include #include -#include -#include -#include -#include -#include -#include -#include +#include "scomplex/chain_calc.hpp" +#include "scomplex/coeff_flow.hpp" +#include "scomplex/path_snapper.hpp" +#include "scomplex/plywriter.hpp" +#include "scomplex/qhull_parsing.hpp" +#include "scomplex/simplicial_complex.hpp" +#include "scomplex/types.hpp" + +#include "tiny_obj_loader.h" +#include "tinyply.h" #include -void call_single_cycle_test(std::string, std::vector, - std::vector>, size_t, bool); +void call_single_cycle_test(std::string, std::vector< size_t >, + std::vector< std::vector< double > >, size_t, bool); -void call_single_cycle_test(std::string, std::vector>, +void call_single_cycle_test(std::string, std::vector< std::vector< double > >, size_t, bool); -void call_single_cycle_test(std::string, std::vector>, +void call_single_cycle_test(std::string, std::vector< std::vector< double > >, size_t); -void call_single_cycle_test(std::string, std::vector, size_t, bool); +void call_single_cycle_test(std::string, std::vector< size_t >, size_t, bool); -void call_single_cycle_test(std::string, std::vector, size_t); +void call_single_cycle_test(std::string, std::vector< size_t >, size_t); int main(int argc, char* argv[]) { YAML::Node input = YAML::LoadFile(argv[1]); std::string paths_file; std::string complex_file; - std::vector> cycle_points; - std::vector cycle_index_points; + std::vector< std::vector< double > > cycle_points; + std::vector< size_t > cycle_index_points; bool in_plane; if (input["single_cycle_test"]) { complex_file = - input["single_cycle_test"]["complex_file"].as(); - in_plane = input["single_cycle_test"]["in_plane"].as(); + input["single_cycle_test"]["complex_file"].as< std::string >(); + in_plane = input["single_cycle_test"]["in_plane"].as< bool >(); size_t null_face; - if (!in_plane) - null_face = input["single_cycle_test"]["null_face"].as(); - else + if (!in_plane) { + null_face = input["single_cycle_test"]["null_face"].as< size_t >(); + } else { null_face = 0; + } if (input["single_cycle_test"]["cycle_index_points"]) { cycle_index_points = - input["single_cycle_test"]["cycle_index_points"].as>(); + input["single_cycle_test"]["cycle_index_points"] + .as< std::vector< size_t > >(); call_single_cycle_test(complex_file, cycle_index_points, null_face, in_plane); } else { cycle_points = input["single_cycle_test"]["cycle_points"] - .as>>(); + .as< std::vector< std::vector< double > > >(); cycle_points.push_back(cycle_points[0]); - call_single_cycle_test(complex_file, cycle_points, null_face, in_plane); + call_single_cycle_test(complex_file, cycle_points, null_face, + in_plane); } } - } void call_single_cycle_test(std::string complex_file, - std::vector> cycle_points, + std::vector< std::vector< double > > cycle_points, size_t null_face, bool in_plane) { call_single_cycle_test(complex_file, {}, cycle_points, null_face, in_plane); } void call_single_cycle_test(std::string complex_file, - std::vector> cycle_points, + std::vector< std::vector< double > > cycle_points, size_t null_face) { call_single_cycle_test(complex_file, {}, cycle_points, null_face, false); } void call_single_cycle_test(std::string complex_file, - std::vector cycle_indices, size_t null_face, - bool in_plane) { + std::vector< size_t > cycle_indices, + size_t null_face, bool in_plane) { call_single_cycle_test(complex_file, cycle_indices, {}, null_face, in_plane); } void call_single_cycle_test(std::string complex_file, - std::vector cycle_indices, + std::vector< size_t > cycle_indices, size_t null_face) { call_single_cycle_test(complex_file, cycle_indices, {}, null_face, false); } +template < typename T > +struct vec3 { + T x, y, z; + + std::vector< double > point() { + std::vector< double > vec_; + vec_.push_back(double(x)); + vec_.push_back(double(y)); + vec_.push_back(double(z)); + return vec_; + } + + std::vector< size_t > simp() { + std::vector< size_t > vec_; + vec_.push_back(size_t(x)); + vec_.push_back(size_t(y)); + vec_.push_back(size_t(z)); + return vec_; + } +}; + // // mother of all single_cycle_tests // void call_single_cycle_test(std::string complex_file, - std::vector cycle_indices, - std::vector> cycle_points, + std::vector< size_t > cycle_indices, + std::vector< std::vector< double > > cycle_points, size_t null_face, bool in_plane) { std::cout << "performing the single cycle test with parameters:\n"; std::cout << "file containing the mesh: " << complex_file << "\n"; @@ -110,40 +133,84 @@ void call_single_cycle_test(std::string complex_file, if (!in_plane) std::cout << "index of face to be null: " << null_face << "\n\n"; - std::vector points_v; - std::vector cells_v; + std::vector< point_t > points_v; + std::vector< cell_t > cells_v; clock_t t0, t1; t0 = clock(); - // PCL shenanigans - - pcl::PolygonMesh::Ptr mesh(new pcl::PolygonMesh{}); std::ifstream file(complex_file); std::string meshtype; - std::getline(file,meshtype); + std::getline(file, meshtype); + if (meshtype == "ply") { - pcl::PLYReader Reader; - Reader.read(complex_file, *mesh); - } else { - pcl::OBJReader Reader; - Reader.read(complex_file, *mesh); - } + tinyply::PlyFile plyMeshFile; + plyMeshFile.parse_header(file); + + std::shared_ptr< tinyply::PlyData > vertices, faces; + try { + vertices = plyMeshFile.request_properties_from_element( + "vertex", {"x", "y", "z"}); + faces = plyMeshFile.request_properties_from_element( + "face", {"vertex_indices"}); + } catch (const std::exception& e) { + std::cerr << "tinyply exception: " << e.what() << std::endl; + } + plyMeshFile.read(file); + { + const size_t numVerticesBytes = vertices->buffer.size_bytes(); + std::vector< vec3< float > > verts(vertices->count); + std::memcpy(verts.data(), vertices->buffer.get(), numVerticesBytes); + for (vec3< float > v : verts) { + points_v.push_back(v.point()); + } + } - pcl::PointCloud::Ptr cloud(new pcl::PointCloud{}); - pcl::fromPCLPointCloud2(mesh->cloud, *cloud); + { + const size_t numFacesBytes = faces->buffer.size_bytes(); + std::vector< vec3< uint32_t > > faces_(faces->count); + std::memcpy(faces_.data(), faces->buffer.get(), numFacesBytes); + for (vec3< uint32_t > f : faces_) { + cells_v.push_back(f.simp()); + } + } - for (auto poly : mesh->polygons) { - cells_v.push_back({poly.vertices[0], poly.vertices[1], poly.vertices[2]}); - } + } else { + tinyobj::attrib_t attrib; + std::vector< tinyobj::shape_t > shapes; + std::vector< tinyobj::material_t > materials; + const char* basepath = {};//= std::string("").c_str(); + std::string err; + const char* filename = complex_file.c_str(); + + bool triangulate = false; + + tinyobj::LoadObj(&attrib, &shapes, &materials, &err, filename, basepath, + triangulate); + + std::vector< double > vert; + for (size_t i = 0; i < attrib.vertices.size(); i++) { + vert.push_back(attrib.vertices.at(i)); + if (i % 3 == 2) { + points_v.push_back(vert); + vert.clear(); + } + } - for (auto pt : cloud->points) { - points_v.push_back({pt.x,pt.y,pt.z}); + std::vector< size_t > face; + tinyobj::shape_t shape = shapes.at(0); + for (size_t i = 0; i < shape.mesh.indices.size(); i++) { + tinyobj::index_t idx = shape.mesh.indices.at(i); + face.push_back(idx.vertex_index); + if (i % 3 == 2) { + cells_v.push_back(face); + face.clear(); + } + } } - // end PCL shenanigans t1 = clock(); std::cout << "mesh has " << points_v.size() << " vertices and " << cells_v.size() << " faces\n"; @@ -169,8 +236,8 @@ void call_single_cycle_test(std::string complex_file, } t0 = clock(); - std::shared_ptr s_comp( - new gsimp::simplicial_complex(points_v, cells_v)); + std::shared_ptr< gsimp::simplicial_complex > s_comp = + std::make_shared< gsimp::simplicial_complex >(points_v, cells_v); t1 = clock(); std::cout << "created complex in " << t1 - t0 << " clock cycles " << float(t1 - t0) / CLOCKS_PER_SEC << " seconds\n"; @@ -178,34 +245,39 @@ void call_single_cycle_test(std::string complex_file, std::cout << "complex comosition:\n"; std::cout << " number of faces: " << s_comp->get_level_size(2) << "\n"; std::cout << " number of edges: " << s_comp->get_level_size(1) << "\n"; - std::cout << " number of vertices: " << s_comp->get_level_size(0) << "\n"; + std::cout << " number of vertices: " << s_comp->get_level_size(0) + << "\n"; // re sort the triangles according to level 2 { - typedef std::pair sortable; - std::vector pairing; - for (auto tri : cells_v) { - // get index of the triangle - size_t ind; - ind = s_comp->cell_to_index(tri); - pairing.push_back({ind,tri}); - } - std::sort(pairing.begin(),pairing.end(),[](const sortable& x,const sortable& y){return x.first < y.first;}); - cells_v.clear(); - for (auto el : pairing) {cells_v.push_back(el.second);} + typedef std::pair< size_t, cell_t > sortable; + std::vector< sortable > pairing; + for (auto tri : cells_v) { + // get index of the triangle + size_t ind; + ind = s_comp->cell_to_index(tri); + pairing.push_back({ind, tri}); + } + std::sort(pairing.begin(), pairing.end(), + [](const sortable& x, const sortable& y) { + return x.first < y.first; + }); + cells_v.clear(); + for (auto el : pairing) { + cells_v.push_back(el.second); + } } - // end sorting t0 = clock(); - std::shared_ptr p_snap( - new gsimp::path_snapper(s_comp)); + std::shared_ptr< gsimp::path_snapper > p_snap = + std::make_shared< gsimp::path_snapper >(s_comp); t1 = clock(); std::cout << "created snapper in " << t1 - t0 << " clock cycles " << float(t1 - t0) / CLOCKS_PER_SEC << " seconds\n"; t0 = clock(); - std::vector snapped; + std::vector< size_t > snapped; if (cycle_indices.size() == 0) snapped = p_snap->snap_path_to_indices(cycle_points); else { @@ -227,17 +299,19 @@ void call_single_cycle_test(std::string complex_file, gsimp::bounding_chain ch_calc(s_comp); t1 = clock(); std::cout << "calculated boundary matrices in " << t1 - t0 - << " clock cycles\ - " << float(t1 - t0) / CLOCKS_PER_SEC << " seconds\n"; + << " clock cycles\ + " << float(t1 - t0) / CLOCKS_PER_SEC + << " seconds\n"; gsimp::chain_t cycle2 = p_snap->index_sequence_to_chain(snapped); t0 = clock(); gsimp::chain_t b_chain_0 = ch_calc.get_bounding_chain(cycle2); t1 = clock(); std::cout << "calculated bounding chain (using Eigen) in " << t1 - t0 << "\ - clock cycles " << float(t1 - t0) / CLOCKS_PER_SEC << " seconds\n"; + clock cycles " + << float(t1 - t0) / CLOCKS_PER_SEC << " seconds\n"; - std::vector> colors; + std::vector< std::vector< int > > colors; for (size_t i = 0; i < gsimp::chain_size(b_chain_0); i++) { if (std::abs(gsimp::chain_val(b_chain_0, i)) > 10e-3) @@ -249,8 +323,8 @@ void call_single_cycle_test(std::string complex_file, std::ofstream my_ply; my_ply.open("my_ply.ply"); - std::vector> edges{}; - std::vector> edge_colors{}; + std::vector< std::vector< size_t > > edges{}; + std::vector< std::vector< int > > edge_colors{}; for (size_t i = 0; i < gsimp::chain_size(cycle); i++) { if (std::abs(gsimp::chain_val(cycle, i)) > 10e-3) @@ -258,16 +332,16 @@ void call_single_cycle_test(std::string complex_file, edge_colors.push_back({0, 0, 255}); } - make_ply(my_ply, s_comp->get_points(), cells_v, colors, edges, - edge_colors); + make_ply(my_ply, s_comp->get_points(), cells_v, colors, edges, edge_colors); my_ply.close(); + gsimp::cell_t null_cell; if (!in_plane) { null_cell = s_comp->index_to_cell(s_comp->dimension(), null_face); std::cout << "null cell: "; std::copy(null_cell.begin(), null_cell.end(), - std::ostream_iterator(std::cout, " ")); + std::ostream_iterator< size_t >(std::cout, " ")); std::cout << "\n"; } diff --git a/tinyobjloader b/tinyobjloader new file mode 160000 index 0000000..d541711 --- /dev/null +++ b/tinyobjloader @@ -0,0 +1 @@ +Subproject commit d541711a794343de4ef5ea76f037c9fb9c127a55 diff --git a/tinyply b/tinyply new file mode 160000 index 0000000..5b1796e --- /dev/null +++ b/tinyply @@ -0,0 +1 @@ +Subproject commit 5b1796e1d3ff32aa55a18e2491c9baf0db3b206e