GCC Code Coverage Report

Directory: ./
File: mesh/triangle_intersection.cpp
Date: 2026-10-17 02:40:23
Exec Total Coverage
Lines: 313 357 87.7%
Functions: 16 20 80.0%
Branches: 209 283 73.9%
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/triangle_intersection.h>
41 #include <geogram/numerics/predicates.h>
42 #include <geogram/numerics/exact_geometry.h>
43 #include <geogram/basic/string.h>
44 #include <geogram/basic/geometry.h>
45 #include <geogram/basic/argused.h>
46 #include <geogram/basic/logger.h>
47 #include <geogram/basic/algorithm.h>
48
49
50 #ifdef GEO_COMPILER_CLANG
51 #pragma GCC diagnostic ignored "-Wswitch-enum"
52 #endif
53
54 // Uncomment to activate debug messages
55 //#define TT_DEBUG
56
57 namespace {
58
59 using namespace GEO;
60
61 /**
62 * \brief Internal implementation class for
63 * triangles_intersections()
64 * \details Keeps the coordinates of the 2x3
65 * vertices of the triangles and a reference
66 * to the result symbolic information
67 */
68 class TriangleTriangleIntersection {
69 public:
70
71 static constexpr int CACHE_UNINITIALIZED = -2;
72
73 ✗ TriangleTriangleIntersection(
74 const vec3& p0, const vec3& p1, const vec3& p2,
75 const vec3& q0, const vec3& q1, const vec3& q2,
76 TriangleIsects* result = nullptr
77 ✗ ) : result_(result) {
78 ✗ p_[0] = p0;
79 ✗ p_[1] = p1;
80 ✗ p_[2] = p2;
81 ✗ p_[3] = q0;
82 ✗ p_[4] = q1;
83 ✗ p_[5] = q2;
84 ✗ for(index_t i=0; i<64; ++i) {
85 ✗ o3d_cache_[i] = CACHE_UNINITIALIZED;
86 }
87 ✗ has_non_degenerate_intersection_ = false;
88 ✗ for(index_t i=0; i<6; ++i) {
89 ✗ p_index_[i] = NO_INDEX;
90 }
91 ✗ has_global_indices_ = false;
92 ✗ }
93
94 TriangleTriangleIntersection(
95 const vec3& p0, const vec3& p1, const vec3& p2,
96 const vec3& q0, const vec3& q1, const vec3& q2,
97 index_t p0_index,
98 index_t p1_index,
99 index_t p2_index,
100 index_t q0_index,
101 index_t q1_index,
102 index_t q2_index,
103 TriangleIsects* result = nullptr
104
2/2
✓ Branch 0 taken 6044538 times.
✓ Branch 1 taken 1007423 times.
 7051961 ) : result_(result) {
105 1007423 p_[0] = p0;
106 1007423 p_[1] = p1;
107 1007423 p_[2] = p2;
108 1007423 p_[3] = q0;
109 1007423 p_[4] = q1;
110 1007423 p_[5] = q2;
111
2/2
✓ Branch 0 taken 64475072 times.
✓ Branch 1 taken 1007423 times.
 65482495 for(index_t i=0; i<64; ++i) {
112 64475072 o3d_cache_[i] = CACHE_UNINITIALIZED;
113 }
114 1007423 has_non_degenerate_intersection_ = false;
115 1007423 p_index_[0] = p0_index;
116 1007423 p_index_[1] = p1_index;
117 1007423 p_index_[2] = p2_index;
118 1007423 p_index_[3] = q0_index;
119 1007423 p_index_[4] = q1_index;
120 1007423 p_index_[5] = q2_index;
121 1007423 has_global_indices_ = true;
122 }
123
124 1007423 void compute() {
125 #ifdef TT_DEBUG
126 Logger::out("TT") << "Call compute()" << std::endl;
127 #endif
128
129
1/2
✓ Branch 0 taken 1007423 times.
✗ Branch 1 not taken.
 1007423 if(result_ != nullptr) {
130 result_->resize(0);
131 }
132
133 // Test for degenerate triangles
134 if(
135
2/4
✓ Branch 1 taken 1007423 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1007423 times.
 2014846 triangle_dim(T1_RGN_P0, T1_RGN_P1, T1_RGN_P2) != 2 ||
136 1007423 triangle_dim(T2_RGN_P0, T2_RGN_P1, T2_RGN_P2) != 2
137 ) {
138 /*
139 Logger::warn("PCK")
140 << "Tri tri intersect: degenerate triangle "
141 << "(not supported)"
142 << std::endl;
143 */
144 ✗ return;
145 }
146
147
1/2
✓ Branch 0 taken 1007423 times.
✗ Branch 1 not taken.
 1007423 if(has_global_indices_) {
148
149 1007423 index_t q_index[3] = {
150 NO_INDEX, NO_INDEX, NO_INDEX
151 };
152
153
2/2
✓ Branch 0 taken 3022269 times.
✓ Branch 1 taken 1007423 times.
 4029692 for(index_t i=0; i<3; ++i) {
154
2/2
✓ Branch 0 taken 9066807 times.
✓ Branch 1 taken 3022269 times.
 12089076 for(index_t j=0; j<3; ++j) {
155
2/2
✓ Branch 0 taken 911777 times.
✓ Branch 1 taken 8155030 times.
 9066807 if(p_index_[i+3] == p_index_[j]) {
156 911777 q_index[i] = j;
157 }
158 }
159 }
160
161 // Early exit tests for configurations where 1 or 2 vertices
162 // are shared.
163 1007423 int nb_shared = (
164 1007423 (q_index[0] != NO_INDEX) +
165 1007423 (q_index[1] != NO_INDEX) +
166 1007423 (q_index[2] != NO_INDEX) ) ;
167 geo_debug_assert(nb_shared != 3);
168
169 // If there is a single shared vertex, early exit if non-shared
170 // edge [q1,q2] is strictly on one side of [p1,p2,p3]
171
2/2
✓ Branch 0 taken 555741 times.
✓ Branch 1 taken 451682 times.
 1007423 if(nb_shared == 1) {
172 555741 int shared_q =
173
2/2
✓ Branch 0 taken 371758 times.
✓ Branch 1 taken 183983 times.
 555741 (q_index[1] != NO_INDEX) + (q_index[2] != NO_INDEX)*2;
174 geo_debug_assert(q_index[shared_q] != NO_INDEX);
175 555741 TriangleRegion q1 =
176 555741 TriangleRegion(T2_RGN_P0 + ((shared_q+1))%3);
177 555741 TriangleRegion q2 =
178 555741 TriangleRegion(T2_RGN_P0 + ((shared_q+2))%3);
179 TriangleRegion p1,p2,p3;
180 555741 get_triangle_vertices(T1_RGN_T, p1,p2,p3);
181 555741 Sign o1 = orient3d(p1,p2,p3,q1);
182 555741 Sign o2 = orient3d(p1,p2,p3,q2);
183
2/2
✓ Branch 0 taken 460962 times.
✓ Branch 1 taken 94779 times.
 555741 if(int(o1)*int(o2) > 0) {
184 460962 return;
185 }
186 }
187
188 // If there are two shared vertices, early exit if non-shared
189 // vertex q is not on support plane of [p1,p2,p3]
190
2/2
✓ Branch 0 taken 178018 times.
✓ Branch 1 taken 368443 times.
 546461 if(nb_shared == 2) {
191 178018 int non_shared_q =
192
2/2
✓ Branch 0 taken 118103 times.
✓ Branch 1 taken 59915 times.
 178018 (q_index[1] == NO_INDEX) + (q_index[2] == NO_INDEX)*2;
193 geo_debug_assert(q_index[non_shared_q] == NO_INDEX);
194 178018 TriangleRegion q = TriangleRegion(T2_RGN_P0 + non_shared_q);
195 TriangleRegion p1,p2,p3;
196 178018 get_triangle_vertices(T1_RGN_T, p1,p2,p3);
197
2/2
✓ Branch 1 taken 166287 times.
✓ Branch 2 taken 11731 times.
 178018 if(orient3d(p1,p2,p3,q) != ZERO) {
198 166287 return;
199 }
200 }
201 }
202
203 // If T1 is strictly on one side of the supporting
204 // plane of T2, then we are sure there is no intersection
205 // and we can stop there.
206 {
207 TriangleRegion p1,p2,p3;
208 TriangleRegion q1,q2,q3;
209 380174 get_triangle_vertices(T1_RGN_T, p1,p2,p3);
210 380174 get_triangle_vertices(T2_RGN_T, q1,q2,q3);
211 380174 Sign o1 = orient3d(q1,q2,q3,p1);
212 380174 Sign o2 = orient3d(q1,q2,q3,p2);
213 380174 Sign o3 = orient3d(q1,q2,q3,p3);
214 380174 if(
215
2/2
✓ Branch 0 taken 135073 times.
✓ Branch 1 taken 245101 times.
 380174 int(o1)*int(o2) == 1 &&
216
2/2
✓ Branch 0 taken 88498 times.
✓ Branch 1 taken 46575 times.
 135073 int(o2)*int(o3) == 1 &&
217
1/2
✓ Branch 0 taken 88498 times.
✗ Branch 1 not taken.
 88498 int(o3)*int(o1) == 1
218 ) {
219 88498 return;
220 }
221 }
222
223 291676 intersect_edge_triangle(T1_RGN_E0, T2_RGN_T);
224 ✗ if(finished()) { return; }
225 291676 intersect_edge_triangle(T1_RGN_E1, T2_RGN_T);
226 ✗ if(finished()) { return; }
227 291676 intersect_edge_triangle(T1_RGN_E2, T2_RGN_T);
228 ✗ if(finished()) { return; }
229
230 291676 intersect_edge_triangle(T2_RGN_E0, T1_RGN_T);
231 ✗ if(finished()) { return; }
232 291676 intersect_edge_triangle(T2_RGN_E1, T1_RGN_T);
233 ✗ if(finished()) { return; }
234 291676 intersect_edge_triangle(T2_RGN_E2, T1_RGN_T);
235 ✗ if(finished()) { return; }
236
237 // The same intersection can appear several times,
238 // remove the duplicates
239
1/2
✓ Branch 0 taken 291676 times.
✗ Branch 1 not taken.
 291676 if(result_ != nullptr) {
240 std::sort(result_->begin(), result_->end());
241 auto p = std::unique(result_->begin(), result_->end());
242 291676 result_->resize(index_t(p - result_->begin()));
243 // sort_unique(*result_); // HERE
244 #ifdef TT_DEBUG
245 std::string message;
246 for(TriangleIsect I: *result_) {
247 message += (" " + String::to_string(I));
248 }
249 Logger::out("II") << "result: " << message << std::endl;
250 #endif
251 }
252 }
253
254 /**
255 * \brief Tests if there as a non-degenerate intersection
256 * \retval true if an intersection different from two colocated
257 * vertices was found.
258 * \retval false otherwise.
259 */
260 bool has_non_degenerate_intersection() const {
261 1007423 return has_non_degenerate_intersection_;
262 }
263
264 protected:
265
266 /**
267 * \brief Tests whether computation is finished
268 * \details If we just want to know whether there is an intersection
269 * then we can stop sooner.
270 */
271 bool finished() const {
272 #ifdef TT_DEBUG
273 return false;
274 #endif
275
11/34
✗ Branch 0 not taken.
✓ Branch 1 taken 8310 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8310 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 562675 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 562675 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 562675 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 291676 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 291676 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 291676 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 291676 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 291676 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✓ Branch 31 taken 291676 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
 3454701 return (result_ == nullptr && has_non_degenerate_intersection_);
276 }
277
278 1750056 void intersect_edge_triangle(TriangleRegion E, TriangleRegion T) {
279
280 #ifdef TT_DEBUG
281 Logger::out("TT") << std::endl;
282 Logger::out("TT") << " ET " << std::make_pair(E,T) << std::endl;
283 #endif
284
285 geo_debug_assert(region_dim(E) == 1);
286 geo_debug_assert(region_dim(T) == 2);
287
288 TriangleRegion R1 = E;
289 TriangleRegion R2 = T;
290
291 TriangleRegion p1,p2,p3;
292 TriangleRegion e1,e2,e3;
293 TriangleRegion q1,q2;
294
295 1750056 get_triangle_vertices(T,p1,p2,p3);
296 1750056 get_triangle_edges(T,e1,e2,e3);
297 1750056 get_edge_vertices(E,q1,q2);
298
299 1750056 Sign o1 = orient3d(p1,p2,p3,q1);
300 1750056 Sign o2 = orient3d(p1,p2,p3,q2);
301
302 // If both extremities of the segment on same side of triangle
303 // plane, then we are sure there is no intersection and we can
304 // stop there.
305
2/2
✓ Branch 0 taken 1402125 times.
✓ Branch 1 taken 347931 times.
 1750056 if(int(o1) * int(o2) == 1) {
306 645241 return;
307 }
308
309
2/2
✓ Branch 0 taken 562675 times.
✓ Branch 1 taken 839450 times.
 1402125 if(o1 == 0 && o2 == 0) {
310 #ifdef TT_DEBUG
311 Logger::out("TT") << " ET coplanar" << std::endl;
312 #endif
313 // Special case: triangle and segment are co-planar
314 1125350 index_t nax = normal_axis(p1,p2,p3);
315
316 // Test whether the extremities of the segment
317 // are in the triangle
318 {
319 562675 Sign a1 = orient2d(q1,p1,p2,nax);
320 562675 Sign a2 = orient2d(q1,p2,p3,nax);
321 562675 Sign a3 = orient2d(q1,p3,p1,nax);
322
323 562675 Sign b1 = orient2d(q2,p1,p2,nax);
324 562675 Sign b2 = orient2d(q2,p2,p3,nax);
325 562675 Sign b3 = orient2d(q2,p3,p1,nax);
326
327 562675 if(
328
2/2
✓ Branch 0 taken 118485 times.
✓ Branch 1 taken 444190 times.
 562675 int(a1)*int(a2) > 0 &&
329
2/2
✓ Branch 0 taken 8310 times.
✓ Branch 1 taken 110175 times.
 118485 int(a2)*int(a3) > 0 &&
330
1/2
✓ Branch 0 taken 8310 times.
✗ Branch 1 not taken.
 8310 int(a3)*int(a1) > 0
331 ) {
332 8310 add_intersection(q1,T);
333 ✗ if(finished()) { return; }
334 }
335
336 562675 if(
337
2/2
✓ Branch 0 taken 117723 times.
✓ Branch 1 taken 444952 times.
 562675 int(b1)*int(b2) > 0 &&
338
2/2
✓ Branch 0 taken 8310 times.
✓ Branch 1 taken 109413 times.
 117723 int(b2)*int(b3) > 0 &&
339
1/2
✓ Branch 0 taken 8310 times.
✗ Branch 1 not taken.
 8310 int(b3)*int(b1) > 0
340 ) {
341 8310 add_intersection(q2,T);
342 ✗ if(finished()) { return; }
343 }
344 }
345
346 562675 intersect_edge_edge_2d(E,e1,nax);
347 ✗ if(finished()) { return; }
348 562675 intersect_edge_edge_2d(E,e2,nax);
349 ✗ if(finished()) { return; }
350 562675 intersect_edge_edge_2d(E,e3,nax);
351 ✗ if(finished()) { return; }
352
353 } else {
354
355
356 // Update symbolic information of segment
357 // if one of the segment vertices is on
358 // the triangle's supporting plane.
359
360
2/2
✓ Branch 0 taken 254544 times.
✓ Branch 1 taken 584906 times.
 839450 if(o1 == ZERO) {
361 254544 R1 = q1;
362
2/2
✓ Branch 0 taken 254544 times.
✓ Branch 1 taken 330362 times.
 584906 } else if(o2 == ZERO) {
363 254544 R1 = q2;
364 }
365
366 839450 Sign oo1 = orient3d(p2,p3,q1,q2);
367 839450 Sign oo2 = orient3d(p3,p1,q1,q2);
368
369
2/2
✓ Branch 0 taken 542140 times.
✓ Branch 1 taken 297310 times.
 839450 if(int(oo1)*int(oo2) == -1) {
370 return;
371 }
372
373 542140 Sign oo3 = orient3d(p1,p2,q1,q2);
374
375 // Update symbolic information of triangle
376 // if intersection is
377 // on a vertex or on an edge
378 542140 int nb_zeros = (oo1 == ZERO) + (oo2 == ZERO) + (oo3 == ZERO);
379 geo_debug_assert(nb_zeros != 3);
380
2/2
✓ Branch 0 taken 45463 times.
✓ Branch 1 taken 496677 times.
 542140 if(nb_zeros == 1) {
381
2/2
✓ Branch 0 taken 22292 times.
✓ Branch 1 taken 23171 times.
 45463 if(oo1 == ZERO) {
382 22292 R2 = e1;
383
2/2
✓ Branch 0 taken 18941 times.
✓ Branch 1 taken 4230 times.
 23171 } else if(oo2 == ZERO) {
384 18941 R2 = e2;
385 } else {
386 4230 R2 = e3;
387 }
388
2/2
✓ Branch 0 taken 289851 times.
✓ Branch 1 taken 206826 times.
 496677 } else if(nb_zeros == 2) {
389
2/2
✓ Branch 0 taken 92368 times.
✓ Branch 1 taken 197483 times.
 289851 if(oo1 != ZERO) {
390 92368 R2 = p1;
391
2/2
✓ Branch 0 taken 95036 times.
✓ Branch 1 taken 102447 times.
 197483 } else if(oo2 != ZERO) {
392 95036 R2 = p2;
393 } else {
394 102447 R2 = p3;
395 }
396 }
397
398 #ifdef TT_DEBUG
399 Logger::out("TT") << o1 << " " << o2
400 << " "
401 << oo1 << " " << oo2 << " " << oo3
402 << std::endl;
403 #endif
404 // Intersection is outside triangle if oo1, oo2 and oo3
405 // do not have the same sign (or zero)
406 bool outside =
407 542140 (int(oo1) * int(oo2) == -1) ||
408
2/2
✓ Branch 0 taken 357076 times.
✓ Branch 1 taken 185064 times.
 542140 (int(oo2) * int(oo3) == -1) ||
409
2/2
✓ Branch 0 taken 343645 times.
✓ Branch 1 taken 13431 times.
 357076 (int(oo3) * int(oo1) == -1) ;
410
411 if(!outside) {
412 343645 add_intersection(R1,R2);
413 }
414 }
415 }
416
417
418 1688025 void intersect_edge_edge_2d(
419 TriangleRegion E1, TriangleRegion E2, index_t nax
420 ) {
421 #ifdef TT_DEBUG
422 Logger::out("TT") << " EE 2d " << std::make_pair(E1,E2)
423 << " axis: " << nax << std::endl;
424 #endif
425 geo_debug_assert(region_dim(E1) == 1);
426 geo_debug_assert(region_dim(E2) == 1);
427
428 TriangleRegion R1 = E1;
429 TriangleRegion R2 = E2;
430
431 TriangleRegion p1,p2;
432 1688025 get_edge_vertices(E1,p1,p2);
433
434 TriangleRegion q1,q2;
435 1688025 get_edge_vertices(E2,q1,q2);
436
437 1688025 Sign a1 = orient2d(q1,q2,p1,nax);
438 1688025 Sign a2 = orient2d(q1,q2,p2,nax);
439
440
2/2
✓ Branch 0 taken 36294 times.
✓ Branch 1 taken 1651731 times.
 1688025 if(a1 == ZERO && a2 == ZERO) {
441 // Special case: 1D
442 // (the 2x2 edge extremities are aligned)
443 36294 intersect_edge_edge_1d(E1,E2);
444 } else {
445 // Update symbolic information if one of E1's vertices is
446 // on the supporting line of E2
447
2/2
✓ Branch 0 taken 168184 times.
✓ Branch 1 taken 1483547 times.
 1651731 if(a1 == ZERO) {
448 168184 R1 = p1;
449
2/2
✓ Branch 0 taken 168244 times.
✓ Branch 1 taken 1315303 times.
 1483547 } else if(a2 == ZERO) {
450 168244 R1 = p2;
451 }
452
453 1651731 Sign b1 = orient2d(p1,p2,q1,nax);
454 1651731 Sign b2 = orient2d(p1,p2,q2,nax);
455
456 // Update symbolic information if one of E2's vertices is
457 // on the supporting line of E1
458
2/2
✓ Branch 0 taken 177614 times.
✓ Branch 1 taken 1474117 times.
 1651731 if(b1 == ZERO) {
459 177614 R2 = q1;
460
2/2
✓ Branch 0 taken 177614 times.
✓ Branch 1 taken 1296503 times.
 1474117 } else if(b2 == ZERO) {
461 177614 R2 = q2;
462 }
463
464
4/4
✓ Branch 0 taken 512419 times.
✓ Branch 1 taken 1139312 times.
✓ Branch 2 taken 396268 times.
✓ Branch 3 taken 116151 times.
 1651731 if( int(a1)*int(a2) != 1 && int(b1)*int(b2) != 1) {
465 396268 add_intersection(R1,R2);
466 }
467 }
468 1688025 }
469
470 36294 void intersect_edge_edge_1d(
471 TriangleRegion E1, TriangleRegion E2
472 ) {
473 #ifdef TT_DEBUG
474 Logger::out("TT") << " EE 1d " << std::make_pair(E1,E2)
475 << std::endl;
476 #endif
477 geo_debug_assert(region_dim(E1) == 1);
478 geo_debug_assert(region_dim(E2) == 1);
479
480 TriangleRegion p1,p2;
481 36294 get_edge_vertices(E1,p1,p2);
482
483 TriangleRegion q1,q2;
484 36294 get_edge_vertices(E2,q1,q2);
485
486 36294 Sign d1 = dot3d(p1,q1,q2);
487 36294 Sign d2 = dot3d(p2,q1,q2);
488 36294 Sign d3 = dot3d(q1,p1,p2);
489 36294 Sign d4 = dot3d(q2,p1,p2);
490
491 // Test for identical vertices
492 // (small optimization, each time there is an identical pair
493 // of points, the two dot3d predicates they participiate to
494 // should return ZERO, so we can filter)
495
496
2/2
✓ Branch 0 taken 24858 times.
✓ Branch 1 taken 11436 times.
 36294 if(d1 == ZERO && d3 == ZERO && points_are_identical(p1,q1)) {
497 1458 add_intersection(p1,q1);
498 }
499
500
2/2
✓ Branch 0 taken 24859 times.
✓ Branch 1 taken 11435 times.
 36294 if(d2 == ZERO && d3 == ZERO && points_are_identical(p2,q1)) {
501 24797 add_intersection(p2,q1);
502 }
503
504
2/2
✓ Branch 0 taken 24859 times.
✓ Branch 1 taken 11435 times.
 36294 if(d1 == ZERO && d4 == ZERO && points_are_identical(p1,q2)) {
505 24797 add_intersection(p1,q2);
506 }
507
508
2/2
✓ Branch 0 taken 24860 times.
✓ Branch 1 taken 11434 times.
 36294 if(d2 == ZERO && d4 == ZERO && points_are_identical(p2,q2)) {
509 1460 add_intersection(p2,q2);
510 }
511
512 // Test for point in segment:
513 // c is in segment [a,b] if (c-a).(c-b) < 0
514
515
2/2
✓ Branch 0 taken 178 times.
✓ Branch 1 taken 36116 times.
 36294 if(d1 == NEGATIVE) {
516 178 add_intersection(p1,E2);
517 }
518
519
2/2
✓ Branch 0 taken 178 times.
✓ Branch 1 taken 36116 times.
 36294 if(d2 == NEGATIVE) {
520 178 add_intersection(p2,E2);
521 }
522
523
2/2
✓ Branch 0 taken 178 times.
✓ Branch 1 taken 36116 times.
 36294 if(d3 == NEGATIVE) {
524 178 add_intersection(E1,q1);
525 }
526
527
2/2
✓ Branch 0 taken 178 times.
✓ Branch 1 taken 36116 times.
 36294 if(d4 == NEGATIVE) {
528 178 add_intersection(E1,q2);
529 }
530 36294 }
531
532 809757 void add_intersection(TriangleRegion R1, TriangleRegion R2) {
533 #ifdef TT_DEBUG
534 Logger::out("TT") << " ==>I: " << std::make_pair(R1,R2)
535 << std::endl;
536 #endif
537
4/4
✓ Branch 1 taken 688367 times.
✓ Branch 2 taken 121390 times.
✓ Branch 4 taken 31586 times.
✓ Branch 5 taken 656781 times.
 809757 if(region_dim(R1) >= 1 || region_dim(R2) >= 1) {
538 152976 has_non_degenerate_intersection_ = true;
539 }
540 if(is_in_T1(R1)) {
541 geo_debug_assert(!is_in_T1(R2));
542
1/2
✓ Branch 0 taken 392834 times.
✗ Branch 1 not taken.
 392834 if(result_ != nullptr) {
543 392834 result_->push_back(std::make_pair(R1,R2));
544 }
545 } else {
546 geo_debug_assert(is_in_T1(R2));
547
1/2
✓ Branch 0 taken 416923 times.
✗ Branch 1 not taken.
 416923 if(result_ != nullptr) {
548 416923 result_->push_back(std::make_pair(R2,R1));
549 }
550 }
551 809757 }
552
553 8151174 Sign orient3d(
554 TriangleRegion i, TriangleRegion j,
555 TriangleRegion k, TriangleRegion l
556 ) const {
557
558 geo_debug_assert(region_dim(i) == 0);
559 geo_debug_assert(region_dim(j) == 0);
560 geo_debug_assert(region_dim(k) == 0);
561 geo_debug_assert(region_dim(l) == 0);
562
563 // Index for the cache (1 bit set for
564 // each vertex)
565
566 8151174 index_t o3d_idx =
567 8151174 (1u << index_t(i)) |
568 8151174 (1u << index_t(j)) |
569 8151174 (1u << index_t(k)) |
570 8151174 (1u << index_t(l)) ;
571
572 geo_debug_assert(o3d_idx < 64);
573
574 // Result of the predicate should be flipped if the
575 // order of the arguments is permutted by an odd permutation
576
577 8151174 bool flip = odd_order(index_t(i),index_t(j),index_t(k),index_t(l));
578
579 // If cache not initialized, set cache value
580
581
2/2
✓ Branch 0 taken 4598651 times.
✓ Branch 1 taken 3552523 times.
 8151174 if(o3d_cache_[o3d_idx] == CACHE_UNINITIALIZED) {
582
2/2
✓ Branch 0 taken 1822886 times.
✓ Branch 1 taken 2775765 times.
 4598651 int o = flip ? -int(PCK::orient_3d(p_[i],p_[j],p_[k],p_[l])) :
583 4598651 int(PCK::orient_3d(p_[i],p_[j],p_[k],p_[l])) ;
584 4598651 o3d_cache_[o3d_idx] = Numeric::int8(o);
585 }
586
587 // Get result from the cache
588
589
2/2
✓ Branch 0 taken 3904083 times.
✓ Branch 1 taken 4247091 times.
 8151174 Sign result = flip ? Sign(-o3d_cache_[o3d_idx])
590 4247091 : Sign( o3d_cache_[o3d_idx]);
591
592 // Sanity check: did our cache return the same result as
593 // directly calling the predicate ?
594 geo_debug_assert(
595 result == PCK::orient_3d(p_[i], p_[j], p_[k], p_[l])
596 );
597
598 8151174 return result;
599 }
600
601 /**
602 * \brief Tests the parity of the permutation of a list of
603 * four distinct indices with respect to the canonical order.
604 */
605 8151174 static bool odd_order(index_t i, index_t j, index_t k, index_t l) {
606 // Implementation: sort the elements (bubble sort is OK for
607 // such a small number), and invert parity each time
608 // two elements are swapped.
609 8151174 index_t tab[4] = { i, j, k, l };
610 constexpr int N = 4;
611 bool result = false;
612
2/2
✓ Branch 0 taken 24453522 times.
✓ Branch 1 taken 8151174 times.
 32604696 for (int I = 0; I < N - 1; ++I) {
613
2/2
✓ Branch 0 taken 48907044 times.
✓ Branch 1 taken 24453522 times.
 73360566 for (int J = 0; J < N - I - 1; ++J) {
614
2/2
✓ Branch 0 taken 14844289 times.
✓ Branch 1 taken 34062755 times.
 48907044 if (tab[J] > tab[J + 1]) {
615 std::swap(tab[J], tab[J + 1]);
616 14844289 result = !result;
617 }
618 }
619 }
620 8151174 return result;
621 }
622
623 10055562 Sign orient2d(
624 TriangleRegion i, TriangleRegion j, TriangleRegion k,
625 index_t normal_axis
626 ) const {
627 // Note: no cache for orient2d (tested, did not bring
628 // any performance gain).
629 geo_debug_assert(region_dim(i) == 0);
630 geo_debug_assert(region_dim(j) == 0);
631 geo_debug_assert(region_dim(k) == 0);
632 double pi[2];
633 double pj[2];
634 double pk[2];
635
2/2
✓ Branch 0 taken 20111124 times.
✓ Branch 1 taken 10055562 times.
 30166686 for(coord_index_t c = 0; c < 2; c++) {
636 20111124 pi[c] = p_[i][index_t((normal_axis + 1 + c) % 3)];
637 20111124 pj[c] = p_[j][index_t((normal_axis + 1 + c) % 3)];
638 20111124 pk[c] = p_[k][index_t((normal_axis + 1 + c) % 3)];
639 }
640 10055562 return PCK::orient_2d(pi,pj,pk);
641 }
642
643 /**
644 * \brief Computes the dot product between two vectors supported
645 * by three points.
646 * \param[in] i , j , k the three points
647 * \return sign(dot(pj-pi,pk-pi))
648 */
649 Sign dot3d(
650 TriangleRegion i, TriangleRegion j, TriangleRegion k
651 ) const {
652 geo_debug_assert(region_dim(i) == 0);
653 geo_debug_assert(region_dim(j) == 0);
654 geo_debug_assert(region_dim(k) == 0);
655 36294 return PCK::dot_3d(p_[i].data(), p_[j].data(), p_[k].data());
656 }
657
658
659 /**
660 * \brief Computes the coordinate along which a triangle can be
661 * projected without introducting degeneracies.
662 * \param[in] v1 , v2 , v3 the three vertices of the triangle
663 * \return the coordinate to be used for 2d computations (0,1 or 2)
664 */
665 coord_index_t normal_axis(
666 TriangleRegion v1, TriangleRegion v2, TriangleRegion v3
667 ) {
668 geo_debug_assert(region_dim(v1) == 0);
669 geo_debug_assert(region_dim(v2) == 0);
670 geo_debug_assert(region_dim(v3) == 0);
671
672 562675 const vec3& p1 = p_[v1];
673 562675 const vec3& p2 = p_[v2];
674 562675 const vec3& p3 = p_[v3];
675
676 562675 return PCK::triangle_normal_axis(p1,p2,p3);
677 }
678
679
680 bool points_are_identical(TriangleRegion i, TriangleRegion j) {
681 geo_debug_assert(region_dim(i) == 0);
682 geo_debug_assert(region_dim(j) == 0);
683 99436 const vec3& p1 = p_[i];
684 99436 const vec3& p2 = p_[j];
685 return
686 161068 (p1[0] == p2[0]) &&
687
16/24
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 5992 times.
✓ Branch 9 taken 18866 times.
✓ Branch 10 taken 2414 times.
✓ Branch 11 taken 3578 times.
✓ Branch 12 taken 24823 times.
✓ Branch 13 taken 36 times.
✓ Branch 14 taken 24801 times.
✓ Branch 15 taken 22 times.
✓ Branch 16 taken 24823 times.
✓ Branch 17 taken 36 times.
✓ Branch 18 taken 24801 times.
✓ Branch 19 taken 22 times.
✓ Branch 20 taken 5994 times.
✓ Branch 21 taken 18866 times.
✓ Branch 22 taken 2416 times.
✓ Branch 23 taken 3578 times.
 99436 (p1[1] == p2[1]) &&
688
8/12
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1458 times.
✓ Branch 5 taken 956 times.
✓ Branch 6 taken 24797 times.
✓ Branch 7 taken 4 times.
✓ Branch 8 taken 24797 times.
✓ Branch 9 taken 4 times.
✓ Branch 10 taken 1460 times.
✓ Branch 11 taken 956 times.
 54432 (p1[2] == p2[2]) ;
689 }
690
691 /**
692 * \brief Detects degenerate triangles
693 * \retval 0 if all the vertices of the triangle are the same point
694 * \retval 1 if the vertices of the triangle are co-linear
695 * \retval 2 otherwise
696 */
697 2014846 index_t triangle_dim(
698 TriangleRegion i, TriangleRegion j, TriangleRegion k
699 ) {
700 geo_debug_assert(region_dim(i) == 0);
701 geo_debug_assert(region_dim(j) == 0);
702 geo_debug_assert(region_dim(k) == 0);
703
704 2014846 const vec3& p1 = p_[i];
705 2014846 const vec3& p2 = p_[j];
706 2014846 const vec3& p3 = p_[k];
707
708
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2014846 times.
 2014846 if(!PCK::aligned_3d(p1.data(), p2.data(), p3.data())) {
709 return 2;
710 }
711 if(points_are_identical(i,j) && points_are_identical(j,k)) {
712 ✗ return 0;
713 }
714 return 1;
715 }
716
717
718 private:
719 vec3 p_[6];
720 TriangleIsects* result_;
721 bool has_non_degenerate_intersection_;
722 mutable Numeric::int8 o3d_cache_[64];
723 index_t p_index_[6];
724 bool has_global_indices_;
725 };
726
727 }
728
729 /****************************************************************************/
730
731 namespace GEO {
732
733 ✗ std::string region_to_string(TriangleRegion rgn) {
734 ✗ const char* strs[T_RGN_NB] = {
735 "T1.P0",
736 "T1.P1",
737 "T1.P2",
738
739 "T2.P0",
740 "T2.P1",
741 "T2.P2",
742
743 "T1.E0",
744 "T1.E1",
745 "T1.E2",
746
747 "T2.E0",
748 "T2.E1",
749 "T2.E2",
750
751 "T1.T",
752 "T2.T"
753 };
754 ✗ geo_assert(int(rgn) < int(T_RGN_NB));
755 ✗ return strs[int(rgn)];
756 }
757
758 // This version returns the symbolic information.
759 ✗ bool triangles_intersections(
760 const vec3& p0, const vec3& p1, const vec3& p2,
761 const vec3& q0, const vec3& q1, const vec3& q2,
762 TriangleIsects& result
763 ) {
764 result.resize(0);
765 TriangleTriangleIntersection I(
766 p0, p1, p2,
767 q0, q1, q2,
768 &result
769 ✗ );
770 ✗ I.compute();
771 ✗ return I.has_non_degenerate_intersection();
772 }
773
774 // This version returns the symbolic information.
775 1007423 bool triangles_intersections(
776 const vec3& p0, const vec3& p1, const vec3& p2,
777 const vec3& q0, const vec3& q1, const vec3& q2,
778 index_t p0_index, index_t p1_index, index_t p2_index,
779 index_t q0_index, index_t q1_index, index_t q2_index,
780 TriangleIsects& result
781 ) {
782 result.resize(0);
783 TriangleTriangleIntersection I(
784 p0, p1, p2,
785 q0, q1, q2,
786 p0_index, p1_index, p2_index,
787 q0_index, q1_index, q2_index,
788 &result
789 );
790 1007423 I.compute();
791 1007423 return I.has_non_degenerate_intersection();
792 }
793
794 // This version is just a predicate (returns true if
795 // there is a non-degenerate intersection, false
796 // otherwise).
797 ✗ bool triangles_intersections(
798 const vec3& p0, const vec3& p1, const vec3& p2,
799 const vec3& q0, const vec3& q1, const vec3& q2
800 ) {
801 TriangleTriangleIntersection I(
802 p0, p1, p2,
803 q0, q1, q2
804 ✗ );
805 ✗ I.compute();
806 ✗ return I.has_non_degenerate_intersection();
807 }
808
809 4899296 coord_index_t region_dim(TriangleRegion r) {
810 geo_debug_assert(index_t(r) < T_RGN_NB);
811 static coord_index_t trgl_rgn_dim[T_RGN_NB] = {
812 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2
813 };
814 4899296 return trgl_rgn_dim[index_t(r)];
815 }
816
817 271532 TriangleRegion swap_T1_T2(TriangleRegion R) {
818 TriangleRegion result=T_RGN_NB;
819
14/16
✓ Branch 0 taken 4432 times.
✓ Branch 1 taken 3754 times.
✓ Branch 2 taken 4741 times.
✓ Branch 3 taken 5947 times.
✓ Branch 4 taken 3630 times.
✓ Branch 5 taken 3986 times.
✓ Branch 6 taken 32363 times.
✓ Branch 7 taken 28804 times.
✓ Branch 8 taken 35706 times.
✓ Branch 9 taken 32857 times.
✓ Branch 10 taken 28585 times.
✓ Branch 11 taken 36297 times.
✓ Branch 12 taken 25966 times.
✓ Branch 13 taken 24464 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 271532 switch(R) {
820 4432 case T1_RGN_P0:
821 result = T2_RGN_P0;
822 4432 break;
823 3754 case T1_RGN_P1:
824 result = T2_RGN_P1;
825 3754 break;
826 4741 case T1_RGN_P2:
827 result = T2_RGN_P2;
828 4741 break;
829 5947 case T2_RGN_P0:
830 result = T1_RGN_P0;
831 5947 break;
832 3630 case T2_RGN_P1:
833 result = T1_RGN_P1;
834 3630 break;
835 3986 case T2_RGN_P2:
836 result = T1_RGN_P2;
837 3986 break;
838 32363 case T1_RGN_E0:
839 result = T2_RGN_E0;
840 32363 break;
841 28804 case T1_RGN_E1:
842 result = T2_RGN_E1;
843 28804 break;
844 35706 case T1_RGN_E2:
845 result = T2_RGN_E2;
846 35706 break;
847 32857 case T2_RGN_E0:
848 result = T1_RGN_E0;
849 32857 break;
850 28585 case T2_RGN_E1:
851 result = T1_RGN_E1;
852 28585 break;
853 36297 case T2_RGN_E2:
854 result = T1_RGN_E2;
855 36297 break;
856 25966 case T1_RGN_T:
857 result = T2_RGN_T;
858 25966 break;
859 24464 case T2_RGN_T:
860 result = T1_RGN_T;
861 24464 break;
862 case T_RGN_NB:
863 ✗ geo_assert_not_reached;
864 }
865 271532 return result;
866 }
867
868 3244163 void get_triangle_vertices(
869 TriangleRegion T,
870 TriangleRegion& p0, TriangleRegion& p1, TriangleRegion& p2
871 ) {
872 geo_debug_assert(region_dim(T) == 2);
873
2/3
✓ Branch 0 taken 1988961 times.
✓ Branch 1 taken 1255202 times.
✗ Branch 2 not taken.
 3244163 switch(T) {
874 1988961 case T1_RGN_T: {
875 1988961 p0 = T1_RGN_P0; p1 = T1_RGN_P1; p2 = T1_RGN_P2;
876 1988961 } break;
877 1255202 case T2_RGN_T: {
878 1255202 p0 = T2_RGN_P0; p1 = T2_RGN_P1; p2 = T2_RGN_P2;
879 1255202 } break;
880 default:
881 ✗ geo_assert_not_reached;
882 }
883 3244163 }
884
885 1750056 void get_triangle_edges(
886 TriangleRegion T,
887 TriangleRegion& e0, TriangleRegion& e1, TriangleRegion& e2
888 ) {
889 geo_debug_assert(region_dim(T) == 2);
890
2/3
✓ Branch 0 taken 875028 times.
✓ Branch 1 taken 875028 times.
✗ Branch 2 not taken.
 1750056 switch(T) {
891 875028 case T1_RGN_T: {
892 875028 e0 = T1_RGN_E0; e1 = T1_RGN_E1; e2 = T1_RGN_E2;
893 875028 } break;
894 875028 case T2_RGN_T: {
895 875028 e0 = T2_RGN_E0; e1 = T2_RGN_E1; e2 = T2_RGN_E2;
896 875028 } break;
897 default:
898 ✗ geo_assert_not_reached;
899 }
900 1750056 }
901
902 5264863 void get_edge_vertices(
903 TriangleRegion E, TriangleRegion& q0, TriangleRegion& q1
904 ) {
905 geo_debug_assert(region_dim(E) == 1);
906
6/7
✓ Branch 0 taken 887421 times.
✓ Branch 1 taken 875660 times.
✓ Branch 2 taken 877623 times.
✓ Branch 3 taken 872769 times.
✓ Branch 4 taken 879037 times.
✓ Branch 5 taken 872353 times.
✗ Branch 6 not taken.
 5264863 switch(E) {
907 887421 case T1_RGN_E0: {
908 887421 q0 = T1_RGN_P1;
909 887421 q1 = T1_RGN_P2;
910 887421 } break;
911 875660 case T1_RGN_E1: {
912 875660 q0 = T1_RGN_P2;
913 875660 q1 = T1_RGN_P0;
914 875660 } break;
915 877623 case T1_RGN_E2: {
916 877623 q0 = T1_RGN_P0;
917 877623 q1 = T1_RGN_P1;
918 877623 } break;
919 872769 case T2_RGN_E0: {
920 872769 q0 = T2_RGN_P1;
921 872769 q1 = T2_RGN_P2;
922 872769 } break;
923 879037 case T2_RGN_E1: {
924 879037 q0 = T2_RGN_P2;
925 879037 q1 = T2_RGN_P0;
926 879037 } break;
927 872353 case T2_RGN_E2: {
928 872353 q0 = T2_RGN_P0;
929 872353 q1 = T2_RGN_P1;
930 872353 } break;
931 default:
932 ✗ geo_assert_not_reached;
933 }
934 5264863 }
935
936 328516 TriangleRegion GEOGRAM_API regions_convex_hull(
937 TriangleRegion R1, TriangleRegion R2
938 ) {
939 geo_debug_assert(is_in_T1(R1) == is_in_T1(R2));
940
2/2
✓ Branch 0 taken 206293 times.
✓ Branch 1 taken 122223 times.
 328516 if(R1 == R2) {
941 return R1;
942 }
943
944 TriangleRegion R = is_in_T1(R1) ? T1_RGN_T : T2_RGN_T;
945
946
4/4
✓ Branch 1 taken 160530 times.
✓ Branch 2 taken 45763 times.
✓ Branch 4 taken 143355 times.
✓ Branch 5 taken 17175 times.
 206293 if(region_dim(R1) == 1 && region_dim(R2) == 0) {
947 TriangleRegion v1,v2;
948 17175 get_edge_vertices(R1,v1,v2);
949
4/4
✓ Branch 0 taken 9493 times.
✓ Branch 1 taken 7682 times.
✓ Branch 2 taken 7584 times.
✓ Branch 3 taken 1909 times.
 17175 if(R2 == v1 || R2 == v2) {
950 R = R1;
951 }
952
4/4
✓ Branch 1 taken 155144 times.
✓ Branch 2 taken 33974 times.
✓ Branch 4 taken 140519 times.
✓ Branch 5 taken 14625 times.
 189118 } else if(region_dim(R2) == 1 && region_dim(R1) == 0) {
953 TriangleRegion v1,v2;
954 14625 get_edge_vertices(R2,v1,v2);
955
4/4
✓ Branch 0 taken 8336 times.
✓ Branch 1 taken 6289 times.
✓ Branch 2 taken 6389 times.
✓ Branch 3 taken 1947 times.
 14625 if(R1 == v1 || R1 == v2) {
956 R = R2;
957 }
958
4/4
✓ Branch 1 taken 10419 times.
✓ Branch 2 taken 164074 times.
✓ Branch 4 taken 10237 times.
✓ Branch 5 taken 182 times.
 174493 } else if(region_dim(R1) == 0 && region_dim(R2) == 0) {
959 24132 for(TriangleRegion E: {
960 T1_RGN_E0, T1_RGN_E1, T1_RGN_E2,
961 T2_RGN_E0, T2_RGN_E1, T2_RGN_E2
962
1/2
✓ Branch 0 taken 34369 times.
✗ Branch 1 not taken.
 34369 }
963 ) {
964 TriangleRegion v1,v2;
965 34369 get_edge_vertices(E,v1,v2);
966
8/8
✓ Branch 0 taken 6742 times.
✓ Branch 1 taken 27627 times.
✓ Branch 2 taken 2336 times.
✓ Branch 3 taken 4406 times.
✓ Branch 4 taken 8929 times.
✓ Branch 5 taken 21034 times.
✓ Branch 6 taken 5831 times.
✓ Branch 7 taken 3098 times.
 34369 if((R1 == v1 && R2 == v2) || (R1 == v2 && R2 == v1)) {
967 R = E;
968 10237 break;
969 }
970 }
971 }
972
973 return R;
974 }
975
976 }
977