GCC Code Coverage Report

Directory: ./
File: lib/geogram/mesh/mesh_surface_intersection_internal.cpp
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 425 523 81.3%
Functions: 27 32 84.4%
Branches: 271 548 49.5%
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 1365 times.
✓ Branch 1 taken 4319 times.
 5684 if(result.w.sign() == ZERO) {
83 return false;
84 }
85
86
1/2
✓ Branch 1 taken 1365 times.
✗ Branch 2 not taken.
 1365 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 1365 times.
✗ Branch 2 not taken.
 1365 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 1365 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 1365 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 67558 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 67558 exact::vec3 D = make_vec3<exact::vec3>(q1,q2);
123 67558 exact::vec3 E1 = make_vec3<exact::vec3>(p1,p2);
124 67558 exact::vec3 E2 = make_vec3<exact::vec3>(p1,p3);
125 67558 exact::vec3 AO = make_vec3<exact::vec3>(p1,q1);
126
1/2
✓ Branch 1 taken 67558 times.
✗ Branch 2 not taken.
 67558 exact::vec3 N = cross(E1,E2);
127
2/4
✓ Branch 1 taken 67558 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 67558 times.
✗ Branch 5 not taken.
 67558 exact::scalar d = -dot(D,N);
128 geo_debug_assert(d.sign() != ZERO);
129
1/4
✓ Branch 1 taken 67558 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 67558 exact::rational t(dot(AO,N),d);
130
1/2
✓ Branch 1 taken 67558 times.
✗ Branch 2 not taken.
 135116 return mix(t,q1,q2);
