GCC Code Coverage Report

Directory: ./
File: lib/geogram/mesh/mesh_AABB.h
Date: 2026-09-07 02:36:43
Exec Total Coverage
Lines: 96 112 85.7%
Functions: 15 23 65.2%
Branches: 61 140 43.6%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #ifndef GEOGRAM_MESH_MESH_AABB
41 #define GEOGRAM_MESH_MESH_AABB
42
43 /**
44 * \file mesh_AABB.h
45 * \brief Axis Aligned Bounding Box trees for accelerating
46 * geometric queries that operate on a Mesh.
47 */
48
49 #include <geogram/basic/common.h>
50 #include <geogram/mesh/mesh.h>
51 #include <geogram/basic/geometry.h>
52 #include <geogram/basic/stopwatch.h>
53
54 namespace GEO {
55
56 /**
57 * \brief Axis Aligned Bounding Box Tree ordering mode
58 * \details One of
59 * - AABB_INDIRECT: leave mesh untouched, store order in separate vector
60 * - AABB_INPLACE: reorder mesh elements in place
61 * - AABB_NOREORDER: use order of mesh elements (was reordered before)
62 */
63 enum AABBReorderMode {
64 AABB_INPLACE, AABB_INDIRECT, AABB_NOREORDER
65 };
66
67 /**
68 * \brief Base class for Axis Aligned Bounding Box trees.
69 * \tparam BOX the box class (Box2d or Box3d).
70 */
71 template <class BOX> class AABB {
72
73 public:
74 /**
75 * \brief Enlarges all the boxes
76 * \param[in] amount the amount that should be subtracted from the lower
77 * bounds and added to the upper bounds of each box.
78 */
79 void enlarge_boxes(double amount) {
80 for(auto& B: bboxes_) {
81 B.enlarge(amount);
82 }
83 }
84
85 protected:
86
87 /**
88 * \brief Initializes this AABB.
89 * \param[in] nb number of items.
90 * \param[in] get_bbox a function(Box&, index_t) that computes the Box
91 * associated with a given index, in [0..nb-1].
92 */
93 166 void initialize(
94 index_t nb, std::function<void(BOX&, index_t)> get_bbox
95 ) {
96 166 nb_ = nb;
97 166 bboxes_.resize(max_node_index(1, 0, nb) + 1);
98 // +1 because size == max_index + 1 !!!
99
2/4
✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 83 times.
✗ Branch 5 not taken.
 166 init_bboxes_recursive(1, 0, nb_, get_bbox);
100 166 }
101
102 /**
103 * \brief Computes all the elements that have a bbox that
104 * intersects a given bbox in a sub-tree of the AABB tree.
105 *
106 * Note that the tree structure is completely implicit,
107 * therefore the bounds of the (continuous) facet indices
108 * sequences that correspond to the elements contained
109 * in the two nodes are sent as well as the node indices.
110 *
111 * \param[in] action a function that takes as argument
112 * an index_t (cell index) invoked for all cells that
113 * has a bounding box that overlaps \p box.
114 * \param[in] box the query box
115 * \param[in] node index of the first node of the AABB tree
116 * \param[in] b index of the first facet in \p node
117 * \param[in] e one position past the index of the last
118 * facet in \p node
119 */
120 void bbox_intersect_recursive(
121 std::function<void(index_t)> action, const BOX& box,
122 index_t node, index_t b, index_t e
123 ) const {
124 geo_debug_assert(e != b);
125
126 // Prune sub-tree that does not have intersection
127 if(!bboxes_overlap(box, bboxes_[node])) {
128 return;
129 }
130
131 // Leaf case
132 if(e == b+1) {
133 action(element_in_leaf(b));
134 return;
135 }
136
137 // Recursion
138 index_t m = b + (e - b) / 2;
139 index_t node_l = 2 * node;
140 index_t node_r = 2 * node + 1;
141
142 bbox_intersect_recursive(action, box, node_l, b, m);
143 bbox_intersect_recursive(action, box, node_r, m, e);
144 }
145
146 /**
147 * \brief Computes all the pairs of intersecting elements
148 * for two sub-trees of the AABB tree.
149 *
150 * Note that the tree structure is completely implicit,
151 * therefore the bounds of the (continuous) facet indices
152 * sequences that correspond to the elementss contained
153 * in the two nodes are sent as well as the node indices.
154 *
155 * \param[in] action a function taking as arguments two
156 * index_t's, invoked of all pairs of elements that have
157 * overlapping bounding boxes.
158 * \param[in] node1 index of the first node of the AABB tree
159 * \param[in] b1 index of the first facet in \p node1
160 * \param[in] e1 one position past the index of the last
161 * facet in \p node1
162 * \param[in] node2 index of the second node of the AABB tree
163 * \param[in] b2 index of the first facet in \p node2
164 * \param[in] e2 one position past the index of the second
165 * facet in \p node2
166 */
167 5596350 void self_intersect_recursive(
168 std::function<void(index_t,index_t)> action,
169 index_t node1, index_t b1, index_t e1,
170 index_t node2, index_t b2, index_t e2
171 ) const {
172
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 5596350 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 5596350 geo_debug_assert(e1 != b1);
173
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 5596350 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 5596350 geo_debug_assert(e2 != b2);
174
175 // Since we are intersecting the AABBTree with *itself*,
176 // we can prune half of the cases by skipping the test
177 // whenever node2's facet index interval is greater than
178 // node1's facet index interval.
179
2/2
✓ Branch 0 taken 114280 times.
✓ Branch 1 taken 5482070 times.
 5596350 if(e2 <= b1) {
180 114280 return;
181 }
182
183 // The acceleration is here:
184 5482070 if(
185
4/4
✓ Branch 0 taken 5253290 times.
✓ Branch 1 taken 228780 times.
✓ Branch 2 taken 1588786 times.
✓ Branch 3 taken 3893284 times.
 10735360 (node1 != node2) &&
186
2/2
✓ Branch 3 taken 1588786 times.
✓ Branch 4 taken 3664504 times.
 5253290 !bboxes_overlap(bboxes_[node1], bboxes_[node2])
187 ) {
188 1588786 return;
189 }
190
191 // Simple case: leaf - leaf intersection.
192
4/4
✓ Branch 0 taken 1879227 times.
✓ Branch 1 taken 2014057 times.
✓ Branch 2 taken 1095851 times.
✓ Branch 3 taken 783376 times.
 3893284 if(b1 + 1 == e1 && b2 + 1 == e2) {
193
2/2
✓ Branch 0 taken 981351 times.
✓ Branch 1 taken 114500 times.
 1095851 if(b1 != b2) {
194 981351 action(element_in_leaf(b1), element_in_leaf(b2));
195 }
196 1095851 return;
197 }
198
199 // If node2 has more elements than node1, then
200 // intersect node2's two children with node1
201 // else
202 // intersect node1's two children with node2
203
2/2
✓ Branch 0 taken 1602762 times.
✓ Branch 1 taken 1194671 times.
 2797433 if(e2 - b2 > e1 - b1) {
204 1602762 index_t m2 = b2 + (e2 - b2) / 2;
205 1602762 index_t node2_l = 2 * node2;
206 1602762 index_t node2_r = 2 * node2 + 1;
207
2/4
✓ Branch 1 taken 1602762 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1602762 times.
✗ Branch 5 not taken.
 1602762 self_intersect_recursive(action, node1, b1, e1, node2_l, b2, m2);
208
2/4
✓ Branch 1 taken 1602762 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1602762 times.
✗ Branch 5 not taken.
 1602762 self_intersect_recursive(action, node1, b1, e1, node2_r, m2, e2);
209 } else {
210 1194671 index_t m1 = b1 + (e1 - b1) / 2;
211 1194671 index_t node1_l = 2 * node1;
212 1194671 index_t node1_r = 2 * node1 + 1;
213
2/4
✓ Branch 1 taken 1194671 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1194671 times.
✗ Branch 5 not taken.
 1194671 self_intersect_recursive(action, node1_l, b1, m1, node2, b2, e2);
214
2/4
✓ Branch 1 taken 1194671 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1194671 times.
✗ Branch 5 not taken.
 1194671 self_intersect_recursive(action, node1_r, m1, e1, node2, b2, e2);
215 }
216 }
217
218 /**
219 * \brief Computes all the pairs of intersecting elements
220 * for two sub-trees of two AABB trees.
221 *
222 * Note that the tree structure is completely implicit,
223 * therefore the bounds of the (continuous) facet indices
224 * sequences that correspond to the elementss contained
225 * in the two nodes are sent as well as the node indices.
226 *
227 * \param[in] action a function taking as arguments two
228 * index_t's, invoked of all pairs of elements that have
229 * overlapping bounding boxes.
230 * \param[in] node1 index of the first node of the AABB tree
231 * \param[in] b1 index of the first facet in \p node1
232 * \param[in] e1 one position past the index of the last
233 * facet in \p node1
234 * \param[in] other the second AABB tree
235 * \param[in] node2 index of the second node of the second AABB tree
236 * \param[in] b2 index of the first facet in \p node2
237 * \param[in] e2 one position past the index of the second
238 * facet in \p node2
239 */
240 void other_intersect_recursive(
241 std::function<void(index_t,index_t)> action,
242 index_t node1, index_t b1, index_t e1,
243 const AABB<BOX>* other,
244 index_t node2, index_t b2, index_t e2
245 ) const {
246 geo_debug_assert(e1 != b1);
247 geo_debug_assert(e2 != b2);
248
249 // The acceleration is here:
250 if(!bboxes_overlap(bboxes_[node1], other->bboxes_[node2])) {
251 return;
252 }
253
254 // Simple case: leaf - leaf intersection.
255 if(b1 + 1 == e1 && b2 + 1 == e2) {
256 action(element_in_leaf(b1), element_in_leaf(b2));
257 return;
258 }
259
260 // If node2 has more elements than node1, then
261 // intersect node2's two children with node1
262 // else
263 // intersect node1's two children with node2
264 if(e2 - b2 > e1 - b1) {
265 index_t m2 = b2 + (e2 - b2) / 2;
266 index_t node2_l = 2 * node2;
267 index_t node2_r = 2 * node2 + 1;
268 other_intersect_recursive(
269 action, node1, b1, e1, other, node2_l, b2, m2
270 );
271 other_intersect_recursive(
272 action, node1, b1, e1, other, node2_r, m2, e2
273 );
274 } else {
275 index_t m1 = b1 + (e1 - b1) / 2;
276 index_t node1_l = 2 * node1;
277 index_t node1_r = 2 * node1 + 1;
278 other_intersect_recursive(
279 action, node1_l, b1, m1, other, node2, b2, e2
280 );
281 other_intersect_recursive(
282 action, node1_r, m1, e1, other, node2, b2, e2
283 );
284 }
285 }
286
287
288 /**
289 * \brief Computes the maximum node index in a subtree
290 * \param[in] node_index node index of the root of the subtree
291 * \param[in] b first facet index in the subtree
292 * \param[in] e one position past the last facet index in the subtree
293 * \return the maximum node index in the subtree rooted at \p node_index
294 */
295 770770 static index_t max_node_index(index_t node_index, index_t b, index_t e) {
296
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 385385 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 770770 geo_debug_assert(e > b);
297
2/2
✓ Branch 0 taken 192734 times.
✓ Branch 1 taken 192651 times.
 770770 if(b + 1 == e) {
298 385468 return node_index;
299 }
300 385302 index_t m = b + (e - b) / 2;
301 385302 index_t childl = 2 * node_index;
302 385302 index_t childr = 2 * node_index + 1;
303 385302 return std::max(
304 385302 max_node_index(childl, b, m),
305
2/4
✓ Branch 1 taken 192651 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 192651 times.
✗ Branch 5 not taken.
 385302 max_node_index(childr, m, e)
306 385302 );
307 }
308
309 /**
310 * \brief Computes the hierarchy of bounding boxes recursively.
311 * \details This function is generic and can be used to compute
312 * a bbox hierarchy of arbitrary elements.
313 * \param[in] node_index the index of the root of the subtree
314 * \param[in] b first element index in the subtree
315 * \param[in] e one position past the last element index in the subtree
316 * \param[in] get_bbox a function that takes a Box3d& and an index_t,
317 * that computes the bbox of an element.
318 */
319 770770 void init_bboxes_recursive(
320 index_t node_index,
321 index_t b, index_t e,
322 std::function<void(BOX&, index_t)> get_bbox
323 ) {
324
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 385385 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 770770 geo_debug_assert(node_index < bboxes_.size());
325
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 385385 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 770770 geo_debug_assert(b != e);
326
2/2
✓ Branch 0 taken 192734 times.
✓ Branch 1 taken 192651 times.
 770770 if(b + 1 == e) {
327 385468 get_bbox(bboxes_[node_index], element_in_leaf(b));
328 385468 return;
329 }
330 385302 index_t m = b + (e - b) / 2;
331 385302 index_t childl = 2 * node_index;
332 385302 index_t childr = 2 * node_index + 1;
333
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 192651 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 385302 geo_debug_assert(childl < bboxes_.size());
334
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 192651 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 385302 geo_debug_assert(childr < bboxes_.size());
335
2/4
✓ Branch 1 taken 192651 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 192651 times.
✗ Branch 5 not taken.
 385302 init_bboxes_recursive(childl, b, m, get_bbox);
336
2/4
✓ Branch 1 taken 192651 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 192651 times.
✗ Branch 5 not taken.
 385302 init_bboxes_recursive(childr, m, e, get_bbox);
337
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 192651 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 385302 geo_debug_assert(childl < bboxes_.size());
338
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 192651 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 385302 geo_debug_assert(childr < bboxes_.size());
339 385302 bbox_union(bboxes_[node_index], bboxes_[childl], bboxes_[childr]);
340 }
341
342 /**
343 * \brief Tests whether this AABB is indirect or in-place
344 * \details An AABB can be indirect (stores a permutation in a vector)
345 * or in-place (elements are reordered in the mesh for instance).
346 * \retval true if the AABB is indirect
347 * \retval false if the AABB uses in-place reordering
348 */
349 11391530 bool indirect() const {
350 11391530 return (reorder_.size() != 0);
351 }
352
353 /**
354 * \brief Gets the element stored in a leaf node
355 * \details If the AABB is indirect, looks-up the element in the
356 * reorder_ permutation, else returns \p n.
357 * \param[in] i the leaf index, between 0 and nb elements - 1
358 * \return the element stored in the leaf node \p n
359 */
360 11391530 index_t element_in_leaf(index_t i) const {
361
1/2
✓ Branch 1 taken 5695765 times.
✗ Branch 2 not taken.
 11391530 return indirect() ? reorder_[i] : i;
362 }
363
364 protected:
365 index_t nb_;
366 vector<BOX> bboxes_;
367 vector<index_t> reorder_; /**< used if indirect, or unused if in-place */
368 };
369
370 typedef AABB<Box2d> AABB2d;
371 typedef AABB<Box3d> AABB3d;
372
373 /**************************************************************/
374
375 /**
376 * \brief Base class for Axis Aligned Bounding Box trees
377 * of mesh elements with 2d boxes.
378 */
379 class GEOGRAM_API MeshAABB2d : public AABB2d {
380 public:
381 /**
382 * \brief MeshAABB2d constructor.
383 */
384 ✗ MeshAABB2d() : mesh_(nullptr) {
385 ✗ }
386
387 /**
388 * \brief Gets the mesh.
389 * \return a const reference to the mesh.
390 */
391 const Mesh* mesh() const {
392 return mesh_;
393 }
394
395 protected:
396 Mesh* mesh_;
397 };
398
399 /**************************************************************/
400
401 /**
402 * \brief Base class for Axis Aligned Bounding Box trees
403 * of mesh elements with 3d boxes.
404 */
405 class GEOGRAM_API MeshAABB3d : public AABB3d {
406 public:
407 /**
408 * \brief MeshAABB3d constructor.
409 */
410 88 MeshAABB3d() : mesh_(nullptr) {
411 88 }
412
413 /**
414 * \brief Gets the mesh.
415 * \return a const reference to the mesh.
416 */
417 1940 const Mesh* mesh() const {
418 1940 return mesh_;
419 }
420
421 protected:
422 Mesh* mesh_;
423 };
424
425 /**************************************************************/
426
427 /**
428 * \brief Axis Aligned Bounding Box tree of mesh facets in 3D.
429 * \details Used to quickly compute facet intersection and
430 * to locate the nearest facet from 3d query points.
431 */
432 class GEOGRAM_API MeshFacetsAABB : public MeshAABB3d {
433 public:
434
435 /**
436 * \brief Stores all the information related with a ray-facet
437 * intersection.
438 */
439 struct Intersection {
440 179996 Intersection() :
441 179996 t(Numeric::max_float64()),
442 179996 f(NO_INDEX),
443 359992 i(NO_INDEX), j(NO_INDEX), k(NO_INDEX)
444 {
445 179996 }
446 vec3 p; /**< the intersection. */
447 double t; /**< the parameter along the intersected ray. */
448 index_t f; /**< the intersected facet. */
449 vec3 N; /**< the normal vector at the intersection. */
450 index_t i,j,k; /**< the vertices of the intersected triangle. */
451 double u,v; /**< the barycentric coordinates in the triangle. */
452 };
453
454
455 /**
456 * \brief MeshFacetsAABB constructor.
457 * \details Creates an uninitialized MeshFacetsAABB.
458 */
459 6 MeshFacetsAABB() {
460 6 }
461
462 /**
463 * \brief Creates an Axis Aligned Bounding Boxes tree for facets.
464 * \param[in] M the input mesh. It can be modified,
465 * and will be triangulated (if
466 * not already a triangular mesh). The facets are
467 * re-ordered depending on \p reorder_mode
468 * \param[in] reorder_mode one of
469 * - AABB_INDIRECT: leave mesh untouched,
470 * store order in separate vector
471 * - AABB_INPLACE: reorder mesh elements in place
472 * - AABB_NOREORDER: use order of mesh elements
473 * (the mesh was reordered before, using mesh_reorder())
474 */
475 60 MeshFacetsAABB(Mesh& M, AABBReorderMode reorder_mode) {
476
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 initialize(M, reorder_mode);
477 60 }
478
479 /**
480 * \brief Creates an Axis Aligned Bounding Boxes tree for facets.
481 * \details Uses AABB_INDIRECT mode (order stored in separate vector).
482 * \param[in] M a const reference to the input mesh.
483 */
484 22 MeshFacetsAABB(const Mesh& M) {
485
1/2
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
 22 initialize(const_cast<Mesh&>(M), AABB_INDIRECT);
486 22 }
487
488 /**
489 * \brief Initializes the Axis Aligned Bounding Boxes tree.
490 * \param[in] M the input mesh. It can be modified,
491 * and will be triangulated (if
492 * not already a triangular mesh). The facets are
493 * re-ordered depending on \p reorder_mode
494 * \param[in] reorder_mode one of
495 * - AABB_INDIRECT: leave mesh untouched,
496 * store order in separate vector
497 * - AABB_INPLACE: reorder mesh elements in place
498 * - AABB_NOREORDER: use order of mesh elements
499 * (the mesh was reordered before, using mesh_reorder())
500 */
501 void initialize(Mesh& M, AABBReorderMode reorder_mode = AABB_INDIRECT);
502
503 #ifndef GOMGEN
504 [[deprecated("use MeshFacetsAABB(Mesh&,AABBReorderMode) instead")]]
505 MeshFacetsAABB(Mesh& M, bool reorder) {
506 initialize(M, reorder ? AABB_INPLACE : AABB_NOREORDER);
507 }
508
509 [[deprecated("use initialize(Mesh&,AABBReorderMode) instead")]]
510 void initialize(Mesh& M, bool reorder) {
511 initialize(M, reorder ? AABB_INPLACE : AABB_NOREORDER);
512 }
513 #endif
514 /**
515 * \brief Computes all the pairs of intersecting facets.
516 * \param[in] action a function that takes two index_t's
517 * and that is invoked of all pairs of facets that have overlapping
518 * bounding boxes. triangles_intersections() needs to be
519 * called to detect the actual intersections.
520 * \param[in] concurrent if set, then action can be called simultaneously
521 * by concurrent threads, else it is only the determination of
522 * overlapping bboxes that is parallelized, then the list of overlapping
523 * bboxes is internally memorized for serializing the calls to action.
524 */
525 60 void compute_facet_bbox_intersections(
526 std::function<void(index_t, index_t)> action,
527 bool concurrent = false
528 ) const {
529 60 if(
530
5/6
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 46 times.
✓ Branch 4 taken 14 times.
✓ Branch 5 taken 46 times.
✓ Branch 6 taken 14 times.
 120 Process::maximum_concurrent_threads() <= 1 ||
531 60 mesh_->facets.nb() <= 1024
532 ) {
533
2/4
✓ Branch 1 taken 46 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 46 times.
✗ Branch 5 not taken.
 46 self_intersect_recursive(
534 action,
535 46 1, 0, mesh_->facets.nb(),
536 46 1, 0, mesh_->facets.nb()
537 );
538 } else {
539
2/4
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
 14 self_bbox_intersections_parallel(action, concurrent);
540 }
541
542 60 }
543
544 /**
545 * \brief Computes all the intersections between a given
546 * box and the bounding boxes of all the facets.
547 * \param[in] action a function that takes an index_t that is
548 * invoked for all facets that have a bounding
549 * box that intersects \p box_in.
550 */
551 void compute_bbox_facet_bbox_intersections(
552 const Box& box_in, std::function<void(index_t)> action
553 ) const {
554 bbox_intersect_recursive(
555 action, box_in, 1, 0, mesh_->facets.nb()
556 );
557 }
558
559 /**
560 * \brief Finds the nearest facet from an arbitrary 3d query point.
561 * \param[in] p query point
562 * \param[out] nearest_point nearest point on the surface
563 * \param[out] sq_dist squared distance between p and the surface or
564 * or Numeric::max_float64() if mesh has no facet
565 * \return the index of the facet nearest to point p or NO_INDEX if
566 * mesh has no facet.
567 */
568 800714 index_t nearest_facet(
569 const vec3& p, vec3& nearest_point, double& sq_dist
570 ) const {
571
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 800714 times.
 800714 if(mesh_->facets.nb() == 0) {
572 ✗ sq_dist = Numeric::max_float64();
573 ✗ return NO_INDEX;
574 }
575 index_t nearest_facet;
576
1/2
✓ Branch 1 taken 800714 times.
✗ Branch 2 not taken.
 800714 get_nearest_facet_hint(p, nearest_facet, nearest_point, sq_dist);
577
1/2
✓ Branch 1 taken 800714 times.
✗ Branch 2 not taken.
 800714 nearest_facet_recursive(
578 p,
579 nearest_facet, nearest_point, sq_dist,
580 800714 1, 0, mesh_->facets.nb()
581 );
582 800714 return nearest_facet;
583 }
584
585 /**
586 * \brief Finds the nearest facet from an arbitrary 3d query point.
587 * \param[in] p query point
588 * \return the index of the facet nearest to point p or NO_INDEX if
589 * mesh has no facet.
590 */
591 index_t nearest_facet(const vec3& p) const {
592 vec3 nearest_point;
593 double sq_dist;
594 return nearest_facet(p, nearest_point, sq_dist);
595 }
596
597
598 /**
599 * \brief Finds the nearest facet from an arbitrary 3d query point taken
600 * into account only certain facets
601 * \param[in] p query point
602 * \param[out] nearest_point nearest point on the surface
603 * \param[out] sq_dist squared distance between p and the surface.
604 * \param[in] filter a function that takes a facet index and that
605 * returns true if the facet should be taken into account or false
606 * otherwise.
607 * \return the index of the facet nearest to point p.
608 */
609 index_t nearest_facet_filtered(
610 const vec3& p, vec3& nearest_point, double& sq_dist,
611 std::function<bool(index_t)> filter
612 ) const {
613 if(mesh_->facets.nb() == 0) {
614 return NO_INDEX;
615 }
616 index_t nearest_facet = NO_INDEX;
617 sq_dist = Numeric::max_float64();
618 nearest_facet_recursive_filtered(
619 p,
620 nearest_facet, nearest_point, sq_dist,
621 1, 0, mesh_->facets.nb(),
622 filter
623 );
624 return nearest_facet;
625 }
626
627 /**
628 * \brief Computes the nearest point and nearest facet from
629 * a query point, using user-specified hint.
630 *
631 * \details The hint is specified as reasonable initial values of
632 * (nearest_facet, nearest_point, sq_dist). If multiple queries
633 * are done on a set of points that has spatial locality,
634 * the hint can be the result of the previous call.
635 *
636 * \param[in] p query point
637 * \param[in,out] nearest_facet the nearest facet so far,
638 * or NO_FACET if not known yet
639 * \param[in,out] nearest_point a point in nearest_facet
640 * \param[in,out] sq_dist squared distance between p and
641 * nearest_point
642 * \note On entry, \p sq_dist needs to be equal to the squared
643 * distance between \p p and \p nearest_point (it is easy to
644 * forget to update it when calling it within a loop).
645 */
646 ✗ void nearest_facet_with_hint(
647 const vec3& p,
648 index_t& nearest_facet, vec3& nearest_point, double& sq_dist
649 ) const {
650 ✗ if(mesh_->facets.nb() == 0) {
651 ✗ nearest_facet = NO_INDEX;
652 ✗ sq_dist = Numeric::max_float64();
653 ✗ nearest_point=vec3(0,0,0);
654 ✗ return;
655 }
656 ✗ if(nearest_facet == NO_FACET) {
657 ✗ get_nearest_facet_hint(
658 p, nearest_facet, nearest_point, sq_dist
659 );
660 }
661 ✗ nearest_facet_recursive(
662 p,
663 nearest_facet, nearest_point, sq_dist,
664 ✗ 1, 0, mesh_->facets.nb()
665 );
666 }
667
668 /**
669 * \brief Computes the distance between an arbitrary 3d query
670 * point and the surface.
671 * \param[in] p query point
672 * \return the squared distance between \p p and the surface.
673 */
674 800714 double squared_distance(const vec3& p) const {
675 800714 vec3 nearest_point;
676 double result;
677
1/2
✓ Branch 1 taken 800714 times.
✗ Branch 2 not taken.
 800714 nearest_facet(p, nearest_point, result);
678 800714 return result;
679 }
680
681 /**
682 * \brief Tests whether there exists an intersection between a ray
683 * and the mesh.
684 * \param[in] R the ray.
685 * \param[in] tmax optional maximum parameter of the intersection along
686 * the ray.
687 * \param[in] ignore_f optional facet to be ignored in intersection
688 * tests.
689 * \retval true if there was an intersection.
690 * \retval false otherwise.
691 */
692 bool ray_intersection(
693 const Ray& R,
694 double tmax = Numeric::max_float64(),
695 index_t ignore_f = NO_INDEX
696 ) const;
697
698
699 /**
700 * \brief Computes the nearest intersection along a ray.
701 * \param[in] R the ray
702 * \param[in,out] I the intersection. If I.t is set, then intersections
703 * further away than I.t are ignored. If I.f is set, then intersection
704 * with facet f is ignored.
705 * \retval true if there was an intersection.
706 * \retval false otherwise.
707 */
708 bool ray_nearest_intersection(const Ray& R, Intersection& I) const;
709
710 /**
711 * \brief Tests whether this surface mesh has an intersection
712 * with a segment.
713 * \param[in] q1 , q2 the two extremities of the segment.
714 * \retval true if there exists an intersection between [q1 , q2]
715 * and a facet of the mesh.
716 * \retval false otherwise.
717 */
718 ✗ bool segment_intersection(const vec3& q1, const vec3& q2) const {
719 ✗ return ray_intersection(Ray(q1, q2-q1), 1.0);
720 }
721
722 /**
723 * \brief Finds the intersection between a segment and a surface that
724 * is nearest to the first extremity of the segment.
725 * \param[in] q1 , q2 the two extremities of the segment.
726 * \param[out] t if there was an intersection, it is t*q2 + (1-t)*q1
727 * \param[out] f the intersected nearest facet or NO_INDEX if there
728 * was no intersection.
729 * \retval true if there exists at least an intersection
730 * between [q1 , q2] and a facet of the mesh.
731 * \retval false otherwise.
732 */
733 bool segment_nearest_intersection(
734 const vec3& q1, const vec3& q2, double& t, index_t& f
735 ) const {
736 Ray R(q1, q2-q1);
737 Intersection I;
738 I.t = 1.0;
739 bool result = ray_nearest_intersection(R,I);
740 t = I.t;
741 f = I.f;
742 return result;
743 }
744
745 /**
746 * \brief Calls a user function for all ray-facet intersection
747 * \param[in] R the ray
748 * \param[in] action the function to be called
749 */
750 void ray_all_intersections(
751 const Ray& R,
752 std::function<void(const Intersection&)> action
753 ) const;
754
755
756 /**
757 * \brief Calls a user function for all line-facet intersections
758 * \param[in] O origin of the line
759 * \param[in] D direction vector of the line
760 * \param[in] action the function to be called
761 */
762 void line_all_intersections(
763 const vec3& O, const vec3& D,
764 std::function<void(const Intersection&)> action
765 ) const;
766
767
768 /**
769 * \brief Tests whether a closed surface contains a point
770 * \pre The surface from which the MeshFacetsAABB was constructed
771 * is closed
772 * \param[in] p the point to be tested
773 * \retval true if the surface contains \p p
774 * \retval false otherwise
775 */
776 bool contains(const vec3& p) const;
777
778 protected:
779
780 /**
781 * \brief Computes a reasonable initialization for
782 * nearest facet search.
783 *
784 * \details A good initialization makes the algorithm faster,
785 * by allowing early pruning of subtrees that provably
786 * do not contain the nearest neighbor.
787 *
788 * \param[in] p query point
789 * \param[out] nearest_facet a facet reasonably near p
790 * \param[out] nearest_point a point in nearest_facet
791 * \param[out] sq_dist squared distance between p and nearest_point
792 */
793 void get_nearest_facet_hint(
794 const vec3& p,
795 index_t& nearest_facet, vec3& nearest_point, double& sq_dist
796 ) const;
797
798 /**
799 * \brief The recursive function used by the implementation
800 * of nearest_facet().
801 *
802 * \details The first call may use get_nearest_facet_hint()
803 * to initialize nearest_facet, nearest_point and sq_dist,
804 * as done in nearest_facet().
805 *
806 * \param[in] p query point
807 * \param[in,out] nearest_facet the nearest facet so far,
808 * \param[in,out] nearest_point a point in nearest_facet
809 * \param[in,out] sq_dist squared distance between p and nearest_point
810 * \param[in] n index of the current node in the AABB tree
811 * \param[in] b index of the first facet in the subtree under node \p n
812 * \param[in] e one position past the index of the last facet in the
813 * subtree under node \p n
814 */
815 void nearest_facet_recursive(
816 const vec3& p,
817 index_t& nearest_facet, vec3& nearest_point, double& sq_dist,
818 index_t n, index_t b, index_t e
819 ) const;
820
821
822 /**
823 * \brief The recursive function used by the implementation
824 * of nearest_facet_filtered().
825 *
826 * \param[in] p query point
827 * \param[in,out] nearest_facet the nearest facet so far,
828 * \param[in,out] nearest_point a point in nearest_facet
829 * \param[in,out] sq_dist squared distance between p and nearest_point
830 * \param[in] n index of the current node in the AABB tree
831 * \param[in] b index of the first facet in the subtree under node \p n
832 * \param[in] e one position past the index of the last facet in the
833 * subtree under node \p n
834 * \param[in] filter a function that takes a facet index and that
835 * returns true if the facet should be taken into account or false
836 * otherwise.
837 */
838 void nearest_facet_recursive_filtered(
839 const vec3& p,
840 index_t& nearest_facet, vec3& nearest_point, double& sq_dist,
841 index_t n, index_t b, index_t e,
842 std::function<bool(index_t)> filter
843 ) const;
844
845 /**
846 * \brief The recursive function used by the implementation
847 * of ray_intersection()
848 * \param[in] R the ray
849 * \param[in] dirinv
850 * precomputed 1/(q2.x-q1.x), 1/(q2.y-q1.y), 1/(q2.z-q1.z)
851 * \param[in] max_t the maximum value of t for an intersection
852 * \param[in] ignore_f facet index to be ignored in tests
853 * \param[in] n index of the current node in the AABB tree
854 * \param[in] b index of the first facet in the subtree under node \p n
855 * \param[in] e one position past the index of the last facet in the
856 * subtree under node \p n
857 * \retval true if their was an intersection
858 * \retval false otherwise
859 */
860 bool ray_intersection_recursive(
861 const Ray& R, const vec3& dirinv, double max_t, index_t ignore_f,
862 index_t n, index_t b, index_t e
863 ) const;
864
865 /**
866 * \brief The recursive function used by the implementation
867 * of ray_nearest_intersection()
868 * \param[in] R the ray
869 * \param[in] dirinv
870 * precomputed 1/(q2.x-q1.x), 1/(q2.y-q1.y), 1/(q2.z-q1.z)
871 * \param[in,out] I the parameters of the nearest intersection
872 * computed so-far. All intersections further away than I.t are
873 * ignored.
874 * \param[in] ignore_f facet index to be ignored in tests
875 * \param[in] n index of the current node in the AABB tree
876 * \param[in] b index of the first facet in the subtree under node \p n
877 * \param[in] e one position past the index of the last facet in the
878 * subtree under node \p n
879 * \param[in] coord the current splitting coordinate, one of 0,1,2
880 */
881 void ray_nearest_intersection_recursive(
882 const Ray& R, const vec3& dirinv, Intersection& I, index_t ignore_f,
883 index_t n, index_t b, index_t e, index_t coord
884 ) const;
885
886
887 /**
888 * \brief The function used to implement ray_all_intersections()
889 * \param[in] R the ray
890 * \param[in] dirinv
891 * precomputed 1/(q2.x-q1.x), 1/(q2.y-q1.y), 1/(q2.z-q1.z)
892 * \param[in] action the function to be called
893 * \param[in] n index of the current node in the AABB tree
894 * \param[in] b index of the first facet in the subtree under node \p n
895 * \param[in] e one position past the index of the last facet in the
896 * subtree under node \p n
897 */
898 void ray_all_intersections_recursive(
899 const Ray& R, const vec3& dirinv,
900 std::function<void(const Intersection&)> action,
901 index_t n, index_t b, index_t e
902 ) const;
903
904
905 /**
906 * \brief The function used to implement line_all_intersections()
907 * \param[in] O origin of the line
908 * \param[in] D direction vector of the line
909 * \param[in] dirinv
910 * precomputed 1/D.x, 1/D.y, 1/D.z
911 * \param[in] action the function to be called
912 * \param[in] n index of the current node in the AABB tree
913 * \param[in] b index of the first facet in the subtree under node \p n
914 * \param[in] e one position past the index of the last facet in the
915 * subtree under node \p n
916 */
917 void line_all_intersections_recursive(
918 const vec3& O, const vec3& D, const vec3& dirinv,
919 std::function<void(const Intersection&)> action,
920 index_t n, index_t b, index_t e
921 ) const;
922
923
924 /**
925 * \brief Computes all the pairs of intersecting elements
926 * in parallel.
927 * \param[in] action a function taking as arguments two
928 * index_t's, invoked of all pairs of elements that have
929 * overlapping bounding boxes.
930 * \param[in] concurrent if set, then action can be called simultaneously
931 * by concurrent threads, else it is only the determination of
932 * overlapping bboxes that is parallelized, then the list of overlapping
933 * bboxes is internally memorized for serializing the calls to action.
934 */
935 void self_bbox_intersections_parallel(
936 std::function<void(index_t, index_t)> action, bool concurrent=false
937 ) const;
938 };
939
940 /***********************************************************************/
941
942 /**
943 * \brief Axis Aligned Bounding Box tree of mesh cells.
944 * \details Used to quickly find the tetrahedron that contains
945 * a given 3d point.
946 */
947 class GEOGRAM_API MeshCellsAABB : public MeshAABB3d {
948 public:
949
950 /**
951 * \brief Symbolic constant for indicating that there
952 * is no containing tetrahedron.
953 * \see containing_tet()
954 */
955 static constexpr index_t NO_TET = NO_INDEX;
956
957 /**
958 * \brief MeshCellsAABB constructor.
959 * \details Creates an uninitialized MeshCellsAABB.
960 */
961 MeshCellsAABB() {
962 }
963
964 /**
965 * \brief Creates an Axis Aligned Bounding Boxes tree for mesh cells.
966 * \param[in] M the input mesh. It can be modified,
967 * The cells are re-ordered depending on \p reorder_mode
968 * \param[in] reorder_mode one of
969 * - AABB_INDIRECT: leave mesh untouched,
970 * store order in separate vector
971 * - AABB_INPLACE: reorder mesh elements in place
972 * - AABB_NOREORDER: use order of mesh elements
973 * (the mesh was reordered before, using mesh_reorder())
974 */
975 MeshCellsAABB(Mesh& M, AABBReorderMode reorder_mode) {
976 initialize(M, reorder_mode);
977 }
978
979 /**
980 * \brief Creates an Axis Aligned Bounding Boxes tree for mesh cells.
981 * \details Uses AABB_INDIRECT mode (order stored in separate vector).
982 * \param[in] M a const reference to the input mesh.
983 */
984 MeshCellsAABB(const Mesh& M) {
985 initialize(const_cast<Mesh&>(M), AABB_INDIRECT);
986 }
987
988 /**
989 * \brief Initializes the Axis Aligned Bounding Boxes tree.
990 * \param[in] M the input mesh. It can be modified,
991 * The cells are re-ordered depending on \p reorder_mode
992 * \param[in] reorder_mode one of
993 * - AABB_INPLACE: reorder mesh elements in place
994 * - AABB_INDIRECT: leave mesh untouched,
995 * store order in separate vector
996 * - AABB_NOREORDER: use order of mesh elements
997 * (the mesh was reordered before, using mesh_reorder())
998 */
999 void initialize(Mesh& M, AABBReorderMode reorder_mode = AABB_INDIRECT);
1000
1001 #ifndef GOMGEN
1002 [[deprecated("use MeshCellsAABB(Mesh&,AABBReorderMode) instead")]]
1003 MeshCellsAABB(Mesh& M, bool reorder) {
1004 initialize(M, reorder ? AABB_INPLACE : AABB_NOREORDER);
1005 }
1006
1007 [[deprecated("use initialize(Mesh&,AABBReorderMode) instead")]]
1008 void initialize(Mesh& M, bool reorder) {
1009 initialize(M, reorder ? AABB_INPLACE : AABB_NOREORDER);
1010 }
1011 #endif
1012
1013 /**
1014 * \brief Finds the index of a tetrahedron that contains a query point
1015 * \param[in] p a const reference to the query point
1016 * \return the index of one of the tetrahedra that contains \p p or
1017 * NO_TET if \p p is outside the mesh.
1018 * \note The input mesh needs to be tetrahedralized. If the mesh has
1019 * arbitrary cells, then one may use instead containing_boxes().
1020 */
1021 index_t containing_tet(const vec3& p) const {
1022 geo_debug_assert(mesh_->cells.are_simplices());
1023 return containing_tet_recursive(
1024 p, 1, 0, mesh_->cells.nb()
1025 );
1026 }
1027
1028 /**
1029 * \brief Computes all the intersections between a given
1030 * box and the bounding boxes of all the cells.
1031 * \param[in] action a function that takes as argument
1032 * an index_t (cell index) invoked for all cells that
1033 * have a bounding box that intersects \p box_in.
1034 */
1035 void compute_bbox_cell_bbox_intersections(
1036 const Box& box_in,
1037 std::function<void(index_t)> action
1038 ) const {
1039 bbox_intersect_recursive(
1040 action, box_in, 1, 0, mesh_->cells.nb()
1041 );
1042 }
1043
1044 /**
1045 * \brief Finds all the cells such that their bounding
1046 * box contain a point.
1047 * \param[in] action a function that takes an index_t
1048 * that is invoked for all cells that have a bounding
1049 * box that contains \p p.
1050 */
1051 void containing_boxes(
1052 const vec3& p, std::function<void(index_t)> action
1053 ) const {
1054 containing_bboxes_recursive(
1055 action, p, 1, 0, mesh_->cells.nb()
1056 );
1057 }
1058
1059 /**
1060 * \brief Computes all the pairs of intersecting cells.
1061 * \param[in] action is a function that takes two index_t's,
1062 * invoked of all pairs of cells that have overlapping
1063 * bounding boxes. Further processing is necessary to
1064 * detect actual cell intersections.
1065 */
1066 void compute_cell_bbox_intersections(
1067 std::function<void(index_t, index_t)> action
1068 ) const {
1069 self_intersect_recursive(
1070 action,
1071 1, 0, mesh_->cells.nb(),
1072 1, 0, mesh_->cells.nb()
1073 );
1074 }
1075
1076 /**
1077 * \brief Computes all the pairs of intersecting cells between this
1078 * AABB and another one.
1079 * \param[in] action is a function that takes two index_t's,
1080 * invoked of all pairs of cells that have overlapping
1081 * bounding boxes. Further processing is necessary to
1082 * detect actual cell intersections.
1083 * \param[in] other the other AABB.
1084 */
1085 void compute_other_cell_bbox_intersections(
1086 MeshCellsAABB* other,
1087 std::function<void(index_t, index_t)> action
1088 ) const {
1089 other_intersect_recursive(
1090 action,
1091 1, 0, mesh_->cells.nb(),
1092 other,
1093 1, 0, other->mesh_->cells.nb()
1094 );
1095 }
1096
1097
1098 protected:
1099
1100 /**
1101 * \brief The recursive function used by the implementation
1102 * of containing_tet().
1103 * \param[in] p a const reference to the query point
1104 * \param[in] n index of the current node in the AABB tree
1105 * \param[in] b index of the first tet in the subtree under node \p n
1106 * \param[in] e one position past the index of the last tet in the
1107 * subtree under node \p n
1108 * \return the index of one of the tetrahedra that contains \p p, or
1109 * NO_TET if \p p is outside the mesh.
1110 */
1111 index_t containing_tet_recursive(
1112 const vec3& p,
1113 index_t n, index_t b, index_t e
1114 ) const;
1115
1116
1117 /**
1118 * \brief Computes all the cells that have a bbox that
1119 * contain a given point in a sub-tree of the AABB tree.
1120 *
1121 * Note that the tree structure is completely implicit,
1122 * therefore the bounds of the (continuous) facet indices
1123 * sequences that correspond to the facets contained
1124 * in the two nodes are sent as well as the node indices.
1125 *
1126 * \param[in] action a function that takes an index_t that is
1127 * invoked for all cells that has a bounding box that
1128 * contains \p p.
1129 * \param[in] p a const reference to the query point
1130 * \param[in] node index of the first node of the AABB tree
1131 * \param[in] b index of the first facet in \p node
1132 * \param[in] e one position past the index of the last
1133 * facet in \p node
1134 */
1135 void containing_bboxes_recursive(
1136 std::function<void(index_t)> action,
1137 const vec3& p,
1138 index_t node, index_t b, index_t e
1139 ) const {
1140 geo_debug_assert(e != b);
1141
1142 // Prune sub-tree that does not have intersection
1143 if(!bboxes_[node].contains(p)) {
1144 return;
1145 }
1146
1147 // Leaf case
1148 if(e == b+1) {
1149 action(element_in_leaf(b));
1150 return;
1151 }
1152
1153 // Recursion
1154 index_t m = b + (e - b) / 2;
1155 index_t node_l = 2 * node;
1156 index_t node_r = 2 * node + 1;
1157
1158 containing_bboxes_recursive(action, p, node_l, b, m);
1159 containing_bboxes_recursive(action, p, node_r, m, e);
1160 }
1161 };
1162
1163 /*******************************************************************/
1164
1165 /**
1166 * \brief Axis Aligned Bounding Box tree of mesh facets in 2D.
1167 * \details Used to quickly find the facet that contains
1168 * a given 2d point.
1169 */
1170 class GEOGRAM_API MeshFacetsAABB2d : public MeshAABB2d {
1171 public:
1172
1173 /**
1174 * \brief Symbolic constant for indicating that there
1175 * is no containing tetrahedron.
1176 * \see containing_tet()
1177 */
1178 static constexpr index_t NO_TRIANGLE = NO_INDEX;
1179
1180 /**
1181 * \brief MeshFacetsAABB2d constructor.
1182 * \details Creates an uninitialized MeshFacetsAABB2d.
1183 */
1184 MeshFacetsAABB2d();
1185
1186 /**
1187 * \brief Creates the Axis Aligned Bounding Boxes tree.
1188 * \param[in] M the input mesh. It can be modified,
1189 * The cells are re-ordered (using Morton's order, see mesh_reorder()).
1190 * \param[in] reorder if not set, Morton re-ordering is
1191 * skipped (but it means that mesh_reorder() was previously
1192 * called else the algorithm will be pretty unefficient).
1193 */
1194 MeshFacetsAABB2d(Mesh& M, bool reorder = true);
1195
1196 /**
1197 * \brief Initializes the Axis Aligned Bounding Boxes tree.
1198 * \param[in] M the input mesh. It can be modified,
1199 * The cells are re-ordered (using Morton's order, see mesh_reorder()).
1200 * \param[in] reorder if not set, Morton re-ordering is
1201 * skipped (but it means that mesh_reorder() was previously
1202 * called else the algorithm will be pretty unefficient).
1203 */
1204 void initialize(Mesh& M, bool reorder = true);
1205
1206 /**
1207 * \brief Finds the index of a facet that contains a query point
1208 * \param[in] p a const reference to the query point
1209 * \return the index of one of the facetthat contains \p p or
1210 * NO_TRIANGLE if \p p is outside the mesh.
1211 * \note The input mesh needs to be triangulated. If the mesh has
1212 * arbitrary cells, then one may use instead containing_boxes().
1213 */
1214 index_t containing_triangle(const vec2& p) const {
1215 geo_debug_assert(mesh_->facets.are_simplices());
1216 return containing_triangle_recursive(
1217 p, 1, 0, mesh_->facets.nb()
1218 );
1219 }
1220
1221 /**
1222 * \brief Computes all the intersections between a given
1223 * box and the bounding boxes of all the facets.
1224 * \param[in] action a function that takes as argument
1225 * an index_t (cell index) invoked for all cells that
1226 * have a bounding box that intersects \p box_in.
1227 */
1228 void compute_bbox_cell_bbox_intersections(
1229 const Box2d& box_in,
1230 std::function<void(index_t)> action
1231 ) const {
1232 bbox_intersect_recursive(
1233 action, box_in, 1, 0, mesh_->facets.nb()
1234 );
1235 }
1236
1237 /**
1238 * \brief Finds all the cells such that their bounding
1239 * box contain a point.
1240 * \param[in] action a function that takes an index_t
1241 * that is invoked for all cells that have a bounding
1242 * box that contains \p p.
1243 */
1244 void containing_boxes(
1245 const vec2& p, std::function<void(index_t)> action
1246 ) const {
1247 containing_bboxes_recursive(
1248 action, p, 1, 0, mesh_->facets.nb()
1249 );
1250 }
1251
1252 /**
1253 * \brief Computes all the pairs of intersecting facets.
1254 * \param[in] action is a function that takes two index_t's,
1255 * invoked of all pairs of cells that have overlapping
1256 * bounding boxes. Further processing is necessary to
1257 * detect actual cell intersections.
1258 */
1259 void compute_facet_bbox_intersections(
1260 std::function<void(index_t, index_t)> action
1261 ) const {
1262 self_intersect_recursive(
1263 action,
1264 1, 0, mesh_->facets.nb(),
1265 1, 0, mesh_->facets.nb()
1266 );
1267 }
1268
1269 /**
1270 * \brief Computes all the pairs of intersecting cells between this
1271 * AABB and another one.
1272 * \param[in] action is a function that takes two index_t's,
1273 * invoked of all pairs of cells that have overlapping
1274 * bounding boxes. Further processing is necessary to
1275 * detect actual cell intersections.
1276 * \param[in] other the other AABB.
1277 */
1278 void compute_other_cell_bbox_intersections(
1279 MeshFacetsAABB2d* other,
1280 std::function<void(index_t, index_t)> action
1281 ) const {
1282 other_intersect_recursive(
1283 action,
1284 1, 0, mesh_->facets.nb(),
1285 other,
1286 1, 0, other->mesh_->facets.nb()
1287 );
1288 }
1289
1290
1291 protected:
1292
1293 /**
1294 * \brief The recursive function used by the implementation
1295 * of containing_triangle().
1296 * \param[in] p a const reference to the query point
1297 * \param[in] n index of the current node in the AABB tree
1298 * \param[in] b index of the first tet in the subtree under node \p n
1299 * \param[in] e one position past the index of the last tet in the
1300 * subtree under node \p n
1301 * \return the index of one of the tetrahedra that contains \p p, or
1302 * NO_TRIANGLE if \p p is outside the mesh.
1303 */
1304 index_t containing_triangle_recursive(
1305 const vec2& p,
1306 index_t n, index_t b, index_t e
1307 ) const;
1308
1309
1310 /**
1311 * \brief Computes all the cells that have a bbox that
1312 * contain a given point in a sub-tree of the AABB tree.
1313 *
1314 * Note that the tree structure is completely implicit,
1315 * therefore the bounds of the (continuous) facet indices
1316 * sequences that correspond to the facets contained
1317 * in the two nodes are sent as well as the node indices.
1318 *
1319 * \param[in] action a function that takes an index_t that is
1320 * invoked for all cells that has a bounding box that
1321 * contains \p p.
1322 * \param[in] p a const reference to the query point
1323 * \param[in] node index of the first node of the AABB tree
1324 * \param[in] b index of the first facet in \p node
1325 * \param[in] e one position past the index of the last
1326 * facet in \p node
1327 */
1328 void containing_bboxes_recursive(
1329 std::function<void(index_t)> action,
1330 const vec2& p,
1331 index_t node, index_t b, index_t e
1332 ) const {
1333 geo_debug_assert(e != b);
1334
1335 // Prune sub-tree that does not have intersection
1336 if(!bboxes_[node].contains(p)) {
1337 return;
1338 }
1339
1340 // Leaf case
1341 if(e == b+1) {
1342 action(element_in_leaf(b));
1343 return;
1344 }
1345
1346 // Recursion
1347 index_t m = b + (e - b) / 2;
1348 index_t node_l = 2 * node;
1349 index_t node_r = 2 * node + 1;
1350
1351 containing_bboxes_recursive(action, p, node_l, b, m);
1352 containing_bboxes_recursive(action, p, node_r, m, e);
1353 }
1354 };
1355
1356
1357 }
1358
1359 #endif
1360