From 8cfc4f42e4ade6960a01026d2eb91b99e6314975 Mon Sep 17 00:00:00 2001 From: Ume no hana <50943607+umenohana13@users.noreply.github.com> Date: Tue, 2 Dec 2025 09:19:15 +0100 Subject: [PATCH] Duality fixes -> changed Complex_duality_data_t: it now contains two shared_ptr -> dualize_complex functions all return a Complex_duality_data_t -> all examples corrected -> added dual_hdvf_cubical example (for cubical data) --- ...f_simplicial.cpp => dual_hdvf_cubical.cpp} | 60 +++++++--- .../HDVF/dual_hdvf_simplicial_cgal.cpp | 17 ++- .../HDVF/dual_hdvf_simplicial_tetgen.cpp | 13 +- HDVF/examples/HDVF/main_dual_hdvf.cpp | 34 +++--- .../CGAL/HDVF/Geometric_chain_complex_tools.h | 113 +++++++++++------- 5 files changed, 137 insertions(+), 100 deletions(-) rename HDVF/examples/HDVF/{dual_hdvf_simplicial.cpp => dual_hdvf_cubical.cpp} (52%) diff --git a/HDVF/examples/HDVF/dual_hdvf_simplicial.cpp b/HDVF/examples/HDVF/dual_hdvf_cubical.cpp similarity index 52% rename from HDVF/examples/HDVF/dual_hdvf_simplicial.cpp rename to HDVF/examples/HDVF/dual_hdvf_cubical.cpp index 689e3411f46..1b7b4772570 100644 --- a/HDVF/examples/HDVF/dual_hdvf_simplicial.cpp +++ b/HDVF/examples/HDVF/dual_hdvf_cubical.cpp @@ -3,8 +3,8 @@ #include #include #include -#include -#include +#include +#include #include #include #include @@ -20,47 +20,71 @@ typedef HDVF::Hdvf_traits_3 Traits; //typedef HDVF::Zp<5,int,true> Coefficient_ring; typedef HDVF::Z2 Coefficient_ring; -typedef HDVF::Simplicial_chain_complex Complex; +typedef HDVF::Cubical_chain_complex Complex; typedef HDVF::Hdvf_duality HDVF_type; -typedef HDVF::Duality_simplicial_complex_tools Tools_type; +typedef HDVF::Duality_cubical_complex_tools Tools_type; typedef HDVF::Sub_chain_complex_mask Sub_chain_complex; +// PRIMAL / DUAL construction + +#define PRIMAL_DUAL true // or false + +// Frame or not the complex (add 1 pixel around) + +#define FRAME true // or false + int main(int argc, char **argv) { std::string filename ; - if (argc > 2) std::cout << "usage: dual_hdvf_simplicial off_file" << std::endl; - else if (argc == 1) filename = "data/mesh_data/two_rings.off"; + if (argc > 2) std::cout << "usage: dual_hdvf_cubical pgm_file" << std::endl; + else if (argc == 1) filename = "data/data_cubical/Eight_10.pgm"; else filename = argv[1]; - // Load cub object - HDVF::Mesh_object_io mesh ; - mesh.read_off(filename); + // Set variables to chose loading mode (PRIMAL or DUAL associated complex, pgm read using Khalimsky coordinates or not) + Complex::Cubical_complex_primal_dual primal_dual; + bool khalimsky; + if (PRIMAL_DUAL) { + primal_dual = Complex::PRIMAL; + khalimsky = true; + } + else { + primal_dual = Complex::DUAL; + khalimsky = false; + } - mesh.print_infos(); + // Load pgm object + HDVF::Cub_object_io cub_object ; + cub_object.read_pgm(filename, khalimsky); + // Frame the object + if (FRAME) + cub_object.frame(); + + cub_object.print_infos(); + // Build K and L - typename Tools_type::Complex_duality_data t(Tools_type::dualize_complex(mesh, 1.5, "tmp/file_K_closed.off", 1)) ; + typename Tools_type::Complex_duality_data t(Tools_type::dualize_complex(cub_object, primal_dual)) ; std::cout << "--- Triangulation built" << std::endl ; - Complex& L(t.L) ; - Sub_chain_complex& K(t.K) ; + std::shared_ptr L(t.L_complex) ; + std::shared_ptr K(t.K_complex) ; std::cout << "--- K,L built" << std::endl ; std::cout << "----> complex informations" << std::endl ; std::cout << "------> complex L" << std::endl ; - std::cout << L; + std::cout << *L; std::cout << "------> subcomplex K" << std::endl ; - std::cout << K << std::endl ; + std::cout << *K << std::endl ; // Create and compute a perfect HDVF - HDVF_type hdvf(L, K, HDVF::OPT_FULL); + HDVF_type hdvf(*L, *K, HDVF::OPT_FULL); hdvf.compute_perfect_hdvf(); // Export K HDVF hdvf.set_mask_K(); - CGAL::IO::write_VTK(hdvf, L, "tmp/res_complex_K", false) ; + CGAL::IO::write_VTK(hdvf, *L, "tmp/res_complex_K", false) ; // Export L-K HDVF hdvf.set_mask_L_K(); - CGAL::IO::write_VTK(hdvf, L, "tmp/res_cocomplex_L_K", false) ; + CGAL::IO::write_VTK(hdvf, *L, "tmp/res_cocomplex_L_K", false) ; // Compute pairing std::vector pairs = hdvf.compute_pairing_hdvf(); // Output pairing diff --git a/HDVF/examples/HDVF/dual_hdvf_simplicial_cgal.cpp b/HDVF/examples/HDVF/dual_hdvf_simplicial_cgal.cpp index 7e9c22e15ad..bdfec5865e4 100644 --- a/HDVF/examples/HDVF/dual_hdvf_simplicial_cgal.cpp +++ b/HDVF/examples/HDVF/dual_hdvf_simplicial_cgal.cpp @@ -47,29 +47,28 @@ int main(int argc, char **argv) // std::cout << sm; // Build K and L - typename Tools_type::Complex_duality_data t; - Tools_type::dualize_complex(sm, t, 1.7, 1 ); + typename Tools_type::Complex_duality_data t(Tools_type::dualize_complex(sm, 1.7, 1)); std::cout << "--- Triangulation built" << std::endl ; - Complex& L(*t.L_complex); - Sub_chain_complex& K(*t.K_complex); + std::shared_ptr L(t.L_complex); + std::shared_ptr K(t.K_complex); std::cout << "--- K,L built" << std::endl ; std::cout << "----> complex informations" << std::endl ; std::cout << "------> complex L" << std::endl ; - std::cout << L; + std::cout << *L; std::cout << "------> subcomplex K" << std::endl ; - std::cout << K << std::endl ; + std::cout << *K << std::endl ; // Create and compute a perfect HDVF - HDVF_type hdvf(L, K, HDVF::OPT_FULL); + HDVF_type hdvf(*L, *K, HDVF::OPT_FULL); hdvf.compute_perfect_hdvf(); // Export K HDVF hdvf.set_mask_K(); - CGAL::IO::write_VTK(hdvf, L, "tmp/res_complex_K", false) ; + CGAL::IO::write_VTK(hdvf, *L, "tmp/res_complex_K", false) ; // Export L-K HDVF hdvf.set_mask_L_K(); - CGAL::IO::write_VTK(hdvf, L, "tmp/res_cocomplex_L_K", false) ; + CGAL::IO::write_VTK(hdvf, *L, "tmp/res_cocomplex_L_K", false) ; // Compute pairing std::vector pairs = hdvf.compute_pairing_hdvf(); // Output pairing diff --git a/HDVF/examples/HDVF/dual_hdvf_simplicial_tetgen.cpp b/HDVF/examples/HDVF/dual_hdvf_simplicial_tetgen.cpp index 6537166a0af..f87be56b667 100644 --- a/HDVF/examples/HDVF/dual_hdvf_simplicial_tetgen.cpp +++ b/HDVF/examples/HDVF/dual_hdvf_simplicial_tetgen.cpp @@ -39,11 +39,10 @@ int main(int argc, char **argv) mesh.print_infos(); // Build K and L - typename Tools_type::Complex_duality_data t; - Tools_type::dualize_complex(mesh, t, 1.5, "tmp/file_K_closed.off", 1) ; + typename Tools_type::Complex_duality_data t = Tools_type::dualize_complex(mesh, 1.5, "tmp/file_K_closed.off", 1) ; std::cout << "--- Triangulation built" << std::endl ; - Complex& L(*t.L_complex) ; - Sub_chain_complex& K(*t.K_complex) ; + std::shared_ptr L(t.L_complex) ; + std::shared_ptr K(t.K_complex) ; std::cout << "--- K,L built" << std::endl ; std::cout << "----> complex informations" << std::endl ; @@ -53,15 +52,15 @@ int main(int argc, char **argv) std::cout << K << std::endl ; // Create and compute a perfect HDVF - HDVF_type hdvf(L, K, HDVF::OPT_FULL); + HDVF_type hdvf(*L, *K, HDVF::OPT_FULL); hdvf.compute_perfect_hdvf(); // Export K HDVF hdvf.set_mask_K(); - CGAL::IO::write_VTK(hdvf, L, "tmp/res_complex_K", false) ; + CGAL::IO::write_VTK(hdvf, *L, "tmp/res_complex_K", false) ; // Export L-K HDVF hdvf.set_mask_L_K(); - CGAL::IO::write_VTK(hdvf, L, "tmp/res_cocomplex_L_K", false) ; + CGAL::IO::write_VTK(hdvf, *L, "tmp/res_cocomplex_L_K", false) ; // Compute pairing std::vector pairs = hdvf.compute_pairing_hdvf(); // Output pairing diff --git a/HDVF/examples/HDVF/main_dual_hdvf.cpp b/HDVF/examples/HDVF/main_dual_hdvf.cpp index 856101b896b..ecbb362b32f 100644 --- a/HDVF/examples/HDVF/main_dual_hdvf.cpp +++ b/HDVF/examples/HDVF/main_dual_hdvf.cpp @@ -177,17 +177,16 @@ void main_code (const Options &options) // Build L (bounding sphere meshed with tetgen), K and L-K - typename ToolsType::Complex_duality_data t; - ToolsType::dualize_complex(mesh,t) ; - Complex& L(*t.L_complex) ; - SubCCType& K(*t.K_complex) ; + typename ToolsType::Complex_duality_data t = ToolsType::dualize_complex(mesh) ; + std::shared_ptr L(t.L_complex) ; + std::shared_ptr K(t.K_complex) ; // Output/export mesh and complex - mesh_complex_output, Complex>(mesh, L, K, options) ; + mesh_complex_output, Complex>(mesh, *L, *K, options) ; // HDVF computation, export, output - HDVF_type& hdvf(dual_HDVF_comput(L, K, options)) ; + HDVF_type& hdvf(dual_HDVF_comput(*L, *K, options)) ; // Export to vtk if (options.with_vtk_export) @@ -196,12 +195,12 @@ void main_code (const Options &options) // K { hdvf.set_mask_K() ; - CGAL::IO::write_VTK(hdvf, L, options.outfile_root+"_complex_K", options.co_faces) ; + CGAL::IO::write_VTK(hdvf, *L, options.outfile_root+"_complex_K", options.co_faces) ; } // L-K { hdvf.set_mask_L_K() ; - CGAL::IO::write_VTK(hdvf, L, options.outfile_root+"_cocomplex_L_K", options.co_faces) ; + CGAL::IO::write_VTK(hdvf, *L, options.outfile_root+"_cocomplex_L_K", options.co_faces) ; } } @@ -240,23 +239,18 @@ void main_code (const Options &options) mesh.frame() ; } - // Complex - Complex* complex = new Complex(mesh, primal_dual); - // Build L, K and L-K - typename ToolsType::Complex_duality_data p; - ToolsType::dualize_complex(*complex,p) ; - delete complex ; - Complex &L(*p.L_complex) ; - SubCCType &K(*p.K_complex) ; + typename ToolsType::Complex_duality_data p = ToolsType::dualize_complex(mesh, primal_dual) ; + std::shared_ptr L(p.L_complex) ; + std::shared_ptr K(p.K_complex) ; // Output/export mesh and complex - mesh_complex_output, Complex>(mesh, L, K, options) ; + mesh_complex_output, Complex>(mesh, *L, *K, options) ; // HDVF computation, export, output - HDVF_type& hdvf(dual_HDVF_comput(L, K, options)) ; + HDVF_type& hdvf(dual_HDVF_comput(*L, *K, options)) ; // Export to vtk if (options.with_vtk_export) @@ -265,12 +259,12 @@ void main_code (const Options &options) // K { hdvf.set_mask_K() ; - CGAL::IO::write_VTK(hdvf, L, options.outfile_root+"_complex_K", options.co_faces) ; + CGAL::IO::write_VTK(hdvf, *L, options.outfile_root+"_complex_K", options.co_faces) ; } // L-K { hdvf.set_mask_L_K() ; - CGAL::IO::write_VTK(hdvf, L, options.outfile_root+"_cocomplex_L_K", options.co_faces) ; + CGAL::IO::write_VTK(hdvf, *L, options.outfile_root+"_cocomplex_L_K", options.co_faces) ; } } diff --git a/HDVF/include/CGAL/HDVF/Geometric_chain_complex_tools.h b/HDVF/include/CGAL/HDVF/Geometric_chain_complex_tools.h index aa9e77c136b..46a895bfeba 100644 --- a/HDVF/include/CGAL/HDVF/Geometric_chain_complex_tools.h +++ b/HDVF/include/CGAL/HDVF/Geometric_chain_complex_tools.h @@ -313,27 +313,19 @@ public: /** \brief Type of sub chain complex mask encoding the sub complex `*K_complex`. */ typedef Sub_chain_complex_mask Sub_chain_complex ; /** \brief Pointer over a complex homeomorphic to \f$\mathbb B^3\f$. */ - Chain_complex* L_complex ; + std::shared_ptr L_complex ; /** \brief Pointer over a sub chain complex mask encoding a sub complex of `*L_complex`. */ - Sub_chain_complex* K_complex ; + std::shared_ptr K_complex ; /** \brief Default constructor*/ - Complex_duality_data_t() : L_complex(NULL), K_complex(NULL) {} + Complex_duality_data_t() : L_complex(nullptr), K_complex(nullptr) {} /** \brief Constructor from a pointer over a chain complex and a sub chain complex mask. * * \warning The `Complex_duality_data_t` destructor deletes its underlying complex and sub chain complex mask. */ - Complex_duality_data_t(Chain_complex* L, Sub_chain_complex* K) : L_complex(L), K_complex(K) {} - Complex_duality_data_t(const Complex_duality_data_t&) = delete; - /** Destructor of the class. - * - * The `Complex_duality_data_t` destructor deletes its underlying complex and sub chain complex mask. - * */ - ~Complex_duality_data_t() { - delete L_complex; - delete K_complex; - } + Complex_duality_data_t(std::shared_ptr L, std::shared_ptr K) : L_complex(L), K_complex(K) {} +// Complex_duality_data_t(const Complex_duality_data_t&) = delete; } ; @@ -370,9 +362,9 @@ public: Duality_simplicial_complex_tools() {} private: - static Sub_chain_complex compute_sub_chain_complex(const Chain_complex& _K, const Chain_complex &L){ + static std::shared_ptr compute_sub_chain_complex(const Chain_complex& _K, const Chain_complex &L){ // Build the Sub_chain_complex_mask encoding K_init inside L - Sub_chain_complex K(Sub_chain_complex(L, false)) ; + std::shared_ptr K(std::make_shared(L, false)); // Visit all cells of K_init and activate the corresponding bit in K for (int q=0; q<=_K.dimension(); ++q) { @@ -380,7 +372,7 @@ private: { const Simplex& simplex(_K.index_to_cell(i,q)) ; const size_t id(L.cell_to_index(simplex)); - K.set_bit_on(q, id) ; + K->set_bit_on(q, id) ; } } return K; @@ -399,12 +391,11 @@ public: * * * \param mesh_io `Mesh_object_io` containing the initial set of simplicies (working mesh). - * @param dualized_complex Output structure containing the dualized complex data. * \param BB_ratio Ratio of the "closing" icosphere diameter with respect to the diameter of the object's bounding box. * \param out_file_prefix Prefix of tetgen intermediate files (default: "file_K_closed.off"). * \param subdiv Subdivision level of the bounding icosphere. */ - static void dualize_complex (Mesh_object_io& mesh_io, Complex_duality_data& dualized_complex, double BB_ratio=1.5, const std::string& out_file_prefix = "file_K_closed.off", unsigned int subdiv = 2) + static Complex_duality_data dualize_complex (Mesh_object_io& mesh_io, double BB_ratio=1.5, const std::string& out_file_prefix = "file_K_closed.off", unsigned int subdiv = 2) { std::cerr << "-- Starting dualize_complex" << std::endl; @@ -439,15 +430,18 @@ public: Tet_object_io tetL(out_file_prefix) ; // Build the output data structure + Complex_duality_data dualized_complex; // Build the associated SimpComplex - dualized_complex.L_complex = new Chain_complex(tetL) ; - Chain_complex& L(*dualized_complex.L_complex); - std::cout << "------ L:" << L; + dualized_complex.L_complex = std::make_shared(tetL) ; + std::shared_ptr L(dualized_complex.L_complex); + std::cout << "------ L:" << *L; // Build the Sub_chain_complex_mask encoding K_init inside L - dualized_complex.K_complex = new Sub_chain_complex(compute_sub_chain_complex(*_K, L)); + dualized_complex.K_complex = compute_sub_chain_complex(*_K, *L); // Remove the temporary complex delete _K; + + return dualized_complex; } /** \brief Generates a subcomplex \f$K\f$K and a complex \f$L\f$ with \f$K\subseteq L\f$ from a `Surface_mesh`. @@ -462,11 +456,10 @@ public: * * * \param mesh `Surface_mesh` containing the initial mesh. - * \param dualized_complex Output structure containing the dualized complex data. * \param BB_ratio Ratio of the "closing" icosphere diameter with respect to the diameter of the object's bounding box. * \param subdiv Subdivision level of the bounding icosphere. */ - static void dualize_complex (Triangle_mesh& mesh, Complex_duality_data& dualized_complex, double BB_ratio=1.5, unsigned int subdiv = 2) + static Complex_duality_data dualize_complex (Triangle_mesh& mesh, double BB_ratio=1.5, unsigned int subdiv = 2) { typedef CGAL::Triangulation_data_structure_3, CGAL::Conforming_constrained_Delaunay_triangulation_cell_base_3> TDS; typedef CGAL::Triangulation_3 Triangulation; @@ -514,39 +507,47 @@ public: // Constrained Delaunay Tetraedrisation preserving plc_facet_map auto ccdt = CGAL::make_conforming_constrained_Delaunay_triangulation_3(mesh, CGAL::parameters::plc_face_id(plc_facet_map)); - -// // Detect refined constrained faces + + // Detect refined constrained faces std::vector > faces_constr(cpt); - for (typename Triangulation::Facet f : ccdt.constrained_facets()) { + for (typename Triangulation::Facet f : ccdt.constrained_facets()) { if (ccdt.is_facet_constrained(f)) { const int i(ccdt.face_constraint_index(f)); faces_constr.at(i).push_back(f); } - } - + } + for (int i=0; i1) + std::cout << i << ": " << faces_constr.at(i).size() << std::endl; + } + Triangulation tri_L = std::move(ccdt).triangulation(); - + + // Build the output object + Complex_duality_data dualized_complex; // Build the associated Triangulation_3_io Triangulation_3_io L_tri_io(tri_L); // Build the associated SimpComplex - dualized_complex.L_complex = new Chain_complex(L_tri_io) ; - Chain_complex& L(*dualized_complex.L_complex); - std::cout << "------ L:" << L; + dualized_complex.L_complex = std::make_shared(L_tri_io) ; + std::shared_ptr L(dualized_complex.L_complex); + std::cout << "------ L:" << *L; - std::vector > indices(L.dimension()+1); - for (int q=0; q<=L.dimension(); ++q){ - for (int i=0; i > indices(L->dimension()+1); + for (int q=0; q<=L->dimension(); ++q){ + for (int i=0; inumber_of_cells(q); ++i){ indices.at(q).push_back(i); } } std::string vtk_file("tmp/test_CDT3.vtk"); - Chain_complex::chain_complex_to_vtk(L, vtk_file, &indices, "int"); + Chain_complex::chain_complex_to_vtk(*L, vtk_file, &indices, "int"); // Build the Sub_chain_complex_mask encoding K_init inside L - dualized_complex.K_complex = new Sub_chain_complex (compute_sub_chain_complex(*_K, L)); + dualized_complex.K_complex = compute_sub_chain_complex(*_K, *L); // Remove the temporary complex delete _K; + + return dualized_complex; } @@ -612,9 +613,8 @@ public: * * * \param K_init Initial cubical chain complex (working mesh). - * \param dualized_complex Output structure containing the dualized complex data. */ - static void dualize_complex (const Chain_complex& K_init, Complex_duality_data& dualized_complex) + static Complex_duality_data dualize_complex (const Chain_complex& K_init) { Cub_object_io tmp_L ; tmp_L.dim = K_init.dimension() ; @@ -628,22 +628,43 @@ public: (tmp_L.ncubs)[dtmp] += 1 ; tmp_L.cubs.push_back(tmpkhal) ; } - dualized_complex.L_complex = new Chain_complex(tmp_L, Chain_complex::PRIMAL); - Chain_complex& L(*dualized_complex.L_complex) ; + + // Build the output object + Complex_duality_data dualized_complex; + dualized_complex.L_complex = std::make_shared(tmp_L, Chain_complex::PRIMAL); + std::shared_ptr L(dualized_complex.L_complex) ; // Build the Sub_chain_complex_mask corresponding to _CC - dualized_complex.K_complex = new Sub_chain_complex (Sub_chain_complex(L, false)) ; - Sub_chain_complex& K(*dualized_complex.K_complex) ; + dualized_complex.K_complex = std::make_shared (Sub_chain_complex(*L, false)) ; + std::shared_ptr K(dualized_complex.K_complex) ; // Visit all cells of _CC and activate the corresponding bit in K for (int q=0; q<=K_init.dimension(); ++q) { for (size_t i=0; i khal(K_init.index_to_cell(i, q)) ; - const size_t j = L.cell_to_index(khal); - K.set_bit_on(q,j) ; + const size_t j = L->cell_to_index(khal); + K->set_bit_on(q,j) ; } } + return dualized_complex; + } + + /** \brief Generates a subcomplex \f$K\f$K and a complex \f$L\f$ with \f$K\subseteq L\f$ from a `Cub_object_io` `K_init`. + * + * `L` is the bounding box of `K_init` (homeomorphic to a ball) and \f$K\f$ is a sub chain complex mask encoding `K_init`. + * + * The following figures shows the resulting complex with an initial simple cubical complex (right - sectional view): + * + * + * + * \param K_init Initial `Cub_object_io` (cubical object). + */ + static Complex_duality_data dualize_complex (const Cub_object_io& K_init, typename Chain_complex::Cubical_complex_primal_dual primal_dual) { + Chain_complex* _K = new Chain_complex(K_init, primal_dual); + Complex_duality_data res(dualize_complex(*_K)); + delete _K; + return res; } } ;