1 // Copyright (c) 2016 GeometryFactory (France). 2 // All rights reserved. 3 // 4 // This file is part of CGAL (www.cgal.org). 5 // 6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Mesh_3/include/CGAL/Mesh_3/experimental/Lipschitz_sizing_experimental.h $ 7 // $Id: Lipschitz_sizing_experimental.h 1faa0e2 2021-04-28T10:55:26+02:00 Sébastien Loriot 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // 11 // Author(s) : Jane Tournois 12 // 13 14 #ifndef CGAL_LIPSCHITZ_SIZING_H 15 #define CGAL_LIPSCHITZ_SIZING_H 16 17 #include <CGAL/license/Mesh_3.h> 18 19 #include <CGAL/AABB_tree.h> 20 #include <CGAL/AABB_traits.h> 21 #include <CGAL/AABB_triangle_primitive.h> 22 23 #include <CGAL/Mesh_3/experimental/AABB_filtered_projection_traits.h> 24 #include <CGAL/Mesh_3/experimental/Facet_patch_id_map.h> 25 26 #include <CGAL/Mesh_3/experimental/Lipschitz_sizing_parameters.h> 27 28 #include <CGAL/Default.h> 29 #include <CGAL/array.h> 30 #include <CGAL/Bbox_3.h> 31 32 #include <memory> 33 34 #include <list> 35 #include <sstream> 36 #include <string> 37 #include <fstream> 38 #include <vector> 39 40 namespace CGAL 41 { 42 namespace Mesh_3 43 { 44 45 template <class Kernel, class MeshDomain 46 , typename AABBTreeTemplate = CGAL::Default 47 #ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 48 , typename Facet_patch_id_map_ = CGAL::Default 49 , typename Patches_ids_ = CGAL::Default 50 #endif 51 > 52 class Lipschitz_sizing 53 { 54 public: 55 typedef Kernel K; 56 typedef typename Kernel::FT FT; 57 typedef typename Kernel::Triangle_3 Triangle; 58 typedef typename Kernel::Point_3 Point_3; 59 60 typedef typename std::list<Triangle>::iterator Tr_iterator; 61 typedef CGAL::AABB_triangle_primitive<K, Tr_iterator> Primitive; 62 typedef CGAL::AABB_traits<K, Primitive> AABB_tr_traits; 63 typedef CGAL::AABB_tree<AABB_tr_traits> AABB_tree; 64 65 typedef typename CGAL::Default::Get<AABBTreeTemplate, AABB_tree>::type Tree; 66 67 typedef typename MeshDomain::Index Index; 68 typedef typename MeshDomain::Corner_index Corner_index; 69 typedef typename MeshDomain::Subdomain_index Subdomain_index; 70 typedef typename MeshDomain::Surface_patch_index Surface_patch_index; 71 72 typedef CGAL::Lipschitz_sizing_parameters<MeshDomain, FT> Parameters; 73 74 #ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 75 typedef typename CGAL::Default::Get<Patches_ids_, 76 typename MeshDomain::Surface_patch_index_set>::type Patches_ids; 77 typedef std::vector<Patches_ids> Patches_ids_map; 78 79 typedef typename CGAL::Default::Get< 80 Facet_patch_id_map_, 81 CGAL::Mesh_3::Facet_patch_id_map<MeshDomain, typename Tree::Primitive> 82 >::type Facet_patch_id_map; 83 84 typedef CGAL::Mesh_3::Filtered_projection_traits< 85 typename Tree::AABB_traits, Facet_patch_id_map> AABB_filtered_traits; 86 87 private: 88 typedef CGAL::Search_traits_3<Kernel> KdTreeTraits; 89 typedef CGAL::Orthogonal_k_neighbor_search<KdTreeTraits> Neighbor_search; 90 typedef typename Neighbor_search::Tree Kd_tree; 91 #endif 92 93 private: 94 //only one of these aabb_trees is needed 95 const Tree* m_ptree; 96 std::shared_ptr<Tree> m_own_ptree; 97 98 const MeshDomain& m_domain; 99 Parameters m_params; 100 101 const std::array<double, 3>& m_vxyz; 102 const CGAL::Bbox_3& m_bbox; 103 const bool m_domain_is_a_box; 104 105 #ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 106 //help to accelerate aabb_tree queries in m_ptree 107 std::shared_ptr<Kd_tree> m_kd_tree; 108 109 Facet_patch_id_map m_facet_patch_id_map; 110 const Patches_ids_map& patches_ids_map; 111 #endif 112 113 public: Lipschitz_sizing(const MeshDomain & domain)114 Lipschitz_sizing(const MeshDomain& domain) 115 : m_ptree(nullptr) 116 , m_own_ptree() 117 , m_domain(domain) 118 , m_params(domain) 119 { 120 } 121 Lipschitz_sizing(const MeshDomain & domain,const Tree * ptree,const std::array<double,3> & vxyz,const CGAL::Bbox_3 & bbox,const bool domain_is_a_box,const Patches_ids_map & patches_ids_map)122 Lipschitz_sizing(const MeshDomain& domain 123 , const Tree* ptree 124 , const std::array<double, 3>& vxyz 125 , const CGAL::Bbox_3& bbox 126 , const bool domain_is_a_box 127 #ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 128 , const Patches_ids_map& patches_ids_map 129 #endif 130 ) 131 : m_ptree(ptree) 132 , m_own_ptree() 133 , m_domain(domain) 134 , m_params(domain) 135 , m_vxyz(vxyz) 136 , m_bbox(bbox) 137 , m_domain_is_a_box(domain_is_a_box) 138 #ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 139 , m_facet_patch_id_map() 140 , patches_ids_map(patches_ids_map) 141 #endif 142 { 143 } 144 operator()145 FT operator()(const Point_3& p, const int dim, const Index& index) const 146 { 147 CGAL_assertion(!m_params.empty()); 148 #ifdef CGAL_MESH_3_LIPSCHITZ_SIZING_VERBOSE 149 std::cout << "D = " << dim << "\t"; 150 #endif 151 if (dim == 3) 152 { 153 return size_in_subdomain(p, m_domain.subdomain_index(index)); 154 } 155 else if (dim == 2) 156 { 157 #ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 158 Surface_patch_index sp_index = m_domain.surface_patch_index(index); 159 160 if(!is_on_cube_boundary(sp_index) 161 && !is_on_cube_boundary(p)) 162 { 163 #ifdef CGAL_MESH_3_LIPSCHITZ_SIZING_VERBOSE 164 std::cout << " \n"; 165 #endif 166 FT size_max; 167 m_params.get_parameters(sp_index, size_max); 168 return size_max; 169 } 170 else 171 { 172 #ifdef CGAL_MESH_3_LIPSCHITZ_SIZING_VERBOSE 173 std::cout << "(on cube) "; 174 #endif 175 const std::pair<Subdomain_index, Subdomain_index>& index 176 = m_params.incident_subdomains(sp_index); 177 178 #ifdef CGAL_MESH_3_IMAGE_EXAMPLE 179 if (index.first == INT_MIN) 180 return size_in_subdomain(p, index.second); 181 else 182 return size_in_subdomain(p, index.first); 183 #else //==POLYHEDRAL EXAMPLE 184 if (!is_in_domain(index.first)) 185 return size_in_subdomain(p, index.second); 186 else 187 return size_in_subdomain(p, index.first); 188 #endif //CGAL_MESH_3_IMAGE_EXAMPLE 189 } 190 191 #else //CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 192 CGAL_assertion(false); 193 return 0.; 194 #endif //CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 195 } 196 197 #ifndef CGAL_MESH_3_IMAGE_EXAMPLE 198 else if (dim == 1) 199 { 200 #ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 201 const typename MeshDomain::Curve_index& curve_id = 202 m_domain.curve_index(index); 203 const Patches_ids& ids = patches_ids_map[curve_id]; 204 205 if (m_domain_is_a_box && ids.size() == 2) 206 { 207 //we are on an edge of the box 208 //same code as when dim == 2 209 Surface_patch_index spi = *(ids.begin()); 210 const std::pair<Subdomain_index, Subdomain_index>& subdomains 211 = m_params.incident_subdomains(spi); 212 if (!is_in_domain(subdomains.first)) 213 return size_in_subdomain(p, subdomains.second); 214 else 215 return size_in_subdomain(p, subdomains.first); 216 } 217 return min_size_in_incident_subdomains(ids); 218 #else 219 CGAL_assertion(false);//should not be used for dimension 1 220 return 0.; 221 #endif 222 } 223 else if (dim == 0) 224 { 225 #ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 226 const Corner_index cid = m_domain.corner_index(index); 227 const Patches_ids& ids = m_domain.corners_incidences_map().find(cid)->second; 228 229 if (m_domain_is_a_box && ids.size() == 3) 230 { 231 //we are on a corner of the box 232 //same code as when dim == 2 233 Surface_patch_index spi = *(ids.begin()); 234 const std::pair<Subdomain_index, Subdomain_index>& subdomains 235 = m_params.incident_subdomains(spi); 236 if (!is_in_domain(subdomains.first)) 237 return size_in_subdomain(p, subdomains.second); 238 else 239 return size_in_subdomain(p, subdomains.first); 240 } 241 242 return min_size_in_incident_subdomains(ids);; 243 #else 244 CGAL_assertion(false);//should not be used for dimension 0 245 return 0; 246 #endif 247 } 248 #endif 249 250 CGAL_assertion(false); 251 return 0.; 252 } 253 254 private: 255 #ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS incident_subdomains(const Patches_ids & ids)256 std::vector<Subdomain_index> incident_subdomains(const Patches_ids& ids) const 257 { 258 std::vector<Subdomain_index> vec; 259 for(Surface_patch_index spi : ids) 260 { 261 const std::pair<Subdomain_index, Subdomain_index>& subdomains 262 = m_params.incident_subdomains(spi); 263 264 if (is_in_domain(subdomains.first)) 265 vec.push_back(subdomains.first); 266 267 if (is_in_domain(subdomains.second)) 268 vec.push_back(subdomains.second); 269 } 270 return vec; 271 } 272 min_size_in_incident_subdomains(const Patches_ids & ids)273 FT min_size_in_incident_subdomains(const Patches_ids& ids) const 274 { 275 FT size = static_cast<FT>((std::numeric_limits<double>::max)()); 276 for(Surface_patch_index spi : ids) 277 { 278 const std::pair<Subdomain_index, Subdomain_index>& subdomains 279 = m_params.incident_subdomains(spi); 280 281 FT k, size_min, size_max; 282 if (is_in_domain(subdomains.first)) 283 { 284 m_params.get_parameters(subdomains.first, k, size_min, size_max); 285 size = (std::min)(size, size_min); 286 } 287 if (is_in_domain(subdomains.second)) 288 { 289 m_params.get_parameters(subdomains.second, k, size_min, size_max); 290 size = (std::min)(size, size_min); 291 } 292 } 293 return size; 294 } 295 #endif 296 297 public: 298 template <typename C3T3> init_aabb_tree_from_c3t3(const C3T3 * p_c3t3)299 void init_aabb_tree_from_c3t3(const C3T3* p_c3t3) 300 { 301 static std::list<Triangle> triangles; 302 for (typename C3T3::Facets_in_complex_iterator 303 fit = p_c3t3->facets_in_complex_begin(); 304 fit != p_c3t3->facets_in_complex_end(); 305 ++fit) 306 { 307 if (!is_on_cube_boundary(*fit)) 308 triangles.push_back(p_c3t3->triangulation().triangle(*fit)); 309 } 310 311 m_own_ptree.reset(new Tree(triangles.begin(), triangles.end())); 312 m_own_ptree->build(); 313 } 314 315 private: 316 template <typename Facet> is_on_cube_boundary(const Facet & f)317 bool is_on_cube_boundary(const Facet& f) const 318 { 319 return is_on_cube_boundary(f.first->surface_patch_index(f.second)); 320 } 321 is_on_cube_boundary(const Surface_patch_index & sp_index)322 bool is_on_cube_boundary(const Surface_patch_index& sp_index) const 323 { 324 const std::pair<Subdomain_index, Subdomain_index>& index 325 = m_params.incident_subdomains(sp_index); 326 327 #ifdef CGAL_MESH_3_IMAGE_EXAMPLE 328 return (index.first == INT_MIN || index.second == INT_MIN); 329 #else //POLYHEDRAL EXAMPLE 330 331 if (m_domain_is_a_box) 332 return !is_in_domain(index.first) || !is_in_domain(index.second); 333 else 334 return false; 335 336 #endif //CGAL_MESH_3_IMAGE_EXAMPLE 337 } 338 is_on_cube_boundary(const Point_3 & p)339 bool is_on_cube_boundary(const Point_3& p) const 340 { 341 #ifndef CGAL_MESH_3_IMAGE_EXAMPLE//POLYHEDRAL EXAMPLE 342 343 if (m_domain_is_a_box) 344 //checks that p is in the outer 'shell' of voxels 345 return p.x() < m_bbox.xmin() + m_vxyz[0] 346 || p.x() > m_bbox.xmax() - m_vxyz[0] 347 || p.y() < m_bbox.ymin() + m_vxyz[1] 348 || p.y() > m_bbox.ymax() - m_vxyz[1] 349 || p.z() < m_bbox.zmin() + m_vxyz[2] 350 || p.z() > m_bbox.zmax() - m_vxyz[2]; 351 else 352 return false; 353 354 #else //IMAGE EXAMPLE 355 CGAL_USE(p); 356 return false; 357 #endif //IMAGE_EXAMPLE 358 } 359 is_in_domain(const Subdomain_index & index)360 bool is_in_domain(const Subdomain_index& index) const 361 { 362 #ifdef CGAL_MESH_3_IMAGE_EXAMPLE 363 return (index != 0 && index != INT_MIN); 364 #else //POLYHEDRAL EXAMPLE 365 return (index != 0); 366 #endif 367 } 368 size_in_subdomain(const Point_3 & p,const Subdomain_index & index)369 FT size_in_subdomain(const Point_3& p, const Subdomain_index& index) const 370 { 371 FT k, size_min, size_max; 372 m_params.get_parameters(index, k, size_min, size_max); 373 374 FT sqdist = 0.; 375 if(m_ptree == nullptr) 376 { 377 sqdist = m_own_ptree->squared_distance(p); 378 } 379 else 380 { 381 Point_3 closest = compute_closest_point(p); 382 sqdist = CGAL::squared_distance(p, closest); 383 } 384 385 FT size = k * CGAL::sqrt(sqdist) + size_min; 386 return (std::min)(size, size_max); 387 } 388 389 #ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS kd_tree()390 void kd_tree() 391 { 392 typedef typename MeshDomain::Polyhedron Polyhedron; 393 if(m_kd_tree.get() == 0) { 394 m_kd_tree.reset(new Kd_tree); 395 for(std::size_t poly_id : m_domain.inside_polyhedra()) { 396 const Polyhedron& poly = m_domain.polyhedra()[poly_id]; 397 for(typename Polyhedron::Vertex_handle v : vertices(poly)) 398 { 399 m_kd_tree->insert(v->point()); 400 } 401 } 402 for(std::size_t poly_id : m_domain.boundary_polyhedra()) { 403 const Polyhedron& poly = m_domain.polyhedra()[poly_id]; 404 for(typename Polyhedron::Vertex_handle v : vertices(poly)) 405 { 406 if(!is_on_cube_boundary(v->point())) 407 m_kd_tree->insert(v->point()); 408 } 409 } 410 m_kd_tree->build(); 411 } 412 } 413 #endif 414 compute_closest_point(const Point_3 & p)415 Point_3 compute_closest_point(const Point_3& p) const 416 { 417 #ifndef CGAL_MESH_3_IMAGE_EXAMPLE //POLYHEDRAL_EXAMPLE 418 419 #ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 420 const std::vector<Surface_patch_index>& boundary_ids = 421 m_domain.boundary_patches(); 422 423 CGAL_STATIC_THREAD_LOCAL_VARIABLE_4(AABB_filtered_traits, 424 projection_traits, 425 boundary_ids.begin(), 426 boundary_ids.end(), 427 m_ptree->traits(), 428 m_facet_patch_id_map); 429 kd_tree();//build it if needed 430 Neighbor_search search(*m_kd_tree, p, 1); 431 projection_traits.reset(search.begin()->first); 432 433 m_ptree->traversal(p, projection_traits); 434 CGAL_assertion(projection_traits.found()); 435 return projection_traits.closest_point(); 436 437 #else //CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 438 return m_ptree->closest_point(p); 439 #endif //CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS 440 441 #else 442 CGAL_assertion(false); 443 return CGAL::ORIGIN;//not used 444 #endif 445 } 446 447 public: add_parameters_for_subdomain(const Subdomain_index & id,const FT & k,const FT & size_min,const FT & size_max)448 void add_parameters_for_subdomain(const Subdomain_index& id 449 , const FT& k 450 , const FT& size_min 451 , const FT& size_max) 452 { 453 m_params.add_subdomain(id, k, size_min, size_max); 454 } 455 456 }; 457 458 }//namespace Mesh_3 459 }//namespace CGAL 460 461 #endif // CGAL_LIPSCHITZ_SIZING_H 462