GCC Code Coverage Report

Directory: ./
File: mesh/mesh_surface_intersection_internal.cpp
Date: 2026-09-27 03:10:11
Exec Total Coverage
Lines: 413 510 81.0%
Functions: 27 32 84.4%
Branches: 263 530 49.6%
Line Branch Exec Source
1 /*
2 * COPYRIGHT (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #include <geogram/mesh/mesh_surface_intersection_internal.h>
41 #include <geogram/mesh/mesh_surface_intersection.h>
42 #include <geogram/basic/debug_stream.h>
43 #include <geogram/basic/boolean_expression.h>
44 #include <stack>
45
46 namespace {
47 using namespace GEO;
48
49 /**
50 * \brief Computes the exact intersection between the support
51 * planes of three triangles
52 * \param[in] p1 , p2 , p3 the three vertices of the first triangle
53 * \param[in] q1 , q2 , q3 the three vertices of the second triangle
54 * \param[in] r1 , r2 , r3 the three vertices of the third triangle
55 * \param[out] result the exact intersection between the three planes
56 * if it exsists
57 * \retval true if the planes have an intersection
58 * \retval false otherwise
59 */
60 5684 bool get_three_planes_intersection(
61 MeshSurfaceIntersection::ExactPoint& result,
62 const vec3& p1, const vec3& p2, const vec3& p3,
63 const vec3& q1, const vec3& q2, const vec3& q3,
64 const vec3& r1, const vec3& r2, const vec3& r3
65 ) {
66 5684 exact::vec3 N1 = triangle_normal<exact::vec3>(p1,p2,p3);
67
1/2
✓ Branch 1 taken 5684 times.
✗ Branch 2 not taken.
 5684 exact::vec3 N2 = triangle_normal<exact::vec3>(q1,q2,q3);
68
1/2
✓ Branch 1 taken 5684 times.
✗ Branch 2 not taken.
 5684 exact::vec3 N3 = triangle_normal<exact::vec3>(r1,r2,r3);
69
70 exact::vec3 B(
71
2/6
✓ Branch 2 taken 5684 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 5684 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 11368 dot(N1,exact::vec3(p1)),
72
2/6
✓ Branch 2 taken 5684 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 5684 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 11368 dot(N2,exact::vec3(q1)),
73
2/4
✓ Branch 2 taken 5684 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 5684 times.
✗ Branch 7 not taken.
 11368 dot(N3,exact::vec3(r1))
74 5684 );
75
76
1/2
✓ Branch 0 taken 5684 times.
✗ Branch 1 not taken.
 5684 result.w = det3x3(
77 N1.x, N1.y, N1.z,
78 N2.x, N2.y, N2.z,
79 N3.x, N3.y, N3.z
80 );
81
82
2/2
✓ Branch 0 taken 1327 times.
✓ Branch 1 taken 4357 times.
 5684 if(result.w.sign() == ZERO) {
83 return false;
84 }
85
86
1/2
✓ Branch 1 taken 1327 times.
✗ Branch 2 not taken.
 1327 result.x = det3x3(
87 B.x, N1.y, N1.z,
88 B.y, N2.y, N2.z,
89 B.z, N3.y, N3.z
90 );
91
92
1/2
✓ Branch 1 taken 1327 times.
✗ Branch 2 not taken.
 1327 result.y = det3x3(
93 N1.x, B.x, N1.z,
94 N2.x, B.y, N2.z,
95 N3.x, B.z, N3.z
96 );
97
98 1327 result.z = det3x3(
99 N1.x, N1.y, B.x,
100 N2.x, N2.y, B.y,
101 N3.x, N3.y, B.z
102 );
103
104 1327 return true;
105 5684 }
106
107 /**
108 * \brief Computes the exact intersection between the support plane
109 * of a triangle and the support line of a segment
110 * \pre The intersection exists
111 * \param[in] p1 , p2 , p3 the three vertices of the triangle
112 * \param[in] q1 , q2 the two vertices of the segment
113 * \return the exact intersection between the plane and the line
114 */
115 67563 MeshSurfaceIntersection::ExactPoint plane_line_intersection(
116 const vec3& p1, const vec3& p2, const vec3& p3,
117 const vec3& q1, const vec3& q2
118 ) {
119 // Moller & Trumbore's algorithm
120 // see: https://stackoverflow.com/questions/42740765/
121 // intersection-between-line-and-triangle-in-3d
122 67563 exact::vec3 D = make_vec3<exact::vec3>(q1,q2);
123 67563 exact::vec3 E1 = make_vec3<exact::vec3>(p1,p2);
124 67563 exact::vec3 E2 = make_vec3<exact::vec3>(p1,p3);
125 67563 exact::vec3 AO = make_vec3<exact::vec3>(p1,q1);
126
1/2
✓ Branch 1 taken 67563 times.
✗ Branch 2 not taken.
 67563 exact::vec3 N = cross(E1,E2);
127
2/4
✓ Branch 1 taken 67563 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 67563 times.
✗ Branch 5 not taken.
 67563 exact::scalar d = -dot(D,N);
128 geo_debug_assert(d.sign() != ZERO);
129
1/4
✓ Branch 1 taken 67563 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 67563 exact::rational t(dot(AO,N),d);
130
1/2
✓ Branch 1 taken 67563 times.
✗ Branch 2 not taken.
 135126 return mix(t,q1,q2);
131 135126 }
132 }
133
134 namespace GEO {
135
136 ✗ void MeshInTriangle::Vertex::print(std::ostream& out) const {
137 ✗ if(sym.f1 != NO_INDEX) {
138 ✗ out << " ( ";
139 ✗ out << sym.f1;
140 ✗ out << region_to_string(sym.R1).substr(2);
141 }
142 ✗ if(sym.f2 != NO_INDEX) {
143 ✗ out << " /\\ ";
144 ✗ out << sym.f2;
145 ✗ out << region_to_string(sym.R2).substr(2);
146 }
147 ✗ if(sym.f1 != NO_INDEX) {
148 ✗ out << " ) ";
149 }
150 ✗ }
151
152 286659 MeshInTriangle::ExactPoint MeshInTriangle::Vertex::compute_geometry() {
153 // Case 1: f1 vertex
154
2/2
✓ Branch 1 taken 46518 times.
✓ Branch 2 taken 240141 times.
 286659 if(region_dim(sym.R1) == 0) {
155
1/2
✓ Branch 0 taken 46518 times.
✗ Branch 1 not taken.
 46518 index_t lv = index_t(sym.R1) - index_t(T1_RGN_P0);
156 geo_debug_assert(lv < 3);
157
1/2
✓ Branch 0 taken 46518 times.
✗ Branch 1 not taken.
 93036 mesh_vertex_index = mesh().facets.vertex(sym.f1,lv);
158 vec3 p = mit->mesh_vertex(mesh_vertex_index);
159 46518 return ExactPoint(p);
160 }
161
162 geo_debug_assert(sym.f1 != NO_INDEX && sym.f2 != NO_INDEX);
163
164 // Case 2: f2 vertex
165
2/2
✓ Branch 1 taken 21655 times.
✓ Branch 2 taken 218486 times.
 240141 if(region_dim(sym.R2) == 0) {
166
1/2
✓ Branch 0 taken 21655 times.
✗ Branch 1 not taken.
 21655 index_t lv = index_t(sym.R2) - index_t(T2_RGN_P0);
167 geo_debug_assert(lv < 3);
168
1/2
✓ Branch 0 taken 21655 times.
✗ Branch 1 not taken.
 43310 mesh_vertex_index = mesh().facets.vertex(sym.f2, lv);
169 vec3 p = mit->mesh_vertex(mesh_vertex_index);
170 21655 return ExactPoint(p);
171 }
172
173 // case 3: f1 /\ f2 edge in 3D or f1 edge /\ f2 edge in 3D
174 if(
175
5/6
✓ Branch 1 taken 191064 times.
✓ Branch 2 taken 27422 times.
✓ Branch 4 taken 191064 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 191064 times.
✓ Branch 7 taken 27422 times.
 436972 (region_dim(sym.R1) == 2 || region_dim(sym.R1) == 1) &&
176 218486 region_dim(sym.R2) == 1
177 ) {
178 191064 vec3 p1 = mit->mesh_facet_vertex(sym.f1, 0);
179 191064 vec3 p2 = mit->mesh_facet_vertex(sym.f1, 1);
180 191064 vec3 p3 = mit->mesh_facet_vertex(sym.f1, 2);
181 191064 index_t e = index_t(sym.R2)-index_t(T2_RGN_E0);
182 geo_debug_assert(e<3);
183 191064 vec3 q1 = mit->mesh_facet_vertex(sym.f2, (e+1)%3);
184 191064 vec3 q2 = mit->mesh_facet_vertex(sym.f2, (e+2)%3);
185
186 bool seg_seg_two_D = (
187
2/2
✓ Branch 1 taken 150982 times.
✓ Branch 2 taken 12660 times.
 354706 region_dim(sym.R1) == 1 &&
188
3/4
✓ Branch 0 taken 163642 times.
✓ Branch 1 taken 27422 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 150982 times.
 342046 PCK::orient_3d(p1,p2,p3,q1) == ZERO &&
189 PCK::orient_3d(p1,p2,p3,q2) == ZERO) ;
190
191 if(!seg_seg_two_D) {
192 40082 return plane_line_intersection(p1,p2,p3,q1,q2);
193 }
194 }
195
196 // case 4: f1 edge /\ f2
197
3/4
✓ Branch 1 taken 178404 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 150982 times.
✓ Branch 5 taken 27422 times.
 178404 if(region_dim(sym.R1) == 1 && region_dim(sym.R2) == 2) {
198 27422 index_t e = index_t(sym.R1)-index_t(T1_RGN_E0);
199 geo_debug_assert(e<3);
200 27422 vec3 p1 = mit->mesh_facet_vertex(sym.f2,0);
201 27422 vec3 p2 = mit->mesh_facet_vertex(sym.f2,1);
202 27422 vec3 p3 = mit->mesh_facet_vertex(sym.f2,2);
203 27422 vec3 q1 = mit->mesh_facet_vertex(sym.f1, (e+1)%3);
204 27422 vec3 q2 = mit->mesh_facet_vertex(sym.f1, (e+2)%3);
205 27422 return plane_line_intersection(p1,p2,p3,q1,q2);
206 }
207
208 // case 5: f1 edge /\ f2 edge in 2D
209
2/4
✓ Branch 1 taken 150982 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 150982 times.
 150982 if(region_dim(sym.R1) == 1 && region_dim(sym.R2) == 1) {
210 150982 index_t e1 = index_t(sym.R1) - index_t(T1_RGN_E0);
211 geo_debug_assert(e1 < 3);
212 150982 index_t e2 = index_t(sym.R2) - index_t(T2_RGN_E0);
213 geo_debug_assert(e2 < 3);
214 150982 vec2 p1 = mit->mesh_facet_vertex_UV(sym.f1, (e1+1)%3);
215 150982 vec2 p2 = mit->mesh_facet_vertex_UV(sym.f1, (e1+2)%3);
216 150982 vec2 q1 = mit->mesh_facet_vertex_UV(sym.f2, (e2+1)%3);
217 150982 vec2 q2 = mit->mesh_facet_vertex_UV(sym.f2, (e2+2)%3);
218 150982 vec3 P1 = mit->mesh_facet_vertex(sym.f1, (e1+1)%3);
219 150982 vec3 P2 = mit->mesh_facet_vertex(sym.f1, (e1+2)%3);
220
221 150982 exact::vec2 D1 = make_vec2<exact::vec2>(p1,p2);
222 150982 exact::vec2 D2 = make_vec2<exact::vec2>(q1,q2);
223
1/2
✓ Branch 1 taken 150982 times.
✗ Branch 2 not taken.
 150982 exact::scalar d = det(D1,D2);
224 geo_debug_assert(d.sign() != ZERO);
225 150982 exact::vec2 AO = make_vec2<exact::vec2>(p1,q1);
226
1/2
✓ Branch 1 taken 150982 times.
✗ Branch 2 not taken.
 150982 exact::rational t(det(AO,D2),d);
227
1/2
✓ Branch 1 taken 150982 times.
✗ Branch 2 not taken.
 150982 return mix(t,P1,P2);
228 301964 }
229
230 // Normally we enumerated all possible cases
231 ✗ geo_assert_not_reached;
232 }
233
234 292343 void MeshInTriangle::Vertex::init_geometry(const ExactPoint& P) {
235 292343 point_exact = P;
236 Numeric::optimize_number_representation(point_exact);
237 #ifndef GEOGRAM_USE_EXACT_NT
238 292343 l = ExactCDT2d::squared_length(P[mit->u_], P[mit->v_], P.w);
239 #endif
240 292343 }
241
242 228 MeshInTriangle::MeshInTriangle(MeshSurfaceIntersection& EM) :
243 228 exact_mesh_(EM),
244 228 mesh_(EM.readonly_mesh()),
245 228 f1_(NO_INDEX),
246 228 dry_run_(false),
247 228 use_pred_cache_insert_buffer_(false)
248 {
249 #ifdef GEOGRAM_USE_EXACT_NT
250 CDTBase2d::exact_incircle_ = true;
251 #else
252 // Since incircle() with expansions computes approximated
253 // lifted coordinate, we need to activate additional
254 // checks for Delaunayization.
255 228 CDTBase2d::exact_incircle_ = false;
256 #endif
257 228 }
258
259 15506 void MeshInTriangle::clear() {
260 15506 vertex_.resize(0);
261 15506 edges_.resize(0);
262 15506 f1_ = NO_INDEX;
263 pred_cache_.clear();
264 15506 pred_cache_insert_buffer_.resize(0);
265 15506 use_pred_cache_insert_buffer_ = false;
266 15506 CDTBase2d::clear();
267 15506 }
268
269 15506 void MeshInTriangle::begin_facet(index_t f) {
270 15506 f1_ = f;
271
272 15506 vec3 p1 = mesh_facet_vertex(f,0);
273 15506 vec3 p2 = mesh_facet_vertex(f,1);
274 15506 vec3 p3 = mesh_facet_vertex(f,2);
275
276 geo_debug_assert(!PCK::aligned_3d(p1,p2,p3));
277
278 15506 f1_normal_axis_ = PCK::triangle_normal_axis(
279 p1,p2,p3
280 );
281
282 15506 u_ = coord_index_t((f1_normal_axis_ + 1) % 3);
283 15506 v_ = coord_index_t((f1_normal_axis_ + 2) % 3);
284
2/2
✓ Branch 0 taken 46518 times.
✓ Branch 1 taken 15506 times.
 62024 for(index_t lv=0; lv<3; ++lv) {
285 46518 vertex_.emplace_back(this, f, lv);
286 }
287
288 15506 CDTBase2d::create_enclosing_triangle(0,1,2);
289
290 15506 edges_.emplace_back(1,2);
291 15506 edges_.emplace_back(2,0);
292 15506 edges_.emplace_back(0,1);
293 15506 }
294
295 262926 index_t MeshInTriangle::add_vertex(
296 index_t f2, TriangleRegion R1, TriangleRegion R2
297 ) {
298 geo_debug_assert(f1_ != NO_INDEX);
299
300 // If vertex is a macro-vertex, return it directly.
301
2/2
✓ Branch 1 taken 22785 times.
✓ Branch 2 taken 240141 times.
 262926 if(region_dim(R1) == 0) {
302 22785 return index_t(R1);
303 }
304
305 // Create the vertex
306 240141 vertex_.emplace_back(this, f1_, f2, R1, R2);
307
308 // Insert it into the triangulation
309
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240141 times.
 240141 index_t v = CDTBase2d::insert(vertex_.size()-1);
310
311 // If it was an existing vertex, return the existing vertex
312
2/2
✓ Branch 0 taken 151848 times.
✓ Branch 1 taken 88293 times.
 240141 if(vertex_.size() > CDTBase2d::nv()) {
313 vertex_.pop_back();
314 }
315 return v;
316 }
317
318 127160 void MeshInTriangle::add_edge(
319 index_t f2,
320 TriangleRegion AR1, TriangleRegion AR2,
321 TriangleRegion BR1, TriangleRegion BR2
322 ) {
323 127160 index_t v1 = add_vertex(f2, AR1, AR2);
324 127160 index_t v2 = add_vertex(f2, BR1, BR2);
325
326 // If both extremities are on the same edge of f1,
327 // we do not add the edge, because it will be generated
328 // when remeshing the edge of f1
329
2/2
✓ Branch 2 taken 49028 times.
✓ Branch 3 taken 78132 times.
 127160 if(region_dim(regions_convex_hull(AR1,BR1)) == 1) {
330 49028 return;
331 }
332
333 // Generate also the combinatorial information of the edge,
334 // that indicates whether both extremities are on the same
335 // edge of f2 (useful later to compute the intersections)
336 78132 edges_.emplace_back(v1,v2,f2,regions_convex_hull(AR2,BR2));
337
338 // Constraints will be added to the triangulation during commit()
339 }
340
341 15506 void MeshInTriangle::commit() {
342
343
2/2
✓ Branch 0 taken 124650 times.
✓ Branch 1 taken 15506 times.
 140156 for(const Edge& E: edges_) {
344 124650 CDTBase2d::insert_constraint(E.v1, E.v2);
345 }
346
347
1/2
✓ Branch 0 taken 15506 times.
✗ Branch 1 not taken.
 15506 if(dry_run_) {
348 return;
349 }
350
351 // Protect global mesh from concurrent accesses
352 15506 exact_mesh_.lock();
353
354 // Create vertices and facets in target mesh
355
2/2
✓ Branch 0 taken 140495 times.
✓ Branch 1 taken 15506 times.
 296496 for(index_t i=0; i<vertex_.size(); ++i) {
356 // Vertex already exists in this MeshInTriangle
357
2/2
✓ Branch 0 taken 49257 times.
✓ Branch 1 taken 91238 times.
 140495 if(vertex_[i].mesh_vertex_index != NO_INDEX) {
358 49257 continue;
359 }
360 91238 vertex_[i].mesh_vertex_index =
361 91238 exact_mesh_.find_or_create_exact_vertex(
362 91238 vertex_[i].point_exact
363 );
364 }
365
366 // Create facets in target mesh
367
2/2
✓ Branch 0 taken 131392 times.
✓ Branch 1 taken 15506 times.
 278290 for(index_t t=0; t<CDTBase2d::nT(); ++t) {
368 index_t i = CDTBase2d::Tv(t,0);
369 index_t j = CDTBase2d::Tv(t,1);
370 index_t k = CDTBase2d::Tv(t,2);
371 131392 i = vertex_[i].mesh_vertex_index;
372 131392 j = vertex_[j].mesh_vertex_index;
373 131392 k = vertex_[k].mesh_vertex_index;
374 131392 index_t new_t = target_mesh().facets.create_triangle(i,j,k);
375 // Copy all attributes from initial facet
376 131392 target_mesh().facets.attributes().copy_item(new_t, f1_);
377 }
378
379 // We are done with modification in the mesh
380 15506 exact_mesh_.unlock();
381 }
382
383 ✗ void MeshInTriangle::get_constraints(Mesh& M, bool with_edges) const {
384 ✗ if(M.vertices.nb() == 0) {
385 ✗ M.vertices.set_dimension(2);
386 ✗ for(index_t v=0; v<vertex_.size(); ++v) {
387 ✗ vec2 p = vertex_[v].get_UV_approx();
388 ✗ M.vertices.create_vertex(p.data());
389 }
390 }
391 ✗ if(with_edges && M.edges.nb() == 0) {
392 ✗ for(const Edge& E: edges_) {
393 ✗ M.edges.create_edge(E.v1, E.v2);
394 }
395 }
396 ✗ }
397
398 /**
399 * \brief Tests the parity of the permutation of a list of
400 * three distinct indices with respect to the canonical order.
401 */
402 5032304 static bool odd_order(index_t i, index_t j, index_t k) {
403 // Implementation: sort the elements (bubble sort is OK for
404 // such a small number), and invert parity each time
405 // two elements are swapped.
406 5032304 index_t tab[3] = { i, j, k};
407 const int N = 3;
408 bool result = false;
409
2/2
✓ Branch 0 taken 10064608 times.
✓ Branch 1 taken 5032304 times.
 15096912 for (int I = 0; I < N - 1; ++I) {
410
2/2
✓ Branch 0 taken 15096912 times.
✓ Branch 1 taken 10064608 times.
 25161520 for (int J = 0; J < N - I - 1; ++J) {
411
2/2
✓ Branch 0 taken 8965107 times.
✓ Branch 1 taken 6131805 times.
 15096912 if (tab[J] > tab[J + 1]) {
412 std::swap(tab[J], tab[J + 1]);
413 8965107 result = !result;
414 }
415 }
416 }
417 5032304 return result;
418 }
419
420 240141 void MeshInTriangle::begin_insert_transaction() {
421 240141 use_pred_cache_insert_buffer_ = true;
422 240141 }
423
424 88293 void MeshInTriangle::commit_insert_transaction() {
425
2/2
✓ Branch 0 taken 879623 times.
✓ Branch 1 taken 88293 times.
 967916 for(const auto& it: pred_cache_insert_buffer_) {
426 879623 pred_cache_[it.first] = it.second;
427 }
428 88293 pred_cache_insert_buffer_.resize(0);
429 88293 use_pred_cache_insert_buffer_ = false;
430 88293 }
431
432 151848 void MeshInTriangle::rollback_insert_transaction() {
433 151848 pred_cache_insert_buffer_.resize(0);
434 151848 use_pred_cache_insert_buffer_ = false;
435 151848 }
436
437
2/2
✓ Branch 0 taken 3220594 times.
✓ Branch 1 taken 1811710 times.
 5032304 Sign MeshInTriangle::orient2d(index_t vx1,index_t vx2,index_t vx3) const {
438
439 trindex K(vx1, vx2, vx3);
440
441
2/2
✓ Branch 0 taken 3531834 times.
✓ Branch 1 taken 1500470 times.
 5032304 if(use_pred_cache_insert_buffer_) {
442 3531834 Sign result = PCK::orient_2d_projected(
443 3531834 vertex_[K.indices[0]].point_exact,
444 3531834 vertex_[K.indices[1]].point_exact,
445 3531834 vertex_[K.indices[2]].point_exact,
446 3531834 f1_normal_axis_
447 );
448 3531834 pred_cache_insert_buffer_.push_back(std::make_pair(K, result));
449
2/2
✓ Branch 0 taken 1743760 times.
✓ Branch 1 taken 1788074 times.
 3531834 if(odd_order(vx1,vx2,vx3)) {
450 1743760 result = Sign(-result);
451 }
452 return result;
453 }
454
455 bool inserted;
456 std::map<trindex, Sign>::iterator it;
457 1500470 std::tie(it,inserted) = pred_cache_.insert(std::make_pair(K,ZERO));
458 Sign result;
459
460
2/2
✓ Branch 0 taken 640845 times.
✓ Branch 1 taken 859625 times.
 1500470 if(inserted) {
461 640845 result = PCK::orient_2d_projected(
462 640845 vertex_[K.indices[0]].point_exact,
463 640845 vertex_[K.indices[1]].point_exact,
464 640845 vertex_[K.indices[2]].point_exact,
465 640845 f1_normal_axis_
466 );
467 640845 it->second = result;
468 } else {
469 859625 result = it->second;
470 }
471
472
2/2
✓ Branch 0 taken 690389 times.
✓ Branch 1 taken 810081 times.
 1500470 if(odd_order(vx1,vx2,vx3)) {
473 690389 result = Sign(-result);
474 }
475
476 return result;
477 }
478
479 185930 Sign MeshInTriangle::incircle(
480 index_t v1,index_t v2,index_t v3,index_t v4
481 ) const {
482 exact::vec2h p1(
483 185930 vertex_[v1].point_exact[u_],
484 185930 vertex_[v1].point_exact[v_],
485 185930 vertex_[v1].point_exact.w
486 185930 );
487 exact::vec2h p2(
488 vertex_[v2].point_exact[u_],
489 vertex_[v2].point_exact[v_],
490 185930 vertex_[v2].point_exact.w
491 185930 );
492 exact::vec2h p3(
493 vertex_[v3].point_exact[u_],
494 vertex_[v3].point_exact[v_],
495 185930 vertex_[v3].point_exact.w
496 185930 );
497 exact::vec2h p4(
498 vertex_[v4].point_exact[u_],
499 vertex_[v4].point_exact[v_],
500 185930 vertex_[v4].point_exact.w
501 185930 );
502 #ifdef GEOGRAM_USE_EXACT_NT
503 return PCK::incircle_2d_SOS(p1,p2,p3,p4);
504 #else
505 185930 return PCK::incircle_2d_SOS_with_lengths(
506 p1,p2,p3,p4,
507 185930 vertex_[v1].l,
508 185930 vertex_[v2].l,
509 185930 vertex_[v3].l,
510
1/2
✓ Branch 1 taken 185930 times.
✗ Branch 2 not taken.
 185930 vertex_[v4].l
511 185930 );
512 #endif
513 185930 }
514
515
1/2
✓ Branch 1 taken 5684 times.
✗ Branch 2 not taken.
 5684 index_t MeshInTriangle::create_intersection(
516 index_t e1, index_t i, index_t j,
517 index_t e2, index_t k, index_t l
518 ) {
519 geo_argused(i);
520 geo_argused(j);
521 geo_argused(k);
522 geo_argused(l);
523
524 ExactPoint I;
525
1/2
✓ Branch 1 taken 5684 times.
✗ Branch 2 not taken.
 5684 get_edge_edge_intersection(e1,e2,I);
526
2/4
✓ Branch 1 taken 5684 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5684 times.
 5684 vertex_.emplace_back(this,I);
527 5684 index_t x = vertex_.size()-1;
528
2/2
✓ Branch 0 taken 5660 times.
✓ Branch 1 taken 24 times.
 5684 CDTBase2d::v2T_.push_back(NO_INDEX);
529 geo_debug_assert(x == CDTBase2d::nv_);
530 5684 ++CDTBase2d::nv_;
531 5684 return x;
532 5684 }
533
534 5684 void MeshInTriangle::get_edge_edge_intersection(
535 index_t e1, index_t e2, ExactPoint& I
536 ) const {
537 5684 index_t f1 = f1_;
538 5684 index_t f2 = edges_[e1].sym.f2;
539 5684 index_t f3 = edges_[e2].sym.f2;
540
541 geo_debug_assert(f1 != NO_INDEX);
542 geo_debug_assert(f2 != NO_INDEX);
543 geo_debug_assert(f3 != NO_INDEX);
544
545 vec3 P[9] = {
546 5684 mesh_facet_vertex(f1,0), mesh_facet_vertex(f1,1),
547 5684 mesh_facet_vertex(f1,2),
548 5684 mesh_facet_vertex(f2,0), mesh_facet_vertex(f2,1),
549 5684 mesh_facet_vertex(f2,2),
550 5684 mesh_facet_vertex(f3,0), mesh_facet_vertex(f3,1),
551 5684 mesh_facet_vertex(f3,2)
552 };
553
554
2/2
✓ Branch 1 taken 4357 times.
✓ Branch 2 taken 1327 times.
 5684 if(!get_three_planes_intersection(
555 I,
556 P[0], P[1], P[2],
557 P[3], P[4], P[5],
558 P[6], P[7], P[8]
559 )) {
560 4357 get_edge_edge_intersection_2D(e1,e2,I);
561 4357 return;
562 }
563 }
564
565 4357 void MeshInTriangle::get_edge_edge_intersection_2D(
566 index_t e1, index_t e2, ExactPoint& I
567 ) const {
568 const Edge& E1 = edges_[e1];
569 const Edge& E2 = edges_[e2];
570
571 if(
572
4/4
✓ Branch 1 taken 4322 times.
✓ Branch 2 taken 35 times.
✓ Branch 3 taken 4298 times.
✓ Branch 4 taken 24 times.
 8679 region_dim(E1.sym.R2) == 1 &&
573 4322 region_dim(E2.sym.R2) == 1
574 ) {
575 4298 index_t le1 = index_t(E1.sym.R2)-index_t(T2_RGN_E0);
576 4298 index_t le2 = index_t(E2.sym.R2)-index_t(T2_RGN_E0);
577 geo_debug_assert(le1 < 3);
578 geo_debug_assert(le2 < 3);
579
580 4298 vec2 p1_uv = mesh_facet_vertex_UV(E1.sym.f2, (le1+1)%3);
581 4298 vec2 p2_uv = mesh_facet_vertex_UV(E1.sym.f2, (le1+2)%3);
582 4298 vec2 q1_uv = mesh_facet_vertex_UV(E2.sym.f2, (le2+1)%3);
583 4298 vec2 q2_uv = mesh_facet_vertex_UV(E2.sym.f2, (le2+2)%3);
584
585 4298 exact::vec2 C1 = make_vec2<exact::vec2>(p1_uv, p2_uv);
586 4298 exact::vec2 C2 = make_vec2<exact::vec2>(q2_uv, q1_uv);
587 4298 exact::vec2 B = make_vec2<exact::vec2>(p1_uv, q1_uv);
588
1/2
✓ Branch 1 taken 4298 times.
✗ Branch 2 not taken.
 4298 exact::scalar d = det(C1,C2);
589 geo_debug_assert(d.sign() != ZERO);
590
1/4
✓ Branch 1 taken 4298 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 4298 exact::rational t(det(B,C2),d);
591 4298 I = mix(
592 t,
593 8596 mesh_facet_vertex(E1.sym.f2,(le1+1)%3),
594
1/2
✓ Branch 1 taken 4298 times.
✗ Branch 2 not taken.
 4298 mesh_facet_vertex(E1.sym.f2,(le1+2)%3)
595 4298 );
596
597 8596 } else {
598 geo_debug_assert(
599 region_dim(E1.sym.R2) == 1 || region_dim(E2.sym.R2) == 1
600 );
601 59 index_t f1 = E1.sym.f2;
602 59 TriangleRegion R1 = E1.sym.R2;
603 59 index_t f2 = E2.sym.f2;
604 59 TriangleRegion R2 = E2.sym.R2;
605
2/2
✓ Branch 1 taken 24 times.
✓ Branch 2 taken 35 times.
 59 if(region_dim(R1) == 1) {
606 std::swap(f1,f2);
607 std::swap(R1,R2);
608 }
609
610 59 index_t e = index_t(R2) - index_t(T2_RGN_E0);
611 geo_debug_assert(e < 3);
612
613 59 I = plane_line_intersection(
614 118 mesh_facet_vertex(f1,0),
615 118 mesh_facet_vertex(f1,1),
616 118 mesh_facet_vertex(f1,2),
617 118 mesh_facet_vertex(f2,(e+1)%3),
618 118 mesh_facet_vertex(f2,(e+2)%3)
619 59 );
620 }
621 4357 }
622
623 ✗ void MeshInTriangle::save(const std::string& filename) const {
624 ✗ Mesh M;
625 ✗ M.vertices.set_dimension(2);
626 ✗ for(index_t v=0; v<CDTBase2d::nv(); ++v) {
627 ✗ vec2 p = vertex_[v].get_UV_approx();
628 ✗ M.vertices.create_vertex(p.data());
629 }
630 ✗ for(index_t t=0; t<CDTBase2d::nT(); ++t) {
631 ✗ M.facets.create_triangle(
632 CDTBase2d::Tv(t,0),
633 CDTBase2d::Tv(t,1),
634 CDTBase2d::Tv(t,2)
635 );
636 }
637
638 Attribute<double> tex_coord;
639 ✗ tex_coord.create_vector_attribute(
640 ✗ M.facet_corners.attributes(), "tex_coord", 2
641 );
642 static double triangle_tex[3][2] = {
643 {0.0, 0.0},
644 {1.0, 0.0},
645 {0.0, 1.0}
646 };
647 ✗ for(index_t c: M.facet_corners) {
648 ✗ tex_coord[2*c] = triangle_tex[c%3][0];
649 ✗ tex_coord[2*c+1] = triangle_tex[c%3][1];
650 }
651 ✗ mesh_save(M, filename);
652 ✗ }
653
654 /**************************************************************************/
655
656 285 CoplanarFacets::CoplanarFacets(
657 MeshSurfaceIntersection& I, bool clear_attributes,
658 double angle_tolerance
659 285 ) :
660
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 I_(I),
661
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 mesh_(I.target_mesh()),
662 285 mesh_copy_(I.readonly_mesh()),
663 285 angle_tolerance_(angle_tolerance),
664
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 facet_group_(I.target_mesh().facets.attributes(),"group"),
665
1/4
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 285 keep_vertex_(I.target_mesh().vertices.attributes(),"keep"),
666 c_is_coplanar_(
667
1/4
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 285 I.target_mesh().facet_corners.attributes(),"is_coplanar"
668 ),
669
3/6
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 57 times.
✓ Branch 4 taken 228 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 570 f_is_flipped_(I.target_mesh().facets.attributes(),"flipped"),
670 halfedges_(*this),
671 285 polylines_(*this)
672 {
673
2/2
✓ Branch 0 taken 57 times.
✓ Branch 1 taken 228 times.
 285 if(clear_attributes) {
674
2/2
✓ Branch 0 taken 123494 times.
✓ Branch 1 taken 57 times.
 123551 for(index_t f: mesh_.facets) {
675 123494 facet_group_[f] = NO_INDEX;
676 }
677
2/2
✓ Branch 0 taken 61753 times.
✓ Branch 1 taken 57 times.
 61810 for(index_t v: mesh_.vertices) {
678 61753 keep_vertex_[v] = false;
679 }
680
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 find_coplanar_facets();
681 }
682
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 f_visited_.assign(mesh_.facets.nb(),false);
683
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 h_visited_.assign(mesh_.facet_corners.nb(),false);
684
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 v_visited_.assign(mesh_.vertices.nb(),false);
685
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 v_idx_.assign(mesh_.vertices.nb(),NO_INDEX);
686 285 }
687
688 57 void CoplanarFacets::find_coplanar_facets() {
689
690 // Positioned by MeshSurfaceIntersection::build_Weiler_model()
691 Attribute<bool> corner_is_on_border(
692
1/2
✓ Branch 2 taken 57 times.
✗ Branch 3 not taken.
 57 mesh_.facet_corners.attributes(), "is_on_border"
693 );
694
695
2/2
✓ Branch 0 taken 370482 times.
✓ Branch 1 taken 57 times.
 370539 for(index_t c: mesh_.facet_corners) {
696 370482 c_is_coplanar_[c] = false;
697 }
698
699 // TODO: when there is an angle tolerance, one should check instead
700 // angle deviation w.r.t. a single seed facet per facet group, because
701 // with the present algorithm, if a large number of tiny facets are
702 // connected (e.g. highly tessellated cylinder), one may group facets
703 // with large angle deviation (without seeing it because each facet has
704 // small angle deviation w.r.t. its neighbors).
705
706
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 parallel_for(
707
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 0, mesh_.facet_corners.nb(),
708
1/4
✓ Branch 0 taken 57 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 57 [&](index_t c1) {
709 370482 index_t f1 = (c1 / 3);
710 370482 index_t le1 = (c1 % 3);
711
1/2
✓ Branch 0 taken 370482 times.
✗ Branch 1 not taken.
 370482 index_t f2 = mesh_.facet_corners.adjacent_facet(c1);
712
713
1/2
✓ Branch 0 taken 370482 times.
✗ Branch 1 not taken.
 370482 if(f2 == NO_INDEX) {
714 216659 return;
715 }
716
717 // do not traverse true borders
718
1/2
✓ Branch 0 taken 370482 times.
✗ Branch 1 not taken.
 370482 if(corner_is_on_border[c1]) {
719 return;
720 }
721
722 index_t v11 = mesh_.facets.vertex(f1,le1);
723 370482 index_t v12 = mesh_.facets.vertex(f1,(le1+1)%3);
724
1/2
✓ Branch 0 taken 370482 times.
✗ Branch 1 not taken.
 370482 index_t le2 = mesh_.facets.find_edge(f2,v12,v11);
725 index_t c2 = mesh_.facets.corner(f2,le2);
726
727 #ifdef GEO_DEBUG
728 index_t v21 = mesh_.facets.vertex(f2,le2);
729 index_t v22 = mesh_.facets.vertex(f2,(le2+1)%3);
730 index_t v13 = mesh_.facets.vertex(f1,(le1+2)%3);
731 index_t v23 = mesh_.facets.vertex(f2,(le2+2)%3);
732
733 geo_debug_assert(v11 == v22);
734 geo_debug_assert(v12 == v21);
735 geo_debug_assert(v11!=v12 && v12!=v13 && v13!=v11);
736 geo_debug_assert(v21!=v22 && v22!=v23 && v23!=v21);
737 #endif
738
739
2/2
✓ Branch 0 taken 185241 times.
✓ Branch 1 taken 185241 times.
 370482 if(c1 > c2) {
740 return;
741 }
742
743 // Small optimization: if both triangles come from same
744 // original facet then they are coplanar
745
2/2
✓ Branch 0 taken 31418 times.
✓ Branch 1 taken 153823 times.
 185241 if(I_.get_initial_facet(f1) == I_.get_initial_facet(f2)) {
746 31418 c_is_coplanar_[c1] = true;
747 31418 c_is_coplanar_[c2] = true;
748 31418 return;
749 }
750
751 // Use original triangles for co-planarity test
752 153823 auto [p1, p2, p3] = I_.get_initial_facet_vertices(f1);
753 153823 auto [q1, q2, q3] = I_.get_initial_facet_vertices(f2);
754
2/2
✓ Branch 1 taken 30362 times.
✓ Branch 2 taken 123461 times.
 153823 if(triangles_are_coplanar(p1,p2,p3,q1,q2,q3)) {
755 30362 c_is_coplanar_[c1] = true;
756 30362 c_is_coplanar_[c2] = true;
757 }
758 }
759 );
760 57 }
761
762 146176 void CoplanarFacets::get(index_t f, index_t group_id) {
763
764 146176 facets_.resize(0);
765 146176 vertices_.resize(0);
766 146176 halfedges_.initialize();
767 146176 polylines_.initialize();
768
769 // Get facets
770 {
771 std::stack<index_t> S;
772
1/2
✓ Branch 0 taken 146176 times.
✗ Branch 1 not taken.
 146176 facet_group_[f] = group_id;
773 f_visited_[f] = true;
774 S.push(f);
775 facets_.push_back(f);
776
2/2
✓ Branch 0 taken 246988 times.
✓ Branch 1 taken 146176 times.
 393164 while(!S.empty()) {
777
2/2
✓ Branch 0 taken 246212 times.
✓ Branch 1 taken 776 times.
 246988 index_t f1 = S.top();
778 S.pop();
779
2/2
✓ Branch 0 taken 740964 times.
✓ Branch 1 taken 246988 times.
 987952 for(index_t le1=0; le1<3; ++le1) {
780
1/2
✓ Branch 0 taken 740964 times.
✗ Branch 1 not taken.
 740964 index_t f2 = mesh_.facets.adjacent(f1,le1);
781 if(
782
5/6
✓ Branch 0 taken 740964 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 594656 times.
✓ Branch 3 taken 146308 times.
✓ Branch 4 taken 100812 times.
✓ Branch 5 taken 493844 times.
 740964 f2 != NO_INDEX && !f_visited_[f2] &&
783
2/2
✓ Branch 0 taken 100812 times.
✓ Branch 1 taken 493844 times.
 594656 c_is_coplanar_[mesh_.facets.corner(f1,le1)]
784 ) {
785
2/2
✓ Branch 0 taken 100036 times.
✓ Branch 1 taken 776 times.
 100812 facet_group_[f2] = facet_group_[f1];
786 f_visited_[f2] = true;
787 S.push(f2);
788 facets_.push_back(f2);
789 }
790 }
791 }
792
2/2
✓ Branch 0 taken 246988 times.
✓ Branch 1 taken 146176 times.
 393164 for(index_t cur_f: facets_) {
793 f_visited_[cur_f] = false;
794 }
795 }
796
797 146176 group_id_ = group_id;
798
799 // Initialize projection coordinates, using original facet
800 {
801 146176 auto [p1, p2, p3] = I_.get_initial_facet_vertices(facets_[0]);
802 146176 coord_index_t projection_axis = PCK::triangle_normal_axis(p1,p2,p3);
803 146176 u_ = coord_index_t((projection_axis+1)%3);
804 146176 v_ = coord_index_t((projection_axis+2)%3);
805 Sign o = PCK::orient_2d(
806 146176 vec2(p1[u_],p1[v_]), vec2(p2[u_],p2[v_]), vec2(p3[u_],p3[v_])
807 );
808 geo_debug_assert(o != ZERO);
809
2/2
✓ Branch 0 taken 71088 times.
✓ Branch 1 taken 75088 times.
 146176 if(o < 0) {
810 std::swap(u_,v_);
811 }
812 }
813
814 // Get vertices and halfedges
815 {
816
2/2
✓ Branch 0 taken 246988 times.
✓ Branch 1 taken 146176 times.
 393164 for(index_t f1: facets_) {
817
2/2
✓ Branch 0 taken 740964 times.
✓ Branch 1 taken 246988 times.
 987952 for(index_t le=0; le<3; ++le) {
818
1/2
✓ Branch 0 taken 740964 times.
✗ Branch 1 not taken.
 740964 index_t f2 = mesh_.facets.adjacent(f1,le);
819 if(
820
3/4
✓ Branch 0 taken 740964 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 247120 times.
✓ Branch 3 taken 493844 times.
 740964 f2 == NO_INDEX ||
821
2/2
✓ Branch 0 taken 247120 times.
✓ Branch 1 taken 493844 times.
 740964 !c_is_coplanar_[mesh_.facets.corner(f1,le)]
822 ) {
823 493844 halfedges_.add(mesh_.facets.corners_begin(f1)+le);
824
1/2
✓ Branch 0 taken 493844 times.
✗ Branch 1 not taken.
 493844 index_t v1 = mesh_.facets.vertex(f1,le);
825
2/2
✓ Branch 0 taken 174907 times.
✓ Branch 1 taken 318937 times.
 493844 index_t v2 = mesh_.facets.vertex(f1,(le+1)%3);
826
2/2
✓ Branch 0 taken 174907 times.
✓ Branch 1 taken 318937 times.
 493844 if(!v_visited_[v1]) {
827
2/2
✓ Branch 0 taken 174285 times.
✓ Branch 1 taken 622 times.
 174907 v_idx_[v1] = vertices_.size();
828 vertices_.push_back(v1);
829 v_visited_[v1] = true;
830 }
831
2/2
✓ Branch 0 taken 307512 times.
✓ Branch 1 taken 186332 times.
 493844 if(!v_visited_[v2]) {
832
2/2
✓ Branch 0 taken 306830 times.
✓ Branch 1 taken 682 times.
 307512 v_idx_[v2] = vertices_.size();
833 vertices_.push_back(v2);
834 v_visited_[v2] = true;
835 }
836 } else {
837 // This one for the particular case of a non-manifold
838 // vertex, such as a cone apex touching a facet
839 // (ThingiCSG/Basic/cube_cone_1.scad)
840 247120 index_t v = mesh_.facets.vertex(f1,le);
841
4/4
✓ Branch 0 taken 71214 times.
✓ Branch 1 taken 175906 times.
✓ Branch 2 taken 11423 times.
✓ Branch 3 taken 59791 times.
 247120 if(keep_vertex_[v] && !v_visited_[v]) {
842 vertices_.push_back(v);
843 v_visited_[v] = true;
844 }
845 }
846 }
847 }
848
2/2
✓ Branch 0 taken 493842 times.
✓ Branch 1 taken 146176 times.
 640018 for(index_t v: vertices_) {
849 v_visited_[v] = false;
850 }
851 }
852
853 // Get polylines
854 {
855 // Get all polylines starting from vertices with more than
856 // 2 incident halfedges.
857
2/2
✓ Branch 0 taken 493842 times.
✓ Branch 1 taken 146176 times.
 640018 for(index_t v: vertices_) {
858
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 493840 times.
 493842 if(halfedges_.nb_halfedges_around_vertex(v) > 1) {
859 for(
860 index_t h=halfedges_.vertex_first_halfedge(v);
861
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 6 h != NO_INDEX; h = halfedges_.next_around_vertex(h)
862 ) {
863
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if(!h_visited_[h]) {
864 polylines_.begin_polyline();
865 index_t h2 = h;
866 do {
867 geo_debug_assert(!h_visited_[h2]);
868 h_visited_[h2] = true;
869 24 polylines_.add_halfedge(h2);
870 24 h2 = halfedges_.next_along_polyline(h2);
871
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 4 times.
 24 } while(h2 != NO_INDEX && h2 != h);
872 polylines_.end_polyline();
873 }
874 }
875 }
876 }
877 // There can be also closed halfedge loops with no irregular vertex
878
4/4
✓ Branch 0 taken 146376 times.
✓ Branch 1 taken 347468 times.
✓ Branch 2 taken 493844 times.
✓ Branch 3 taken 146176 times.
 640020 for(index_t h: halfedges_) {
879
2/2
✓ Branch 0 taken 146376 times.
✓ Branch 1 taken 347468 times.
 493844 if(!h_visited_[h]) {
880 polylines_.begin_polyline();
881 index_t h2 = h;
882 do {
883 geo_debug_assert(!h_visited_[h2]);
884 h_visited_[h2] = true;
885 493820 polylines_.add_halfedge(h2);
886 493820 h2 = halfedges_.next_along_polyline(h2);
887 geo_debug_assert(h2 != NO_INDEX);
888
2/2
✓ Branch 0 taken 347444 times.
✓ Branch 1 taken 146376 times.
 493820 } while(h2 != h);
889 polylines_.end_polyline();
890 }
891 }
892
2/2
✓ Branch 0 taken 493844 times.
✓ Branch 1 taken 146176 times.
 640020 for(index_t h: halfedges_) {
893 h_visited_[h] = false;
894 }
895 }
896 146176 }
897
898
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 73088 times.
 73088 void CoplanarFacets::mark_vertices_to_keep() {
899
2/2
✓ Branch 0 taken 73190 times.
✓ Branch 1 taken 73088 times.
 146278 for(index_t P: polylines_) {
900 73190 index_t first_v = polylines_.first_vertex(P);
901 73190 index_t last_v = polylines_.last_vertex(P);
902
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 73190 times.
 73190 if(first_v != last_v) {
903 ✗ keep_vertex_[first_v] = true;
904 keep_vertex_[last_v] = true;
905 }
906 73190 index_t v1 = polylines_.prev_first_vertex(P);
907 index_t v2 = NO_INDEX;
908 index_t v3 = NO_INDEX;
909
2/2
✓ Branch 0 taken 246922 times.
✓ Branch 1 taken 73190 times.
 320112 for(index_t h: polylines_.halfedges(P)) {
910
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 246922 times.
 246922 if(v1 == NO_INDEX) {
911 ✗ continue;
912 }
913 v2 = halfedges_.vertex(h,0);
914 v3 = halfedges_.vertex(h,1);
915
1/2
✓ Branch 1 taken 246922 times.
✗ Branch 2 not taken.
 246922 ExactPoint p1 = I_.exact_vertex(v1);
916
1/2
✓ Branch 1 taken 246922 times.
✗ Branch 2 not taken.
 246922 ExactPoint p2 = I_.exact_vertex(v2);
917
1/2
✓ Branch 1 taken 246922 times.
✗ Branch 2 not taken.
 246922 ExactPoint p3 = I_.exact_vertex(v3);
918
3/4
✓ Branch 1 taken 246922 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 238479 times.
✓ Branch 4 taken 8443 times.
 246922 if(!edges_are_colinear(p1,p2,p3)) {
919 238479 keep_vertex_[v2] = true;
920 }
921 v1 = v2;
922 246922 }
923 }
924 73088 }
925
926 ✗ void CoplanarFacets::save_borders(const std::string& filename) {
927 ✗ Mesh borders;
928 ✗ borders.vertices.set_dimension(2);
929 index_t cur_idx = 0;
930 ✗ for(index_t v: vertices_) {
931 ✗ vec3 p = mesh_.vertices.point(v);
932 ✗ vec2 q(p[u_],p[v_]);
933 ✗ borders.vertices.create_vertex(q.data());
934 ✗ v_idx_[v] = cur_idx;
935 ✗ ++cur_idx;
936 }
937
938 ✗ for(index_t h: halfedges_) {
939 index_t v1 = halfedges_.vertex(h,0);
940 index_t v2 = halfedges_.vertex(h,1);
941 ✗ v1 = v_idx_[v1];
942 ✗ v2 = v_idx_[v2];
943 geo_debug_assert(v1 != NO_INDEX);
944 geo_debug_assert(v2 != NO_INDEX);
945 borders.edges.create_edge(v1,v2);
946 }
947
948 ✗ Attribute<bool> selection(borders.vertices.attributes(), "selection");
949 ✗ for(index_t v: vertices_) {
950 geo_debug_assert(v_idx_[v] != NO_INDEX);
951 ✗ selection[v_idx_[v]] = keep_vertex_[v];
952 }
953 ✗ mesh_save(borders,filename);
954 ✗ }
955
956 ✗ void CoplanarFacets::save_facet_group(const std::string& filename) {
957 ✗ Mesh M;
958 ✗ Attribute<bool> keep_vertex(M.vertices.attributes(),"keep");
959 ✗ M.vertices.set_dimension(2);
960 ✗ for(index_t f: facets_) {
961 ✗ for(index_t lv=0; lv<3; ++lv) {
962 ✗ index_t v = mesh_.facets.vertex(f,lv);
963 ✗ v_idx_[v] = NO_INDEX;
964 }
965 }
966 ✗ for(index_t f: facets_) {
967 ✗ for(index_t lv=0; lv<3; ++lv) {
968 ✗ index_t v = mesh_.facets.vertex(f,lv);
969 ✗ if(v_idx_[v] == NO_INDEX) {
970 ✗ vec3 p = mesh_.vertices.point(v);
971 ✗ vec2 q(p[u_], p[v_]);
972 ✗ v_idx_[v] = M.vertices.create_vertex(q.data());
973 ✗ keep_vertex[v_idx_[v]] = keep_vertex_[v];
974 }
975 }
976 ✗ M.facets.create_triangle(
977 v_idx_[mesh_.facets.vertex(f,0)],
978 v_idx_[mesh_.facets.vertex(f,1)],
979 ✗ v_idx_[mesh_.facets.vertex(f,2)]
980 );
981 }
982
983 ✗ for(index_t f: facets_) {
984 ✗ for(index_t lv=0; lv<3; ++lv) {
985 ✗ index_t v = mesh_.facets.vertex(f,lv);
986 ✗ v_idx_[v] = NO_INDEX;
987 }
988 }
989
990 ✗ M.facets.connect();
991 ✗ mesh_save(M,filename);
992 ✗ }
993
994 8069 void CoplanarFacets::triangulate() {
995
996 // Compute 2D projected BBOX
997 8069 double umin = Numeric::max_float64();
998 8069 double vmin = Numeric::max_float64();
999 8069 double umax = -Numeric::max_float64();
1000 8069 double vmax = -Numeric::max_float64();
1001
2/2
✓ Branch 0 taken 58475 times.
✓ Branch 1 taken 8069 times.
 66544 for(index_t f: facets_) {
1002
2/2
✓ Branch 0 taken 175425 times.
✓ Branch 1 taken 58475 times.
 233900 for(index_t lv=0; lv<3; ++lv) {
1003
1/2
✓ Branch 0 taken 175425 times.
✗ Branch 1 not taken.
 175425 index_t vx = mesh_.facets.vertex(f,lv);
1004 175425 double u = mesh_.vertices.point(vx)[u_];
1005 175425 double v = mesh_.vertices.point(vx)[v_];
1006 175425 umin = std::min(umin, u);
1007 175425 umax = std::max(umax, u);
1008 175425 vmin = std::min(vmin, v);
1009 175425 vmax = std::max(vmax, v);
1010 }
1011 }
1012 8069 double d = std::max(umax-umin, vmax-vmin);
1013 8069 d *= 10.0;
1014 d = std::max(d, 1.0);
1015 8069 umin-=d; vmin-=d; umax+=d; vmax+=d;
1016
1017 // Create CDT
1018 8069 CDT.clear();
1019 8069 CDT.create_enclosing_rectangle(umin, vmin, umax, vmax);
1020
1021
2/2
✓ Branch 0 taken 51864 times.
✓ Branch 1 taken 8069 times.
 59933 for(index_t v: vertices_) {
1022
2/2
✓ Branch 0 taken 43480 times.
✓ Branch 1 taken 8384 times.
 51864 if(keep_vertex_[v]) {
1023
1/2
✓ Branch 1 taken 43480 times.
✗ Branch 2 not taken.
 43480 ExactPoint P = I_.exact_vertex(v);
1024
1/2
✓ Branch 2 taken 43480 times.
✗ Branch 3 not taken.
 43480 v_idx_[v] = CDT.insert(exact::vec2h(P[u_], P[v_], P.w), v);
1025 43480 } else {
1026 8384 v_idx_[v] = NO_INDEX;
1027 }
1028 }
1029
1030 // Insert constraints
1031
2/2
✓ Branch 0 taken 8171 times.
✓ Branch 1 taken 8069 times.
 16240 for(index_t P: polylines_) {
1032 vector<index_t> Pvertices;
1033 8171 index_t v = polylines_.first_vertex(P);
1034
2/2
✓ Branch 0 taken 7214 times.
✓ Branch 1 taken 957 times.
 8171 if(keep_vertex_[v]) {
1035 Pvertices.push_back(v);
1036 }
1037
3/4
✓ Branch 0 taken 51865 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 51865 times.
✓ Branch 3 taken 8171 times.
 60036 for(index_t h: polylines_.halfedges(P)) {
1038
2/2
✓ Branch 0 taken 43481 times.
✓ Branch 1 taken 8384 times.
 51865 v = halfedges_.vertex(h,1);
1039
2/2
✓ Branch 0 taken 43481 times.
✓ Branch 1 taken 8384 times.
 51865 if(keep_vertex_[v]) {
1040 Pvertices.push_back(v);
1041 }
1042 }
1043 if(
1044
2/4
✓ Branch 0 taken 8171 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8171 times.
✗ Branch 3 not taken.
 16342 polylines_.first_vertex(P) == polylines_.last_vertex(P) &&
1045 Pvertices.size() != 0
1046 ) {
1047 Pvertices.push_back(Pvertices[0]);
1048 }
1049
1050
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 58866 times.
✓ Branch 2 taken 50695 times.
✓ Branch 3 taken 8171 times.
 58866 for(index_t i=0; i+1<Pvertices.size(); ++i) {
1051 50695 index_t v1 = v_idx_[Pvertices[i]];
1052 50695 index_t v2 = v_idx_[Pvertices[i+1]];
1053 geo_debug_assert(v1 != NO_INDEX);
1054 geo_debug_assert(v2 != NO_INDEX);
1055
1/2
✓ Branch 1 taken 50695 times.
✗ Branch 2 not taken.
 50695 CDT.insert_constraint(v1,v2,NO_INDEX);
1056 }
1057 }
1058
1059 8069 CDT.remove_external_triangles(true);
1060 8069 }
1061
1062 153823 bool CoplanarFacets::triangles_are_coplanar(
1063 const vec3& p1, const vec3& p2, const vec3& p3,
1064 const vec3& q1, const vec3& q2, const vec3& q3
1065 ) const {
1066 exact::vec3 N1 = cross(
1067 307646 make_vec3<exact::vec3>(p1,p2), make_vec3<exact::vec3>(p1,p3)
1068
1/2
✓ Branch 1 taken 153823 times.
✗ Branch 2 not taken.
 153823 );
1069 exact::vec3 N2 = cross(
1070
1/2
✓ Branch 4 taken 153823 times.
✗ Branch 5 not taken.
 307646 make_vec3<exact::vec3>(q1,q2), make_vec3<exact::vec3>(q1,q3)
1071
1/2
✓ Branch 1 taken 153823 times.
✗ Branch 2 not taken.
 153823 );
1072
1073
5/6
✓ Branch 0 taken 35742 times.
✓ Branch 1 taken 118081 times.
✓ Branch 2 taken 30325 times.
✓ Branch 3 taken 5417 times.
✓ Branch 4 taken 30325 times.
✗ Branch 5 not taken.
 219890 if(N1.x.sign() == ZERO && N1.y.sign() == ZERO && N1.z.sign() == ZERO) {
1074 std::cerr << std::endl;
1075 std::cerr << "degenerate triangle" << std::endl;
1076 std::cerr << "aligned: " << PCK::aligned_3d(p1,p2,p3) << std::endl;
1077 return false;
1078 }
1079
1080
5/6
✓ Branch 0 taken 36457 times.
✓ Branch 1 taken 117366 times.
✓ Branch 2 taken 30988 times.
✓ Branch 3 taken 5469 times.
✓ Branch 4 taken 30988 times.
✗ Branch 5 not taken.
 221268 if(N2.x.sign() == ZERO && N2.y.sign() == ZERO && N2.z.sign() == ZERO) {
1081 std::cerr << std::endl;
1082 std::cerr << "degenerate triangle" << std::endl;
1083 std::cerr << "aligned: " << PCK::aligned_3d(q1,q2,q3) << std::endl;
1084 return false;
1085 }
1086
1087 // Tolerance for co-planarity test
1088
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 153823 times.
 153823 if(angle_tolerance_ != 0.0) {
1089 ✗ double threshold = cos(angle_tolerance_ * M_PI / 180.0);
1090 ✗ exact::scalar left = geo_sqr(dot(N1,N2));
1091 exact::scalar right =
1092 ✗ exact::scalar(threshold*threshold)*length2(N1)*length2(N2);
1093 return left > right;
1094 }
1095
1096 // Exact version
1097
1/2
✓ Branch 1 taken 153823 times.
✗ Branch 2 not taken.
 153823 exact::vec3 N12 = cross(N1,N2);
1098
6/6
✓ Branch 0 taken 38767 times.
✓ Branch 1 taken 115056 times.
✓ Branch 2 taken 32852 times.
✓ Branch 3 taken 5915 times.
✓ Branch 4 taken 2490 times.
✓ Branch 5 taken 30362 times.
 225442 if((N12.x.sign()!=ZERO) || (N12.y.sign()!=ZERO) ||(N12.z.sign()!=ZERO)) {
1099 123461 return false;
1100 }
1101
1102 return true;
1103 153823 }
1104
1105 /**************************************************************************/
1106
1107 246922 bool CoplanarFacets::edges_are_colinear(
1108 const ExactPoint& P1, const ExactPoint& P2, const ExactPoint& P3
1109 ) const {
1110
1111
1/2
✓ Branch 0 taken 246922 times.
✗ Branch 1 not taken.
 246922 if(angle_tolerance_ == 0.0) {
1112 246922 return PCK::on_segment_3d(P2,P1,P3);
1113 }
1114
1115 ✗ ExactPoint UU = P1-P2;
1116 ✗ exact::vec3 U(UU.x, UU.y, UU.z);
1117 ✗ if(UU.w.sign() == NEGATIVE) {
1118 U.x.negate(); U.y.negate(); U.z.negate();
1119 }
1120 ✗ ExactPoint VV = P3-P2;
1121 ✗ exact::vec3 V(VV.x, VV.y, VV.z);
1122 ✗ if(VV.w.sign() == NEGATIVE) {
1123 V.x.negate(); V.y.negate(); V.z.negate();
1124 }
1125
1126 ✗ double threshold = cos(angle_tolerance_ * M_PI / 180.0);
1127
1128 ✗ exact::scalar left = dot(U,V);
1129
1130 ✗ if(left.sign() == POSITIVE) {
1131 return false;
1132 }
1133
1134 ✗ left = geo_sqr(left);
1135 exact::scalar right =
1136 ✗ exact::scalar(threshold*threshold)*length2(U)*length2(V);
1137
1138 return left > right;
1139 ✗ }
1140
1141 /**************************************************************************/
1142
1143 }
1144