diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/Adaptive_remeshing_sizing_field.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/Adaptive_remeshing_sizing_field.h index ee0e5e4bcdd..f6e7479d22f 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/Adaptive_remeshing_sizing_field.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/Adaptive_remeshing_sizing_field.h @@ -28,6 +28,8 @@ #include #include +#include + #include #include @@ -84,6 +86,8 @@ private: using Distance = typename Neighbor_search::Distance; using Splitter = typename Neighbor_search::Splitter; + using Medial_axis_kd_tree = internal::Medial_axis_kd_tree; + using Triangle_vec = std::vector; using Triangle_iter = typename Triangle_vec::iterator; using Triangle_primitive = CGAL::AABB_triangle_primitive; @@ -97,6 +101,8 @@ public: Adaptive_remeshing_sizing_field(const Tr& tr) : m_gt(tr.geom_traits()) , m_kd_tree(points_with_info(tr), Splitter(), Kd_traits(Point_property_map())) + , m_k_lipschitz(0.5) + , m_medial_axis_kd_tree(tr) { m_kd_tree.build(); build_aabb_trees(tr); @@ -123,6 +129,9 @@ public: template FT operator()(const Bare_point& p, const int& dim, const Index& i) const { +// if(dim < 3) +// return FT(0);//automatically adapt to 3D sizing around the point + // Find nearest vertex and local size before remeshing Point_property_map pp_map; Distance dist(pp_map); @@ -134,9 +143,6 @@ public: dist); const auto [pi, size, dimension] = search.begin()->first; - if (dim < 3) - return size; - // measure distance to input surfaces Bare_point closest_point = p; FT shortest_distance = (std::numeric_limits::max)(); @@ -153,10 +159,16 @@ public: } } shortest_distance = CGAL::approximate_sqrt(shortest_distance); - FT div_max_distance = 1. / 10.; - //?? (10 as a magic number for distance_field.getMaxDist();) - return size - size / (shortest_distance * div_max_distance); + FT lfs = m_medial_axis_kd_tree.distance_to_medial_axis(p);// closest_point); + FT res = m_k_lipschitz * shortest_distance + lfs; + +// std::cout << "res = " << res +// << "\tshortest_distance = " << shortest_distance +// << "\tlfs = " << lfs << std::endl; + + const FT min_size = 0.5; + return (std::max)(res, min_size); } private: @@ -190,6 +202,9 @@ private: std::vector m_aabb_trees; const GT& m_gt; + const FT m_k_lipschitz; + + Medial_axis_kd_tree m_medial_axis_kd_tree; }; diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/Medial_axis_kd_tree.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/Medial_axis_kd_tree.h new file mode 100644 index 00000000000..b045a0490d4 --- /dev/null +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/Medial_axis_kd_tree.h @@ -0,0 +1,214 @@ +// Copyright (c) 2024 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// +// Author(s) : Jane Tournois +// +//****************************************************************************** +// File Description : Defines a Kd-tree for the medial axis of a 3D triangulation +//****************************************************************************** + +#ifndef CGAL_TETRAHEDRAL_REMESHING_MEDIAL_AXIS_KD_TREE_H +#define CGAL_TETRAHEDRAL_REMESHING_MEDIAL_AXIS_KD_TREE_H + +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include +#include + +namespace CGAL +{ +namespace Tetrahedral_remeshing +{ +namespace internal +{ +/** + * @class Medial_axis_kd_tree + * @tparam Tr a triangulation + * + * @todo add template parameter to take an input aabb_tree + */ +template +class Medial_axis_kd_tree +{ + // Types + typedef typename Tr::Geom_traits GT; + typedef typename Tr::Geom_traits::Point_3 Bare_point; + typedef typename Tr::Point Tr_point; + typedef typename GT::FT FT; + + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Cell_handle Cell_handle; + + struct Point_with_info + { + Bare_point p; + }; + +private: + struct Point_property_map + { + using Self = Point_property_map; + using value_type = Bare_point; + using reference = value_type; //TODO : why can't that be value_type& ? + using key_type = Point_with_info; + using category = boost::readable_property_map_tag; + + const value_type operator[](const key_type& pwi) const { return pwi.p; } + friend const value_type get(const Self&, const key_type& pwi) { return pwi.p; } + }; + +private: + using Kd_traits = CGAL::Search_traits_adapter >; + using Neighbor_search = CGAL::Orthogonal_k_neighbor_search; + using Kd_tree = typename Neighbor_search::Tree; + using Distance = typename Neighbor_search::Distance; + using Splitter = typename Neighbor_search::Splitter; + +public: + /** + * Constructor + */ + Medial_axis_kd_tree(const Tr & tr) + : m_kd_tree(poles(tr), Splitter(), Kd_traits(Point_property_map())) + { + m_kd_tree.build(); + } + +private: + std::vector poles(const Tr& tr) const + { + CGAL::Delaunay_triangulation_3 dt; + + auto tr_cp = tr.geom_traits().construct_point_3_object(); + auto dt_cp = dt.geom_traits().construct_point_3_object(); + + for(auto v : tr.finite_vertex_handles()) + { + if(v->in_dimension() < 3) + dt.insert(dt_cp(v->point())); + } + + for(auto f : tr.finite_facets()) + { + if(f.first->is_facet_on_surface(f.second)) + { + auto t = tr.triangle(f); + dt.insert(dt_cp(CGAL::centroid(tr_cp(t[0]), tr_cp(t[1]), tr_cp(t[2])))); + dt.insert(dt_cp(CGAL::midpoint(tr_cp(t[0]), tr_cp(t[1])))); + dt.insert(dt_cp(CGAL::midpoint(tr_cp(t[0]), tr_cp(t[2])))); + dt.insert(dt_cp(CGAL::midpoint(tr_cp(t[1]), tr_cp(t[2])))); + } + } + + std::ofstream os("dt.xyz"); + for(auto v : dt.finite_vertex_handles()) + os << dt_cp(v->point()) << std::endl; + os.close(); + + std::vector points; + for (auto c : dt.finite_cell_handles()) + { + auto dual_pt_in_dt = dt_cp(dt.dual(c)); + Bare_point p = tr_cp(dual_pt_in_dt); + points.push_back(Point_with_info{p}); + } + + std::ofstream ofs("medial_axis.xyz"); + for (auto p : points) + ofs << p.p << std::endl; + ofs.close(); + + return points; + } + + std::vector skeleton_points(const Tr& tr) const + { + using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3; + using SM = CGAL::Surface_mesh; + using Skeletonization = CGAL::Mean_curvature_flow_skeletonization; + using Skeleton = typename Skeletonization::Skeleton; + using Skeleton_vertex = typename Skeleton::vertex_descriptor; + + C3t3 c3t3; + c3t3.triangulation() = tr; + c3t3.rescan_after_load_of_triangulation(); + + SM sm; + CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, sm); + Skeleton skeleton; + CGAL::extract_mean_curvature_flow_skeleton(sm, skeleton); + + std::vector points; + for (Skeleton_vertex v : CGAL::make_range(vertices(skeleton))) + { + for (auto vd : skeleton[v].vertices) + { + points.push_back(Point_with_info{skeleton[v].point}); +// output << "2 " << skeleton[v].point << " " +// << get(CGAL::vertex_point, tmesh, vd) << "\n"; + } + } + + std::ofstream ofs("skeleton.xyz"); + for (auto p : points) + ofs << p.p << std::endl; + ofs.close(); + + return points; + } + +public: + FT distance_to_medial_axis(const Bare_point& p) const + { + const int nb_nearest_neighbors = 5; + Point_property_map pp_map; + Distance dist(pp_map); + Neighbor_search search(m_kd_tree, + p, //query point + nb_nearest_neighbors, //nb nearest neighbors + 0, //epsilon + true, //search nearest + dist); + + typename GT::Vector_3 sump = CGAL::NULL_VECTOR; + for (auto it = search.begin(); it != search.end(); ++it) + sump += typename GT::Vector_3(CGAL::ORIGIN, it->first.p); + sump = sump / nb_nearest_neighbors; + + Bare_point psump = CGAL::ORIGIN + sump; + return CGAL::approximate_sqrt(CGAL::squared_distance(p, psump)); + +// const auto pi = search.begin()->first; +// return CGAL::approximate_sqrt(CGAL::squared_distance(p, pi.p)); + } + +private: + Kd_tree m_kd_tree; +}; +//end class Medial_axis_kd_tree + +} // end namespace internal +} // end namespace Tetrahedral_remeshing +} //namespace CGAL + +#endif // CGAL_TETRAHEDRAL_REMESHING_MEDIAL_AXIS_KD_TREE_H