| 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 |
1/2✓ Branch 2 taken 83 times.
✗ Branch 3 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 | 5600588 | 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 | geo_debug_assert(e1 != b1); | ||
| 173 | 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 5486308 times.
✓ Branch 1 taken 114280 times.
|
5600588 | if(e2 <= b1) { |
| 180 | return; | ||
| 181 | } | ||
| 182 | |||
| 183 | // The acceleration is here: | ||
| 184 | if( | ||
| 185 |
4/4✓ Branch 0 taken 5257528 times.
✓ Branch 1 taken 228780 times.
✓ Branch 2 taken 3666623 times.
✓ Branch 3 taken 1590905 times.
|
10743836 | (node1 != node2) && |
| 186 | !bboxes_overlap(bboxes_[node1], bboxes_[node2]) | ||
| 187 | ) { | ||
| 188 | return; | ||
| 189 | } | ||
| 190 | |||
| 191 | // Simple case: leaf - leaf intersection. | ||
| 192 |
4/4✓ Branch 0 taken 1879599 times.
✓ Branch 1 taken 2015804 times.
✓ Branch 2 taken 1095851 times.
✓ Branch 3 taken 783748 times.
|
3895403 | 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 1603961 times.
✓ Branch 1 taken 1195591 times.
|
2799552 | if(e2 - b2 > e1 - b1) { |
| 204 | 1603961 | index_t m2 = b2 + (e2 - b2) / 2; | |
| 205 | 1603961 | index_t node2_l = 2 * node2; | |
| 206 | 1603961 | index_t node2_r = 2 * node2 + 1; | |
| 207 |
1/2✓ Branch 2 taken 1603961 times.
✗ Branch 3 not taken.
|
1603961 | self_intersect_recursive(action, node1, b1, e1, node2_l, b2, m2); |
| 208 |
1/2✓ Branch 2 taken 1603961 times.
✗ Branch 3 not taken.
|
3207922 | self_intersect_recursive(action, node1, b1, e1, node2_r, m2, e2); |
| 209 | } else { | ||
| 210 | 1195591 | index_t m1 = b1 + (e1 - b1) / 2; | |
| 211 | 1195591 | index_t node1_l = 2 * node1; | |
| 212 | 1195591 | index_t node1_r = 2 * node1 + 1; | |
| 213 |
1/2✓ Branch 2 taken 1195591 times.
✗ Branch 3 not taken.
|
1195591 | self_intersect_recursive(action, node1_l, b1, m1, node2, b2, e2); |
| 214 |
1/2✓ Branch 2 taken 1195591 times.
✗ Branch 3 not taken.
|
2391182 | 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 | geo_debug_assert(e > b); | ||
| 297 |
2/2✓ Branch 0 taken 192651 times.
✓ Branch 1 taken 192734 times.
|
770770 | if(b + 1 == e) { |
| 298 | 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 | return std::max( | ||
| 304 | 385302 | max_node_index(childl, b, m), | |
| 305 | 385302 | max_node_index(childr, m, e) | |
| 306 | ); | ||
| 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 | geo_debug_assert(node_index < bboxes_.size()); | ||
| 325 | 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 | geo_debug_assert(childl < bboxes_.size()); | ||
| 334 | geo_debug_assert(childr < bboxes_.size()); | ||
| 335 |
1/2✓ Branch 2 taken 192651 times.
✗ Branch 3 not taken.
|
385302 | init_bboxes_recursive(childl, b, m, get_bbox); |
| 336 |
1/2✓ Branch 2 taken 192651 times.
✗ Branch 3 not taken.
|
385302 | init_bboxes_recursive(childr, m, e, get_bbox); |
| 337 | geo_debug_assert(childl < bboxes_.size()); | ||
| 338 | geo_debug_assert(childr < bboxes_.size()); | ||
| 339 | 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 | bool indirect() const { | ||
| 350 | 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 | index_t element_in_leaf(index_t i) const { | ||
| 361 |
6/30✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 192734 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 981351 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 981351 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 491836 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 2247907 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 800714 times.
✗ Branch 29 not taken.
|
4714542 | 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 | 88 | class GEOGRAM_API MeshAABB3d : public AABB3d { | |
| 406 | public: | ||
| 407 | /** | ||
| 408 | * \brief MeshAABB3d constructor. | ||
| 409 | */ | ||
| 410 |
2/3✓ Branch 0 taken 1 times.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
82 | MeshAABB3d() : mesh_(nullptr) { |
| 411 | } | ||
| 412 | |||
| 413 | /** | ||
| 414 | * \brief Gets the mesh. | ||
| 415 | * \return a const reference to the mesh. | ||
| 416 | */ | ||
| 417 | const Mesh* mesh() const { | ||
| 418 |
1/2✓ Branch 0 taken 388 times.
✗ Branch 1 not taken.
|
388 | 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 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
88 | 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 |
0/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
89998 | Intersection() : |
| 441 | 89998 | t(Numeric::max_float64()), | |
| 442 | 89998 | f(NO_INDEX), | |
| 443 |
0/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
89998 | i(NO_INDEX), j(NO_INDEX), k(NO_INDEX) |
| 444 | { | ||
| 445 | } | ||
| 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 | MeshFacetsAABB() { | ||
| 460 | } | ||
| 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 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
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 |
1/2✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
|
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 | if( | ||
| 530 |
3/4✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 46 times.
✓ Branch 4 taken 14 times.
|
60 | Process::maximum_concurrent_threads() <= 1 || |
| 531 |
2/2✓ Branch 0 taken 46 times.
✓ Branch 1 taken 14 times.
|
60 | mesh_->facets.nb() <= 1024 |
| 532 | ) { | ||
| 533 |
1/2✓ Branch 1 taken 46 times.
✗ Branch 2 not taken.
|
92 | self_intersect_recursive( |
| 534 | action, | ||
| 535 | 1, 0, mesh_->facets.nb(), | ||
| 536 | 46 | 1, 0, mesh_->facets.nb() | |
| 537 | ); | ||
| 538 | } else { | ||
| 539 |
1/2✓ Branch 2 taken 14 times.
✗ Branch 3 not taken.
|
28 | 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 0 not taken.
✓ Branch 1 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 | 800714 | get_nearest_facet_hint(p, nearest_facet, nearest_point, sq_dist); | |
| 577 | 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 | double squared_distance(const vec3& p) const { | ||
| 675 | vec3 nearest_point; | ||
| 676 | double result; | ||
| 677 | 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 |