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