131 135116 }
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 286652 MeshInTriangle::ExactPoint MeshInTriangle::Vertex::compute_geometry() {
153 // Case 1: f1 vertex
154
2/2
✓ Branch 1 taken 46515 times.
✓ Branch 2 taken 240137 times.
 286652 if(region_dim(sym.R1) == 0) {
155
1/2
✓ Branch 0 taken 46515 times.
✗ Branch 1 not taken.
 46515 index_t lv = index_t(sym.R1) - index_t(T1_RGN_P0);
156 geo_debug_assert(lv < 3);
157
1/2
✓ Branch 0 taken 46515 times.
✗ Branch 1 not taken.
 93030 mesh_vertex_index = mesh().facets.vertex(sym.f1,lv);
158 vec3 p = mit->mesh_vertex(mesh_vertex_index);
159 46515 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 218482 times.
 240137 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 191062 times.
✓ Branch 2 taken 27420 times.
✓ Branch 4 taken 191062 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 191062 times.
✓ Branch 7 taken 27420 times.
 436964 (region_dim(sym.R1) == 2 || region_dim(sym.R1) == 1) &&
176 218482 region_dim(sym.R2) == 1
177 ) {
178 191062 vec3 p1 = mit->mesh_facet_vertex(sym.f1, 0);
179 191062 vec3 p2 = mit->mesh_facet_vertex(sym.f1, 1);
180 191062 vec3 p3 = mit->mesh_facet_vertex(sym.f1, 2);
181 191062 index_t e = index_t(sym.R2)-index_t(T2_RGN_E0);
182 geo_debug_assert(e<3);
183 191062 vec3 q1 = mit->mesh_facet_vertex(sym.f2, (e+1)%3);
184 191062 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.
 354704 region_dim(sym.R1) == 1 &&
188
3/4
✓ Branch 0 taken 163642 times.
✓ Branch 1 taken 27420 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 150982 times.
 342044 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 40080 return plane_line_intersection(p1,p2,p3,q1,q2);
193 }
194 }
195
196 // case 4: f1 edge /\ f2
197
3/4
✓ Branch 1 taken 178402 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 150982 times.
✓ Branch 5 taken 27420 times.
 178402 if(region_dim(sym.R1) == 1 && region_dim(sym.R2) == 2) {
198 27420 index_t e = index_t(sym.R1)-index_t(T1_RGN_E0);
199 geo_debug_assert(e<3);
200 27420 vec3 p1 = mit->mesh_facet_vertex(sym.f2,0);
201 27420 vec3 p2 = mit->mesh_facet_vertex(sym.f2,1);
202 27420 vec3 p3 = mit->mesh_facet_vertex(sym.f2,2);
203 27420 vec3 q1 = mit->mesh_facet_vertex(sym.f1, (e+1)%3);
204 27420 vec3 q2 = mit->mesh_facet_vertex(sym.f1, (e+2)%3);
205 27420 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 292336 void MeshInTriangle::Vertex::init_geometry(const ExactPoint& P) {
235 292336 point_exact = P;
236 Numeric::optimize_number_representation(point_exact);
237 #ifndef GEOGRAM_USE_EXACT_NT
238
3/8
✓ Branch 3 taken 292336 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 292336 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 292336 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 1461680 l = (geo_sqr(P[mit->u_]) + geo_sqr(P[mit->v_])).estimate() /
239
1/4
✓ Branch 1 taken 292336 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 584672 geo_sqr(P.w).estimate() ;
240 #endif
241 292336 }
242
243 229 MeshInTriangle::MeshInTriangle(MeshSurfaceIntersection& EM) :
244 229 exact_mesh_(EM),
245 229 mesh_(EM.readonly_mesh()),
246 229 f1_(NO_INDEX),
247 229 dry_run_(false),
248 229 use_pred_cache_insert_buffer_(false)
249 {
250 #ifdef GEOGRAM_USE_EXACT_NT
251 CDTBase2d::exact_incircle_ = true;
252 #else
253 // Since incircle() with expansions computes approximated
254 // lifted coordinate, we need to activate additional
255 // checks for Delaunayization.
256 229 CDTBase2d::exact_incircle_ = false;
257 #endif
258 229 }
259
260 15505 void MeshInTriangle::clear() {
261 15505 vertex_.resize(0);
262 15505 edges_.resize(0);
263 15505 f1_ = NO_INDEX;
264 pred_cache_.clear();
265 15505 pred_cache_insert_buffer_.resize(0);
266 15505 use_pred_cache_insert_buffer_ = false;
267 15505 CDTBase2d::clear();
268 15505 }
269
270 15505 void MeshInTriangle::begin_facet(index_t f) {
271 15505 f1_ = f;
272
273 15505 latest_f2_ = NO_INDEX;
274 15505 latest_f2_count_ = 0;
275
276 15505 vec3 p1 = mesh_facet_vertex(f,0);
277 15505 vec3 p2 = mesh_facet_vertex(f,1);
278 15505 vec3 p3 = mesh_facet_vertex(f,2);
279
280 geo_debug_assert(!PCK::aligned_3d(p1,p2,p3));
281
282 15505 f1_normal_axis_ = PCK::triangle_normal_axis(
283 p1,p2,p3
284 );
285
286 15505 u_ = coord_index_t((f1_normal_axis_ + 1) % 3);
287 15505 v_ = coord_index_t((f1_normal_axis_ + 2) % 3);
288
2/2
✓ Branch 0 taken 46515 times.
✓ Branch 1 taken 15505 times.
 62020 for(index_t lv=0; lv<3; ++lv) {
289 93030 vertex_.push_back(Vertex(this, f, lv));
290 }
291
292 15505 CDTBase2d::create_enclosing_triangle(0,1,2);
293
294 15505 edges_.push_back(Edge(1,2));
295 15505 edges_.push_back(Edge(2,0));
296 15505 edges_.push_back(Edge(0,1));
297
298 15505 has_planar_isect_ = false;
299 15505 }
300
301 262922 index_t MeshInTriangle::add_vertex(
302 index_t f2, TriangleRegion R1, TriangleRegion R2
303 ) {
304 geo_debug_assert(f1_ != NO_INDEX);
305
306 // If the same f2 comes more than twice, then
307 // we got a planar facet /\ facet intersection
308 // (and it is good to know it, see get_constraints())
309
3/4
✓ Branch 0 taken 262922 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 153495 times.
✓ Branch 3 taken 109427 times.
 262922 if(f2 != NO_INDEX && f2 == latest_f2_) {
310 153495 ++latest_f2_count_;
311
2/2
✓ Branch 0 taken 36894 times.
✓ Branch 1 taken 116601 times.
 153495 if(latest_f2_count_ > 2) {
312 36894 has_planar_isect_ = true;
313 }
314 } else {
315 109427 latest_f2_ = f2;
316 109427 latest_f2_count_ = 0;
317 }
318
319 // If vertex is a macro-vertex, return it directly.
320
2/2
✓ Branch 1 taken 240137 times.
✓ Branch 2 taken 22785 times.
 262922 if(region_dim(R1) == 0) {
321 return index_t(R1);
322 }
323
324 // Create the vertex
325 240137 vertex_.push_back(Vertex(this, f1_, f2, R1, R2));
326
327 // Insert it into the triangulation
328
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240137 times.
 240137 index_t v = CDTBase2d::insert(vertex_.size()-1);
329
330 // If it was an existing vertex, return the existing vertex
331
2/2
✓ Branch 0 taken 151847 times.
✓ Branch 1 taken 88290 times.
 240137 if(vertex_.size() > CDTBase2d::nv()) {
332 vertex_.pop_back();
333 }
334 return v;
335 }
336
337 127158 void MeshInTriangle::add_edge(
338 index_t f2,
339 TriangleRegion AR1, TriangleRegion AR2,
340 TriangleRegion BR1, TriangleRegion BR2
341 ) {
342 127158 index_t v1 = add_vertex(f2, AR1, AR2);
343 127158 index_t v2 = add_vertex(f2, BR1, BR2);
344
345 // If both extremities are on the same edge of f1,
346 // we do not add the edge, because it will be generated
347 // when remeshing the edge of f1
348
2/2
✓ Branch 2 taken 78130 times.
✓ Branch 3 taken 49028 times.
 127158 if(region_dim(regions_convex_hull(AR1,BR1)) == 1) {
349 return;
350 }
351
352 // Generate also the combinatorial information of the edge,
353 // that indicates whether both extremities are on the same
354 // edge of f2 (useful later to compute the intersections)
355 156260 edges_.push_back(Edge(v1,v2,f2,regions_convex_hull(AR2,BR2)));
356
357 // Constraints will be added to the triangulation during commit()
358 }
359
360 15505 void MeshInTriangle::commit() {
361
362
2/2
✓ Branch 0 taken 124645 times.
✓ Branch 1 taken 15505 times.
 140150 for(const Edge& E: edges_) {
363 124645 CDTBase2d::insert_constraint(E.v1, E.v2);
364 }
365
366
1/2
✓ Branch 0 taken 15505 times.
✗ Branch 1 not taken.
 15505 if(dry_run_) {
367 return;
368 }
369
370 // Protect global mesh from concurrent accesses
371 15505 exact_mesh_.lock();
372
373 // Create vertices and facets in target mesh
374
2/2
✓ Branch 0 taken 140489 times.
✓ Branch 1 taken 15505 times.
 296483 for(index_t i=0; i<vertex_.size(); ++i) {
375 // Vertex already exists in this MeshInTriangle
376
2/2
✓ Branch 0 taken 49252 times.
✓ Branch 1 taken 91237 times.
 140489 if(vertex_[i].mesh_vertex_index != NO_INDEX) {
377 49252 continue;
378 }
379 91237 vertex_[i].mesh_vertex_index =
380 91237 exact_mesh_.find_or_create_exact_vertex(
381 91237 vertex_[i].point_exact
382 );
383 }
384
385 // Create facets in target mesh
386
2/2
✓ Branch 0 taken 131387 times.
✓ Branch 1 taken 15505 times.
 278279 for(index_t t=0; t<CDTBase2d::nT(); ++t) {
387 index_t i = CDTBase2d::Tv(t,0);
388 index_t j = CDTBase2d::Tv(t,1);
389 index_t k = CDTBase2d::Tv(t,2);
390 131387 i = vertex_[i].mesh_vertex_index;
391 131387 j = vertex_[j].mesh_vertex_index;
392 131387 k = vertex_[k].mesh_vertex_index;
393 131387 index_t new_t = target_mesh().facets.create_triangle(i,j,k);
394 // Copy all attributes from initial facet
395 131387 target_mesh().facets.attributes().copy_item(new_t, f1_);
396 }
397
398 // We are done with modification in the mesh
399 15505 exact_mesh_.unlock();
400 }
401
402 ✗ void MeshInTriangle::get_constraints(Mesh& M, bool with_edges) const {
403 ✗ if(M.vertices.nb() == 0) {
404 ✗ M.vertices.set_dimension(2);
405 ✗ for(index_t v=0; v<vertex_.size(); ++v) {
406 ✗ vec2 p = vertex_[v].get_UV_approx();
407 ✗ M.vertices.create_vertex(p.data());
408 }
409 }
410 ✗ if(with_edges && M.edges.nb() == 0) {
411 ✗ for(const Edge& E: edges_) {
412 ✗ M.edges.create_edge(E.v1, E.v2);
413 }
414 }
415 ✗ }
416
417 /**
418 * \brief Tests the parity of the permutation of a list of
419 * three distinct indices with respect to the canonical order.
420 */
421 5044317 static bool odd_order(index_t i, index_t j, index_t k) {
422 // Implementation: sort the elements (bubble sort is OK for
423 // such a small number), and invert parity each time
424 // two elements are swapped.
425 5044317 index_t tab[3] = { i, j, k};
426 const int N = 3;
427 bool result = false;
428
2/2
✓ Branch 0 taken 10088634 times.
✓ Branch 1 taken 5044317 times.
 15132951 for (int I = 0; I < N - 1; ++I) {
429
2/2
✓ Branch 0 taken 15132951 times.
✓ Branch 1 taken 10088634 times.
 25221585 for (int J = 0; J < N - I - 1; ++J) {
430
2/2
✓ Branch 0 taken 8981432 times.
✓ Branch 1 taken 6151519 times.
 15132951 if (tab[J] > tab[J + 1]) {
431 std::swap(tab[J], tab[J + 1]);
432 8981432 result = !result;
433 }
434 }
435 }
436 5044317 return result;
437 }
438
439 240137 void MeshInTriangle::begin_insert_transaction() {
440 240137 use_pred_cache_insert_buffer_ = true;
441 240137 }
442
443 88290 void MeshInTriangle::commit_insert_transaction() {
444
2/2
✓ Branch 0 taken 893540 times.
✓ Branch 1 taken 88290 times.
 981830 for(const auto& it: pred_cache_insert_buffer_) {
445 893540 pred_cache_[it.first] = it.second;
446 }
447 88290 pred_cache_insert_buffer_.resize(0);
448 88290 use_pred_cache_insert_buffer_ = false;
449 88290 }
450
451 151847 void MeshInTriangle::rollback_insert_transaction() {
452 151847 pred_cache_insert_buffer_.resize(0);
453 151847 use_pred_cache_insert_buffer_ = false;
454 151847 }
455
456
2/2
✓ Branch 0 taken 3227616 times.
✓ Branch 1 taken 1816701 times.
 5044317 Sign MeshInTriangle::orient2d(index_t vx1,index_t vx2,index_t vx3) const {
457
458 trindex K(vx1, vx2, vx3);
459
460
2/2
✓ Branch 0 taken 3538575 times.
✓ Branch 1 taken 1505742 times.
 5044317 if(use_pred_cache_insert_buffer_) {
461 3538575 Sign result = PCK::orient_2d_projected(
462 3538575 vertex_[K.indices[0]].point_exact,
463 3538575 vertex_[K.indices[1]].point_exact,
464 3538575 vertex_[K.indices[2]].point_exact,
465 3538575 f1_normal_axis_
466 );
467 3538575 pred_cache_insert_buffer_.push_back(std::make_pair(K, result));
468
2/2
✓ Branch 0 taken 1743667 times.
✓ Branch 1 taken 1794908 times.
 3538575 if(odd_order(vx1,vx2,vx3)) {
469 1743667 result = Sign(-result);
470 }
471 return result;
472 }
473
474 bool inserted;
475 std::map<trindex, Sign>::iterator it;
476 1505742 std::tie(it,inserted) = pred_cache_.insert(std::make_pair(K,ZERO));
477 Sign result;
478
479
2/2
✓ Branch 0 taken 642212 times.
✓ Branch 1 taken 863530 times.
 1505742 if(inserted) {
480 642212 result = PCK::orient_2d_projected(
481 642212 vertex_[K.indices[0]].point_exact,
482 642212 vertex_[K.indices[1]].point_exact,
483 642212 vertex_[K.indices[2]].point_exact,
484 642212 f1_normal_axis_
485 );
486 642212 it->second = result;
487 } else {
488 863530 result = it->second;
489 }
490
491
2/2
✓ Branch 0 taken 693255 times.
✓ Branch 1 taken 812487 times.
 1505742 if(odd_order(vx1,vx2,vx3)) {
492 693255 result = Sign(-result);
493 }
494
495 return result;
496 }
497
498 187443 Sign MeshInTriangle::incircle(
499 index_t v1,index_t v2,index_t v3,index_t v4
500 ) const {
501 exact::vec2h p1(
502 187443 vertex_[v1].point_exact[u_],
503 187443 vertex_[v1].point_exact[v_],
504 187443 vertex_[v1].point_exact.w
505 187443 );
506 exact::vec2h p2(
507 vertex_[v2].point_exact[u_],
508 vertex_[v2].point_exact[v_],
509 187443 vertex_[v2].point_exact.w
510 187443 );
511 exact::vec2h p3(
512 vertex_[v3].point_exact[u_],
513 vertex_[v3].point_exact[v_],
514 187443 vertex_[v3].point_exact.w
515 187443 );
516 exact::vec2h p4(
517 vertex_[v4].point_exact[u_],
518 vertex_[v4].point_exact[v_],
519 187443 vertex_[v4].point_exact.w
520 187443 );
521 #ifdef GEOGRAM_USE_EXACT_NT
522 return PCK::incircle_2d_SOS(p1,p2,p3,p4);
523 #else
524 187443 return PCK::incircle_2d_SOS_with_lengths(
525 p1,p2,p3,p4,
526 187443 vertex_[v1].l,
527 187443 vertex_[v2].l,
528 187443 vertex_[v3].l,
529
1/2
✓ Branch 1 taken 187443 times.
✗ Branch 2 not taken.
 187443 vertex_[v4].l
530 187443 );
531 #endif
532 187443 }
533
534
1/2
✓ Branch 1 taken 5684 times.
✗ Branch 2 not taken.
 5684 index_t MeshInTriangle::create_intersection(
535 index_t e1, index_t i, index_t j,
536 index_t e2, index_t k, index_t l
537 ) {
538 geo_argused(i);
539 geo_argused(j);
540 geo_argused(k);
541 geo_argused(l);
542
543 ExactPoint I;
544
1/2
✓ Branch 1 taken 5684 times.
✗ Branch 2 not taken.
 5684 get_edge_edge_intersection(e1,e2,I);
545
2/4
✓ Branch 1 taken 5684 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5684 times.
 11368 vertex_.push_back(Vertex(this,I));
546 5684 index_t x = vertex_.size()-1;
547
2/2
✓ Branch 0 taken 5660 times.
✓ Branch 1 taken 24 times.
 5684 CDTBase2d::v2T_.push_back(NO_INDEX);
548 geo_debug_assert(x == CDTBase2d::nv_);
549 5684 ++CDTBase2d::nv_;
550 5684 return x;
551 5684 }
552
553 5684 void MeshInTriangle::get_edge_edge_intersection(
554 index_t e1, index_t e2, ExactPoint& I
555 ) const {
556 5684 index_t f1 = f1_;
557 5684 index_t f2 = edges_[e1].sym.f2;
558 5684 index_t f3 = edges_[e2].sym.f2;
559
560 geo_debug_assert(f1 != NO_INDEX);
561 geo_debug_assert(f2 != NO_INDEX);
562 geo_debug_assert(f3 != NO_INDEX);
563
564 vec3 P[9] = {
565 5684 mesh_facet_vertex(f1,0), mesh_facet_vertex(f1,1),
566 5684 mesh_facet_vertex(f1,2),
567 5684 mesh_facet_vertex(f2,0), mesh_facet_vertex(f2,1),
568 5684 mesh_facet_vertex(f2,2),
569 5684 mesh_facet_vertex(f3,0), mesh_facet_vertex(f3,1),
570 5684 mesh_facet_vertex(f3,2)
571 };
572
573
2/2
✓ Branch 1 taken 4319 times.
✓ Branch 2 taken 1365 times.
 5684 if(!get_three_planes_intersection(
574 I,
575 P[0], P[1], P[2],
576 P[3], P[4], P[5],
577 P[6], P[7], P[8]
578 )) {
579 4319 get_edge_edge_intersection_2D(e1,e2,I);
580 4319 return;
581 }
582 }
583
584 4319 void MeshInTriangle::get_edge_edge_intersection_2D(
585 index_t e1, index_t e2, ExactPoint& I
586 ) const {
587 const Edge& E1 = edges_[e1];
588 const Edge& E2 = edges_[e2];
589
590 if(
591
4/4
✓ Branch 1 taken 4289 times.
✓ Branch 2 taken 30 times.
✓ Branch 3 taken 4261 times.
✓ Branch 4 taken 28 times.
 8608 region_dim(E1.sym.R2) == 1 &&
592 4289 region_dim(E2.sym.R2) == 1
593 ) {
594 4261 index_t le1 = index_t(E1.sym.R2)-index_t(T2_RGN_E0);
595 4261 index_t le2 = index_t(E2.sym.R2)-index_t(T2_RGN_E0);
596 geo_debug_assert(le1 < 3);
597 geo_debug_assert(le2 < 3);
598
599 4261 vec2 p1_uv = mesh_facet_vertex_UV(E1.sym.f2, (le1+1)%3);
600 4261 vec2 p2_uv = mesh_facet_vertex_UV(E1.sym.f2, (le1+2)%3);
601 4261 vec2 q1_uv = mesh_facet_vertex_UV(E2.sym.f2, (le2+1)%3);
602 4261 vec2 q2_uv = mesh_facet_vertex_UV(E2.sym.f2, (le2+2)%3);
603
604 4261 exact::vec2 C1 = make_vec2<exact::vec2>(p1_uv, p2_uv);
605 4261 exact::vec2 C2 = make_vec2<exact::vec2>(q2_uv, q1_uv);
606 4261 exact::vec2 B = make_vec2<exact::vec2>(p1_uv, q1_uv);
607
1/2
✓ Branch 1 taken 4261 times.
✗ Branch 2 not taken.
 4261 exact::scalar d = det(C1,C2);
608 geo_debug_assert(d.sign() != ZERO);
609
1/4
✓ Branch 1 taken 4261 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 4261 exact::rational t(det(B,C2),d);
610 4261 I = mix(
611 t,
612 8522 mesh_facet_vertex(E1.sym.f2,(le1+1)%3),
613
1/2
✓ Branch 1 taken 4261 times.
✗ Branch 2 not taken.
 4261 mesh_facet_vertex(E1.sym.f2,(le1+2)%3)
614 4261 );
615
616 8522 } else {
617 geo_debug_assert(
618 region_dim(E1.sym.R2) == 1 || region_dim(E2.sym.R2) == 1
619 );
620 58 index_t f1 = E1.sym.f2;
621 58 TriangleRegion R1 = E1.sym.R2;
622 58 index_t f2 = E2.sym.f2;
623 58 TriangleRegion R2 = E2.sym.R2;
624
2/2
✓ Branch 1 taken 28 times.
✓ Branch 2 taken 30 times.
 58 if(region_dim(R1) == 1) {
625 std::swap(f1,f2);
626 std::swap(R1,R2);
627 }
628
629 index_t e = index_t(R2) - index_t(T2_RGN_E0);
630 geo_debug_assert(e < 3);
631
632 58 I = plane_line_intersection(
633 116 mesh_facet_vertex(f1,0),
634 116 mesh_facet_vertex(f1,1),
635 116 mesh_facet_vertex(f1,2),
636 116 mesh_facet_vertex(f2,(e+1)%3),
637 116 mesh_facet_vertex(f2,(e+2)%3)
638 58 );
639 }
640 4319 }
641
642 ✗ void MeshInTriangle::save(const std::string& filename) const {
643 ✗ Mesh M;
644 ✗ M.vertices.set_dimension(2);
645 ✗ for(index_t v=0; v<CDTBase2d::nv(); ++v) {
646 ✗ vec2 p = vertex_[v].get_UV_approx();
647 ✗ M.vertices.create_vertex(p.data());
648 }
649 ✗ for(index_t t=0; t<CDTBase2d::nT(); ++t) {
650 ✗ M.facets.create_triangle(
651 CDTBase2d::Tv(t,0),
652 CDTBase2d::Tv(t,1),
653 CDTBase2d::Tv(t,2)
654 );
655 }
656
657 Attribute<double> tex_coord;
658 ✗ tex_coord.create_vector_attribute(
659 ✗ M.facet_corners.attributes(), "tex_coord", 2
660 );
661 static double triangle_tex[3][2] = {
662 {0.0, 0.0},
663 {1.0, 0.0},
664 {0.0, 1.0}
665 };
666 ✗ for(index_t c: M.facet_corners) {
667 ✗ tex_coord[2*c] = triangle_tex[c%3][0];
668 ✗ tex_coord[2*c+1] = triangle_tex[c%3][1];
669 }
670 ✗ mesh_save(M, filename);
671 ✗ }
672
673 /**************************************************************************/
674
675 285 CoplanarFacets::CoplanarFacets(
676 MeshSurfaceIntersection& I, bool clear_attributes,
677 double angle_tolerance
678 285 ) :
679
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 I_(I),
680
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 mesh_(I.target_mesh()),
681 285 mesh_copy_(I.readonly_mesh()),
682 285 angle_tolerance_(angle_tolerance),
683
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 facet_group_(I.target_mesh().facets.attributes(),"group"),
684
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"),
685 c_is_coplanar_(
686
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"
687 ),
688
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"),
689 halfedges_(*this),
690 285 polylines_(*this)
691 {
692
2/2
✓ Branch 0 taken 57 times.
✓ Branch 1 taken 228 times.
 285 if(clear_attributes) {
693
2/2
✓ Branch 0 taken 123492 times.
✓ Branch 1 taken 57 times.
 123549 for(index_t f: mesh_.facets) {
694 123492 facet_group_[f] = NO_INDEX;
695 }
696
2/2
✓ Branch 0 taken 61752 times.
✓ Branch 1 taken 57 times.
 61809 for(index_t v: mesh_.vertices) {
697 61752 keep_vertex_[v] = false;
698 }
699
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 find_coplanar_facets();
700 }
701
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 f_visited_.assign(mesh_.facets.nb(),false);
702
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 h_visited_.assign(mesh_.facet_corners.nb(),false);
703
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 v_visited_.assign(mesh_.vertices.nb(),false);
704
1/2
✓ Branch 1 taken 285 times.
✗ Branch 2 not taken.
 285 v_idx_.assign(mesh_.vertices.nb(),NO_INDEX);
705 285 }
706
707 57 void CoplanarFacets::find_coplanar_facets() {
708
709 // Positioned by MeshSurfaceIntersection::build_Weiler_model()
710 Attribute<bool> corner_is_on_border(
711
1/2
✓ Branch 2 taken 57 times.
✗ Branch 3 not taken.
 57 mesh_.facet_corners.attributes(), "is_on_border"
712 );
713
714
2/2
✓ Branch 0 taken 370476 times.
✓ Branch 1 taken 57 times.
 370533 for(index_t c: mesh_.facet_corners) {
715 370476 c_is_coplanar_[c] = false;
716 }
717
718 // TODO: when there is an angle tolerance, one should check instead
719 // angle deviation w.r.t. a single seed facet per facet group, because
720 // with the present algorithm, if a large number of tiny facets are
721 // connected (e.g. highly tessellated cylinder), one may group facets
722 // with large angle deviation (without seeing it because each facet has
723 // small angle deviation w.r.t. its neighbors).
724
725
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 parallel_for(
726
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 0, mesh_.facet_corners.nb(),
727
1/4
✓ Branch 0 taken 57 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 57 [&](index_t c1) {
728 370476 index_t f1 = (c1 / 3);
729 370476 index_t le1 = (c1 % 3);
730
1/2
✓ Branch 0 taken 370476 times.
✗ Branch 1 not taken.
 370476 index_t f2 = mesh_.facet_corners.adjacent_facet(c1);
731
732
1/2
✓ Branch 0 taken 370476 times.
✗ Branch 1 not taken.
 370476 if(f2 == NO_INDEX) {
733 216544 return;
734 }
735
736 // do not traverse true borders
737
1/2
✓ Branch 0 taken 370476 times.
✗ Branch 1 not taken.
 370476 if(corner_is_on_border[c1]) {
738 return;
739 }
740
741 index_t v11 = mesh_.facets.vertex(f1,le1);
742 370476 index_t v12 = mesh_.facets.vertex(f1,(le1+1)%3);
743
1/2
✓ Branch 0 taken 370476 times.
✗ Branch 1 not taken.
 370476 index_t le2 = mesh_.facets.find_edge(f2,v12,v11);
744 index_t c2 = mesh_.facets.corner(f2,le2);
745
746 #ifdef GEO_DEBUG
747 index_t v21 = mesh_.facets.vertex(f2,le2);
748 index_t v22 = mesh_.facets.vertex(f2,(le2+1)%3);
749 index_t v13 = mesh_.facets.vertex(f1,(le1+2)%3);
750 index_t v23 = mesh_.facets.vertex(f2,(le2+2)%3);
751
752 geo_debug_assert(v11 == v22);
753 geo_debug_assert(v12 == v21);
754 geo_debug_assert(v11!=v12 && v12!=v13 && v13!=v11);
755 geo_debug_assert(v21!=v22 && v22!=v23 && v23!=v21);
756 #endif
757
758
2/2
✓ Branch 0 taken 185238 times.
✓ Branch 1 taken 185238 times.
 370476 if(c1 > c2) {
759 return;
760 }
761
762 // Small optimization: if both triangles come from same
763 // original facet then they are coplanar
764
2/2
✓ Branch 0 taken 31306 times.
✓ Branch 1 taken 153932 times.
 185238 if(I_.get_initial_facet(f1) == I_.get_initial_facet(f2)) {
765 31306 c_is_coplanar_[c1] = true;
766 31306 c_is_coplanar_[c2] = true;
767 31306 return;
768 }
769
770 // Use original triangles for co-planarity test
771 153932 auto [p1, p2, p3] = I_.get_initial_facet_vertices(f1);
772 153932 auto [q1, q2, q3] = I_.get_initial_facet_vertices(f2);
773
2/2
✓ Branch 1 taken 30473 times.
✓ Branch 2 taken 123459 times.
 153932 if(triangles_are_coplanar(p1,p2,p3,q1,q2,q3)) {
774 30473 c_is_coplanar_[c1] = true;
775 30473 c_is_coplanar_[c2] = true;
776 }
777 }
778 );
779 57 }
780
781 146174 void CoplanarFacets::get(index_t f, index_t group_id) {
782
783 146174 facets_.resize(0);
784 146174 vertices_.resize(0);
785 146174 halfedges_.initialize();
786 146174 polylines_.initialize();
787
788 // Get facets
789 {
790 std::stack<index_t> S;
791
1/2
✓ Branch 0 taken 146174 times.
✗ Branch 1 not taken.
 146174 facet_group_[f] = group_id;
792 f_visited_[f] = true;
793 S.push(f);
794 facets_.push_back(f);
795
2/2
✓ Branch 0 taken 246984 times.
✓ Branch 1 taken 146174 times.
 393158 while(!S.empty()) {
796
2/2
✓ Branch 0 taken 246642 times.
✓ Branch 1 taken 342 times.
 246984 index_t f1 = S.top();
797 S.pop();
798
2/2
✓ Branch 0 taken 740952 times.
✓ Branch 1 taken 246984 times.
 987936 for(index_t le1=0; le1<3; ++le1) {
799
1/2
✓ Branch 0 taken 740952 times.
✗ Branch 1 not taken.
 740952 index_t f2 = mesh_.facets.adjacent(f1,le1);
800 if(
801
5/6
✓ Branch 0 taken 740952 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 594646 times.
✓ Branch 3 taken 146306 times.
✓ Branch 4 taken 100810 times.
✓ Branch 5 taken 493836 times.
 740952 f2 != NO_INDEX && !f_visited_[f2] &&
802
2/2
✓ Branch 0 taken 100810 times.
✓ Branch 1 taken 493836 times.
 594646 c_is_coplanar_[mesh_.facets.corner(f1,le1)]
803 ) {
804
2/2
✓ Branch 0 taken 100468 times.
✓ Branch 1 taken 342 times.
 100810 facet_group_[f2] = facet_group_[f1];
805 f_visited_[f2] = true;
806 S.push(f2);
807 facets_.push_back(f2);
808 }
809 }
810 }
811
2/2
✓ Branch 0 taken 246984 times.
✓ Branch 1 taken 146174 times.
 393158 for(index_t cur_f: facets_) {
812 246984 f_visited_[cur_f] = false;
813 }
814 }
815
816 146174 group_id_ = group_id;
817
818 // Initialize projection coordinates, using original facet
819 {
820 146174 auto [p1, p2, p3] = I_.get_initial_facet_vertices(facets_[0]);
821 146174 coord_index_t projection_axis = PCK::triangle_normal_axis(p1,p2,p3);
822 146174 u_ = coord_index_t((projection_axis+1)%3);
823 146174 v_ = coord_index_t((projection_axis+2)%3);
824 Sign o = PCK::orient_2d(
825 146174 vec2(p1[u_],p1[v_]), vec2(p2[u_],p2[v_]), vec2(p3[u_],p3[v_])
826 );
827 geo_debug_assert(o != ZERO);
828
2/2
✓ Branch 0 taken 71086 times.
✓ Branch 1 taken 75088 times.
 146174 if(o < 0) {
829 std::swap(u_,v_);
830 }
831 }
832
833 // Get vertices and halfedges
834 {
835
2/2
✓ Branch 0 taken 246984 times.
✓ Branch 1 taken 146174 times.
 393158 for(index_t f1: facets_) {
836
2/2
✓ Branch 0 taken 740952 times.
✓ Branch 1 taken 246984 times.
 987936 for(index_t le=0; le<3; ++le) {
837
1/2
✓ Branch 0 taken 740952 times.
✗ Branch 1 not taken.
 740952 index_t f2 = mesh_.facets.adjacent(f1,le);
838 if(
839
3/4
✓ Branch 0 taken 740952 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 247116 times.
✓ Branch 3 taken 493836 times.
 740952 f2 == NO_INDEX ||
840
2/2
✓ Branch 0 taken 247116 times.
✓ Branch 1 taken 493836 times.
 740952 !c_is_coplanar_[mesh_.facets.corner(f1,le)]
841 ) {
842 493836 halfedges_.add(mesh_.facets.corners_begin(f1)+le);
843
1/2
✓ Branch 0 taken 493836 times.
✗ Branch 1 not taken.
 493836 index_t v1 = mesh_.facets.vertex(f1,le);
844 493836 index_t v2 = mesh_.facets.vertex(f1,(le+1)%3);
845
2/2
✓ Branch 0 taken 175015 times.
✓ Branch 1 taken 318821 times.
 493836 if(!v_visited_[v1]) {
846
2/2
✓ Branch 0 taken 174419 times.
✓ Branch 1 taken 596 times.
 175015 v_idx_[v1] = vertices_.size();
847 vertices_.push_back(v1);
848 v_visited_[v1] = true;
849 }
850
2/2
✓ Branch 0 taken 307277 times.
✓ Branch 1 taken 186559 times.
 493836 if(!v_visited_[v2]) {
851
2/2
✓ Branch 0 taken 306593 times.
✓ Branch 1 taken 684 times.
 307277 v_idx_[v2] = vertices_.size();
852 vertices_.push_back(v2);
853 v_visited_[v2] = true;
854 }
855 } else {
856 // This one for the particular case of a non-manifold
857 // vertex, such as a cone apex touching a facet
858 // (ThingiCSG/Basic/cube_cone_1.scad)
859 247116 index_t v = mesh_.facets.vertex(f1,le);
860
4/4
✓ Branch 0 taken 71199 times.
✓ Branch 1 taken 175917 times.
✓ Branch 2 taken 11540 times.
✓ Branch 3 taken 59659 times.
 247116 if(keep_vertex_[v] && !v_visited_[v]) {
861 vertices_.push_back(v);
862 v_visited_[v] = true;
863 }
864 }
865 }
866 }
867
2/2
✓ Branch 0 taken 493832 times.
✓ Branch 1 taken 146174 times.
 640006 for(index_t v: vertices_) {
868 493832 v_visited_[v] = false;
869 }
870 }
871
872 // Get polylines
873 {
874 // Get all polylines starting from vertices with more than
875 // 2 incident halfedges.
876
2/2
✓ Branch 0 taken 493832 times.
✓ Branch 1 taken 146174 times.
 640006 for(index_t v: vertices_) {
877
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 493828 times.
 493832 if(halfedges_.nb_halfedges_around_vertex(v) > 1) {
878 for(
879 index_t h=halfedges_.vertex_first_halfedge(v);
880
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 4 times.
 12 h != NO_INDEX; h = halfedges_.next_around_vertex(h)
881 ) {
882
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
 8 if(!h_visited_[h]) {
883 polylines_.begin_polyline();
884 index_t h2 = h;
885 do {
886 geo_debug_assert(!h_visited_[h2]);
887
1/2
✓ Branch 0 taken 76 times.
✗ Branch 1 not taken.
 76 h_visited_[h2] = true;
888 76 polylines_.add_halfedge(h2);
889 76 h2 = halfedges_.next_along_polyline(h2);
890
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 8 times.
 76 } while(h2 != NO_INDEX && h2 != h);
891 polylines_.end_polyline();
892 }
893 }
894 }
895 }
896 // There can be also closed halfedge loops with no irregular vertex
897
2/2
✓ Branch 0 taken 493836 times.
✓ Branch 1 taken 146174 times.
 640010 for(index_t h: halfedges_) {
898
2/2
✓ Branch 0 taken 146372 times.
✓ Branch 1 taken 347464 times.
 493836 if(!h_visited_[h]) {
899 polylines_.begin_polyline();
900 index_t h2 = h;
901 do {
902 geo_debug_assert(!h_visited_[h2]);
903
2/2
✓ Branch 0 taken 492318 times.
✓ Branch 1 taken 1442 times.
 493760 h_visited_[h2] = true;
904 493760 polylines_.add_halfedge(h2);
905 493760 h2 = halfedges_.next_along_polyline(h2);
906 geo_debug_assert(h2 != NO_INDEX);
907
2/2
✓ Branch 0 taken 347388 times.
✓ Branch 1 taken 146372 times.
 493760 } while(h2 != h);
908 polylines_.end_polyline();
909 }
910 }
911
2/2
✓ Branch 0 taken 493836 times.
✓ Branch 1 taken 146174 times.
 640010 for(index_t h: halfedges_) {
912 493836 h_visited_[h] = false;
913 }
914 }
915 146174 }
916
917
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 73087 times.
 73087 void CoplanarFacets::mark_vertices_to_keep() {
918
2/2
✓ Branch 0 taken 73190 times.
✓ Branch 1 taken 73087 times.
 146277 for(index_t P: polylines_) {
919 73190 index_t first_v = polylines_.first_vertex(P);
920 73190 index_t last_v = polylines_.last_vertex(P);
921
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 73190 times.
 73190 if(first_v != last_v) {
922 ✗ keep_vertex_[first_v] = true;
923 keep_vertex_[last_v] = true;
924 }
925 73190 index_t v1 = polylines_.prev_first_vertex(P);
926 index_t v2 = NO_INDEX;
927 index_t v3 = NO_INDEX;
928
2/2
✓ Branch 0 taken 246918 times.
✓ Branch 1 taken 73190 times.
 320108 for(index_t h: polylines_.halfedges(P)) {
929
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 246918 times.
 246918 if(v1 == NO_INDEX) {
930 ✗ continue;
931 }
932 v2 = halfedges_.vertex(h,0);
933 v3 = halfedges_.vertex(h,1);
934
1/2
✓ Branch 1 taken 246918 times.
✗ Branch 2 not taken.
 246918 ExactPoint p1 = I_.exact_vertex(v1);
935
1/2
✓ Branch 1 taken 246918 times.
✗ Branch 2 not taken.
 246918 ExactPoint p2 = I_.exact_vertex(v2);
936
1/2
✓ Branch 1 taken 246918 times.
✗ Branch 2 not taken.
 246918 ExactPoint p3 = I_.exact_vertex(v3);
937
3/4
✓ Branch 1 taken 246918 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 238479 times.
✓ Branch 4 taken 8439 times.
 246918 if(!edges_are_colinear(p1,p2,p3)) {
938 238479 keep_vertex_[v2] = true;
939 }
940 v1 = v2;
941 246918 }
942 }
943 73087 }
944
945 ✗ void CoplanarFacets::save_borders(const std::string& filename) {
946 ✗ Mesh borders;
947 ✗ borders.vertices.set_dimension(2);
948 index_t cur_idx = 0;
949 ✗ for(index_t v: vertices_) {
950 ✗ vec3 p = mesh_.vertices.point(v);
951 ✗ vec2 q(p[u_],p[v_]);
952 ✗ borders.vertices.create_vertex(q.data());
953 ✗ v_idx_[v] = cur_idx;
954 ✗ ++cur_idx;
955 }
956
957 ✗ for(index_t h: halfedges_) {
958 index_t v1 = halfedges_.vertex(h,0);
959 index_t v2 = halfedges_.vertex(h,1);
960 ✗ v1 = v_idx_[v1];
961 ✗ v2 = v_idx_[v2];
962 geo_debug_assert(v1 != NO_INDEX);
963 geo_debug_assert(v2 != NO_INDEX);
964 ✗ borders.edges.create_edge(v1,v2);
965 }
966
967 ✗ Attribute<bool> selection(borders.vertices.attributes(), "selection");
968 ✗ for(index_t v: vertices_) {
969 geo_debug_assert(v_idx_[v] != NO_INDEX);
970 ✗ selection[v_idx_[v]] = keep_vertex_[v];
971 }
972 ✗ mesh_save(borders,filename);
973 ✗ }
974
975 ✗ void CoplanarFacets::save_facet_group(const std::string& filename) {
976 ✗ Mesh M;
977 ✗ Attribute<bool> keep_vertex(M.vertices.attributes(),"keep");
978 ✗ M.vertices.set_dimension(2);
979 ✗ for(index_t f: facets_) {
980 ✗ for(index_t lv=0; lv<3; ++lv) {
981 ✗ index_t v = mesh_.facets.vertex(f,lv);
982 ✗ v_idx_[v] = NO_INDEX;
983 }
984 }
985 ✗ for(index_t f: facets_) {
986 ✗ for(index_t lv=0; lv<3; ++lv) {
987 ✗ index_t v = mesh_.facets.vertex(f,lv);
988 ✗ if(v_idx_[v] == NO_INDEX) {
989 ✗ vec3 p = mesh_.vertices.point(v);
990 ✗ vec2 q(p[u_], p[v_]);
991 ✗ v_idx_[v] = M.vertices.create_vertex(q.data());
992 ✗ keep_vertex[v_idx_[v]] = keep_vertex_[v];
993 }
994 }
995 ✗ M.facets.create_triangle(
996 v_idx_[mesh_.facets.vertex(f,0)],
997 v_idx_[mesh_.facets.vertex(f,1)],
998 ✗ v_idx_[mesh_.facets.vertex(f,2)]
999 );
1000 }
1001
1002 ✗ for(index_t f: facets_) {
1003 ✗ for(index_t lv=0; lv<3; ++lv) {
1004 ✗ index_t v = mesh_.facets.vertex(f,lv);
1005 ✗ v_idx_[v] = NO_INDEX;
1006 }
1007 }
1008
1009 ✗ M.facets.connect();
1010 ✗ mesh_save(M,filename);
1011 ✗ }
1012
1013 8067 void CoplanarFacets::triangulate() {
1014
1015 // Compute 2D projected BBOX
1016 8067 double umin = Numeric::max_float64();
1017 8067 double vmin = Numeric::max_float64();
1018 8067 double umax = -Numeric::max_float64();
1019 8067 double vmax = -Numeric::max_float64();
1020
2/2
✓ Branch 0 taken 58472 times.
✓ Branch 1 taken 8067 times.
 66539 for(index_t f: facets_) {
1021
2/2
✓ Branch 0 taken 175416 times.
✓ Branch 1 taken 58472 times.
 233888 for(index_t lv=0; lv<3; ++lv) {
1022
1/2
✓ Branch 0 taken 175416 times.
✗ Branch 1 not taken.
 175416 index_t vx = mesh_.facets.vertex(f,lv);
1023 175416 double u = mesh_.vertices.point(vx)[u_];
1024 175416 double v = mesh_.vertices.point(vx)[v_];
1025 175416 umin = std::min(umin, u);
1026 175416 umax = std::max(umax, u);
1027 175416 vmin = std::min(vmin, v);
1028 175416 vmax = std::max(vmax, v);
1029 }
1030 }
1031 8067 double d = std::max(umax-umin, vmax-vmin);
1032 8067 d *= 10.0;
1033 d = std::max(d, 1.0);
1034 8067 umin-=d; vmin-=d; umax+=d; vmax+=d;
1035
1036 // Create CDT
1037 8067 CDT.clear();
1038 8067 CDT.create_enclosing_rectangle(umin, vmin, umax, vmax);
1039
1040
2/2
✓ Branch 0 taken 51856 times.
✓ Branch 1 taken 8067 times.
 59923 for(index_t v: vertices_) {
1041
2/2
✓ Branch 0 taken 43476 times.
✓ Branch 1 taken 8380 times.
 51856 if(keep_vertex_[v]) {
1042
1/2
✓ Branch 1 taken 43476 times.
✗ Branch 2 not taken.
 43476 ExactPoint P = I_.exact_vertex(v);
1043
1/2
✓ Branch 2 taken 43476 times.
✗ Branch 3 not taken.
 43476 v_idx_[v] = CDT.insert(exact::vec2h(P[u_], P[v_], P.w), v);
1044 43476 } else {
1045 8380 v_idx_[v] = NO_INDEX;
1046 }
1047 }
1048
1049 // Insert constraints
1050
2/2
✓ Branch 0 taken 8170 times.
✓ Branch 1 taken 8067 times.
 16237 for(index_t P: polylines_) {
1051 vector<index_t> Pvertices;
1052 8170 index_t v = polylines_.first_vertex(P);
1053
2/2
✓ Branch 0 taken 7214 times.
✓ Branch 1 taken 956 times.
 8170 if(keep_vertex_[v]) {
1054 Pvertices.push_back(v);
1055 }
1056
3/4
✓ Branch 0 taken 51858 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 51858 times.
✓ Branch 3 taken 8170 times.
 60028 for(index_t h: polylines_.halfedges(P)) {
1057
2/2
✓ Branch 0 taken 43478 times.
✓ Branch 1 taken 8380 times.
 51858 v = halfedges_.vertex(h,1);
1058
2/2
✓ Branch 0 taken 43478 times.
✓ Branch 1 taken 8380 times.
 51858 if(keep_vertex_[v]) {
1059 Pvertices.push_back(v);
1060 }
1061 }
1062 if(
1063
2/4
✓ Branch 0 taken 8170 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8170 times.
✗ Branch 3 not taken.
 16340 polylines_.first_vertex(P) == polylines_.last_vertex(P) &&
1064 Pvertices.size() != 0
1065 ) {
1066 Pvertices.push_back(Pvertices[0]);
1067 }
1068
1069
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 58862 times.
✓ Branch 2 taken 50692 times.
✓ Branch 3 taken 8170 times.
 58862 for(index_t i=0; i+1<Pvertices.size(); ++i) {
1070 50692 index_t v1 = v_idx_[Pvertices[i]];
1071 50692 index_t v2 = v_idx_[Pvertices[i+1]];
1072 geo_debug_assert(v1 != NO_INDEX);
1073 geo_debug_assert(v2 != NO_INDEX);
1074
1/2
✓ Branch 1 taken 50692 times.
✗ Branch 2 not taken.
 50692 CDT.insert_constraint(v1,v2,NO_INDEX);
1075 }
1076 }
1077
1078 8067 CDT.remove_external_triangles(true);
1079 8067 }
1080
1081 153932 bool CoplanarFacets::triangles_are_coplanar(
1082 const vec3& p1, const vec3& p2, const vec3& p3,
1083 const vec3& q1, const vec3& q2, const vec3& q3
1084 ) const {
1085 exact::vec3 N1 = cross(
1086 307864 make_vec3<exact::vec3>(p1,p2), make_vec3<exact::vec3>(p1,p3)
1087
1/2
✓ Branch 1 taken 153932 times.
✗ Branch 2 not taken.
 153932 );
1088 exact::vec3 N2 = cross(
1089
1/2
✓ Branch 4 taken 153932 times.
✗ Branch 5 not taken.
 307864 make_vec3<exact::vec3>(q1,q2), make_vec3<exact::vec3>(q1,q3)
1090
1/2
✓ Branch 1 taken 153932 times.
✗ Branch 2 not taken.
 153932 );
1091
1092
5/6
✓ Branch 0 taken 35851 times.
✓ Branch 1 taken 118081 times.
✓ Branch 2 taken 30418 times.
✓ Branch 3 taken 5433 times.
✓ Branch 4 taken 30418 times.
✗ Branch 5 not taken.
 220201 if(N1.x.sign() == ZERO && N1.y.sign() == ZERO && N1.z.sign() == ZERO) {
1093 std::cerr << std::endl;
1094 std::cerr << "degenerate triangle" << std::endl;
1095 std::cerr << "aligned: " << PCK::aligned_3d(p1,p2,p3) << std::endl;
1096 return false;
1097 }
1098
1099
5/6
✓ Branch 0 taken 36556 times.
✓ Branch 1 taken 117376 times.
✓ Branch 2 taken 31129 times.
✓ Branch 3 taken 5427 times.
✓ Branch 4 taken 31129 times.
✗ Branch 5 not taken.
 221617 if(N2.x.sign() == ZERO && N2.y.sign() == ZERO && N2.z.sign() == ZERO) {
1100 std::cerr << std::endl;
1101 std::cerr << "degenerate triangle" << std::endl;
1102 std::cerr << "aligned: " << PCK::aligned_3d(q1,q2,q3) << std::endl;
1103 return false;
1104 }
1105
1106 // Tolerance for co-planarity test
1107
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 153932 times.
 153932 if(angle_tolerance_ != 0.0) {
1108 ✗ double threshold = cos(angle_tolerance_ * M_PI / 180.0);
1109 ✗ exact::scalar left = geo_sqr(dot(N1,N2));
1110 exact::scalar right =
1111 ✗ exact::scalar(threshold*threshold)*length2(N1)*length2(N2);
1112 return left > right;
1113 }
1114
1115 // Exact version
1116
1/2
✓ Branch 1 taken 153932 times.
✗ Branch 2 not taken.
 153932 exact::vec3 N12 = cross(N1,N2);
1117
6/6
✓ Branch 0 taken 38880 times.
✓ Branch 1 taken 115052 times.
✓ Branch 2 taken 32965 times.
✓ Branch 3 taken 5915 times.
✓ Branch 4 taken 2492 times.
✓ Branch 5 taken 30473 times.
 225777 if((N12.x.sign()!=ZERO) || (N12.y.sign()!=ZERO) ||(N12.z.sign()!=ZERO)) {
1118 123459 return false;
1119 }
1120
1121 return true;
1122 153932 }
1123
1124 /**************************************************************************/
1125
1126 246918 bool CoplanarFacets::edges_are_colinear(
1127 const ExactPoint& P1, const ExactPoint& P2, const ExactPoint& P3
1128 ) const {
1129
1130
1/2
✓ Branch 0 taken 246918 times.
✗ Branch 1 not taken.
 246918 if(angle_tolerance_ == 0.0) {
1131 246918 return PCK::on_segment_3d(P2,P1,P3);
1132 }
1133
1134 ✗ ExactPoint UU = P1-P2;
1135 ✗ exact::vec3 U(UU.x, UU.y, UU.z);
1136 ✗ if(UU.w.sign() == NEGATIVE) {
1137 U.x.negate(); U.y.negate(); U.z.negate();
1138 }
1139 ✗ ExactPoint VV = P3-P2;
1140 ✗ exact::vec3 V(VV.x, VV.y, VV.z);
1141 ✗ if(VV.w.sign() == NEGATIVE) {
1142 V.x.negate(); V.y.negate(); V.z.negate();
1143 }
1144
1145 ✗ double threshold = cos(angle_tolerance_ * M_PI / 180.0);
1146
1147 ✗ exact::scalar left = dot(U,V);
1148
1149 ✗ if(left.sign() == POSITIVE) {
1150 return false;
1151 }
1152
1153 ✗ left = geo_sqr(left);
1154 exact::scalar right =
1155 ✗ exact::scalar(threshold*threshold)*length2(U)*length2(V);
1156
1157 return left > right;
1158 ✗ }
1159
1160 /**************************************************************************/
1161
1162 }
1163