| 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 | #include <geogram/mesh/mesh_AABB.h> | ||
| 41 | #include <geogram/mesh/mesh_reorder.h> | ||
| 42 | #include <geogram/mesh/mesh_geometry.h> | ||
| 43 | #include <geogram/mesh/mesh_repair.h> | ||
| 44 | #include <geogram/numerics/predicates.h> | ||
| 45 | #include <geogram/basic/geometry_nd.h> | ||
| 46 | #include <geogram/basic/algorithm.h> | ||
| 47 | |||
| 48 | #include <stack> | ||
| 49 | |||
| 50 | namespace { | ||
| 51 | |||
| 52 | using namespace GEO; | ||
| 53 | |||
| 54 | /** | ||
| 55 | * \brief Finds the nearest point in a mesh facet from a query point. | ||
| 56 | * \param[in] M the mesh | ||
| 57 | * \param[in] p the query point | ||
| 58 | * \param[in] f index of the facet in \p M | ||
| 59 | * \param[out] nearest_p the point of facet \p f nearest to \p p | ||
| 60 | * \param[out] squared_dist the squared distance between | ||
| 61 | * \p p and \p nearest_p | ||
| 62 | */ | ||
| 63 | 2247907 | void get_point_facet_nearest_point( | |
| 64 | const Mesh& M, | ||
| 65 | const vec3& p, | ||
| 66 | index_t f, | ||
| 67 | vec3& nearest_p, | ||
| 68 | double& squared_dist | ||
| 69 | ) { | ||
| 70 | 2247907 | squared_dist = Numeric::max_float64(); | |
| 71 |
8/14✓ Branch 1 taken 2247907 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2247907 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2247907 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2247907 times.
✗ Branch 11 not taken.
✓ Branch 16 taken 2247907 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 4495814 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 2247907 times.
✓ Branch 22 taken 2247907 times.
|
6743721 | for(auto [ p1, p2, p3] : M.facets.triangle_points(f)) { |
| 72 | double lambda1, lambda2, lambda3; // barycentric coords,unused. | ||
| 73 | 2247907 | vec3 cur_nearest_p; | |
| 74 |
1/2✓ Branch 1 taken 2247907 times.
✗ Branch 2 not taken.
|
2247907 | double cur_squared_dist = Geom::point_triangle_squared_distance( |
| 75 | p, p1, p2, p3, cur_nearest_p, lambda1, lambda2, lambda3 | ||
| 76 | ); | ||
| 77 |
1/2✓ Branch 0 taken 2247907 times.
✗ Branch 1 not taken.
|
2247907 | if(cur_squared_dist < squared_dist) { |
| 78 | 2247907 | squared_dist = cur_squared_dist; | |
| 79 | 2247907 | nearest_p = cur_nearest_p; | |
| 80 | } | ||
| 81 | } | ||
| 82 | 2247907 | } | |
| 83 | |||
| 84 | /** | ||
| 85 | * \brief Computes the squared distance between a point and a Box. | ||
| 86 | * \param[in] p the point | ||
| 87 | * \param[in] B the box | ||
| 88 | * \return the squared distance between \p p and \p B | ||
| 89 | * \pre p is inside B | ||
| 90 | */ | ||
| 91 | 13295474 | double inner_point_box_squared_distance( | |
| 92 | const vec3& p, | ||
| 93 | const Box& B | ||
| 94 | ) { | ||
| 95 |
2/8✓ Branch 1 taken 13295474 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 13295474 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
13295474 | geo_debug_assert(B.contains(p)); |
| 96 |
1/2✓ Branch 1 taken 13295474 times.
✗ Branch 2 not taken.
|
13295474 | double result = geo_sqr(p[0] - B.xyz_min[0]); |
| 97 |
1/2✓ Branch 1 taken 13295474 times.
✗ Branch 2 not taken.
|
13295474 | result = std::min(result, geo_sqr(p[0] - B.xyz_max[0])); |
| 98 |
2/2✓ Branch 0 taken 26590948 times.
✓ Branch 1 taken 13295474 times.
|
39886422 | for(coord_index_t c = 1; c < 3; ++c) { |
| 99 |
1/2✓ Branch 1 taken 26590948 times.
✗ Branch 2 not taken.
|
26590948 | result = std::min(result, geo_sqr(p[c] - B.xyz_min[c])); |
| 100 |
1/2✓ Branch 1 taken 26590948 times.
✗ Branch 2 not taken.
|
26590948 | result = std::min(result, geo_sqr(p[c] - B.xyz_max[c])); |
| 101 | } | ||
| 102 | 13295474 | return result; | |
| 103 | } | ||
| 104 | |||
| 105 | /** | ||
| 106 | * \brief Computes the squared distance between a point and a Box | ||
| 107 | * with negative sign if the point is inside the Box. | ||
| 108 | * \param[in] p the point | ||
| 109 | * \param[in] B the box | ||
| 110 | * \return the signed squared distance between \p p and \p B | ||
| 111 | */ | ||
| 112 | 29434698 | double point_box_signed_squared_distance( | |
| 113 | const vec3& p, | ||
| 114 | const Box& B | ||
| 115 | ) { | ||
| 116 | 29434698 | bool inside = true; | |
| 117 | 29434698 | double result = 0.0; | |
| 118 |
2/2✓ Branch 0 taken 88304094 times.
✓ Branch 1 taken 29434698 times.
|
117738792 | for(coord_index_t c = 0; c < 3; c++) { |
| 119 |
2/2✓ Branch 1 taken 11496525 times.
✓ Branch 2 taken 76807569 times.
|
88304094 | if(p[c] < B.xyz_min[c]) { |
| 120 | 11496525 | inside = false; | |
| 121 | 11496525 | result += geo_sqr(p[c] - B.xyz_min[c]); | |
| 122 |
2/2✓ Branch 1 taken 10746041 times.
✓ Branch 2 taken 66061528 times.
|
76807569 | } else if(p[c] > B.xyz_max[c]) { |
| 123 | 10746041 | inside = false; | |
| 124 | 10746041 | result += geo_sqr(p[c] - B.xyz_max[c]); | |
| 125 | } | ||
| 126 | } | ||
| 127 |
2/2✓ Branch 0 taken 13295474 times.
✓ Branch 1 taken 16139224 times.
|
29434698 | if(inside) { |
| 128 | 13295474 | result = -inner_point_box_squared_distance(p, B); | |
| 129 | } | ||
| 130 | 29434698 | return result; | |
| 131 | } | ||
| 132 | |||
| 133 | /** | ||
| 134 | * \brief Computes the squared distance between a point and the | ||
| 135 | * center of a box. | ||
| 136 | * \param[in] p the point | ||
| 137 | * \param[in] B the box | ||
| 138 | * \return the squared distance between \p p and the center of \p B | ||
| 139 | */ | ||
| 140 | 15850044 | double point_box_center_squared_distance( | |
| 141 | const vec3& p, const Box& B | ||
| 142 | ) { | ||
| 143 | 15850044 | double result = 0.0; | |
| 144 |
2/2✓ Branch 0 taken 47550132 times.
✓ Branch 1 taken 15850044 times.
|
63400176 | for(coord_index_t c = 0; c < 3; ++c) { |
| 145 | 47550132 | double d = p[c] - 0.5 * (B.xyz_min[c] + B.xyz_max[c]); | |
| 146 | 47550132 | result += geo_sqr(d); | |
| 147 | } | ||
| 148 | 15850044 | return result; | |
| 149 | } | ||
| 150 | |||
| 151 | /** | ||
| 152 | * \brief Tests whether a mesh tetrahedron contains a given point | ||
| 153 | * \param[in] M a const reference to the mesh | ||
| 154 | * \param[in] t the index of the tetrahedron in \p M | ||
| 155 | * \param[in] p a const reference to the point | ||
| 156 | * \retval true if the tetrahedron \p t or its boundary contains | ||
| 157 | * the point \p p | ||
| 158 | * \retval false otherwise | ||
| 159 | */ | ||
| 160 | ✗ | bool mesh_tet_contains_point( | |
| 161 | const Mesh& M, index_t t, const vec3& p | ||
| 162 | ) { | ||
| 163 | ✗ | const vec3& p0 = M.cells.point(t,0); | |
| 164 | ✗ | const vec3& p1 = M.cells.point(t,1); | |
| 165 | ✗ | const vec3& p2 = M.cells.point(t,2); | |
| 166 | ✗ | const vec3& p3 = M.cells.point(t,3); | |
| 167 | |||
| 168 | Sign s[4]; | ||
| 169 | ✗ | s[0] = PCK::orient_3d(p, p1, p2, p3); | |
| 170 | ✗ | s[1] = PCK::orient_3d(p0, p, p2, p3); | |
| 171 | ✗ | s[2] = PCK::orient_3d(p0, p1, p, p3); | |
| 172 | ✗ | s[3] = PCK::orient_3d(p0, p1, p2, p); | |
| 173 | |||
| 174 | return ( | ||
| 175 | ✗ | (s[0] >= 0 && s[1] >= 0 && s[2] >= 0 && s[3] >= 0) || | |
| 176 | ✗ | (s[0] <= 0 && s[1] <= 0 && s[2] <= 0 && s[3] <= 0) | |
| 177 | ✗ | ); | |
| 178 | } | ||
| 179 | |||
| 180 | |||
| 181 | /** | ||
| 182 | * \brief Tests whether a mesh triangle contains a given point | ||
| 183 | * \param[in] M a const reference to the mesh | ||
| 184 | * \param[in] t the index of the triangle in \p M | ||
| 185 | * \param[in] p a const reference to the point | ||
| 186 | * \retval true if the triangle \p t or its boundary contains | ||
| 187 | * the point \p p | ||
| 188 | * \retval false otherwise | ||
| 189 | */ | ||
| 190 | ✗ | bool mesh_triangle_contains_point( | |
| 191 | const Mesh& M, index_t t, const vec2& p | ||
| 192 | ) { | ||
| 193 | ✗ | index_t i = M.facets.vertex(t,0); | |
| 194 | ✗ | index_t j = M.facets.vertex(t,1); | |
| 195 | ✗ | index_t k = M.facets.vertex(t,2); | |
| 196 | ✗ | vec2 p0 = M.vertices.point<2>(i); | |
| 197 | ✗ | vec2 p1 = M.vertices.point<2>(j); | |
| 198 | ✗ | vec2 p2 = M.vertices.point<2>(k); | |
| 199 | |||
| 200 | Sign s[3]; | ||
| 201 | ✗ | s[0] = PCK::orient_2d(p, p1, p2); | |
| 202 | ✗ | s[1] = PCK::orient_2d(p0, p, p2); | |
| 203 | ✗ | s[2] = PCK::orient_2d(p0, p1, p ); | |
| 204 | |||
| 205 | return ( | ||
| 206 | ✗ | (s[0] >= 0 && s[1] >= 0 && s[2] >= 0 ) || | |
| 207 | ✗ | (s[0] <= 0 && s[1] <= 0 && s[2] <= 0 ) | |
| 208 | ✗ | ); | |
| 209 | } | ||
| 210 | |||
| 211 | /** | ||
| 212 | * \brief Computes the intersection between a ray and a triangle. | ||
| 213 | * \param[in] O origin of the ray. | ||
| 214 | * \param[in] D direction vector of the ray. | ||
| 215 | * \param[in] A , B , C the three vertices of the triangle. | ||
| 216 | * \param[out] t the intersection point is O + t D when it exists. | ||
| 217 | * \param[out] u , v the intersection point is A + u (B-A) + v (C-A). | ||
| 218 | * when it exists. | ||
| 219 | * \param[out] N the normal to the triangle. | ||
| 220 | * \param[in] bidirectional if set, computes a line-triangle intersection | ||
| 221 | * (intersections before \p O are also reported) | ||
| 222 | * \retval true if there is an intersection point | ||
| 223 | * \retval false otherwise | ||
| 224 | */ | ||
| 225 | 492438 | bool ray_triangle_intersection( | |
| 226 | const vec3& O, const vec3& D, | ||
| 227 | const vec3& A, const vec3& B, const vec3& C, | ||
| 228 | double& t, double& u, double& v, vec3& N, | ||
| 229 | bool bidirectional = false | ||
| 230 | ) { | ||
| 231 | // M\"oller and Trumbore, | ||
| 232 | // Fast, Minimum Storage Ray-Triangle Intersection, | ||
| 233 | // Journal of Graphics Tools, vol. 2, 1997, p. 21-28 | ||
| 234 | // (with small adaptations: branchless, and reusing the normal vector) | ||
| 235 | // | ||
| 236 | // Let E1 = B-A; E2 = C-A, write ray eqn (1) and triangle eqn (2), then | ||
| 237 | // write equality between (1) and (2) at intersection point (3): | ||
| 238 | // | ||
| 239 | // (1) O + tD = A + uE1 + vE2 | ||
| 240 | // (2) uE1 + vE2 -tD = O-A | ||
| 241 | // | ||
| 242 | // [u] | ||
| 243 | // (3) [E1|E2|-D] [v] = O-A | ||
| 244 | // [t] | ||
| 245 | // | ||
| 246 | // (where [E1|E2|-D] is the 3x3 matrix with E1,E2,-D as its columns) | ||
| 247 | // | ||
| 248 | // Using Cramer's formula for the solution of: | ||
| 249 | // | ||
| 250 | // [a11 a12 a13][x1] [b1] | ||
| 251 | // [a12 a22 a23][x2] = [b2] | ||
| 252 | // [a31 a32 a33][x3] [b3] | ||
| 253 | // | ||
| 254 | // gives: | ||
| 255 | // | ||
| 256 | // |b1 a12 a13| |a11 a12 a13| | ||
| 257 | // x1 = |b2 a22 a23| / |a21 a22 a23| | ||
| 258 | // |b3 a32 a33| |a31 a32 a33| | ||
| 259 | // | ||
| 260 | // |a11 b1 a13| |a11 a12 a13| | ||
| 261 | // x2 = |a21 b2 a23| / |a21 a22 a23| | ||
| 262 | // |a31 b3 a33| |a31 a32 a33| | ||
| 263 | // | ||
| 264 | // |a11 a12 b1| |a11 a12 a13| | ||
| 265 | // x3 = |a21 a22 b2| / |a21 a22 a23| | ||
| 266 | // |a31 a32 b3| |a31 a32 a33| | ||
| 267 | // | ||
| 268 | // Now we get: | ||
| 269 | // | ||
| 270 | // u = (O-A,E2,-D) / (E1,E2,-D) | ||
| 271 | // v = (E1,O-A,-D) / (E1,E2,-D) | ||
| 272 | // t = (E1,E2,O-A) / (E1,E2,-D) | ||
| 273 | // | ||
| 274 | // where (A,B,C) denotes the determinant of the 3x3 matrix | ||
| 275 | // with A,B,C as its column vectors. | ||
| 276 | // | ||
| 277 | // Now we use the following identities: | ||
| 278 | // (A,B,C) = dot(A,cross(B,C)) (develop the det w.r.t. first column) | ||
| 279 | // (B,A,C) = -(A,B,C) (swapping two cols changes sign) | ||
| 280 | // (B,C,A) = (A,B,C) (circular perm does not change sign) | ||
| 281 | // | ||
| 282 | // Now we get: | ||
| 283 | // | ||
| 284 | // u = -(E2,O-A,D) / (D,E1,E2) | ||
| 285 | // v = (E1,O-A,D) / (D,E1,E2) | ||
| 286 | // t = -(O-A,E1,E2) / (D,E1,E2) | ||
| 287 | // | ||
| 288 | // Using N=cross(E1,E2); AO = O-A; DAO = cross(D,AO) | ||
| 289 | 492438 | vec3 E1(B-A); | |
| 290 | 492438 | vec3 E2(C-A); | |
| 291 | 492438 | N = cross(E1,E2); | |
| 292 | 492438 | double det = -dot(D,N); | |
| 293 | 492438 | double invdet = 1.0/det; | |
| 294 | 492438 | vec3 AO = O - A; | |
| 295 | 492438 | vec3 DAO = cross(AO,D); | |
| 296 | 492438 | u = dot(E2,DAO) * invdet; | |
| 297 | 492438 | v = -dot(E1,DAO) * invdet; | |
| 298 | 492438 | t = dot(AO,N) * invdet; | |
| 299 | return ( | ||
| 300 |
1/2✓ Branch 0 taken 492438 times.
✗ Branch 1 not taken.
|
492438 | (fabs(det) >= 1e-20) && |
| 301 |
2/2✓ Branch 0 taken 353033 times.
✓ Branch 1 taken 139405 times.
|
492438 | (bidirectional || t >= 0.0) && |
| 302 |
2/2✓ Branch 0 taken 294903 times.
✓ Branch 1 taken 58130 times.
|
353033 | (u >= 0.0) && |
| 303 |
3/4✓ Branch 0 taken 492438 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 232514 times.
✓ Branch 3 taken 62389 times.
|
1217390 | (v >= 0.0) && |
| 304 |
2/2✓ Branch 0 taken 166675 times.
✓ Branch 1 taken 65839 times.
|
232514 | ((u+v) <= 1.0) |
| 305 | 492438 | ); | |
| 306 | } | ||
| 307 | |||
| 308 | 7851352 | inline double max3(double x1, double x2, double x3) { | |
| 309 | 7851352 | return std::max(x1,std::max(x2,x3)); | |
| 310 | } | ||
| 311 | |||
| 312 | 7851352 | inline double min3(double x1, double x2, double x3) { | |
| 313 | 7851352 | return std::min(x1,std::min(x2,x3)); | |
| 314 | } | ||
| 315 | |||
| 316 | |||
| 317 | // https://tavianator.com/fast-branchless-raybounding-box-intersections/ | ||
| 318 | // https://tavianator.com/fast-branchless-raybounding-box-intersections-part-2-nans/ | ||
| 319 | // http://www.flipcode.com/archives/SSE_RayBox_Intersection_Test.shtml | ||
| 320 | // http://psgraphics.blogspot.com/2016/02/ray-box-intersection-and-fmin.html | ||
| 321 | |||
| 322 | /** | ||
| 323 | * \brief Tests whether a segment intersects a box. | ||
| 324 | * \param[in] q1 the first extremity of the segment. | ||
| 325 | * \param[in] dirinv precomputed 1/(q2.x-q1.x), 1/(q2.y-q1.y), 1/(q2.z-q1.z) | ||
| 326 | * where q2 denotes the second extremity of the segment. | ||
| 327 | * \param[in] box the box. | ||
| 328 | * \param[in] T the maximum acceptable value for the intersection parameter. | ||
| 329 | * Can be used to early-prune boxes while traversing the tree. | ||
| 330 | * \param[in] bidirectional if set, computes a line-box intersection | ||
| 331 | * (intersections before \p q1 are also reported) | ||
| 332 | * \retval true if [q1,q2] intersects the box. | ||
| 333 | * \retval false otherwise. | ||
| 334 | */ | ||
| 335 | 7851352 | bool ray_box_intersection( | |
| 336 | const vec3& q1, const vec3& dirinv, const Box& box, double T = 1.0, | ||
| 337 | bool bidirectional = false | ||
| 338 | ) { | ||
| 339 | // This version: slab method. | ||
| 340 | // Step 1: compute | ||
| 341 | // (tx1, tx2) : parameters of intersection with slab {xmin <= x <= xmax} | ||
| 342 | // (ty1, ty2) : parameters of intersection with slab {ymin <= y <= ymax} | ||
| 343 | // (tz1, tz2) : parameters of intersection with slab {zmin <= z <= zmax} | ||
| 344 | // (note: they are unordered, it is possible that tx1 > tx2) | ||
| 345 | // This defines three intervals: | ||
| 346 | // Ix = [ min(tx1,tx2) ... max(tx1,tx2) ] | ||
| 347 | // Iy = [ min(ty1,ty2) ... max(ty1,ty2) ] | ||
| 348 | // Iz = [ min(tz1,tz2) ... max(tz1,tz2) ] | ||
| 349 | // The intersection between [q1,q2] and the slab {xmin <= x <= xmax} is | ||
| 350 | // the set of points {q1 + t(q2-q1)} where t in Ix | ||
| 351 | |||
| 352 | // Q: what does it do if one of the fracs is zero ? | ||
| 353 | // normally the tests with inf do what they should | ||
| 354 | // (to be tested) | ||
| 355 | |||
| 356 | 7851352 | double tx1 = dirinv.x*(box.xyz_min[0] - q1.x); | |
| 357 | 7851352 | double tx2 = dirinv.x*(box.xyz_max[0] - q1.x); | |
| 358 | |||
| 359 | 7851352 | double ty1 = dirinv.y*(box.xyz_min[1] - q1.y); | |
| 360 | 7851352 | double ty2 = dirinv.y*(box.xyz_max[1] - q1.y); | |
| 361 | |||
| 362 | 7851352 | double tz1 = dirinv.z*(box.xyz_min[2] - q1.z); | |
| 363 | 7851352 | double tz2 = dirinv.z*(box.xyz_max[2] - q1.z); | |
| 364 | |||
| 365 | // now compute the intersection of the three intervals | ||
| 366 | // Ix /\ Iy /\ Iz | ||
| 367 | // this gives us the range of t that corresponds to points in the | ||
| 368 | // box (because the box is the intersection of the 3 slabs) | ||
| 369 | // it starts at the maximum of the left bounds of the 3 intervals | ||
| 370 | // it stops at the minimum of the right bounds of the 3 intervals | ||
| 371 | |||
| 372 | double tmin = | ||
| 373 |
1/2✓ Branch 4 taken 7851352 times.
✗ Branch 5 not taken.
|
7851352 | max3(std::min(tx1,tx2), std::min(ty1,ty2), std::min(tz1,tz2)); |
| 374 | |||
| 375 | double tmax = | ||
| 376 |
1/2✓ Branch 4 taken 7851352 times.
✗ Branch 5 not taken.
|
7851352 | min3(std::max(tx1,tx2), std::max(ty1,ty2), std::max(tz1,tz2)); |
| 377 | |||
| 378 | // There is no intersection if the interval is empty (tmin > tmax) | ||
| 379 | // or if the interval is outside [0,T] | ||
| 380 | // Note: the test is tmin <= tmax, because a bbox can be infinitely | ||
| 381 | // thin (for instance, the bbox of a triangle orthogonal to one | ||
| 382 | // of the axes). | ||
| 383 | |||
| 384 |
7/8✓ Branch 0 taken 7851352 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5881216 times.
✓ Branch 3 taken 1970136 times.
✓ Branch 4 taken 4461872 times.
✓ Branch 5 taken 1419344 times.
✓ Branch 6 taken 4327386 times.
✓ Branch 7 taken 134486 times.
|
15702704 | return (bidirectional || tmax >= 0.0) && (tmin <= tmax) && (tmin <= T); |
| 385 | } | ||
| 386 | } | ||
| 387 | |||
| 388 | /****************************************************************************/ | ||
| 389 | |||
| 390 | namespace GEO { | ||
| 391 | |||
| 392 | /************************************************************************/ | ||
| 393 | |||
| 394 | 83 | void MeshFacetsAABB::initialize(Mesh& M, AABBReorderMode reorder_mode) { | |
| 395 | 83 | mesh_ = &M; | |
| 396 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 83 times.
|
83 | if(mesh_->facets.nb() == 0) { |
| 397 | ✗ | return; | |
| 398 | } | ||
| 399 |
1/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 83 times.
✗ Branch 3 not taken.
|
83 | switch(reorder_mode) { |
| 400 | ✗ | case AABB_NOREORDER: | |
| 401 | ✗ | break; | |
| 402 | ✗ | case AABB_INPLACE: | |
| 403 | ✗ | reorder_.clear(); | |
| 404 | ✗ | mesh_reorder(*mesh_, MESH_ORDER_MORTON, MESH_FACETS); | |
| 405 | ✗ | break; | |
| 406 | 83 | case AABB_INDIRECT: | |
| 407 | 83 | compute_mesh_elements_spatial_order( | |
| 408 | 83 | *mesh_, MESH_FACETS, reorder_, MESH_ORDER_MORTON | |
| 409 | ); | ||
| 410 | 83 | break; | |
| 411 | } | ||
| 412 |
1/2✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
|
83 | AABB::initialize( |
| 413 | 83 | mesh_->facets.nb(), | |
| 414 | 166 | [this](Box& B, index_t f) { | |
| 415 | // Get facet bbox | ||
| 416 |
2/2✓ Branch 0 taken 578202 times.
✓ Branch 1 taken 192734 times.
|
770936 | for(coord_index_t coord = 0; coord < 3; ++coord) { |
| 417 | 578202 | B.xyz_min[coord] = Numeric::max_float64(); | |
| 418 | 578202 | B.xyz_max[coord] = -Numeric::max_float64(); | |
| 419 | } | ||
| 420 |
6/10✓ Branch 1 taken 192734 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 192734 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 192734 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 578380 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 578380 times.
✓ Branch 15 taken 192734 times.
|
771114 | for(const vec3& p: mesh_->facets.points(f)) { |
| 421 |
2/2✓ Branch 0 taken 1735140 times.
✓ Branch 1 taken 578380 times.
|
2313520 | for(coord_index_t coord = 0; coord < 3; ++coord) { |
| 422 |
1/2✓ Branch 1 taken 1735140 times.
✗ Branch 2 not taken.
|
1735140 | B.xyz_min[coord] = std::min(B.xyz_min[coord], p[coord]); |
| 423 |
1/2✓ Branch 1 taken 1735140 times.
✗ Branch 2 not taken.
|
1735140 | B.xyz_max[coord] = std::max(B.xyz_max[coord], p[coord]); |
| 424 | } | ||
| 425 | } | ||
| 426 | 192734 | } | |
| 427 | ); | ||
| 428 | } | ||
| 429 | |||
| 430 | 800714 | void MeshFacetsAABB::get_nearest_facet_hint( | |
| 431 | const vec3& p, | ||
| 432 | index_t& nearest_f, vec3& nearest_point, double& sq_dist | ||
| 433 | ) const { | ||
| 434 | |||
| 435 | // Find a good initial value for nearest_f by traversing | ||
| 436 | // the boxes and selecting the child such that the center | ||
| 437 | // of its bounding box is nearer to the query point. | ||
| 438 | // For a large mesh (20M facets) this gains up to 10% | ||
| 439 | // performance as compared to picking nearest_f randomly. | ||
| 440 | 800714 | index_t b = 0; | |
| 441 | 800714 | index_t e = mesh_->facets.nb(); | |
| 442 | 800714 | index_t n = 1; | |
| 443 | |||
| 444 |
2/2✓ Branch 0 taken 7925022 times.
✓ Branch 1 taken 800714 times.
|
8725736 | while(e != b + 1) { |
| 445 | 7925022 | index_t m = b + (e - b) / 2; | |
| 446 | 7925022 | index_t childl = 2 * n; | |
| 447 | 7925022 | index_t childr = 2 * n + 1; | |
| 448 | 7925022 | if( | |
| 449 | 7925022 | point_box_center_squared_distance(p, bboxes_[childl]) < | |
| 450 |
2/2✓ Branch 2 taken 3773136 times.
✓ Branch 3 taken 4151886 times.
|
7925022 | point_box_center_squared_distance(p, bboxes_[childr]) |
| 451 | ) { | ||
| 452 | 3773136 | e = m; | |
| 453 | 3773136 | n = childl; | |
| 454 | } else { | ||
| 455 | 4151886 | b = m; | |
| 456 | 4151886 | n = childr; | |
| 457 | } | ||
| 458 | } | ||
| 459 | 800714 | nearest_f = element_in_leaf(b); | |
| 460 | |||
| 461 | 800714 | index_t v = mesh_->facets.vertex(nearest_f, 0); | |
| 462 | 800714 | nearest_point = mesh_->vertices.point(v); | |
| 463 | 800714 | sq_dist = Geom::distance2(p, nearest_point); | |
| 464 | 800714 | } | |
| 465 | |||
| 466 | 16965256 | void MeshFacetsAABB::nearest_facet_recursive( | |
| 467 | const vec3& p, | ||
| 468 | index_t& nearest_f, vec3& nearest_point, double& sq_dist, | ||
| 469 | index_t n, index_t b, index_t e | ||
| 470 | ) const { | ||
| 471 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 16965256 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
16965256 | geo_debug_assert(e > b); |
| 472 | |||
| 473 | // If node is a leaf: compute point-facet distance | ||
| 474 | // and replace current if nearer | ||
| 475 |
2/2✓ Branch 0 taken 2247907 times.
✓ Branch 1 taken 14717349 times.
|
16965256 | if(b + 1 == e) { |
| 476 |
1/2✓ Branch 1 taken 2247907 times.
✗ Branch 2 not taken.
|
2247907 | index_t f = element_in_leaf(b); |
| 477 | 2247907 | vec3 cur_nearest_point; | |
| 478 | double cur_sq_dist; | ||
| 479 | 2247907 | get_point_facet_nearest_point( | |
| 480 |
1/2✓ Branch 1 taken 2247907 times.
✗ Branch 2 not taken.
|
2247907 | *mesh_, p, f, cur_nearest_point, cur_sq_dist |
| 481 | ); | ||
| 482 |
2/2✓ Branch 0 taken 1578471 times.
✓ Branch 1 taken 669436 times.
|
2247907 | if(cur_sq_dist < sq_dist) { |
| 483 | 1578471 | nearest_f = f; | |
| 484 | 1578471 | nearest_point = cur_nearest_point; | |
| 485 | 1578471 | sq_dist = cur_sq_dist; | |
| 486 | } | ||
| 487 | 2247907 | return; | |
| 488 | } | ||
| 489 | 14717349 | index_t m = b + (e - b) / 2; | |
| 490 | 14717349 | index_t childl = 2 * n; | |
| 491 | 14717349 | index_t childr = 2 * n + 1; | |
| 492 | |||
| 493 | 14717349 | double dl = point_box_signed_squared_distance(p, bboxes_[childl]); | |
| 494 | 14717349 | double dr = point_box_signed_squared_distance(p, bboxes_[childr]); | |
| 495 | |||
| 496 | // Traverse the "nearest" child first, so that it has more chances | ||
| 497 | // to prune the traversal of the other child. | ||
| 498 |
2/2✓ Branch 0 taken 6432094 times.
✓ Branch 1 taken 8285255 times.
|
14717349 | if(dl < dr) { |
| 499 |
2/2✓ Branch 0 taken 5919078 times.
✓ Branch 1 taken 513016 times.
|
6432094 | if(dl < sq_dist) { |
| 500 | 5919078 | nearest_facet_recursive( | |
| 501 | p, | ||
| 502 | nearest_f, nearest_point, sq_dist, | ||
| 503 | childl, b, m | ||
| 504 | ); | ||
| 505 | } | ||
| 506 |
2/2✓ Branch 0 taken 813000 times.
✓ Branch 1 taken 5619094 times.
|
6432094 | if(dr < sq_dist) { |
| 507 | 813000 | nearest_facet_recursive( | |
| 508 | p, | ||
| 509 | nearest_f, nearest_point, sq_dist, | ||
| 510 | childr, m, e | ||
| 511 | ); | ||
| 512 | } | ||
| 513 | } else { | ||
| 514 |
2/2✓ Branch 0 taken 7722642 times.
✓ Branch 1 taken 562613 times.
|
8285255 | if(dr < sq_dist) { |
| 515 | 7722642 | nearest_facet_recursive( | |
| 516 | p, | ||
| 517 | nearest_f, nearest_point, sq_dist, | ||
| 518 | childr, m, e | ||
| 519 | ); | ||
| 520 | } | ||
| 521 |
2/2✓ Branch 0 taken 1709822 times.
✓ Branch 1 taken 6575433 times.
|
8285255 | if(dl < sq_dist) { |
| 522 | 1709822 | nearest_facet_recursive( | |
| 523 | p, | ||
| 524 | nearest_f, nearest_point, sq_dist, | ||
| 525 | childl, b, m | ||
| 526 | ); | ||
| 527 | } | ||
| 528 | } | ||
| 529 | } | ||
| 530 | |||
| 531 | |||
| 532 | ✗ | void MeshFacetsAABB::nearest_facet_recursive_filtered( | |
| 533 | const vec3& p, | ||
| 534 | index_t& nearest_f, vec3& nearest_point, double& sq_dist, | ||
| 535 | index_t n, index_t b, index_t e, | ||
| 536 | std::function<bool(index_t)> filter | ||
| 537 | ) const { | ||
| 538 | ✗ | geo_debug_assert(e > b); | |
| 539 | |||
| 540 | // If node is a leaf: compute point-facet distance | ||
| 541 | // and replace current if nearer | ||
| 542 | ✗ | if(b + 1 == e) { | |
| 543 | ✗ | index_t f = element_in_leaf(b); | |
| 544 | ✗ | if(filter(f)) { | |
| 545 | ✗ | vec3 cur_nearest_point; | |
| 546 | double cur_sq_dist; | ||
| 547 | ✗ | get_point_facet_nearest_point( | |
| 548 | ✗ | *mesh_, p, f, cur_nearest_point, cur_sq_dist | |
| 549 | ); | ||
| 550 | ✗ | if(cur_sq_dist < sq_dist) { | |
| 551 | ✗ | nearest_f = f; | |
| 552 | ✗ | nearest_point = cur_nearest_point; | |
| 553 | ✗ | sq_dist = cur_sq_dist; | |
| 554 | } | ||
| 555 | } | ||
| 556 | ✗ | return; | |
| 557 | } | ||
| 558 | ✗ | index_t m = b + (e - b) / 2; | |
| 559 | ✗ | index_t childl = 2 * n; | |
| 560 | ✗ | index_t childr = 2 * n + 1; | |
| 561 | |||
| 562 | ✗ | double dl = point_box_signed_squared_distance(p, bboxes_[childl]); | |
| 563 | ✗ | double dr = point_box_signed_squared_distance(p, bboxes_[childr]); | |
| 564 | |||
| 565 | // Traverse the "nearest" child first, so that it has more chances | ||
| 566 | // to prune the traversal of the other child. | ||
| 567 | ✗ | if(dl < dr) { | |
| 568 | ✗ | if(dl < sq_dist) { | |
| 569 | ✗ | nearest_facet_recursive_filtered( | |
| 570 | p, | ||
| 571 | nearest_f, nearest_point, sq_dist, | ||
| 572 | childl, b, m, | ||
| 573 | filter | ||
| 574 | ); | ||
| 575 | } | ||
| 576 | ✗ | if(dr < sq_dist) { | |
| 577 | ✗ | nearest_facet_recursive_filtered( | |
| 578 | p, | ||
| 579 | nearest_f, nearest_point, sq_dist, | ||
| 580 | childr, m, e, | ||
| 581 | filter | ||
| 582 | ); | ||
| 583 | } | ||
| 584 | } else { | ||
| 585 | ✗ | if(dr < sq_dist) { | |
| 586 | ✗ | nearest_facet_recursive_filtered( | |
| 587 | p, | ||
| 588 | nearest_f, nearest_point, sq_dist, | ||
| 589 | childr, m, e, | ||
| 590 | filter | ||
| 591 | ); | ||
| 592 | } | ||
| 593 | ✗ | if(dl < sq_dist) { | |
| 594 | ✗ | nearest_facet_recursive_filtered( | |
| 595 | p, | ||
| 596 | nearest_f, nearest_point, sq_dist, | ||
| 597 | childl, b, m, | ||
| 598 | filter | ||
| 599 | ); | ||
| 600 | } | ||
| 601 | } | ||
| 602 | } | ||
| 603 | |||
| 604 | ✗ | bool MeshFacetsAABB::ray_intersection( | |
| 605 | const Ray& R, double tmax, index_t ignore_f | ||
| 606 | ) const { | ||
| 607 | vec3 dirinv( | ||
| 608 | ✗ | 1.0/R.direction.x, | |
| 609 | ✗ | 1.0/R.direction.y, | |
| 610 | ✗ | 1.0/R.direction.z | |
| 611 | ✗ | ); | |
| 612 | ✗ | return ray_intersection_recursive( | |
| 613 | ✗ | R, dirinv, tmax, ignore_f, 1, 0, mesh_->facets.nb() | |
| 614 | ✗ | ); | |
| 615 | } | ||
| 616 | |||
| 617 | 179996 | bool MeshFacetsAABB::ray_nearest_intersection( | |
| 618 | const Ray& R, Intersection& I | ||
| 619 | ) const { | ||
| 620 | 179996 | index_t f = I.f; | |
| 621 | vec3 dirinv( | ||
| 622 | 359992 | 1.0/R.direction.x, | |
| 623 | 359992 | 1.0/R.direction.y, | |
| 624 | 359992 | 1.0/R.direction.z | |
| 625 | 179996 | ); | |
| 626 |
1/2✓ Branch 1 taken 179996 times.
✗ Branch 2 not taken.
|
179996 | ray_nearest_intersection_recursive( |
| 627 | 179996 | R, dirinv, I, f, 1, 0, mesh_->facets.nb(), 0 | |
| 628 | ); | ||
| 629 |
2/2✓ Branch 0 taken 156083 times.
✓ Branch 1 taken 23913 times.
|
179996 | if(I.f != f) { |
| 630 | 156083 | I.p = R.origin + I.t * R.direction; | |
| 631 | 156083 | return true; | |
| 632 | } | ||
| 633 | 23913 | return false; | |
| 634 | } | ||
| 635 | |||
| 636 | ✗ | void MeshFacetsAABB::ray_all_intersections( | |
| 637 | const Ray& R, | ||
| 638 | std::function<void(const Intersection&)> action | ||
| 639 | ) const { | ||
| 640 | vec3 dirinv( | ||
| 641 | ✗ | 1.0/R.direction.x, | |
| 642 | ✗ | 1.0/R.direction.y, | |
| 643 | ✗ | 1.0/R.direction.z | |
| 644 | ✗ | ); | |
| 645 | ✗ | ray_all_intersections_recursive( | |
| 646 | R, dirinv, action, | ||
| 647 | ✗ | 1, 0, mesh_->facets.nb() | |
| 648 | ); | ||
| 649 | ✗ | } | |
| 650 | |||
| 651 | ✗ | void MeshFacetsAABB::line_all_intersections( | |
| 652 | const vec3& O, const vec3& D, | ||
| 653 | std::function<void(const Intersection&)> action | ||
| 654 | ) const { | ||
| 655 | vec3 dirinv( | ||
| 656 | ✗ | 1.0/D.x, | |
| 657 | ✗ | 1.0/D.y, | |
| 658 | ✗ | 1.0/D.z | |
| 659 | ✗ | ); | |
| 660 | ✗ | line_all_intersections_recursive( | |
| 661 | O, D, dirinv, action, | ||
| 662 | ✗ | 1, 0, mesh_->facets.nb() | |
| 663 | ); | ||
| 664 | ✗ | } | |
| 665 | |||
| 666 | ✗ | bool MeshFacetsAABB::ray_intersection_recursive( | |
| 667 | const Ray& R, const vec3& dirinv, double tmax, index_t ignore_f, | ||
| 668 | index_t n, index_t b, index_t e | ||
| 669 | ) const { | ||
| 670 | ✗ | if(!ray_box_intersection(R.origin, dirinv, bboxes_[n], tmax)) { | |
| 671 | ✗ | return false; | |
| 672 | } | ||
| 673 | ✗ | if(b + 1 == e) { | |
| 674 | ✗ | index_t f = element_in_leaf(b); | |
| 675 | ✗ | if(f == ignore_f) { | |
| 676 | ✗ | return false; | |
| 677 | } | ||
| 678 | ✗ | for(auto [ p1, p2, p3] : mesh_->facets.triangle_points(f)) { | |
| 679 | ✗ | vec3 N; | |
| 680 | double t,u,v; | ||
| 681 | ✗ | if( | |
| 682 | ✗ | ray_triangle_intersection( | |
| 683 | ✗ | R.origin, R.direction, p1, p2, p3, t, u, v, N | |
| 684 | ✗ | ) && t < tmax | |
| 685 | ) { | ||
| 686 | ✗ | return true; | |
| 687 | } | ||
| 688 | } | ||
| 689 | ✗ | return false; | |
| 690 | } | ||
| 691 | ✗ | index_t m = b + (e - b) / 2; | |
| 692 | ✗ | index_t childl = 2 * n; | |
| 693 | ✗ | index_t childr = 2 * n + 1; | |
| 694 | return ( | ||
| 695 | ✗ | ray_intersection_recursive(R,dirinv,tmax,ignore_f, childl,b,m) || | |
| 696 | ✗ | ray_intersection_recursive(R,dirinv,tmax,ignore_f, childr,m,e) | |
| 697 | ✗ | ); | |
| 698 | } | ||
| 699 | |||
| 700 | ✗ | bool MeshFacetsAABB::contains(const vec3& p) const { | |
| 701 | ✗ | vec3 D = normalize( | |
| 702 | ✗ | vec3( | |
| 703 | ✗ | Numeric::random_float64(), | |
| 704 | ✗ | Numeric::random_float64(), | |
| 705 | ✗ | Numeric::random_float64() | |
| 706 | ✗ | ) | |
| 707 | ); | ||
| 708 | ✗ | index_t nb_intersections = 0; | |
| 709 | ✗ | ray_all_intersections( | |
| 710 | ✗ | Ray(p, 1e6*D), | |
| 711 | ✗ | [&nb_intersections](const MeshFacetsAABB::Intersection&) { | |
| 712 | ✗ | ++nb_intersections; | |
| 713 | ✗ | } | |
| 714 | ); | ||
| 715 | ✗ | return (nb_intersections & 1) != 0; | |
| 716 | } | ||
| 717 | |||
| 718 | 7851352 | void MeshFacetsAABB::ray_nearest_intersection_recursive( | |
| 719 | const Ray& R, const vec3& dirinv, Intersection& I, index_t ignore_f, | ||
| 720 | index_t n, index_t b, index_t e, index_t coord | ||
| 721 | ) const { | ||
| 722 |
2/2✓ Branch 2 taken 3523966 times.
✓ Branch 3 taken 4327386 times.
|
7851352 | if(!ray_box_intersection(R.origin, dirinv, bboxes_[n], I.t)) { |
| 723 | 3523966 | return; | |
| 724 | } | ||
| 725 |
2/2✓ Branch 0 taken 491708 times.
✓ Branch 1 taken 3835678 times.
|
4327386 | if(b + 1 == e) { |
| 726 | 491708 | index_t f = element_in_leaf(b); | |
| 727 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 491708 times.
|
491708 | if(f == ignore_f) { |
| 728 | ✗ | return; | |
| 729 | } | ||
| 730 |
6/10✓ Branch 1 taken 491708 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 491708 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 491708 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 492438 times.
✗ Branch 11 not taken.
✓ Branch 17 taken 492438 times.
✓ Branch 18 taken 491708 times.
|
984146 | for( auto [ v1, v2, v3] : mesh_->facets.triangles(f)) { |
| 731 |
1/2✓ Branch 1 taken 492438 times.
✗ Branch 2 not taken.
|
492438 | const vec3& p1 = mesh_->vertices.point(v1); |
| 732 |
1/2✓ Branch 1 taken 492438 times.
✗ Branch 2 not taken.
|
492438 | const vec3& p2 = mesh_->vertices.point(v2); |
| 733 |
1/2✓ Branch 1 taken 492438 times.
✗ Branch 2 not taken.
|
492438 | const vec3& p3 = mesh_->vertices.point(v3); |
| 734 | 492438 | vec3 N; | |
| 735 | double t,u,v; | ||
| 736 | 492438 | if( | |
| 737 | 984876 | ray_triangle_intersection( | |
| 738 |
1/2✓ Branch 1 taken 492438 times.
✗ Branch 2 not taken.
|
492438 | R.origin,R.direction,p1,p2,p3,t,u,v,N |
| 739 |
6/6✓ Branch 0 taken 166675 times.
✓ Branch 1 taken 325763 times.
✓ Branch 2 taken 166672 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 166672 times.
✓ Branch 5 taken 325766 times.
|
492438 | ) && t<I.t |
| 740 | ) { | ||
| 741 | 166672 | I.t = t; | |
| 742 | 166672 | I.u = u; | |
| 743 | 166672 | I.v = v; | |
| 744 | 166672 | I.N = N; | |
| 745 | 166672 | I.i = v1; | |
| 746 | 166672 | I.j = v2; | |
| 747 | 166672 | I.k = v3; | |
| 748 | 166672 | I.f = f; | |
| 749 | } | ||
| 750 | } | ||
| 751 | 491708 | return; | |
| 752 | } | ||
| 753 | 3835678 | index_t m = b + (e - b) / 2; | |
| 754 | 3835678 | index_t childl = 2 * n; | |
| 755 | 3835678 | index_t childr = 2 * n + 1; | |
| 756 |
2/2✓ Branch 1 taken 1938080 times.
✓ Branch 2 taken 1897598 times.
|
3835678 | if(dirinv[coord] < 0.0) { |
| 757 | 1938080 | ray_nearest_intersection_recursive( | |
| 758 | 1938080 | R, dirinv, I, ignore_f, childr, m, e, (coord+1)%3 | |
| 759 | ); | ||
| 760 | 1938080 | ray_nearest_intersection_recursive( | |
| 761 | 1938080 | R, dirinv, I, ignore_f, childl, b, m, (coord+1)%3 | |
| 762 | ); | ||
| 763 | } else { | ||
| 764 | 1897598 | ray_nearest_intersection_recursive( | |
| 765 | 1897598 | R, dirinv, I, ignore_f, childl, b, m, (coord+1)%3 | |
| 766 | ); | ||
| 767 | 1897598 | ray_nearest_intersection_recursive( | |
| 768 | 1897598 | R, dirinv, I, ignore_f, childr, m, e, (coord+1)%3 | |
| 769 | ); | ||
| 770 | } | ||
| 771 | } | ||
| 772 | |||
| 773 | |||
| 774 | ✗ | void MeshFacetsAABB::ray_all_intersections_recursive( | |
| 775 | const Ray& R, const vec3& dirinv, | ||
| 776 | std::function<void(const Intersection&)> action, | ||
| 777 | index_t n, index_t b, index_t e | ||
| 778 | ) const { | ||
| 779 | ✗ | Intersection I; | |
| 780 | ✗ | if(!ray_box_intersection(R.origin, dirinv, bboxes_[n], I.t)) { | |
| 781 | ✗ | return; | |
| 782 | } | ||
| 783 | ✗ | if(b + 1 == e) { | |
| 784 | ✗ | index_t f = element_in_leaf(b); | |
| 785 | ✗ | for( auto [ v1, v2, v3] : mesh_->facets.triangles(f)) { | |
| 786 | ✗ | const vec3& p1 = mesh_->vertices.point(v1); | |
| 787 | ✗ | const vec3& p2 = mesh_->vertices.point(v2); | |
| 788 | ✗ | const vec3& p3 = mesh_->vertices.point(v3); | |
| 789 | ✗ | vec3 N; | |
| 790 | ✗ | if( | |
| 791 | ✗ | ray_triangle_intersection( | |
| 792 | ✗ | R.origin,R.direction,p1,p2,p3,I.t,I.u,I.v,I.N | |
| 793 | ) | ||
| 794 | ) { | ||
| 795 | ✗ | I.i = v1; | |
| 796 | ✗ | I.j = v2; | |
| 797 | ✗ | I.k = v3; | |
| 798 | ✗ | I.f = f; | |
| 799 | ✗ | action(I); | |
| 800 | } | ||
| 801 | } | ||
| 802 | ✗ | return; | |
| 803 | } | ||
| 804 | ✗ | index_t m = b + (e - b) / 2; | |
| 805 | ✗ | index_t childl = 2 * n; | |
| 806 | ✗ | index_t childr = 2 * n + 1; | |
| 807 | ✗ | ray_all_intersections_recursive(R, dirinv, action, childr, m, e); | |
| 808 | ✗ | ray_all_intersections_recursive(R, dirinv, action, childl, b, m); | |
| 809 | } | ||
| 810 | |||
| 811 | ✗ | void MeshFacetsAABB::line_all_intersections_recursive( | |
| 812 | const vec3& O, const vec3& D, const vec3& dirinv, | ||
| 813 | std::function<void(const Intersection&)> action, | ||
| 814 | index_t n, index_t b, index_t e | ||
| 815 | ) const { | ||
| 816 | static constexpr bool bidirectional = true; | ||
| 817 | ✗ | Intersection I; | |
| 818 | ✗ | if(!ray_box_intersection(O, dirinv, bboxes_[n], I.t, bidirectional)) { | |
| 819 | ✗ | return; | |
| 820 | } | ||
| 821 | ✗ | if(b + 1 == e) { | |
| 822 | ✗ | index_t f = element_in_leaf(b); | |
| 823 | ✗ | for( auto [ v1, v2, v3] : mesh_->facets.triangles(f)) { | |
| 824 | ✗ | const vec3& p1 = mesh_->vertices.point(v1); | |
| 825 | ✗ | const vec3& p2 = mesh_->vertices.point(v2); | |
| 826 | ✗ | const vec3& p3 = mesh_->vertices.point(v3); | |
| 827 | ✗ | vec3 N; | |
| 828 | ✗ | if( | |
| 829 | ✗ | ray_triangle_intersection( | |
| 830 | O,D,p1,p2,p3,I.t,I.u,I.v,I.N,bidirectional | ||
| 831 | ) | ||
| 832 | ) { | ||
| 833 | ✗ | I.i = v1; | |
| 834 | ✗ | I.j = v2; | |
| 835 | ✗ | I.k = v3; | |
| 836 | ✗ | I.f = f; | |
| 837 | ✗ | action(I); | |
| 838 | } | ||
| 839 | } | ||
| 840 | ✗ | return; | |
| 841 | } | ||
| 842 | ✗ | index_t m = b + (e - b) / 2; | |
| 843 | ✗ | index_t childl = 2 * n; | |
| 844 | ✗ | index_t childr = 2 * n + 1; | |
| 845 | ✗ | line_all_intersections_recursive(O, D, dirinv, action, childr, m, e); | |
| 846 | ✗ | line_all_intersections_recursive(O, D, dirinv, action, childl, b, m); | |
| 847 | } | ||
| 848 | |||
| 849 | 14 | void MeshFacetsAABB::self_bbox_intersections_parallel( | |
| 850 | std::function<void(index_t, index_t)> action, bool concurrent | ||
| 851 | ) const { | ||
| 852 | |||
| 853 | // the parameter of a job that computes the intersection | ||
| 854 | // between two subtrees. | ||
| 855 | struct Job { | ||
| 856 | index_t node1; index_t b1; index_t e1; | ||
| 857 | index_t node2; index_t b2; index_t e2; | ||
| 858 | }; | ||
| 859 | |||
| 860 | // maximum number of facets in a job that will be | ||
| 861 | // run in parallel | ||
| 862 | static constexpr index_t max_job_size = 1024; | ||
| 863 | |||
| 864 |
1/2✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
|
14 | std::stack<Job> S; |
| 865 | 14 | std::vector<Job> jobs; | |
| 866 |
1/2✓ Branch 3 taken 14 times.
✗ Branch 4 not taken.
|
14 | S.push({1,0,mesh_->facets.nb(),1,0,mesh_->facets.nb()}); |
| 867 | |||
| 868 | // De-recursified version of the algorithm in self_intersect_recursive() | ||
| 869 |
2/2✓ Branch 1 taken 4566 times.
✓ Branch 2 taken 14 times.
|
4580 | while(!S.empty()) { |
| 870 | 4566 | Job J = S.top(); | |
| 871 | 4566 | S.pop(); | |
| 872 | |||
| 873 | // Since we are intersecting the AABBTree with *itself*, | ||
| 874 | // we can prune half of the cases by skipping the test | ||
| 875 | // whenever node2's facet index interval is greater than | ||
| 876 | // node1's facet index interval. | ||
| 877 |
2/2✓ Branch 0 taken 160 times.
✓ Branch 1 taken 4406 times.
|
4566 | if(J.e2 <= J.b1) { |
| 878 | 2290 | continue; | |
| 879 | } | ||
| 880 | |||
| 881 | // The acceleration is here: | ||
| 882 | 5098 | if( | |
| 883 |
4/4✓ Branch 0 taken 4072 times.
✓ Branch 1 taken 334 times.
✓ Branch 2 taken 692 times.
✓ Branch 3 taken 3714 times.
|
8478 | (J.node1 != J.node2) && |
| 884 |
4/6✓ Branch 1 taken 4072 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4072 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 692 times.
✓ Branch 8 taken 3380 times.
|
4072 | !bboxes_overlap(bboxes_[J.node1], bboxes_[J.node2]) |
| 885 | ) { | ||
| 886 | 692 | continue; | |
| 887 | } | ||
| 888 | |||
| 889 | // Simple case: leaf - leaf intersection. | ||
| 890 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 3714 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
3714 | if(J.b1 + 1 == J.e1 && J.b2 + 1 == J.e2) { |
| 891 | ✗ | if(J.b1 != J.b2) { | |
| 892 | ✗ | action(element_in_leaf(J.b1), element_in_leaf(J.b2)); | |
| 893 | } | ||
| 894 | ✗ | continue; | |
| 895 | } | ||
| 896 | |||
| 897 | // If job is small enough, push it to the list | ||
| 898 | // of jobs to be run in parallel | ||
| 899 |
4/4✓ Branch 0 taken 2174 times.
✓ Branch 1 taken 1540 times.
✓ Branch 2 taken 1438 times.
✓ Branch 3 taken 736 times.
|
3714 | if(J.e1 - J.b1 <= max_job_size && J.e2 - J.b2 <= max_job_size) { |
| 900 |
1/2✓ Branch 1 taken 1438 times.
✗ Branch 2 not taken.
|
1438 | jobs.push_back(J); |
| 901 | 1438 | continue; | |
| 902 | } | ||
| 903 | |||
| 904 | // If node2 has more elements than node1, then | ||
| 905 | // intersect node2's two children with node1 | ||
| 906 | // else | ||
| 907 | // intersect node1's two children with node2 | ||
| 908 |
2/2✓ Branch 0 taken 1323 times.
✓ Branch 1 taken 953 times.
|
2276 | if(J.e2 - J.b2 > J.e1 - J.b1) { |
| 909 | 1323 | index_t m2 = J.b2 + (J.e2 - J.b2) / 2; | |
| 910 | 1323 | index_t node2_l = 2 * J.node2; | |
| 911 | 1323 | index_t node2_r = 2 * J.node2 + 1; | |
| 912 |
1/2✓ Branch 1 taken 1323 times.
✗ Branch 2 not taken.
|
1323 | S.push({J.node1, J.b1, J.e1, node2_l, J.b2, m2}); |
| 913 |
1/2✓ Branch 1 taken 1323 times.
✗ Branch 2 not taken.
|
1323 | S.push({J.node1, J.b1, J.e1, node2_r, m2, J.e2}); |
| 914 | } else { | ||
| 915 | 953 | index_t m1 = J.b1 + (J.e1 - J.b1) / 2; | |
| 916 | 953 | index_t node1_l = 2 * J.node1; | |
| 917 | 953 | index_t node1_r = 2 * J.node1 + 1; | |
| 918 |
1/2✓ Branch 1 taken 953 times.
✗ Branch 2 not taken.
|
953 | S.push({node1_l, J.b1, m1, J.node2, J.b2, J.e2}); |
| 919 |
1/2✓ Branch 1 taken 953 times.
✗ Branch 2 not taken.
|
953 | S.push({node1_r, m1, J.e1, J.node2, J.b2, J.e2}); |
| 920 | } | ||
| 921 | } | ||
| 922 | |||
| 923 | // Random shuffling avoids configurations with some cores that | ||
| 924 | // have all the hard work to do while other ones finish early | ||
| 925 |
1/2✓ Branch 3 taken 14 times.
✗ Branch 4 not taken.
|
14 | GEO::random_shuffle(jobs.begin(), jobs.end()); |
| 926 | |||
| 927 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 14 times.
|
14 | if(concurrent) { |
| 928 | ✗ | parallel_for( | |
| 929 | ✗ | 0, index_t(jobs.size()), | |
| 930 | ✗ | [&](index_t i) { | |
| 931 | ✗ | const Job& J = jobs[i]; | |
| 932 | ✗ | self_intersect_recursive( | |
| 933 | action, | ||
| 934 | ✗ | J.node1, J.b1, J.e1, | |
| 935 | ✗ | J.node2, J.b2, J.e2 | |
| 936 | ); | ||
| 937 | ✗ | } | |
| 938 | ); | ||
| 939 | } else { | ||
| 940 | // Temporarily memorize intersecting pairs in a per-thread vector | ||
| 941 | // (inserting in same vector generates too much contention) | ||
| 942 | vector<vector<std::pair<index_t,index_t>>> all_candidates( | ||
| 943 | Process::maximum_concurrent_threads() | ||
| 944 |
2/4✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
|
14 | ); |
| 945 |
1/2✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
|
14 | parallel_for_slice( |
| 946 | 14 | 0, index_t(jobs.size()), | |
| 947 |
1/2✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
|
28 | [&](index_t b, index_t e) { |
| 948 |
1/2✓ Branch 1 taken 49 times.
✗ Branch 2 not taken.
|
49 | index_t t = Thread::current()->id(); |
| 949 |
2/2✓ Branch 0 taken 1438 times.
✓ Branch 1 taken 49 times.
|
1487 | for(index_t i=b; i<e; ++i) { |
| 950 | 1438 | const Job& J = jobs[i]; | |
| 951 |
1/2✓ Branch 1 taken 1438 times.
✗ Branch 2 not taken.
|
1438 | self_intersect_recursive( |
| 952 | 1438 | [&](index_t f1, index_t f2) { | |
| 953 |
1/2✓ Branch 3 taken 865317 times.
✗ Branch 4 not taken.
|
865317 | all_candidates[t].push_back({f1,f2}); |
| 954 | 865317 | }, | |
| 955 | 1438 | J.node1, J.b1, J.e1, | |
| 956 | 1438 | J.node2, J.b2, J.e2 | |
| 957 | ); | ||
| 958 | } | ||
| 959 | 49 | } | |
| 960 | ); | ||
| 961 | |||
| 962 | // Call the action for each pair of intersecting boxes | ||
| 963 |
2/2✓ Branch 2 taken 56 times.
✓ Branch 3 taken 14 times.
|
84 | for(const auto& candidates: all_candidates) { |
| 964 |
2/2✓ Branch 2 taken 865317 times.
✓ Branch 3 taken 56 times.
|
865429 | for(const auto& F1F2: candidates) { |
| 965 |
1/2✓ Branch 1 taken 865317 times.
✗ Branch 2 not taken.
|
865317 | action(F1F2.first, F1F2.second); |
| 966 | } | ||
| 967 | } | ||
| 968 | 14 | } | |
| 969 | 14 | } | |
| 970 | |||
| 971 | /****************************************************************************/ | ||
| 972 | |||
| 973 | ✗ | void MeshCellsAABB::initialize(Mesh& M, AABBReorderMode reorder_mode) { | |
| 974 | ✗ | mesh_ = &M; | |
| 975 | ✗ | switch(reorder_mode) { | |
| 976 | ✗ | case AABB_NOREORDER: | |
| 977 | ✗ | break; | |
| 978 | ✗ | case AABB_INPLACE: | |
| 979 | ✗ | reorder_.clear(); | |
| 980 | ✗ | mesh_reorder(*mesh_, MESH_ORDER_MORTON); | |
| 981 | ✗ | break; | |
| 982 | ✗ | case AABB_INDIRECT: | |
| 983 | ✗ | compute_mesh_elements_spatial_order( | |
| 984 | ✗ | *mesh_, MESH_CELLS, reorder_, MESH_ORDER_MORTON | |
| 985 | ); | ||
| 986 | ✗ | break; | |
| 987 | } | ||
| 988 | ✗ | AABB::initialize( | |
| 989 | ✗ | mesh_->cells.nb(), | |
| 990 | ✗ | [this](Box& B, index_t c) { | |
| 991 | // Get cell bbox | ||
| 992 | ✗ | for(coord_index_t coord = 0; coord < 3; ++coord) { | |
| 993 | ✗ | B.xyz_min[coord] = Numeric::max_float64(); | |
| 994 | ✗ | B.xyz_max[coord] = -Numeric::max_float64(); | |
| 995 | } | ||
| 996 | ✗ | for(const vec3& p: mesh_->cells.points(c)) { | |
| 997 | ✗ | for(coord_index_t coord = 0; coord < 3; ++coord) { | |
| 998 | ✗ | B.xyz_min[coord] = std::min( | |
| 999 | ✗ | B.xyz_min[coord], p[coord] | |
| 1000 | ); | ||
| 1001 | ✗ | B.xyz_max[coord] = std::max( | |
| 1002 | ✗ | B.xyz_max[coord], p[coord] | |
| 1003 | ); | ||
| 1004 | } | ||
| 1005 | } | ||
| 1006 | ✗ | } | |
| 1007 | ); | ||
| 1008 | ✗ | } | |
| 1009 | |||
| 1010 | ✗ | index_t MeshCellsAABB::containing_tet_recursive( | |
| 1011 | const vec3& p, | ||
| 1012 | index_t n, index_t b, index_t e | ||
| 1013 | ) const { | ||
| 1014 | |||
| 1015 | ✗ | if(!bboxes_[n].contains(p)) { | |
| 1016 | ✗ | return NO_TET; | |
| 1017 | } | ||
| 1018 | |||
| 1019 | ✗ | if(e==b+1) { | |
| 1020 | ✗ | index_t t = element_in_leaf(b); | |
| 1021 | ✗ | if(mesh_tet_contains_point(*mesh_, t, p)) { | |
| 1022 | ✗ | return t; | |
| 1023 | } else { | ||
| 1024 | ✗ | return NO_TET; | |
| 1025 | } | ||
| 1026 | } | ||
| 1027 | |||
| 1028 | ✗ | index_t m = b + (e - b) / 2; | |
| 1029 | ✗ | index_t childl = 2 * n; | |
| 1030 | ✗ | index_t childr = 2 * n + 1; | |
| 1031 | |||
| 1032 | ✗ | index_t result = containing_tet_recursive( | |
| 1033 | p, childl, b, m | ||
| 1034 | ); | ||
| 1035 | ✗ | if(result == NO_TET) { | |
| 1036 | ✗ | result = containing_tet_recursive(p, childr, m, e); | |
| 1037 | } | ||
| 1038 | ✗ | return result; | |
| 1039 | } | ||
| 1040 | |||
| 1041 | /****************************************************************************/ | ||
| 1042 | |||
| 1043 | ✗ | MeshFacetsAABB2d::MeshFacetsAABB2d() { | |
| 1044 | ✗ | } | |
| 1045 | |||
| 1046 | ✗ | MeshFacetsAABB2d::MeshFacetsAABB2d(Mesh& M, bool reorder) { | |
| 1047 | ✗ | initialize(M, reorder); | |
| 1048 | ✗ | } | |
| 1049 | |||
| 1050 | ✗ | void MeshFacetsAABB2d::initialize(Mesh& M, bool reorder) { | |
| 1051 | ✗ | bool was_2d = (M.vertices.dimension() == 2); | |
| 1052 | ✗ | if(was_2d) { | |
| 1053 | // It is a bit stupid, spatial sort is just implemented | ||
| 1054 | // in 3D for now, so if input mesh was 2D, we temporarily | ||
| 1055 | // create z=0 coordinates for all vertices. | ||
| 1056 | ✗ | M.vertices.set_dimension(3); | |
| 1057 | } | ||
| 1058 | ✗ | mesh_ = &M; | |
| 1059 | ✗ | if(reorder) { | |
| 1060 | ✗ | mesh_reorder(*mesh_, MESH_ORDER_MORTON); | |
| 1061 | } | ||
| 1062 | ✗ | AABB::initialize( | |
| 1063 | ✗ | mesh_->facets.nb(), | |
| 1064 | ✗ | [this](Box2d& B, index_t f) { | |
| 1065 | // Get facet bbox | ||
| 1066 | ✗ | for(coord_index_t coord = 0; coord < 2; ++coord) { | |
| 1067 | ✗ | B.xy_min[coord] = Numeric::max_float64(); | |
| 1068 | ✗ | B.xy_max[coord] = -Numeric::max_float64(); | |
| 1069 | } | ||
| 1070 | ✗ | for(const vec2& p: mesh_->facets.points<2>(f)) { | |
| 1071 | ✗ | for(coord_index_t coord = 0; coord < 2; ++coord) { | |
| 1072 | ✗ | B.xy_min[coord] = std::min( | |
| 1073 | ✗ | B.xy_min[coord], p[coord] | |
| 1074 | ); | ||
| 1075 | ✗ | B.xy_max[coord] = std::max( | |
| 1076 | ✗ | B.xy_max[coord], p[coord] | |
| 1077 | ); | ||
| 1078 | } | ||
| 1079 | } | ||
| 1080 | ✗ | } | |
| 1081 | ); | ||
| 1082 | ✗ | if(was_2d) { | |
| 1083 | ✗ | M.vertices.set_dimension(2); | |
| 1084 | } | ||
| 1085 | ✗ | } | |
| 1086 | |||
| 1087 | ✗ | index_t MeshFacetsAABB2d::containing_triangle_recursive( | |
| 1088 | const vec2& p, | ||
| 1089 | index_t n, index_t b, index_t e | ||
| 1090 | ) const { | ||
| 1091 | |||
| 1092 | ✗ | if(!bboxes_[n].contains(p)) { | |
| 1093 | ✗ | return NO_TRIANGLE; | |
| 1094 | } | ||
| 1095 | |||
| 1096 | ✗ | if(e==b+1) { | |
| 1097 | ✗ | index_t f = element_in_leaf(b); | |
| 1098 | ✗ | if(mesh_triangle_contains_point(*mesh_, f, p)) { | |
| 1099 | ✗ | return b; | |
| 1100 | } else { | ||
| 1101 | ✗ | return NO_TRIANGLE; | |
| 1102 | } | ||
| 1103 | } | ||
| 1104 | |||
| 1105 | ✗ | index_t m = b + (e - b) / 2; | |
| 1106 | ✗ | index_t childl = 2 * n; | |
| 1107 | ✗ | index_t childr = 2 * n + 1; | |
| 1108 | |||
| 1109 | ✗ | index_t result = containing_triangle_recursive( | |
| 1110 | p, childl, b, m | ||
| 1111 | ); | ||
| 1112 | ✗ | if(result == NO_TRIANGLE) { | |
| 1113 | ✗ | result = containing_triangle_recursive(p, childr, m, e); | |
| 1114 | } | ||
| 1115 | ✗ | return result; | |
| 1116 | } | ||
| 1117 | |||
| 1118 | /****************************************************************************/ | ||
| 1119 | |||
| 1120 | } | ||
| 1121 |