GCC Code Coverage Report

Directory: ./
File: lib/geogram/mesh/mesh_surface_intersection.cpp
Date: 2026-09-07 02:37:58
Exec Total Coverage
Lines: 847 1102 76.9%
Functions: 43 53 81.1%
Branches: 924 2208 41.8%
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.h>
41 #include <geogram/mesh/mesh_surface_intersection_internal.h>
42 #include <geogram/mesh/mesh_AABB.h>
43 #include <geogram/mesh/mesh_repair.h>
44 #include <geogram/mesh/mesh_fill_holes.h>
45 #include <geogram/mesh/mesh_geometry.h>
46 #include <geogram/mesh/mesh_topology.h>
47 #include <geogram/mesh/mesh_io.h>
48 #include <geogram/mesh/boxes_intersections.h>
49 #include <geogram/mesh/index.h>
50 #include <geogram/delaunay/CDT_2d.h>
51 #include <geogram/numerics/predicates.h>
52 #include <geogram/numerics/expansion_nt.h>
53 #include <geogram/basic/stopwatch.h>
54 #include <geogram/basic/permutation.h>
55 #include <geogram/basic/boolean_expression.h>
56 #include <geogram/basic/debug_stream.h>
57 #include <geogram/basic/algorithm.h>
58
59 #include <sstream>
60 #include <stack>
61
62 #ifdef GEO_COMPILER_CLANG
63 // I'm using long long
64 #pragma GCC diagnostic ignored "-Wc++98-compat-pedantic"
65 #endif
66
67 // If defined, displays status messages and saves files whenever some
68 // error conditions are met.
69 // #define MESH_SURFACE_INTERSECTION_DEBUG
70
71 namespace {
72 using namespace GEO;
73 /**
74 * \brief Removes all facets that have their three vertices aligned
75 */
76 61 void remove_degenerate_triangles(Mesh& M) {
77
1/2
✓ Branch 2 taken 61 times.
✗ Branch 3 not taken.
 61 vector<index_t> remove_f(M.facets.nb());
78
1/2
✓ Branch 3 taken 61 times.
✗ Branch 4 not taken.
 61 parallel_for(0, M.facets.nb(), [&M, &remove_f](index_t f) {
79 118536 const vec3& p1 = M.facets.point(f,0);
80 118536 const vec3& p2 = M.facets.point(f,1);
81 118536 const vec3& p3 = M.facets.point(f,2);
82 118536 remove_f[f] = PCK::aligned_3d(p1,p2,p3);
83 118536 });
84
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 M.facets.delete_elements(remove_f);
85 61 }
86
87 /**
88 * \brief Enumerates the connected components in a facet attribute
89 * \param[in] M a reference to the mesh
90 * \param[in] attribute the name of the facet attribute
91 * \return the number of found connected components
92 */
93 131 index_t get_surface_connected_components(
94 Mesh& M, const std::string& attribute = "chart"
95 ) {
96
1/2
✓ Branch 2 taken 131 times.
✗ Branch 3 not taken.
 131 Attribute<index_t> chart(M.facets.attributes(), attribute);
97
1/2
✓ Branch 1 taken 131 times.
✗ Branch 2 not taken.
 262 return GEO::get_connected_components(M, chart);
98 131 }
99
100 /**
101 * \brief Computes the intersection between two mesh triangular facets
102 * \details This function is just a wrapper around triangles_intersections()
103 * for Mesh facets.
104 * \param[in] M the mesh
105 * \param[in] f1 , f2 the two facets
106 * \param[out] I a vector of triangle intersections
107 * \retval true if there was an intersection
108 * \retval false otherwise
109 */
110 1007411 bool mesh_facets_intersect(
111 Mesh& M, index_t f1, index_t f2, TriangleIsects& I
112 ) {
113
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1007411 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1007411 geo_debug_assert(M.facets.nb_vertices(f1) == 3);
114
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1007411 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1007411 geo_debug_assert(M.facets.nb_vertices(f2) == 3);
115 1007411 index_t v1 = M.facets.vertex(f1,0);
116 1007411 index_t v2 = M.facets.vertex(f1,1);
117 1007411 index_t v3 = M.facets.vertex(f1,2);
118 1007411 index_t w1 = M.facets.vertex(f2,0);
119 1007411 index_t w2 = M.facets.vertex(f2,1);
120 1007411 index_t w3 = M.facets.vertex(f2,2);
121 2014822 return triangles_intersections(
122 1007411 M.vertices.point(v1),
123 1007411 M.vertices.point(v2),
124 1007411 M.vertices.point(v3),
125 1007411 M.vertices.point(w1),
126 1007411 M.vertices.point(w2),
127 1007411 M.vertices.point(w3),
128 v1,v2,v3,
129 w1,w2,w3,
130 I
131 1007411 );
132 }
133
134 /**
135 * \brief Tests whether a segment intersects a triangle
136 * \param[in] q1 , q2 the two extremities of the segment
137 * \param[in] p1 , p2 , p3 the three verties of the triangle
138 * \retval true if the segment has an intersection with the
139 * interior of the triangle
140 * \retval false otherwise
141 * \details Degenerate configurations (segment passing through vertex,
142 * edge, or co-planar with triangle) are symbolically perturbed.
143 */
144 2104 template <class POINT> bool segment_triangle_intersection_SOS(
145 const POINT& q1, const POINT& q2,
146 const POINT& p1, const POINT& p2, const POINT& p3
147 ) {
148 2104 Sign o1 = PCK::orient_3d_SOS(q1,p1,p2,p3);
149 2104 Sign o2 = PCK::orient_3d_SOS(q2,p1,p2,p3);
150
151 // There is no intersection if q1 and q2 are on the
152 // same side of the supporting plane of (p1,p2,p3)
153
2/2
✓ Branch 0 taken 713 times.
✓ Branch 1 taken 339 times.
 2104 if(o1 == o2) {
154 1426 return false;
155 }
156
157 // There is an intersection if the three tetrahedra
158 // formed by [q1,q2] and the three edges of the triangle
159 // have the same orientation
160 678 Sign s1 = PCK::orient_3d_SOS(q1,q2,p1,p2);
161 678 Sign s2 = PCK::orient_3d_SOS(q1,q2,p2,p3);
162
2/2
✓ Branch 0 taken 229 times.
✓ Branch 1 taken 110 times.
 678 if(s1*s2 < 0) {
163 458 return false;
164 }
165 220 Sign s3 = PCK::orient_3d_SOS(q1,q2,p3,p1);
166
3/4
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 107 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 220 return(s2*s3 > 0 && s3*s1 > 0);
167 }
168 }
169
170
171 namespace GEO {
172
173 61 MeshSurfaceIntersection::MeshSurfaceIntersection(Mesh& M) :
174 61 lock_(GEOGRAM_SPINLOCK_INIT),
175 61 mesh_(M),
176
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 61 times.
✗ Branch 6 not taken.
 122 vertex_to_exact_point_(M.vertices.attributes(), "exact_point"),
177 61 dry_run_(false),
178
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 halfedges_(*this),
179
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 radial_bundles_(*this),
180
3/6
✓ Branch 2 taken 61 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 61 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 61 times.
✗ Branch 10 not taken.
 183 radial_polylines_(*this)
181 {
182
2/2
✓ Branch 4 taken 59900 times.
✓ Branch 5 taken 61 times.
 59961 for(index_t v: mesh_.vertices) {
183
1/2
✓ Branch 1 taken 59900 times.
✗ Branch 2 not taken.
 59900 vertex_to_exact_point_[v] = nullptr;
184 }
185 61 verbose_ = false;
186 61 fine_verbose_ = false;
187 61 delaunay_ = true;
188 61 detect_intersecting_neighbors_ = true;
189 61 use_radial_sort_ = true;
190 61 monster_threshold_ = NO_INDEX;
191 61 has_operand_bits_ = false;
192 61 skeleton_ = nullptr;
193 61 skeleton_trim_fins_ = false;
194 61 interpolate_attributes_ = false;
195 61 }
196
197 61 MeshSurfaceIntersection::~MeshSurfaceIntersection() {
198 61 vertex_to_exact_point_.destroy();
199 61 }
200
201 ✗ void MeshSurfaceIntersection::remove_external_shell() {
202 ✗ vector<index_t> remove_f;
203 ✗ mark_external_shell(remove_f);
204 ✗ mesh_.facets.delete_elements(remove_f);
205 ✗ mesh_.facets.connect();
206 ✗ }
207
208 28 void MeshSurfaceIntersection::remove_internal_shells() {
209
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 28 times.
 28 if(mesh_.facets.nb() == 0) {
210 ✗ return;
211 }
212 28 vector<index_t> remove_f;
213
1/2
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
 28 mark_external_shell(remove_f);
214
2/2
✓ Branch 2 taken 225066 times.
✓ Branch 3 taken 28 times.
 225122 for(index_t& i: remove_f) {
215 225066 i=1-i;
216 }
217
1/2
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
 28 mesh_.facets.delete_elements(remove_f);
218 28 }
219
220 ✗ void MeshSurfaceIntersection::remove_fins() {
221 // TODO
222 ✗ }
223
224 61 void MeshSurfaceIntersection::intersect_prologue() {
225
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
 122 Stopwatch W("Prologue", verbose_);
226
2/2
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 56 times.
 61 if(!mesh_.facets.are_simplices()) {
227
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 tessellate_facets(mesh_,3);
228 }
229
230
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 Attribute<index_t> operand_bit;
231
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 operand_bit.bind_if_is_defined(
232
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 122 mesh_.facets.attributes(), "operand_bit"
233 );
234 61 has_operand_bits_ = operand_bit.is_bound();
235
2/2
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 46 times.
 61 if(!operand_bit.is_bound()) {
236 // TODO: not good, there might be more than 32 components
237 // (and why do we need to do that BTW?)
238
2/4
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
 30 get_surface_connected_components(mesh_,"operand_bit");
239
2/4
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 15 times.
✗ Branch 6 not taken.
 15 operand_bit.bind(mesh_.facets.attributes(), "operand_bit");
240
2/2
✓ Branch 4 taken 54108 times.
✓ Branch 5 taken 15 times.
 54123 for(index_t f: mesh_.facets) {
241
2/4
✓ Branch 1 taken 54108 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 54108 times.
✗ Branch 5 not taken.
 54108 operand_bit[f] = (index_t(1) << operand_bit[f]) ;
242 }
243 }
244
245
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 remove_degenerate_triangles(mesh_);
246
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 mesh_colocate_vertices_no_check(mesh_);
247
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 mesh_remove_bad_facets_no_check(mesh_);
248
249 // Set symbolic perturbation mode to lexicographic order
250 // on point coordinates instead of point indices only,
251 // Needed to get compatible triangulations on coplanar faces
252 // (example, cubes that touch on a facet).
253
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 SOS_bkp_ = PCK::get_SOS_mode();
254
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 PCK::set_SOS_mode(PCK::SOS_LEXICO);
255 61 }
256
257 61 void MeshSurfaceIntersection::intersect_get_intersections(
258 vector<IsectInfo>& intersections
259 ) {
260
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
 122 Stopwatch Wtot("Find isects", verbose_);
261 {
262 61 vector<std::pair<index_t, index_t> > FF;
263 1007411 auto report_BB = [this,&FF](index_t f1, index_t f2) {
264 // Optionally skip facet pairs that share a vertex or an edge
265 1007411 if(
266
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1007411 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1007411 times.
 1007411 !detect_intersecting_neighbors_ && (
267 ✗ (mesh_.facets.find_adjacent(f1,f2)!=NO_INDEX) ||
268 ✗ (mesh_.facets.find_common_vertex(f1,f2)!=NO_INDEX)
269 )
270 ) {
271 ✗ return;
272 }
273
2/4
✓ Branch 1 taken 1007411 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1007411 times.
✗ Branch 5 not taken.
 1007411 FF.push_back(std::make_pair(f1,f2));
274 61 };
275
276 // For now, keeping parallel AABB construction and evaluation,
277 // ZE seems to be slower, see stats in:
278 // https://github.com/BrunoLevy/geogram/issues/362
279 static constexpr bool use_ZE = false; // Zoromodian-Edelsbrunner
280 static constexpr bool ZE_groups = false;
281 static constexpr bool ZE_groups_auto_isect = true;
282
283 // Implementation using AABB
284 if(!use_ZE) {
285
4/10
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 61 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 61 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 183 Stopwatch* W = new Stopwatch("AABB build", verbose_);
286
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 MeshFacetsAABB AABB(mesh_, AABB_INDIRECT);
287
1/2
✓ Branch 0 taken 61 times.
✗ Branch 1 not taken.
 61 delete W;
288 // Get candidate pairs of intersecting facets
289
4/10
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 61 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 61 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 183 W = new Stopwatch("AABB box-box", verbose_);
290
1/2
✓ Branch 2 taken 61 times.
✗ Branch 3 not taken.
 61 AABB.compute_facet_bbox_intersections(report_BB);
291
1/2
✓ Branch 0 taken 61 times.
✗ Branch 1 not taken.
 61 delete W;
292 61 }
293
294 // Alternative implementation using Zoromodian-Edelsbrunner
295 if(use_ZE) {
296 Stopwatch W("AABB box-box", verbose_);
297 vector<Box3d> boxes(mesh_.facets.nb());
298 parallel_for(
299 0, mesh_.facets.nb(), [&](index_t f) {
300 vec3 p0 = mesh_.facets.point(f,0);
301 vec3 p1 = mesh_.facets.point(f,1);
302 vec3 p2 = mesh_.facets.point(f,2);
303 boxes[f].xyz_min[0] = std::min(std::min(p0.x,p1.x),p2.x);
304 boxes[f].xyz_min[1] = std::min(std::min(p0.y,p1.y),p2.y);
305 boxes[f].xyz_min[2] = std::min(std::min(p0.z,p1.z),p2.z);
306 boxes[f].xyz_max[0] = std::max(std::max(p0.x,p1.x),p2.x);
307 boxes[f].xyz_max[1] = std::max(std::max(p0.y,p1.y),p2.y);
308 boxes[f].xyz_max[2] = std::max(std::max(p0.z,p1.z),p2.z);
309 }
310 );
311
312
313 index_t nb_groups = 0;
314 vector<index_t> indices;
315 vector<index_t> group_ptr;
316 bool has_facets_in_several_groups = false;
317
318 auto has_one_bit = [](index_t x)->bool {
319 return ((x != 0) && ((x & (x - 1))) == 0);
320 };
321
322 // Compute groups
323 if(ZE_groups && has_operand_bits_) {
324 Attribute<index_t> operand_bit;
325 operand_bit.bind_if_is_defined(
326 mesh_.facets.attributes(), "operand_bit"
327 );
328 if(operand_bit.is_bound()) {
329 vector<index_t> group_size(32,0);
330 for(index_t f: mesh_.facets) {
331 for(index_t g=0; g<32; ++g) {
332 if(operand_bit[f] == (1u << g)) {
333 ++group_size[g];
334 } else if (!has_one_bit(operand_bit[f])) {
335 has_facets_in_several_groups = true;
336 break;
337 }
338 }
339 }
340 // If some facets belong to several groups (this
341 // can happen when some input facets are exactly
342 // co-planar), then we do not use group-based
343 // computation (it would be possible, but it
344 // would require to split the set of facets in
345 // several groups according to the groups they
346 // belong to). It is also possible to consider
347 // all the facets in multiple groups as a single
348 // group, but doing so is slow for high-arity
349 // operators (example: dragon_bas.scad).
350 if(!has_facets_in_several_groups) {
351 for(index_t g=0; g<32; ++g) {
352 if(group_size[g] != 0) {
353 nb_groups = std::max(nb_groups,g);
354 }
355 }
356 ++nb_groups;
357 group_ptr.resize(nb_groups+1);
358 group_ptr[0] = 0;
359 for(index_t g=0; g<nb_groups; ++g) {
360 group_ptr[g+1] = group_ptr[g] + group_size[g];
361 }
362 indices.resize(group_ptr[nb_groups]);
363 vector<index_t> group_insert = group_ptr;
364 for(index_t f: mesh_.facets) {
365 for(index_t g=0; g<32; ++g) {
366 if(operand_bit[f] == (1u << g)) {
367 indices[group_insert[g]] = f;
368 ++group_insert[g];
369 }
370 }
371 }
372 }
373 }
374 }
375
376 if(nb_groups == 0) {
377 boxes_intersections(boxes, report_BB);
378 } else {
379 for(index_t g1 = 0; g1 < nb_groups; ++g1) {
380 for(index_t g2 = 0; g2 < nb_groups; ++g2) {
381 if(ZE_groups_auto_isect || g1 != g2) {
382 index_t b1 = group_ptr[g1];
383 index_t e1 = group_ptr[g1+1];
384 index_t b2 = group_ptr[g2];
385 index_t e2 = group_ptr[g2+1];
386 boxes_intersections_hybrid_impl(
387 boxes.data(),
388 indices.data()+b1, indices.data()+e1,
389 boxes.data(),
390 indices.data()+b2, indices.data()+e2,
391 report_BB
392 );
393 }
394 }
395 }
396 }
397 }
398
399 // Compute facet-facet intersections in parallel
400
4/10
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 61 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 61 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 183 Stopwatch* W = new Stopwatch("AABB tri-tri", verbose_);
401 61 Process::spinlock lock = GEOGRAM_SPINLOCK_INIT;
402
1/2
✓ Branch 2 taken 61 times.
✗ Branch 3 not taken.
 61 parallel_for(
403
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 122 0, FF.size(), [&](index_t i) {
404
1/2
✓ Branch 1 taken 1007411 times.
✗ Branch 2 not taken.
 1007411 TriangleIsects I;
405
1/2
✓ Branch 1 taken 1007411 times.
✗ Branch 2 not taken.
 1007411 index_t f1 = FF[i].first;
406
1/2
✓ Branch 1 taken 1007411 times.
✗ Branch 2 not taken.
 1007411 index_t f2 = FF[i].second;
407
3/4
✓ Branch 1 taken 1007411 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 970610 times.
✓ Branch 4 taken 36801 times.
 1007411 if(!mesh_facets_intersect(mesh_,f1, f2, I)) {
408 970610 return;
409 }
410
411 36801 Process::acquire_spinlock(lock);
412
413
2/2
✓ Branch 1 taken 10907 times.
✓ Branch 2 taken 25894 times.
 36801 if(I.size() > 2) {
414 // Coplanar intersection: to generate the edges,
415 // test validity of all possible pairs of vertices.
416
2/2
✓ Branch 1 taken 41988 times.
✓ Branch 2 taken 10907 times.
 52895 for(index_t i1=0; i1< I.size(); ++i1) {
417
2/2
✓ Branch 0 taken 61612 times.
✓ Branch 1 taken 41988 times.
 103600 for(index_t i2=0; i2<i1; ++i2) {
418 IsectInfo II = {
419 f1, f2,
420 184836 I[i1].first, I[i1].second,
421 184836 I[i2].first, I[i2].second
422
4/8
✓ Branch 1 taken 61612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 61612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 61612 times.
✗ Branch 11 not taken.
 61612 };
423
424 // Valid edges are the ones where both
425 // extremities are on the same edge
426 // of f1 or on the same edge of f2
427 // (note: it is a *combinatorial*
428 // convex hull).
429
1/2
✓ Branch 1 taken 61612 times.
✗ Branch 2 not taken.
 61612 TriangleRegion AB1=regions_convex_hull(
430 II.A_rgn_f1,II.B_rgn_f1
431 );
432
433
1/2
✓ Branch 1 taken 61612 times.
✗ Branch 2 not taken.
 61612 TriangleRegion AB2=regions_convex_hull(
434 II.A_rgn_f2, II.B_rgn_f2
435 );
436
437
8/10
✓ Branch 1 taken 61612 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 40292 times.
✓ Branch 4 taken 21320 times.
✓ Branch 6 taken 40292 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 20668 times.
✓ Branch 9 taken 19624 times.
✓ Branch 10 taken 41988 times.
✓ Branch 11 taken 19624 times.
 61612 if(region_dim(AB1)==1 || region_dim(AB2)==1) {
438
1/2
✓ Branch 1 taken 41988 times.
✗ Branch 2 not taken.
 41988 intersections.push_back(II);
439
1/2
✓ Branch 1 taken 41988 times.
✗ Branch 2 not taken.
 41988 II.flip();
440
1/2
✓ Branch 1 taken 41988 times.
✗ Branch 2 not taken.
 41988 intersections.push_back(II);
441 }
442 }
443 }
444 } else {
445 // Intersection is either a segment
446 // or a vertex of f2.
447
1/2
✓ Branch 1 taken 25894 times.
✗ Branch 2 not taken.
 25894 TriangleRegion A_rgn_f1 = I[0].first;
448
1/2
✓ Branch 1 taken 25894 times.
✗ Branch 2 not taken.
 25894 TriangleRegion A_rgn_f2 = I[0].second;
449
450 25894 TriangleRegion B_rgn_f1 = A_rgn_f1;
451 25894 TriangleRegion B_rgn_f2 = A_rgn_f2;
452
453
2/2
✓ Branch 1 taken 21591 times.
✓ Branch 2 taken 4303 times.
 25894 if(I.size() == 2) {
454
1/2
✓ Branch 1 taken 21591 times.
✗ Branch 2 not taken.
 21591 B_rgn_f1 = I[1].first;
455
1/2
✓ Branch 1 taken 21591 times.
✗ Branch 2 not taken.
 21591 B_rgn_f2 = I[1].second;
456 }
457
458 IsectInfo II = {
459 f1, f2,
460 A_rgn_f1, A_rgn_f2,
461 B_rgn_f1, B_rgn_f2
462 25894 };
463
1/2
✓ Branch 1 taken 25894 times.
✗ Branch 2 not taken.
 25894 intersections.push_back(II);
464
1/2
✓ Branch 1 taken 25894 times.
✗ Branch 2 not taken.
 25894 II.flip();
465
1/2
✓ Branch 1 taken 25894 times.
✗ Branch 2 not taken.
 25894 intersections.push_back(II);
466 }
467 36801 Process::release_spinlock(lock);
468 }
469 );
470
1/2
✓ Branch 0 taken 61 times.
✗ Branch 1 not taken.
 61 delete W;
471 61 }
472 61 }
473
474 61 void MeshSurfaceIntersection::intersect_remesh_intersections(
475 vector<IsectInfo>& intersections
476 ) {
477
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
 122 Stopwatch W("CDT",verbose_);
478
479 // Keep track of original facet ids: they will be copied to the subfacets
480 // whenever a facet is split.
481 // NOTE: of course, facet ids are no longer valids once coplanar facets
482 // are merged.
483
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 61 times.
✗ Branch 6 not taken.
 122 original_facet_id_.bind(mesh_.facets.attributes(), "original_facet_id");
484
2/2
✓ Branch 4 taken 118500 times.
✓ Branch 5 taken 61 times.
 118561 for(index_t f: mesh_.facets) {
485
1/2
✓ Branch 1 taken 118500 times.
✗ Branch 2 not taken.
 118500 original_facet_id_[f] = f;
486 }
487
488 // We need to copy the initial mesh, because MeshInTriangle needs
489 // to access it in parallel threads, and without a copy, the internal
490 // arrays of the mesh can be modified whenever there is a
491 // reallocation. Without copying, we would need to insert many
492 // locks (each time the mesh is accessed).
493
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 mesh_copy_.copy(mesh_);
494
495 {
496 // Sort intersections by f1, so that all intersections between f1
497 // and another facet appear as a contiguous sequence.
498
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 GEO::sort(
499 61 intersections.begin(), intersections.end(),
500 ✗ [](const IsectInfo& a, const IsectInfo& b) -> bool {
501 1831588 return (a.f1 < b.f1);
502 }
503 );
504
505 // Now iterate on all intersections, and identify
506 // the [b,e[ intervals that correspond to the same f1 facet.
507 // Get starting indices of intersections in same facet.
508 61 vector<index_t> start;
509 {
510 61 index_t b=0;
511
2/2
✓ Branch 1 taken 15505 times.
✓ Branch 2 taken 61 times.
 15566 while(b < intersections.size()) {
512
1/2
✓ Branch 1 taken 15505 times.
✗ Branch 2 not taken.
 15505 start.push_back(b);
513 15505 index_t e = b;
514 15505 while(
515
4/4
✓ Branch 1 taken 151213 times.
✓ Branch 2 taken 56 times.
✓ Branch 3 taken 135764 times.
✓ Branch 4 taken 15505 times.
 302482 e < intersections.size() &&
516
4/6
✓ Branch 1 taken 151213 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 151213 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 135764 times.
✓ Branch 7 taken 15449 times.
 151213 intersections[e].f1 == intersections[b].f1
517 ) {
518 135764 ++e;
519 }
520 15505 b = e;
521 }
522
1/2
✓ Branch 2 taken 61 times.
✗ Branch 3 not taken.
 61 start.push_back(intersections.size());
523 }
524
525 // Display intersection stats
526
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 47 times.
 61 if(verbose_) {
527 14 index_t nb_intersections = intersections.size()/2;
528 14 index_t nb_intersected_triangles = (start.size()-1)/2;
529 14 index_t max_intersections_in_triangle = 0;
530
2/2
✓ Branch 1 taken 4527 times.
✓ Branch 2 taken 14 times.
 4541 for(index_t i=0; i+1<start.size(); ++i) {
531 9054 max_intersections_in_triangle = std::max(
532 max_intersections_in_triangle,
533
2/4
✓ Branch 1 taken 4527 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4527 times.
✗ Branch 5 not taken.
 4527 index_t(start[i+1]-start[i])
534 );
535 }
536
3/6
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14 times.
✗ Branch 8 not taken.
 42 Logger::out("CDT") << "Intersections: "
537
2/4
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
 14 << nb_intersections << std::endl;
538
3/6
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14 times.
✗ Branch 8 not taken.
 42 Logger::out("CDT") << "Intersected triangles: "
539
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 << nb_intersected_triangles
540
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 << std::endl;
541
3/6
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14 times.
✗ Branch 8 not taken.
 42 Logger::out("CDT") << "Max intersections in triangle: "
542
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 << max_intersections_in_triangle
543
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 << std::endl;
544 }
545
546
547 #define TRIANGULATE_IN_PARALLEL
548
549 61 index_t f_done = 0;
550 61 index_t f_tot = (start.size()-1);
551
552 #ifdef TRIANGULATE_IN_PARALLEL
553
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 parallel_for_slice(
554
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 122 0,start.size()-1, [&](index_t k1, index_t k2) {
555 #else
556 index_t k1 = 0;
557 index_t k2 = start.size()-1;
558 #endif
559
560
1/2
✓ Branch 1 taken 229 times.
✗ Branch 2 not taken.
 229 MeshInTriangle MIT(*this);
561 229 MIT.set_delaunay(delaunay_);
562 229 MIT.set_dry_run(dry_run_);
563
564
1/2
✓ Branch 1 taken 229 times.
✗ Branch 2 not taken.
 229 index_t tid = Thread::current_id();
565
566
2/2
✓ Branch 0 taken 15505 times.
✓ Branch 1 taken 229 times.
 15734 for(index_t k=k1; k<k2; ++k) {
567
1/2
✓ Branch 1 taken 15505 times.
✗ Branch 2 not taken.
 15505 index_t b = start[k];
568
1/2
✓ Branch 1 taken 15505 times.
✗ Branch 2 not taken.
 15505 index_t e = start[k+1];
569
570
2/2
✓ Branch 0 taken 4527 times.
✓ Branch 1 taken 10978 times.
 15505 if(fine_verbose_) {
571 4527 ++f_done;
572
2/4
✓ Branch 1 taken 4527 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4527 times.
✗ Branch 5 not taken.
 9054 Logger::out("Isect")
573 4527 << String::format(
574 "[%2d] %5d/%5d %6d:%3d",
575 4527 int(tid), int(f_done), int(f_tot),
576
2/4
✓ Branch 1 taken 4527 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4527 times.
✗ Branch 5 not taken.
 4527 int(intersections[b].f1), int(e-b)
577
1/2
✓ Branch 1 taken 4527 times.
✗ Branch 2 not taken.
 4527 )
578
1/2
✓ Branch 1 taken 4527 times.
✗ Branch 2 not taken.
 4527 << std::endl;
579 }
580
581
2/4
✓ Branch 1 taken 15505 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15505 times.
✗ Branch 5 not taken.
 15505 MIT.begin_facet(intersections[b].f1);
582
2/2
✓ Branch 0 taken 135764 times.
✓ Branch 1 taken 15505 times.
 151269 for(index_t i=b; i<e; ++i) {
583
1/2
✓ Branch 1 taken 135764 times.
✗ Branch 2 not taken.
 135764 const IsectInfo& II = intersections[i];
584
585 // Each IsectInfo is either an individual vertex
586 // or a segment with two vertices.
587 // Each vertex is represented combinatorially.
588 // The MeshInTriangle knows how to compute the
589 // geometry from the combinatorial information.
590
591
2/2
✓ Branch 1 taken 8606 times.
✓ Branch 2 taken 127158 times.
 135764 if(II.is_point()) {
592 8606 MIT.add_vertex(
593 8606 II.f2,
594
1/2
✓ Branch 1 taken 8606 times.
✗ Branch 2 not taken.
 8606 II.A_rgn_f1, II.A_rgn_f2
595 );
596 } else {
597 127158 MIT.add_edge(
598 127158 II.f2,
599 127158 II.A_rgn_f1, II.A_rgn_f2,
600
1/2
✓ Branch 1 taken 127158 times.
✗ Branch 2 not taken.
 127158 II.B_rgn_f1, II.B_rgn_f2
601 );
602 }
603 }
604
605
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15505 times.
 15505 if(e-b >= monster_threshold_) {
606 ✗ index_t f = intersections[b].f1;
607 ✗ MIT.save_constraints(
608 ✗ "constraints_"+String::to_string(f)+".geogram"
609 );
610 }
611
612 // Inserts constraints
613 // and creates new vertices in shared mesh
614
1/2
✓ Branch 1 taken 15505 times.
✗ Branch 2 not taken.
 15505 MIT.commit();
615
616 // For debugging, optionally save "monsters"
617 // (that is, triangles that have a huge number
618 // of intersections).
619
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15505 times.
 15505 if(e-b >= monster_threshold_) {
620 ✗ index_t f = intersections[b].f1;
621 ✗ MIT.save(
622 ✗ "triangulation_"+String::to_string(f)+".geogram"
623 );
624 //MIT.save_constraints(
625 // "constraints_"+String::to_string(f)+".geogram"
626 //);
627 std::ofstream out(
628 ✗ "facet_"+String::to_string(f)+".obj"
629 ✗ );
630 ✗ for(index_t v: mesh_copy_.facets.vertices(f)) {
631 ✗ vec3 p = mesh_copy_.vertices.point(v);
632 ✗ out << "v " << p << std::endl;
633 }
634 ✗ out << "f ";
635 ✗ for(
636 ✗ index_t lv=0;
637 ✗ lv<mesh_copy_.facets.nb_vertices(f); ++lv
638 ) {
639 ✗ out << lv+1 << " ";
640 }
641 ✗ out << std::endl;
642 ✗ }
643
644 // Clear it so that it is clean for next triangle.
645
1/2
✓ Branch 1 taken 15505 times.
✗ Branch 2 not taken.
 15505 MIT.clear();
646 }
647
2/2
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 177 times.
 229 if(fine_verbose_) {
648
2/4
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 52 times.
✗ Branch 5 not taken.
 104 Logger::out("Isect")
649
2/4
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 52 times.
✗ Branch 5 not taken.
 104 << String::format("[%2d] done",int(tid))
650
1/2
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
 52 << std::endl;
651 }
652 #ifdef TRIANGULATE_IN_PARALLEL
653 229 });
654 #endif
655 61 }
656 61 }
657
658 61 void MeshSurfaceIntersection::intersect_epilogue(
659 const vector<IsectInfo>& intersections
660 ) {
661
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
 122 Stopwatch W_epilogue("Epilogue", verbose_);
662 // Vertices coming from intersections may land exactly
663 // on an existing vertex (see #111)
664 {
665
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
 122 Stopwatch W("I on v", verbose_);
666
1/2
✓ Branch 2 taken 61 times.
✗ Branch 3 not taken.
 61 vector<index_t> v2v(mesh_.vertices.nb());
667
2/2
✓ Branch 4 taken 84490 times.
✓ Branch 5 taken 61 times.
 84551 for(index_t v : mesh_.vertices) {
668
1/2
✓ Branch 1 taken 84490 times.
✗ Branch 2 not taken.
 84490 v2v[v] = v;
669 }
670
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 parallel_for(
671 61 0, mesh_.vertices.nb(),
672 122 [&](index_t v) {
673 // If the vertex is an original vertex, not
674 // coming from an intersection, check whether
675 // it also exists as an intersection
676
2/2
✓ Branch 1 taken 59595 times.
✓ Branch 2 taken 24895 times.
 84490 if(vertex_to_exact_point_[v] == nullptr) {
677
1/2
✓ Branch 1 taken 59595 times.
✗ Branch 2 not taken.
 59595 const vec3& xyz = mesh_.vertices.point(v);
678
4/8
✓ Branch 1 taken 59595 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 59595 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 59595 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 59595 times.
✗ Branch 11 not taken.
 59595 ExactPoint K(xyz[0], xyz[1], xyz[2], 1.0);
679
1/2
✓ Branch 1 taken 59595 times.
✗ Branch 2 not taken.
 59595 auto it = exact_point_to_vertex_.find(K);
680
2/2
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 59591 times.
 59595 if(it != exact_point_to_vertex_.end()) {
681
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 v2v[v] = it->second;
682 }
683 59595 }
684 84490 }
685 );
686
2/2
✓ Branch 5 taken 749661 times.
✓ Branch 6 taken 61 times.
 749722 for(index_t c : mesh_.facet_corners) {
687
2/4
✓ Branch 1 taken 749661 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 749661 times.
✗ Branch 5 not taken.
 749661 index_t v = v2v[mesh_.facet_corners.vertex(c)];
688
1/2
✓ Branch 1 taken 749661 times.
✗ Branch 2 not taken.
 749661 mesh_.facet_corners.set_vertex(c, v);
689 }
690 61 }
691
692
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 61 times.
 61 if(interpolate_attributes_) {
693 ✗ Stopwatch W("Attrib", verbose_);
694 Attribute<index_t> original_facet_id(
695 ✗ mesh_.facets.attributes(), "original_facet_id"
696 ✗ );
697
698 ✗ for(index_t c: mesh_.facet_corners) {
699 ✗ index_t f = c/3;
700 ✗ index_t f0 = original_facet_id[f];
701 // If this is an original facet, there is nothing to do !
702 ✗ if(f == f0) {
703 ✗ continue;
704 }
705
706 ✗ index_t v = mesh_.facet_corners.vertex(c);
707 ✗ index_t v1 = mesh_.facets.vertex(f0,0);
708 ✗ index_t v2 = mesh_.facets.vertex(f0,1);
709 ✗ index_t v3 = mesh_.facets.vertex(f0,2);
710
711 ✗ index_t c1 = mesh_.facets.corner(f0,0);
712 ✗ index_t c2 = mesh_.facets.corner(f0,1);
713 ✗ index_t c3 = mesh_.facets.corner(f0,2);
714
715 ✗ if(v == v1) {
716 ✗ mesh_.facet_corners.attributes().copy_item(c,c1);
717 ✗ continue;
718 }
719
720 ✗ if(v == v2) {
721 ✗ mesh_.facet_corners.attributes().copy_item(c,c2);
722 ✗ continue;
723 }
724
725 ✗ if(v == v3) {
726 ✗ mesh_.facet_corners.attributes().copy_item(c,c3);
727 ✗ continue;
728 }
729
730 ✗ vec3 p1=mesh_.vertices.point(v1);
731 ✗ vec3 p2=mesh_.vertices.point(v2);
732 ✗ vec3 p3=mesh_.vertices.point(v3);
733 ✗ vec3 p =mesh_.vertices.point(v);
734
735 ✗ double a = Geom::triangle_area(p1,p2,p3);
736 ✗ double a1 = Geom::triangle_area(p ,p2,p3);
737 ✗ double a2 = Geom::triangle_area(p1,p ,p3);
738 ✗ double a3 = Geom::triangle_area(p1,p2,p );
739
740 ✗ mesh_.facet_corners.attributes().zero_item(c);
741 ✗ mesh_.facet_corners.attributes().madd_item(c, a1/a, c1);
742 ✗ mesh_.facet_corners.attributes().madd_item(c, a2/a, c2);
743 ✗ mesh_.facet_corners.attributes().madd_item(c, a3/a, c3);
744 }
745 ✗ }
746
747
748 // Remove original facets that have intersections.
749 {
750
1/2
✓ Branch 2 taken 61 times.
✗ Branch 3 not taken.
 61 vector<index_t> has_intersections(mesh_.facets.nb(), 0);
751
2/2
✓ Branch 2 taken 135764 times.
✓ Branch 3 taken 61 times.
 135886 for(const IsectInfo& II: intersections) {
752
1/2
✓ Branch 1 taken 135764 times.
✗ Branch 2 not taken.
 135764 has_intersections[II.f1] = 1;
753
1/2
✓ Branch 1 taken 135764 times.
✗ Branch 2 not taken.
 135764 has_intersections[II.f2] = 1;
754 }
755
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 mesh_.facets.delete_elements(has_intersections);
756 61 }
757
758 // There can be duplicated facets coming from
759 // tesselated co-planar facets.
760 // Note: this updates operand_bit attribute
761
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 mesh_remove_bad_facets_no_check(mesh_);
762
763 #ifdef MESH_SURFACE_INTERSECTION_DEBUG
764 std::cerr << "Sanity check: verify that there is no degenerate triangle"
765 << std::endl;
766 // Sanity check: do we have facets with their three vertices
767 // aligned ? Normally cannot happen since we have eliminated
768 // them during intersection, but who knows ?
769 // Actually this happens sometimes...
770 {
771 Attribute<bool> selected(mesh_.facets.attributes(), "selection");
772 for(index_t t: mesh_.facets) {
773 if(PCK::aligned_3d(
774 exact_vertex(mesh_.facets.vertex(t,0)),
775 exact_vertex(mesh_.facets.vertex(t,1)),
776 exact_vertex(mesh_.facets.vertex(t,2))
777 )) {
778 selected[t] = true;
779 std::cerr << "FACET HAS 3 ALIGNED VERTICES" << std::endl;
780 } else {
781 selected[t] = false;
782 }
783
784 }
785 }
786 std::cerr << "There is no degenerate triangle" << std::endl;
787 #endif
788
789
2/2
✓ Branch 0 taken 58 times.
✓ Branch 1 taken 3 times.
 61 if(use_radial_sort_) {
790
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 build_Weiler_model();
791 }
792
793
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 PCK::set_SOS_mode(SOS_bkp_);
794 61 }
795
796 61 void MeshSurfaceIntersection::intersect() {
797
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 61 times.
 61 if(mesh_.facets.nb() == 0) {
798 ✗ return;
799 }
800
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
 122 Stopwatch W("Intersect", verbose_);
801
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 intersect_prologue();
802 61 vector<IsectInfo> intersections;
803
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 intersect_get_intersections(intersections);
804
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 intersect_remesh_intersections(intersections);
805
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 intersect_epilogue(intersections);
806 61 }
807
808 807842 MeshSurfaceIntersection::ExactPoint MeshSurfaceIntersection::exact_vertex(
809 index_t v
810 ) const {
811
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 807842 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 807842 geo_debug_assert(v < mesh_.vertices.nb());
812 807842 const ExactPoint* p = vertex_to_exact_point_[v];
813
2/2
✓ Branch 0 taken 150767 times.
✓ Branch 1 taken 657075 times.
 807842 if(p != nullptr) {
814 150767 return *p;
815 }
816 657075 const vec3& xyz = mesh_.vertices.point(v);
817 657075 return ExactPoint(xyz[0], xyz[1], xyz[2], 1.0);
818 }
819
820 91237 index_t MeshSurfaceIntersection::find_or_create_exact_vertex(
821 const ExactPoint& p
822 ) {
823 91237 std::map<ExactPoint,index_t,ExactPointCompare>::iterator it;
824 bool inserted;
825
1/2
✓ Branch 1 taken 91237 times.
✗ Branch 2 not taken.
 91237 std::tie(it, inserted) = exact_point_to_vertex_.insert(
826
1/2
✓ Branch 1 taken 91237 times.
✗ Branch 2 not taken.
 182474 std::make_pair(p,NO_INDEX)
827 91237 );
828
2/2
✓ Branch 0 taken 66342 times.
✓ Branch 1 taken 24895 times.
 91237 if(!inserted) {
829 66342 return it->second;
830 }
831
1/2
✓ Branch 1 taken 24895 times.
✗ Branch 2 not taken.
 24895 vec3 p_inexact = PCK::approximate(p);
832
1/2
✓ Branch 2 taken 24895 times.
✗ Branch 3 not taken.
 24895 index_t v = mesh_.vertices.create_vertex(p_inexact.data());
833 24895 it->second = v;
834
1/2
✓ Branch 2 taken 24895 times.
✗ Branch 3 not taken.
 24895 vertex_to_exact_point_[v] = &(it->first);
835 24895 return v;
836 }
837
838 /************************ Radial sort ***********************************/
839
840 3820 void MeshSurfaceIntersection::RadialSort::init(index_t h_ref) {
841 3820 degenerate_ = false;
842 3820 h_ref_ = NO_INDEX; // so that normal() computes N_ref_
843
2/4
✓ Branch 1 taken 3820 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3820 times.
✗ Branch 5 not taken.
 3820 N_ref_ = normal(h_ref);
844 3820 h_ref_ = h_ref;
845 3820 }
846
847 24479 bool MeshSurfaceIntersection::RadialSort::operator()(
848 index_t h1, index_t h2
849 ) const {
850
851 // To order facets around radial edges, we use two different things:
852 //
853 // 1) a 2D coordinate system. For each halfedge h, this gives
854 // two signs. The 2D coordinate system is defined by:
855 // - a reference halfedge h_ref passed to init(h_ref)
856 // - u_sign (su): location of opposite vertex relative to facet
857 // incident to h_ref: h_orient(h_ref_, h)
858 // - v_sign (sv): dot product between the normals to the facets
859 // incident to h and incident to h_ref: h_refNorient(h)
860 // u_sign and v_sign determine 4 quadrants (in fact 8 regions if one
861 // counts when they vanish). If two halfedges land in two different
862 // regions then their order is determined. The order of the regions
863 // is stored in su_sv_to_linear_index[3][3].
864 //
865 // 2) now if two halfedges h1,h2 land in the same region, their order is
866 // given by orient3d(p1,p2,q1,q2) = horient(h1,h2) where p1,p2 are
867 // the vertices of the radial edge and q1 (resp q2) the opposite
868 // vertices of h1 (resp h2).
869
870 24479 Sign su1 = h_orient(h_ref_, h1);
871 24479 Sign su2 = h_orient(h_ref_, h2);
872
873 // Optimization: quick exit, h1 appears first if it is the
874 // only one in positive u half.
875
2/2
✓ Branch 0 taken 2122 times.
✓ Branch 1 taken 22357 times.
 24479 if(su1 * su2 < 0) {
876 2122 return (su1 > 0);
877 }
878
879 22357 Sign sv1 = h_refNorient(h1);
880 22357 Sign sv2 = h_refNorient(h2);
881
882 // Map (su,sv) sign pair into a region linear index ("pseudo-angle")
883 //
884 // su
885 // ^
886 // | + 2
887 // | 3 | 1
888 // | 0 4---+---0
889 // | 5 | 7
890 // | - 6
891 // |
892 // | - 0 +
893 // *--------------> sv
894 static int su_sv_to_linear_index[3][3] = {
895 // -, 0, + <--------- su = h_orient(h_ref, h)
896 {5, 4, 3}, // - <--- sv = h_Norient(h)
897 {6,-1, 2}, // 0
898 {7, 0, 1} // +
899 };
900
901 22357 int theta1 = su_sv_to_linear_index[sv1+1][su1+1];
902 22357 int theta2 = su_sv_to_linear_index[sv2+1][su2+1];
903
904
2/4
✓ Branch 0 taken 22357 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 22357 times.
 22357 if(theta1 == -1 || theta2 == -1) {
905 ✗ report_problem("Triangle with both zero orient and zero Norient");
906 ✗ return false;
907 }
908
909
2/2
✓ Branch 0 taken 22355 times.
✓ Branch 1 taken 2 times.
 22357 if(theta1 != theta2) { // Different thetas: then we know the order
910 22355 return (theta2 > theta1);
911 }
912
913
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if((theta1 & 1) == 0) { // Same thetas, = to 0,2,4,6 (should not happen)
914 ✗ report_problem("Both triangles in same reference half-plane");
915 ✗ return false;
916 }
917
918 // Same thetas, both = to 1,3,5 or 7, measure relative orientation
919 2 Sign o_12 = h_orient(h1,h2);
920
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(o_12 == ZERO) {
921 ✗ report_problem("Both triangles in same half-plane");
922 ✗ return false;
923 }
924 2 return o_12 > 0;
925 }
926
927 48960 Sign MeshSurfaceIntersection::RadialSort::h_orient(
928 index_t h1, index_t h2
929 ) const {
930
931 // Normally we compute -orient3d(p1,p2,q1,q2) where
932 // p1,p2 denote the radial edge and q1 (resp q2) the opposite
933 // vertex to h1 (resp h2).
934 // Replacing a facet of a tetrahedron with a co-planar triangle
935 // of same orientation does not change signed volume, so instead
936 // of the facet adjacent to h1, we use the facet of the original
937 // mesh (that contains the facet adjacent to h1), that has simpler
938 // coordinates.
939
940
2/2
✓ Branch 0 taken 16550 times.
✓ Branch 1 taken 32410 times.
 48960 if(h1 == h2) {
941 16550 return ZERO;
942 }
943 32410 index_t f1 = I_.halfedges_.facet(h1);
944 32410 index_t f2 = I_.halfedges_.facet(h2);
945
1/2
✓ Branch 1 taken 32410 times.
✗ Branch 2 not taken.
 32410 index_t w1 = I_.halfedges_.vertex(h1,2);
946
1/2
✓ Branch 1 taken 32410 times.
✗ Branch 2 not taken.
 32410 index_t w2 = I_.halfedges_.vertex(h2,2);
947
948 // Optimization (using original points as often as possible)
949 // If w1 is an original vertex, use w1 and original facet of h2
950 // (there is no minus sign, because args are swapped, we have h2 first)
951
3/4
✓ Branch 1 taken 32410 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5173 times.
✓ Branch 4 taken 27237 times.
 32410 if(I_.is_original_vertex(w1)) {
952
1/2
✓ Branch 1 taken 5173 times.
✗ Branch 2 not taken.
 5173 vec3 p0 = mesh_.vertices.point(w1);
953
1/2
✓ Branch 1 taken 5173 times.
✗ Branch 2 not taken.
 5173 auto [q0, q1, q2] = I_.get_initial_facet_vertices(f2);
954
1/2
✓ Branch 1 taken 5173 times.
✗ Branch 2 not taken.
 5173 return Sign(PCK::orient_3d(q0,q1,q2,p0));
955 }
956
957 // Optimization (using original points as often as possible)
958 // If w2 is an original vertex, use w2 and original facet of h1
959
3/4
✓ Branch 1 taken 27237 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3604 times.
✓ Branch 4 taken 23633 times.
 27237 if(I_.is_original_vertex(w2)) {
960
1/2
✓ Branch 1 taken 3604 times.
✗ Branch 2 not taken.
 3604 vec3 q0 = mesh_.vertices.point(w2);
961
1/2
✓ Branch 1 taken 3604 times.
✗ Branch 2 not taken.
 3604 auto [p0, p1, p2] = I_.get_initial_facet_vertices(f1);
962
1/2
✓ Branch 1 taken 3604 times.
✗ Branch 2 not taken.
 3604 return Sign(-PCK::orient_3d(p0,p1,p2,q0));
963 }
964
965 // General case: use w2 (it's an intersection point) and original
966 // facet of h1 converted to exact points.
967
1/2
✓ Branch 1 taken 23633 times.
✗ Branch 2 not taken.
 23633 auto [p0, p1, p2] = I_.get_initial_facet_vertices(f1);
968
969
1/2
✓ Branch 1 taken 23633 times.
✗ Branch 2 not taken.
 23633 ExactPoint pp0(p0.x, p0.y, p0.z, 1.0);
970
1/2
✓ Branch 1 taken 23633 times.
✗ Branch 2 not taken.
 23633 ExactPoint pp1(p1.x, p1.y, p1.z, 1.0);
971
1/2
✓ Branch 1 taken 23633 times.
✗ Branch 2 not taken.
 23633 ExactPoint pp2(p2.x, p2.y, p2.z, 1.0);
972
1/2
✓ Branch 1 taken 23633 times.
✗ Branch 2 not taken.
 23633 const ExactPoint& q2 = I_.exact_vertex(w2);
973
974
1/2
✓ Branch 1 taken 23633 times.
✗ Branch 2 not taken.
 23633 return Sign(-PCK::orient_3d(pp0,pp1,pp2,q2));
975
976 23633 }
977
978 44714 Sign MeshSurfaceIntersection::RadialSort::h_refNorient(index_t h2) const {
979
2/2
✓ Branch 0 taken 16550 times.
✓ Branch 1 taken 28164 times.
 44714 if(h2 == h_ref_) {
980 16550 return POSITIVE;
981 }
982
5/8
✓ Branch 0 taken 44 times.
✓ Branch 1 taken 28120 times.
✓ Branch 3 taken 33 times.
✓ Branch 4 taken 11 times.
✓ Branch 6 taken 33 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 28164 static PCK::PredicateStats stats("h_refNorient");
983 28164 stats.log_invoke();
984
1/2
✓ Branch 1 taken 28164 times.
✗ Branch 2 not taken.
 28164 exact::vec3 N2 = normal(h2);
985
2/4
✓ Branch 1 taken 28164 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28164 times.
✗ Branch 5 not taken.
 28164 return dot(N_ref_,N2).sign();
986 28164 }
987
988 31984 exact::vec3 MeshSurfaceIntersection::RadialSort::normal(index_t h) const {
989
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 31984 times.
 31984 if(h == h_ref_) {
990 ✗ return N_ref_;
991 }
992 31984 auto [p1, p2, p3] = I_.get_initial_facet_vertices(
993 31984 I_.halfedges_.facet(h)
994
1/2
✓ Branch 2 taken 31984 times.
✗ Branch 3 not taken.
 31984 );
995 exact::vec3 N = cross(
996
1/2
✓ Branch 1 taken 31984 times.
✗ Branch 2 not taken.
 63968 make_vec3<exact::vec3>(p1,p2),
997
1/2
✓ Branch 1 taken 31984 times.
✗ Branch 2 not taken.
 63968 make_vec3<exact::vec3>(p1,p3)
998
1/2
✓ Branch 1 taken 31984 times.
✗ Branch 2 not taken.
 31984 );
999
1/2
✓ Branch 1 taken 31984 times.
✗ Branch 2 not taken.
 31984 Numeric::optimize_number_representation(N);
1000 31984 return N;
1001 31984 }
1002
1003 ✗ void MeshSurfaceIntersection::RadialSort::report_problem(
1004 const char* message
1005 ) const {
1006 ✗ degenerate_ = true;
1007 ✗ Logger::err("RadialSort") << message << std::endl;
1008 ✗ }
1009
1010 /*****************************************************************/
1011
1012 28 void MeshSurfaceIntersection::mark_external_shell(
1013 vector<index_t>& on_external_shell
1014 ) {
1015 // Chart attribute corresponds to volumetric regions
1016
2/4
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 28 times.
✗ Branch 6 not taken.
 28 Attribute<index_t> chart(mesh_.facets.attributes(), "chart");
1017
1018 // Get nb charts
1019 28 index_t nb_charts = 0;
1020
2/2
✓ Branch 4 taken 225066 times.
✓ Branch 5 taken 28 times.
 225094 for(index_t f: mesh_.facets) {
1021
1/2
✓ Branch 1 taken 225066 times.
✗ Branch 2 not taken.
 225066 nb_charts = std::max(nb_charts, chart[f]+1);
1022 }
1023
1024 // Get connected components by traversing both alpha2 and alpha3 links,
1025 // and orient facets coherently
1026 28 index_t nb_components = 0;
1027
1/2
✓ Branch 2 taken 28 times.
✗ Branch 3 not taken.
 28 vector<index_t> component(mesh_.facets.nb(), NO_INDEX);
1028 {
1029
2/2
✓ Branch 4 taken 225066 times.
✓ Branch 5 taken 28 times.
 225094 for(index_t f:mesh_.facets) {
1030
3/4
✓ Branch 1 taken 225066 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 31 times.
✓ Branch 4 taken 225035 times.
 225066 if(component[f] == NO_INDEX) {
1031
1/2
✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
 31 std::stack<index_t> S;
1032
1/2
✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
 31 component[f] = nb_components;
1033
1/2
✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
 31 S.push(f);
1034
2/2
✓ Branch 1 taken 225066 times.
✓ Branch 2 taken 31 times.
 225097 while(!S.empty()) {
1035 225066 index_t f1 = S.top();
1036 225066 S.pop();
1037
1038 {
1039
1/2
✓ Branch 1 taken 225066 times.
✗ Branch 2 not taken.
 225066 index_t f2 = halfedges_.facet_alpha3(f1);
1040
3/4
✓ Branch 1 taken 225066 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 72442 times.
✓ Branch 4 taken 152624 times.
 225066 if(component[f2] == NO_INDEX) {
1041
2/4
✓ Branch 1 taken 72442 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72442 times.
✗ Branch 5 not taken.
 72442 component[f2]=component[f1];
1042
1/2
✓ Branch 1 taken 72442 times.
✗ Branch 2 not taken.
 72442 S.push(f2);
1043 }
1044 }
1045
1046
2/2
✓ Branch 0 taken 675198 times.
✓ Branch 1 taken 225066 times.
 900264 for(index_t le1=0; le1<3; ++le1) {
1047
1/2
✓ Branch 1 taken 675198 times.
✗ Branch 2 not taken.
 675198 index_t f2 = mesh_.facets.adjacent(f1,le1);
1048
6/8
✓ Branch 0 taken 675198 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 675198 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 152593 times.
✓ Branch 6 taken 522605 times.
✓ Branch 7 taken 152593 times.
✓ Branch 8 taken 522605 times.
 675198 if(f2!=NO_INDEX && component[f2]==NO_INDEX) {
1049 #ifdef GEO_DEBUG
1050
1/2
✓ Branch 1 taken 152593 times.
✗ Branch 2 not taken.
 152593 index_t le2 = mesh_.facets.find_adjacent(f2,f1);
1051
3/10
✓ Branch 1 taken 152593 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 152593 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 152593 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 152593 geo_debug_assert(
1052 mesh_.facets.vertex(f1,le1) !=
1053 mesh_.facets.vertex(f2,le2)
1054 );
1055 #endif
1056
2/4
✓ Branch 1 taken 152593 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 152593 times.
✗ Branch 5 not taken.
 152593 component[f2] = component[f1];
1057
1/2
✓ Branch 1 taken 152593 times.
✗ Branch 2 not taken.
 152593 S.push(f2);
1058 }
1059 }
1060 }
1061 31 }
1062 225066 ++nb_components;
1063 }
1064 }
1065
1066 // Compute the volume enclosed by each chart
1067
1068
1/2
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
 28 vector<double> chart_volume(nb_charts,0.0);
1069
2/2
✓ Branch 4 taken 225066 times.
✓ Branch 5 taken 28 times.
 225094 for(index_t f: mesh_.facets) {
1070
1/2
✓ Branch 1 taken 225066 times.
✗ Branch 2 not taken.
 225066 vec3 p1 = mesh_.facets.point(f,0);
1071
1/2
✓ Branch 1 taken 225066 times.
✗ Branch 2 not taken.
 225066 vec3 p2 = mesh_.facets.point(f,1);
1072
1/2
✓ Branch 1 taken 225066 times.
✗ Branch 2 not taken.
 225066 vec3 p3 = mesh_.facets.point(f,2);
1073
2/4
✓ Branch 3 taken 225066 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 225066 times.
✗ Branch 7 not taken.
 225066 chart_volume[chart[f]] += dot(p1,cross(p2,p3)) / 6.0;
1074 }
1075
1076
2/2
✓ Branch 1 taken 1124 times.
✓ Branch 2 taken 28 times.
 1152 for(index_t c=0; c<chart_volume.size(); ++c) {
1077
2/4
✓ Branch 1 taken 1124 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1124 times.
✗ Branch 5 not taken.
 1124 chart_volume[c] = ::fabs(chart_volume[c]);
1078 }
1079
1080 // For each component, find the chart that encloses the largest
1081 // volume (it is the external boundary of the component)
1082
1083
1/2
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
 28 vector<double> max_chart_volume_in_component(nb_components, 0.0);
1084 vector<index_t> chart_with_max_volume_in_component(
1085 nb_components, NO_INDEX
1086
1/2
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
 28 );
1087
1088
2/2
✓ Branch 5 taken 225066 times.
✓ Branch 6 taken 28 times.
 225094 for(index_t f: mesh_.facets) {
1089
2/4
✓ Branch 1 taken 225066 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 225066 times.
✗ Branch 5 not taken.
 225066 double V = chart_volume[chart[f]];
1090
4/6
✓ Branch 1 taken 225066 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 225066 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 108866 times.
✓ Branch 7 taken 116200 times.
 225066 if( V >= max_chart_volume_in_component[component[f]]) {
1091
2/4
✓ Branch 1 taken 108866 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 108866 times.
✗ Branch 5 not taken.
 108866 max_chart_volume_in_component[component[f]] = V;
1092
3/6
✓ Branch 1 taken 108866 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 108866 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 108866 times.
✗ Branch 8 not taken.
 108866 chart_with_max_volume_in_component[component[f]] = chart[f];
1093 }
1094 }
1095
1096
1/2
✓ Branch 2 taken 28 times.
✗ Branch 3 not taken.
 28 on_external_shell.resize(mesh_.facets.nb());
1097
2/2
✓ Branch 4 taken 225066 times.
✓ Branch 5 taken 28 times.
 225094 for(index_t f: mesh_.facets) {
1098 225066 on_external_shell[f] = (
1099
4/8
✓ Branch 1 taken 225066 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 225066 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 225066 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 225066 times.
✗ Branch 11 not taken.
 225066 chart[f] == chart_with_max_volume_in_component[component[f]]
1100 );
1101 }
1102 28 }
1103
1104 /*************************************************************************/
1105
1106 58 void MeshSurfaceIntersection::build_Weiler_model() {
1107
1108
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
 174 Stopwatch Wweiler("Weiler",verbose_);
1109
1110 Attribute<bool> corner_is_on_border(
1111
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 116 mesh_.facet_corners.attributes(), "is_on_border"
1112
1/2
✓ Branch 2 taken 58 times.
✗ Branch 3 not taken.
 58 );
1113
2/2
✓ Branch 4 taken 601587 times.
✓ Branch 5 taken 58 times.
 601645 for(index_t c: mesh_.facet_corners) {
1114
1/2
✓ Branch 2 taken 601587 times.
✗ Branch 3 not taken.
 601587 corner_is_on_border[c] = false;
1115 }
1116
1117
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 halfedges_.initialize();
1118
1119 // Step 1: Duplicate all surfaces and create alpha3 links
1120 {
1121 58 index_t nf = mesh_.facets.nb();
1122
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 mesh_.facets.create_triangles(nf);
1123
2/2
✓ Branch 0 taken 200529 times.
✓ Branch 1 taken 58 times.
 200587 for(index_t f1=0; f1<nf; ++f1) {
1124 200529 index_t f2 = f1+nf;
1125
2/4
✓ Branch 1 taken 200529 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 200529 times.
✗ Branch 5 not taken.
 200529 mesh_.facets.set_vertex(f2,0,mesh_.facets.vertex(f1,2));
1126
2/4
✓ Branch 1 taken 200529 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 200529 times.
✗ Branch 5 not taken.
 200529 mesh_.facets.set_vertex(f2,1,mesh_.facets.vertex(f1,1));
1127
2/4
✓ Branch 1 taken 200529 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 200529 times.
✗ Branch 5 not taken.
 200529 mesh_.facets.set_vertex(f2,2,mesh_.facets.vertex(f1,0));
1128
1129 // Copy attributes
1130
1/2
✓ Branch 2 taken 200529 times.
✗ Branch 3 not taken.
 200529 mesh_.facets.attributes().copy_item(f2,f1);
1131
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 200529 times.
 200529 if(interpolate_attributes_) {
1132 ✗ mesh_.facet_corners.attributes().copy_item(3*f2, 3*f1+2);
1133 ✗ mesh_.facet_corners.attributes().copy_item(3*f2+1,3*f1+1);
1134 ✗ mesh_.facet_corners.attributes().copy_item(3*f2+2,3*f1 );
1135 }
1136
1137 // Sew halfedges (need to be done *after* copy attributes,
1138 // since alpha3 is stored as an .. **attribute** !!!
1139
1/2
✓ Branch 1 taken 200529 times.
✗ Branch 2 not taken.
 200529 halfedges_.sew3(3*f1, 3*f2+1);
1140
1/2
✓ Branch 1 taken 200529 times.
✗ Branch 2 not taken.
 200529 halfedges_.sew3(3*f1+1,3*f2 );
1141
1/2
✓ Branch 1 taken 200529 times.
✗ Branch 2 not taken.
 200529 halfedges_.sew3(3*f1+2,3*f2+2);
1142 }
1143 }
1144
1145 // Memorize flipped status for facets before classification
1146 // (that changes total number of facets). Used by get_initial_facet()
1147
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 58 times.
✗ Branch 6 not taken.
 116 f_is_flipped_.bind(mesh_.facets.attributes(), "flipped");
1148
2/2
✓ Branch 4 taken 401058 times.
✓ Branch 5 taken 58 times.
 401116 for(index_t f: mesh_.facets) {
1149
1/2
✓ Branch 3 taken 401058 times.
✗ Branch 4 not taken.
 401058 f_is_flipped_[f] = (f >= mesh_.facets.nb() / 2);
1150 }
1151
1152 // Step 2: Clear all facet-facet links
1153
2/2
✓ Branch 5 taken 1203174 times.
✓ Branch 6 taken 58 times.
 1203232 for(index_t c: mesh_.facet_corners) {
1154
1/2
✓ Branch 1 taken 1203174 times.
✗ Branch 2 not taken.
 1203174 mesh_.facet_corners.set_adjacent_facet(c,NO_INDEX);
1155 }
1156
1157 // Step 3: Compute halfedges bundles and radial polylines
1158 {
1159
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
 116 Stopwatch W("Bundles",verbose_);
1160
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 radial_bundles_.initialize();
1161 58 }
1162
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 radial_polylines_.initialize();
1163
1164
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 58 times.
 58 if(skeleton_ != nullptr) {
1165 ✗ geo_assert(!dry_run_);
1166 ✗ radial_polylines_.get_skeleton(*skeleton_, skeleton_trim_fins_);
1167 }
1168
1169 // Step 4: Connect manifold edges
1170
3/4
✓ Branch 2 taken 58 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 563370 times.
✓ Branch 8 taken 58 times.
 563428 for(index_t bndl: radial_bundles_) {
1171
2/4
✓ Branch 1 taken 563370 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 563370 times.
 563370 if(radial_bundles_.nb_halfedges(bndl) == 1) {
1172 // If there is only one halfedge in the radial bundle,
1173 // this is the original surface border, then we
1174 // create the "hem" by connecting the halfedge on the
1175 // border to its mate in the duplicated facet
1176 //
1177 // Note: if an object has a tiny hole, it has the effect
1178 // of creating a single shell instead of two shells with
1179 // opposite orientation, it may be a problem sometimes...
1180 ✗ index_t h = radial_bundles_.halfedge(bndl,0);
1181 ✗ halfedges_.sew2(h,halfedges_.alpha3(h));
1182 ✗ corner_is_on_border[h] = true;
1183 ✗ corner_is_on_border[halfedges_.alpha3(h)] = true;
1184
3/4
✓ Branch 1 taken 563370 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 519698 times.
✓ Branch 4 taken 43672 times.
 563370 } else if(radial_bundles_.nb_halfedges(bndl) == 2) {
1185 // If there are two halfedges in the radial bundle,
1186 // then it is a standard manifold edge
1187
1/2
✓ Branch 1 taken 519698 times.
✗ Branch 2 not taken.
 519698 index_t h1 = radial_bundles_.halfedge(bndl,0);
1188
1/2
✓ Branch 1 taken 519698 times.
✗ Branch 2 not taken.
 519698 index_t h2 = radial_bundles_.halfedge(bndl,1);
1189
2/4
✓ Branch 1 taken 519698 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 519698 times.
✗ Branch 5 not taken.
 519698 halfedges_.sew2(h1,halfedges_.alpha3(h2));
1190 }
1191 }
1192
1193 // Step 5: get charts
1194 // After this step, charts are interconnected triangles with coherent
1195 // orientation bordered by non-manifold radial edges
1196 58 index_t nb_charts = 0;
1197 {
1198
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
 174 Stopwatch W("Charts",verbose_);
1199
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
 58 nb_charts = get_surface_connected_components(mesh_, "chart");
1200 58 }
1201
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 44 times.
 58 if(verbose_) {
1202
2/4
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
 42 Logger::out("Weiler")
1203
4/8
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 14 times.
✗ Branch 11 not taken.
 14 << "Found " << nb_charts << " charts" << std::endl;
1204 }
1205
1206 // I had before an assertion check that the number of charts is even
1207 // but I removed it:
1208 // the number of charts is not necessarily even, since there
1209 // can be "fins" (see e.g. "saturn" example that creates 5 charts).
1210
1211 // Step 6: Radial sort
1212 {
1213
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
 116 Stopwatch("Radial sort", verbose_);
1214
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 radial_polylines_.radial_sort();
1215 }
1216
1217 // Step 7: Create alpha2 links
1218
3/4
✓ Branch 2 taken 58 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 3820 times.
✓ Branch 8 taken 58 times.
 3878 for(index_t P: radial_polylines_) {
1219
4/6
✓ Branch 1 taken 3820 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 21836 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 21836 times.
✓ Branch 11 taken 3820 times.
 25656 for(index_t bndl: radial_polylines_.bundles(P)) {
1220
2/4
✓ Branch 1 taken 21836 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 21836 times.
 21836 if(!radial_bundles_.is_sorted(bndl)) {
1221 ✗ continue;
1222 }
1223
1/2
✓ Branch 1 taken 21836 times.
✗ Branch 2 not taken.
 21836 index_t N = radial_bundles_.nb_halfedges(bndl);
1224
2/2
✓ Branch 0 taken 81889 times.
✓ Branch 1 taken 21836 times.
 103725 for(index_t i=0; i<N; ++i) {
1225
2/2
✓ Branch 0 taken 60053 times.
✓ Branch 1 taken 21836 times.
 81889 index_t i_next = (i == N-1) ? 0 : i+1;
1226
2/2
✓ Branch 0 taken 21836 times.
✓ Branch 1 taken 60053 times.
 81889 index_t i_prev = (i == 0 ) ? N-1 : i-1;
1227
1/2
✓ Branch 1 taken 81889 times.
✗ Branch 2 not taken.
 81889 index_t h = radial_bundles_.halfedge(bndl,i);
1228
1/2
✓ Branch 1 taken 81889 times.
✗ Branch 2 not taken.
 81889 index_t h_next = radial_bundles_.halfedge(bndl,i_next);
1229
1/2
✓ Branch 1 taken 81889 times.
✗ Branch 2 not taken.
 81889 index_t h_prev = radial_bundles_.halfedge(bndl,i_prev);
1230
2/4
✓ Branch 1 taken 81889 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81889 times.
✗ Branch 5 not taken.
 81889 halfedges_.sew2(h,halfedges_.alpha3(h_next));
1231
2/4
✓ Branch 1 taken 81889 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81889 times.
✗ Branch 5 not taken.
 81889 halfedges_.sew2(h_prev,halfedges_.alpha3(h));
1232 }
1233 }
1234 }
1235
1236 // Step 8: Identify regions
1237 {
1238
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
 174 Stopwatch W("Regions",verbose_);
1239
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
 58 index_t nb_regions = get_surface_connected_components(mesh_,"chart");
1240
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 44 times.
 58 if(verbose_) {
1241
2/4
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
 42 Logger::out("Weiler")
1242
4/8
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 14 times.
✗ Branch 11 not taken.
 14 << "Found " << nb_regions << " regions" << std::endl;
1243 }
1244 58 }
1245
1246 58 }
1247
1248
1249 /***********************************************************************/
1250
1251 58 void MeshSurfaceIntersection::RadialBundles::initialize() {
1252 // Sorted vector of halfedges
1253 // In this step, we only insert the ones such that v2 < v1
1254 // For the ones such that v1 > v2, we can deduce the information
1255 // from the first half.
1256 // Note: some of them come from the original halfedges and some of
1257 // them from the ones we just created by duplicating the faces.
1258 58 H_.reserve(mesh_.facet_corners.nb());
1259 58 H_.resize(0);
1260
2/2
✓ Branch 5 taken 1203174 times.
✓ Branch 6 taken 58 times.
 1203232 for(index_t h: mesh_.facet_corners) {
1261
4/6
✓ Branch 1 taken 1203174 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203174 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 601587 times.
✓ Branch 7 taken 601587 times.
 1203174 if(I_.halfedges_.vertex(h,0) < I_.halfedges_.vertex(h,1)) {
1262
1/2
✓ Branch 1 taken 601587 times.
✗ Branch 2 not taken.
 601587 H_.push_back(h);
1263 }
1264 }
1265
1266
1267 // Step 2: Lexicographic sort of the H array, so that "bundles" will be
1268 // contiguous. By "bundle", I mean the set of halfedges having the same
1269 // extremities in the same order.
1270
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 GEO::sort(
1271 58 H_.begin(), H_.end(),
1272 58 [&](index_t h1, index_t h2)->bool {
1273
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 13785350 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 13785350 geo_debug_assert(h1 < mesh_.facet_corners.nb());
1274
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 13785350 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 13785350 geo_debug_assert(h2 < mesh_.facet_corners.nb());
1275 13785350 index_t v10 = I_.halfedges_.vertex(h1,0);
1276 13785350 index_t v20 = I_.halfedges_.vertex(h2,0);
1277
2/2
✓ Branch 0 taken 5245252 times.
✓ Branch 1 taken 8540098 times.
 13785350 if(v10 < v20) {
1278 5245252 return true;
1279 }
1280
2/2
✓ Branch 0 taken 4686729 times.
✓ Branch 1 taken 3853369 times.
 8540098 if(v10 > v20) {
1281 4686729 return false;
1282 }
1283 3853369 return (I_.halfedges_.vertex(h1,1) < I_.halfedges_.vertex(h2,1));
1284 }
1285 );
1286
1287 // Step 3: find "bundles" of halfedges in the sorted array. Two halfedges
1288 // are in the same bundle if they connect the same pair of vertices (in
1289 // the same order).
1290 // Halfedges of bundle B are H[b],H[b+1]...H[e-1]
1291 // where b = bndl_start[B] and e = bndl_start[B+1]
1292 // (there are bndl_start.size()-1 bundles)
1293 58 bndl_start_.resize(0);
1294
2/2
✓ Branch 1 taken 281685 times.
✓ Branch 2 taken 58 times.
 281743 for(index_t b=0, e=0; b<H_.size(); b=e) {
1295
1/2
✓ Branch 1 taken 281685 times.
✗ Branch 2 not taken.
 281685 bndl_start_.push_back(b);
1296
2/4
✓ Branch 1 taken 281685 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 281685 times.
✗ Branch 5 not taken.
 281685 index_t v0 = I_.halfedges_.vertex(H_[b],0);
1297
2/4
✓ Branch 1 taken 281685 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 281685 times.
✗ Branch 5 not taken.
 281685 index_t v1 = I_.halfedges_.vertex(H_[b],1);
1298 281685 e = b+1;
1299 281685 while(
1300 601587 e < H_.size() &&
1301
8/10
✓ Branch 0 taken 601529 times.
✓ Branch 1 taken 58 times.
✓ Branch 3 taken 601529 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 601529 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 524407 times.
✓ Branch 9 taken 77122 times.
✓ Branch 10 taken 319902 times.
✓ Branch 11 taken 281685 times.
 1125994 I_.halfedges_.vertex(H_[e],0) == v0 &&
1302
4/6
✓ Branch 1 taken 524407 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 524407 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 319902 times.
✓ Branch 7 taken 204505 times.
 524407 I_.halfedges_.vertex(H_[e],1) == v1
1303 ) {
1304 319902 ++e;
1305 }
1306 }
1307
1/2
✓ Branch 2 taken 58 times.
✗ Branch 3 not taken.
 58 bndl_start_.push_back(H_.size());
1308
1309 // Step 4: construct second half of the sorted halfedges array from
1310 // the first half. Remember, at this point, H contains a mixture of
1311 // original halfedges and new ones, but they are all such that v1 < v2.
1312 // We now store all the bundles such that v1 > v2 (by "mirroring" the
1313 // initial ones, that is, traversing alpha3).
1314 {
1315 58 index_t bndl_start_size = bndl_start_.size();
1316
2/2
✓ Branch 0 taken 281685 times.
✓ Branch 1 taken 58 times.
 281743 for(index_t bndl=0; bndl+1<bndl_start_size; ++bndl) {
1317 281685 index_t b = bndl_start_[bndl];
1318 281685 index_t e = bndl_start_[bndl+1];
1319
2/2
✓ Branch 0 taken 601587 times.
✓ Branch 1 taken 281685 times.
 883272 for(index_t i=b; i<e; ++i) {
1320
3/6
✓ Branch 1 taken 601587 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 601587 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 601587 times.
✗ Branch 8 not taken.
 601587 H_.push_back(I_.halfedges_.alpha3(H_[i]));
1321 }
1322
1/2
✓ Branch 2 taken 281685 times.
✗ Branch 3 not taken.
 281685 bndl_start_.push_back(H_.size());
1323 }
1324 }
1325
1326 // Step 5: chain bundles around vertices
1327 58 v_first_bndl_.assign(mesh_.vertices.nb(), NO_INDEX);
1328 58 bndl_next_around_v_.assign(bndl_start_.size()-1, NO_INDEX);
1329
2/2
✓ Branch 1 taken 563370 times.
✓ Branch 2 taken 58 times.
 563428 for(index_t bndl = 0; bndl < nb(); ++bndl) {
1330 // Skip regular bundles (that are inside charts),
1331 // we only need chaining along non-manifold radial polylines
1332
2/2
✓ Branch 1 taken 519698 times.
✓ Branch 2 taken 43672 times.
 563370 if(nb_halfedges(bndl) == 2) {
1333 519698 continue;
1334 }
1335 43672 index_t v1 = vertex(bndl,0);
1336 43672 bndl_next_around_v_[bndl] = v_first_bndl_[v1];
1337 43672 v_first_bndl_[v1] = bndl;
1338 }
1339
1340
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
 58 bndl_is_sorted_.assign(nb(), false);
1341
1342
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 if(!facet_chart_.is_bound()) {
1343
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 58 times.
✗ Branch 6 not taken.
 174 facet_chart_.bind(mesh_.facets.attributes(), "chart");
1344 }
1345
1346
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 44 times.
 58 if(I_.verbose_) {
1347
6/12
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 14 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 14 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 14 times.
✗ Branch 17 not taken.
 42 Logger::out("Weiler") << "Found " << nb() << " bundles"
1348
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 << std::endl;
1349 }
1350 58 }
1351
1352 21836 void MeshSurfaceIntersection::RadialBundles::get_sorted_incident_charts(
1353 index_t bndl, vector<ChartPos>& chart_pos
1354 ) {
1355 21836 chart_pos.resize(0);
1356
4/6
✓ Branch 1 taken 21836 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 81889 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 81889 times.
✓ Branch 10 taken 21836 times.
 103725 for(index_t h: halfedges(bndl)) {
1357
1/2
✓ Branch 1 taken 81889 times.
✗ Branch 2 not taken.
 81889 chart_pos.push_back(
1358
2/4
✓ Branch 1 taken 81889 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81889 times.
✗ Branch 5 not taken.
 81889 std::make_pair(
1359 163778 facet_chart_[I_.halfedges_.facet(h)], chart_pos.size()
1360 )
1361 );
1362 }
1363 21836 std::sort(
1364 chart_pos.begin(), chart_pos.end(),
1365 117126 [](const ChartPos& a, const ChartPos& b)->bool {
1366 117126 return a.first < b.first;
1367 }
1368 );
1369 21836 }
1370
1371 58 void MeshSurfaceIntersection::RadialPolylines::initialize() {
1372
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 B_.resize(0);
1373
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 polyline_start_.resize(0);
1374
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 polyline_start_.push_back(0);
1375
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
 58 vector<bool> bndl_visited(I_.radial_bundles_.nb(),false);
1376
3/4
✓ Branch 2 taken 58 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 563370 times.
✓ Branch 8 taken 58 times.
 563428 for(index_t bndl: I_.radial_bundles_) {
1377
1378 // Bundle already visited, or regular internal edgge
1379
7/8
✓ Branch 2 taken 523518 times.
✓ Branch 3 taken 39852 times.
✓ Branch 5 taken 523518 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 519698 times.
✓ Branch 8 taken 3820 times.
✓ Branch 9 taken 559550 times.
✓ Branch 10 taken 3820 times.
 563370 if(bndl_visited[bndl] || I_.radial_bundles_.nb_halfedges(bndl)==2) {
1380 559550 continue;
1381 }
1382
1383 // Find first bundle of the poly line:
1384 // Traverse bundles on border backwards until a
1385 // non-manifold vertex is reached or until we have
1386 // looped to our starting point.
1387 3820 index_t bndl_first = bndl;
1388 for(;;) {
1389 index_t bndl_p =
1390
1/2
✓ Branch 1 taken 14574 times.
✗ Branch 2 not taken.
 14574 I_.radial_bundles_.prev_along_polyline(bndl_first);
1391
4/4
✓ Branch 0 taken 10889 times.
✓ Branch 1 taken 3685 times.
✓ Branch 2 taken 10754 times.
✓ Branch 3 taken 135 times.
 14574 if(bndl_p == NO_INDEX || bndl_p == bndl) {
1392 break;
1393 }
1394 10754 bndl_first = bndl_p;
1395 10754 }
1396
1397 // Traverse polyline by traversing bundles on border forward
1398 // until a non-manifold vertex is reached or until we have looped
1399 // to our starting point.
1400 3820 index_t bndl_cur = bndl_first;
1401 for(;;) {
1402
1/2
✓ Branch 1 taken 21836 times.
✗ Branch 2 not taken.
 21836 B_.push_back(bndl_cur);
1403 21836 bndl_visited[bndl_cur] = true;
1404
1/2
✓ Branch 1 taken 21836 times.
✗ Branch 2 not taken.
 21836 bndl_visited[I_.radial_bundles_.opposite(bndl_cur)] = true;
1405
1/2
✓ Branch 1 taken 21836 times.
✗ Branch 2 not taken.
 21836 index_t bndl_n=I_.radial_bundles_.next_along_polyline(bndl_cur);
1406
4/4
✓ Branch 0 taken 18151 times.
✓ Branch 1 taken 3685 times.
✓ Branch 2 taken 18016 times.
✓ Branch 3 taken 135 times.
 21836 if(bndl_n == NO_INDEX || bndl_n == bndl_first) {
1407 break;
1408 }
1409 18016 bndl_cur = bndl_n;
1410 18016 }
1411
1/2
✓ Branch 2 taken 3820 times.
✗ Branch 3 not taken.
 3820 polyline_start_.push_back(B_.size());
1412 }
1413
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 44 times.
 58 if(I_.verbose_) {
1414
6/12
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 14 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 14 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 14 times.
✗ Branch 17 not taken.
 42 Logger::out("Weiler") << "Found " <<nb()<< " polylines"
1415
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 << std::endl;
1416 }
1417 58 }
1418
1419 ✗ void MeshSurfaceIntersection::RadialPolylines::get_skeleton(
1420 Mesh& skeleton, bool trim_fins
1421 ) {
1422 ✗ skeleton.clear();
1423 ✗ skeleton.vertices.set_dimension(3);
1424 Attribute<bool> new_v_selection(
1425 ✗ skeleton.vertices.attributes(), "selection"
1426 ✗ );
1427 Attribute<bool> v_selection(
1428 ✗ mesh_.vertices.attributes(), "selection"
1429 ✗ );
1430 ✗ vector<index_t> v_id(mesh_.vertices.nb(), NO_INDEX);
1431 ✗ for(index_t bndl: B_) {
1432 ✗ for(index_t lv=0; lv<2; ++lv) {
1433 ✗ index_t v = I_.radial_bundles_.vertex(bndl,lv);
1434 ✗ if(v_id[v] == NO_INDEX) {
1435 ✗ v_id[v] = skeleton.vertices.create_vertex(
1436 ✗ mesh_.vertices.point(v)
1437 );
1438 ✗ if(I_.radial_bundles_.nb_bundles_around_vertex(v) != 2) {
1439 ✗ v_selection[v] = true;
1440 ✗ new_v_selection[v_id[v]] = true;
1441 }
1442 }
1443 }
1444 }
1445 ✗ for(index_t polyline: *this) {
1446 ✗ for(index_t bndl: bundles(polyline)) {
1447 ✗ if(trim_fins && I_.radial_bundles_.nb_halfedges(bndl) < 3) {
1448 ✗ continue;
1449 }
1450 ✗ index_t v1 = I_.radial_bundles_.vertex(bndl,0);
1451 ✗ index_t v2 = I_.radial_bundles_.vertex(bndl,1);
1452 ✗ geo_debug_assert(v_id[v1] != NO_INDEX);
1453 ✗ geo_debug_assert(v_id[v2] != NO_INDEX);
1454 ✗ skeleton.edges.create_edge(v_id[v1], v_id[v2]);
1455 }
1456 }
1457 ✗ skeleton.vertices.remove_isolated();
1458 ✗ }
1459
1460 58 void MeshSurfaceIntersection::RadialPolylines::radial_sort() {
1461
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 44 times.
 58 if(I_.verbose_) {
1462
2/4
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
 42 Logger::out("Radial sort")
1463
4/8
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 14 times.
✗ Branch 11 not taken.
 14 << "Nb radial polylines:" << nb() << std::endl;
1464 }
1465
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
 116 Stopwatch W("Radial sort",I_.verbose_);
1466
1467 58 index_t nb_sorted = 0;
1468
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 index_t nb_to_sort = nb();
1469
1470 // For each polyline in parallel
1471
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
 58 parallel_for_slice(
1472 0, nb(),
1473
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 116 [&](index_t b, index_t e) {
1474
1/2
✓ Branch 1 taken 190 times.
✗ Branch 2 not taken.
 190 index_t tid = Thread::current_id();
1475
1/2
✓ Branch 1 taken 190 times.
✗ Branch 2 not taken.
 190 RadialSort RS(I_);
1476
1477 190 vector<RadialBundles::ChartPos> ref_chart_to_radial_id;
1478 190 vector<RadialBundles::ChartPos> cur_chart_to_radial_id;
1479 190 vector<index_t> bndl_h;
1480
1481 // P: current polyline
1482
2/2
✓ Branch 0 taken 3820 times.
✓ Branch 1 taken 190 times.
 4010 for(index_t P = b; P < e; ++P) {
1483
1484 3820 index_t bndl_ref = NO_INDEX; // reference bundle
1485 3820 index_t N = NO_INDEX; // nb halfedges in ref bundle
1486
1487
3/6
✓ Branch 1 taken 3820 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 3820 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 3820 times.
✗ Branch 11 not taken.
 3820 for(index_t bndl: bundles(P)) {
1488
1/2
✓ Branch 1 taken 3820 times.
✗ Branch 2 not taken.
 3820 bool OK = I_.radial_bundles_.radial_sort(bndl, RS);
1489
1/2
✓ Branch 0 taken 3820 times.
✗ Branch 1 not taken.
 3820 if(OK) {
1490
1/2
✓ Branch 1 taken 3820 times.
✗ Branch 2 not taken.
 3820 I_.radial_bundles_.get_sorted_incident_charts(
1491 bndl, ref_chart_to_radial_id
1492 );
1493 13909 for(
1494 3820 index_t i=0;
1495
2/2
✓ Branch 1 taken 10089 times.
✓ Branch 2 taken 3820 times.
 13909 i+1<ref_chart_to_radial_id.size(); ++i
1496 ) {
1497
1/2
✓ Branch 0 taken 10089 times.
✗ Branch 1 not taken.
 20178 OK = OK && (
1498
1/2
✓ Branch 1 taken 10089 times.
✗ Branch 2 not taken.
 10089 ref_chart_to_radial_id[i].first !=
1499
2/4
✓ Branch 1 taken 10089 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10089 times.
✗ Branch 4 not taken.
 10089 ref_chart_to_radial_id[i+1].first
1500 );
1501 }
1502 }
1503
1/2
✓ Branch 0 taken 3820 times.
✗ Branch 1 not taken.
 3820 if(OK) {
1504 3820 bndl_ref = bndl;
1505
1/2
✓ Branch 1 taken 3820 times.
✗ Branch 2 not taken.
 3820 N = I_.radial_bundles_.nb_halfedges(bndl);
1506 3820 break;
1507 }
1508 }
1509
1510 // Copy order to all other bundles in polyline
1511
4/6
✓ Branch 1 taken 3820 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 21836 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 21836 times.
✓ Branch 11 taken 3820 times.
 25656 for(index_t bndl: bundles(P)) {
1512
2/2
✓ Branch 0 taken 3820 times.
✓ Branch 1 taken 18016 times.
 21836 if(bndl == bndl_ref) {
1513 3820 continue;
1514 }
1515
1516
2/4
✓ Branch 1 taken 18016 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 18016 times.
 18016 if(I_.radial_bundles_.nb_halfedges(bndl)==1) {
1517 ✗ continue;
1518 }
1519
1520 // Necessary condition: reference bndl should be OK
1521 // and current bndl should have same nbr of halfedges
1522
1/2
✓ Branch 0 taken 18016 times.
✗ Branch 1 not taken.
 36032 bool OK = (bndl_ref != NO_INDEX) && (
1523
2/4
✓ Branch 1 taken 18016 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 18016 times.
✗ Branch 4 not taken.
 18016 I_.radial_bundles_.nb_halfedges(bndl)==N
1524 18016 );
1525
1526 // Now check that charts around current bndl are the
1527 // same as charts around reference bndl. It can happen
1528 // that they differ even if N matches (example21.csg)
1529
1/2
✓ Branch 0 taken 18016 times.
✗ Branch 1 not taken.
 18016 if(OK) {
1530
1/2
✓ Branch 1 taken 18016 times.
✗ Branch 2 not taken.
 18016 I_.radial_bundles_.get_sorted_incident_charts(
1531 bndl, cur_chart_to_radial_id
1532 );
1533
2/2
✓ Branch 0 taken 67980 times.
✓ Branch 1 taken 18016 times.
 85996 for(index_t i=0; i<N; ++i) {
1534
1/2
✓ Branch 0 taken 67980 times.
✗ Branch 1 not taken.
 135960 OK = OK && (
1535
1/2
✓ Branch 1 taken 67980 times.
✗ Branch 2 not taken.
 67980 cur_chart_to_radial_id[i].first ==
1536
2/4
✓ Branch 1 taken 67980 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 67980 times.
✗ Branch 4 not taken.
 67980 ref_chart_to_radial_id[i].first
1537 );
1538 }
1539 }
1540
1541
1/2
✓ Branch 0 taken 18016 times.
✗ Branch 1 not taken.
 18016 if(OK) {
1542 // Here ref_bnfl has a valid order, and bndl has
1543 // the same surrounding charts as ref_bndl,
1544 // so we can reuse the radial order computed
1545 // in ref_bndl
1546
1/2
✓ Branch 1 taken 18016 times.
✗ Branch 2 not taken.
 18016 bndl_h.assign(N, NO_INDEX);
1547
2/2
✓ Branch 0 taken 67980 times.
✓ Branch 1 taken 18016 times.
 85996 for(index_t i=0; i<N; ++i) {
1548 // This halfedge ...
1549 203940 index_t h = I_.radial_bundles_.halfedge(
1550
2/4
✓ Branch 1 taken 67980 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 67980 times.
✗ Branch 5 not taken.
 67980 bndl, cur_chart_to_radial_id[i].second
1551 );
1552 // ... goes here.
1553
2/4
✓ Branch 1 taken 67980 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 67980 times.
✗ Branch 5 not taken.
 67980 bndl_h[ref_chart_to_radial_id[i].second] = h;
1554 }
1555
1/2
✓ Branch 1 taken 18016 times.
✗ Branch 2 not taken.
 18016 I_.radial_bundles_.set_sorted_halfedges(bndl,bndl_h);
1556 } else {
1557 // Else compute the radial sort geometrically
1558 ✗ OK = I_.radial_bundles_.radial_sort(bndl, RS);
1559 // May return !OK (if not using geogram_plus) when it
1560 // cannot sort (example_022.csg and example_024.csg)
1561 ✗ if(!OK) {
1562 ✗ std::cerr << std::endl
1563 << "FATAL ERROR: "
1564 << "Did not manage to sort a bundle"
1565 ✗ << " in Polyline of length "
1566 ✗ << nb_bundles(P)
1567 ✗ << std::endl;
1568 std::cerr
1569 << "(if you reached this point, you may"
1570 ✗ << " need geogramplus, contact TESSAEL)"
1571 ✗ << std::endl;
1572 ✗ geo_assert_not_reached;
1573 }
1574 }
1575 }
1576
2/2
✓ Branch 0 taken 1500 times.
✓ Branch 1 taken 2320 times.
 3820 if(I_.fine_verbose_) {
1577 1500 ++nb_sorted;
1578
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 1488 times.
 1500 if(!(nb_sorted%100)) {
1579
2/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
 24 Logger::out("Radial sort")
1580 12 << String::format(
1581 "[%2d] %6d/%6d",
1582
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 int(tid),int(nb_sorted),int(nb_to_sort)
1583
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 )
1584
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 << std::endl;
1585 }
1586 }
1587 }
1588
2/2
✓ Branch 0 taken 49 times.
✓ Branch 1 taken 141 times.
 190 if(I_.fine_verbose_) {
1589
2/4
✓ Branch 1 taken 49 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49 times.
✗ Branch 5 not taken.
 98 Logger::out("Radial sort")
1590
2/4
✓ Branch 1 taken 49 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49 times.
✗ Branch 5 not taken.
 98 << String::format("[%2d] done",int(tid))
1591
1/2
✓ Branch 1 taken 49 times.
✗ Branch 2 not taken.
 49 << std::endl;
1592 }
1593 190 }
1594 );
1595 58 }
1596
1597 /***********************************************************************/
1598 }
1599
1600
1601 namespace {
1602 using namespace GEO;
1603
1604 /**
1605 * \brief Gets the position of the leftmost
1606 * bit set in a 32 bits integer
1607 * \param[in] x the integer
1608 * \return the position of the leftmost bit
1609 * set, or NO_INDEX if the specified integer
1610 * is zero.
1611 */
1612 30 inline index_t leftmost_bit_set(index_t x) {
1613 30 index_t result = NO_INDEX;
1614
2/2
✓ Branch 0 taken 960 times.
✓ Branch 1 taken 30 times.
 990 for(index_t i=0; i<32; ++i) {
1615
2/2
✓ Branch 0 taken 67 times.
✓ Branch 1 taken 893 times.
 960 if((x&1) != 0) {
1616 67 result = i;
1617 }
1618 960 x = x >> 1;
1619 }
1620 30 return result;
1621 }
1622 }
1623
1624 /***************************************************/
1625
1626 namespace GEO {
1627
1628 30 void MeshSurfaceIntersection::classify(const std::string& expr) {
1629
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 30 times.
 30 if(mesh_.facets.nb() == 0) {
1630 ✗ return;
1631 }
1632
2/4
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
 90 Stopwatch W("Classify", verbose_);
1633
1634 // Takes as input a Weiler model, with duplicated interfaces,
1635 // operand bit (that indices for each triangle the set of operands
1636 // it corresponds to), volumetric alpha3 links and correct facet
1637 // adjacency links. It computes the operand_inclusion_bits attribute,
1638 // that indicates for each triangle the set of operands that contains
1639 // it, then evalues the boolean expression \p expr on all facets,
1640 // and keeps the facets on the boundary of the region where \p expr
1641 // evaluates to true.
1642
1643 // Chart attribute corresponds to volumetric regions
1644
2/4
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 30 times.
✗ Branch 6 not taken.
 60 Attribute<index_t> chart(mesh_.facets.attributes(), "chart");
1645
1646 // For each facet, bit n set if facet belongs to the boundary
1647 // of operand n. There can be several bit sets if two operands
1648 // are tangent (and share facets).
1649 // For each facet pair (f, g=halfedges_.facet_alpha3(f)),
1650 // we have operand_bit[f] = operand_bit[g]
1651
2/4
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 30 times.
✗ Branch 6 not taken.
 60 Attribute<index_t> operand_bit(mesh_.facets.attributes(), "operand_bit");
1652
1653 // For each facet, bit n set if facet is inside operand n.
1654 // For each facet pair (f, g=halfedges_.facet_alpha3(f)),
1655 // we have nth-bit(operand_bit[f]) != nth-bit(operand_bit[g]),
1656 // that is, among a facet pair (f,g) on the boundary of operand n,
1657 // one of f,g is considered to be "outside" the operand, and the
1658 // other one is considered to be "inside".
1659 // The one that is considered to be "outside" is the one that belongs
1660 // to the chart that encloses the largest volume.
1661 Attribute<index_t> operand_inclusion_bits(
1662
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 60 mesh_.facets.attributes(), "operand_inclusion_bits"
1663
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 );
1664
1665 // Get nb charts and nb operands
1666 30 index_t nb_charts = 0;
1667 30 index_t nb_operands = 0;
1668
2/2
✓ Branch 4 taken 175992 times.
✓ Branch 5 taken 30 times.
 176022 for(index_t f: mesh_.facets) {
1669
1/2
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
 175992 nb_charts = std::max(nb_charts, chart[f]+1);
1670
1/2
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
 175992 nb_operands = nb_operands | operand_bit[f];
1671 }
1672
1673
1/2
✓ Branch 0 taken 30 times.
✗ Branch 1 not taken.
 30 if(nb_operands != 0) {
1674 30 nb_operands = leftmost_bit_set(nb_operands) + 1;
1675 }
1676
1677
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
 30 if(verbose_) {
1678 ✗ Logger::out("Weiler") << "nb operands=" << nb_operands << std::endl;
1679 }
1680
1681 // Get connected components, obtained by traversing all facet
1682 // adjacency links and volumetric alpha3 links
1683
1684 30 index_t nb_components = 0;
1685
1686 Attribute<index_t> facet_component(
1687
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 60 mesh_.facets.attributes(), "component"
1688
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 );
1689
1690
2/2
✓ Branch 4 taken 175992 times.
✓ Branch 5 taken 30 times.
 176022 for(index_t f: mesh_.facets) {
1691
1/2
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
 175992 facet_component[f] = NO_INDEX;
1692 }
1693 // TODO: Switch to vector when I will not need to visualize anymore
1694 // vector<index_t> facet_component(mesh_.facets.nb(), NO_INDEX);
1695
1696 // one vertex per component
1697 30 vector<index_t> component_vertex;
1698
1699 // ith bit is 1 if component is inside ith operand
1700 30 vector<index_t> component_inclusion_bits;
1701
1702 // get the connected components obtained by traversing both alpha2
1703 // and alpha3 links
1704 {
1705
2/2
✓ Branch 5 taken 175992 times.
✓ Branch 6 taken 30 times.
 176022 for(index_t f:mesh_.facets) {
1706
3/4
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 36 times.
✓ Branch 4 taken 175956 times.
 175992 if(facet_component[f] == NO_INDEX) {
1707
1708
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 36 component_vertex.push_back(mesh_.facets.vertex(f,0));
1709
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 component_inclusion_bits.push_back(0);
1710
1711
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 std::stack<index_t> S;
1712
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 facet_component[f] = nb_components;
1713
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 S.push(f);
1714
2/2
✓ Branch 1 taken 175992 times.
✓ Branch 2 taken 36 times.
 176028 while(!S.empty()) {
1715 175992 index_t f1 = S.top();
1716 175992 S.pop();
1717 {
1718
1/2
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
 175992 index_t f2 = halfedges_.facet_alpha3(f1);
1719
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 175992 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 175992 geo_debug_assert(f2 != NO_INDEX);
1720
3/4
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 55809 times.
✓ Branch 4 taken 120183 times.
 175992 if(facet_component[f2] == NO_INDEX) {
1721
2/4
✓ Branch 1 taken 55809 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55809 times.
✗ Branch 5 not taken.
 55809 facet_component[f2]=facet_component[f1];
1722
1/2
✓ Branch 1 taken 55809 times.
✗ Branch 2 not taken.
 55809 S.push(f2);
1723 }
1724 }
1725
1726
2/2
✓ Branch 0 taken 527976 times.
✓ Branch 1 taken 175992 times.
 703968 for(index_t le1=0; le1<3; ++le1) {
1727
1/2
✓ Branch 1 taken 527976 times.
✗ Branch 2 not taken.
 527976 index_t f2 = mesh_.facets.adjacent(f1,le1);
1728 527976 if(
1729
3/4
✓ Branch 0 taken 527976 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 120147 times.
✓ Branch 3 taken 407829 times.
 1055952 f2 != NO_INDEX &&
1730
3/4
✓ Branch 1 taken 527976 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 120147 times.
✓ Branch 4 taken 407829 times.
 527976 facet_component[f2] == NO_INDEX
1731 ) {
1732 #ifdef GEO_DEBUG
1733
1/2
✓ Branch 1 taken 120147 times.
✗ Branch 2 not taken.
 120147 index_t le2 = mesh_.facets.find_adjacent(f2,f1);
1734
3/10
✓ Branch 1 taken 120147 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 120147 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 120147 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 120147 geo_debug_assert(
1735 mesh_.facets.vertex(f1,le1) !=
1736 mesh_.facets.vertex(f2,le2)
1737 );
1738 #endif
1739
2/4
✓ Branch 1 taken 120147 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 120147 times.
✗ Branch 5 not taken.
 120147 facet_component[f2] = facet_component[f1];
1740
1/2
✓ Branch 1 taken 120147 times.
✗ Branch 2 not taken.
 120147 S.push(f2);
1741 }
1742 }
1743 }
1744 36 ++nb_components;
1745 36 }
1746 }
1747
1748 // Prefer original vertices for starting raytracing
1749
2/2
✓ Branch 5 taken 175992 times.
✓ Branch 6 taken 30 times.
 176022 for(index_t f: mesh_.facets) {
1750
1/2
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
 175992 index_t component = facet_component[f];
1751
1/2
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
 175992 index_t v = component_vertex[component];
1752 // If component's vertex is already an original vertex
1753 // we are done.
1754
2/4
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 175992 times.
✗ Branch 4 not taken.
 175992 if(is_original_vertex(v)) {
1755 175992 continue;
1756 }
1757 // See if we can find an original vertex in facet's vertices,
1758 // if yes, replace component's vertex.
1759 ✗ for(index_t w: mesh_.facets.vertices(f)) {
1760 ✗ if(is_original_vertex(w)) {
1761 ✗ component_vertex[component] = w;
1762 ✗ break;
1763 }
1764 }
1765 }
1766 }
1767
1768 // Compute the volume enclosed by each chart
1769
1770
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 vector<double> chart_volume(nb_charts,0.0);
1771
2/2
✓ Branch 4 taken 175992 times.
✓ Branch 5 taken 30 times.
 176022 for(index_t f: mesh_.facets) {
1772
1/2
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
 175992 vec3 p1 = mesh_.facets.point(f,0);
1773
1/2
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
 175992 vec3 p2 = mesh_.facets.point(f,1);
1774
1/2
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
 175992 vec3 p3 = mesh_.facets.point(f,2);
1775
2/4
✓ Branch 3 taken 175992 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 175992 times.
✗ Branch 7 not taken.
 175992 chart_volume[chart[f]] += dot(p1,cross(p2,p3)) / 6.0;
1776 }
1777
1778 // For each component, find the chart that encloses the largest
1779 // volume (it is the external boundary of the component)
1780
1781
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 vector<double> max_chart_volume_in_component(nb_components, 0.0);
1782 vector<index_t> chart_with_max_volume_in_component(
1783 nb_components, NO_INDEX
1784
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 );
1785
1786
2/2
✓ Branch 5 taken 175992 times.
✓ Branch 6 taken 30 times.
 176022 for(index_t f: mesh_.facets) {
1787
2/4
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 175992 times.
✗ Branch 5 not taken.
 175992 double V = chart_volume[chart[f]];
1788 351984 if( ::fabs(V) >=
1789
4/6
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 175992 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 76263 times.
✓ Branch 7 taken 99729 times.
 175992 ::fabs(max_chart_volume_in_component[facet_component[f]])
1790 ) {
1791
2/4
✓ Branch 1 taken 76263 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 76263 times.
✗ Branch 5 not taken.
 76263 max_chart_volume_in_component[facet_component[f]] = V;
1792
1/2
✓ Branch 1 taken 76263 times.
✗ Branch 2 not taken.
 76263 chart_with_max_volume_in_component[facet_component[f]] =
1793
2/4
✓ Branch 1 taken 76263 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 76263 times.
✗ Branch 5 not taken.
 76263 chart[f];
1794 }
1795 }
1796
1797
1798 // If there is more than one component, one needs to check
1799 // whether some components "float" inside other ones.
1800 // To do that, we determine the component inclusion bits
1801 // by launching a ray from a vertex of the component, and
1802 // checking parity of the number of intersections for each operand.
1803
1804
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 25 times.
 30 if(nb_components > 1) {
1805
1806 // Copy facet_component[] from mesh_ to mesh_copy_,
1807 // following original_facet_id_ links.
1808 // Note: when some input facet are overlapping and co-planar,
1809 // some original facets may become orphan (facet_component_copy_
1810 // remains NO_INDEX for those facets).
1811 // They are ignored in raytracing. This may be a problem in very
1812 // specific configurations (to be investigated).
1813 // TODO: find a way of keeping original facet <-> intersection facet
1814 // relations even when facets are merged.
1815 {
1816 Attribute<index_t> facet_component_copy(
1817
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 10 mesh_copy_.facets.attributes(), "component"
1818
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 );
1819
2/2
✓ Branch 4 taken 830 times.
✓ Branch 5 taken 5 times.
 835 for(index_t f: mesh_copy_.facets) {
1820
1/2
✓ Branch 1 taken 830 times.
✗ Branch 2 not taken.
 830 facet_component_copy[f] = NO_INDEX;
1821 }
1822
2/2
✓ Branch 4 taken 3836 times.
✓ Branch 5 taken 5 times.
 3841 for(index_t f: mesh_.facets) {
1823
1/2
✓ Branch 1 taken 3836 times.
✗ Branch 2 not taken.
 3836 index_t original_f = original_facet_id_[f];
1824
1/2
✓ Branch 1 taken 3836 times.
✗ Branch 2 not taken.
 3836 index_t component = facet_component[f];
1825
1/2
✓ Branch 1 taken 3836 times.
✗ Branch 2 not taken.
 3836 facet_component_copy[original_f] = component;
1826 }
1827 5 }
1828
1829
1830
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if(verbose_) {
1831 ✗ Logger::out("Weiler") << "Classifying " << nb_components
1832 ✗ << " components using ray tracing"
1833 ✗ << std::endl;
1834 }
1835
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 parallel_for(
1836
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 10 0, nb_components, [&](index_t component) {
1837 22 component_inclusion_bits[component] =
1838 11 compute_component_inclusion_bits(
1839 11 component, component_vertex[component]
1840 );
1841 11 }
1842 );
1843
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if(verbose_) {
1844 ✗ Logger::out("Weiler") << "Done." << std::endl;
1845 }
1846 }
1847
1848 // Compute operand inclusion bits for each facet,
1849 // by propagating component's inclusion bits
1850 // from component's external shell
1851
1852 {
1853
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 vector<index_t> visited(mesh_.facets.nb(), false);
1854
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 std::stack<index_t> S;
1855
1856
2/2
✓ Branch 5 taken 175992 times.
✓ Branch 6 taken 30 times.
 176022 for(index_t f: mesh_.facets) {
1857
1/2
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
 175992 if(chart[f] ==
1858
4/6
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 175992 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 49804 times.
✓ Branch 7 taken 126188 times.
 175992 chart_with_max_volume_in_component[facet_component[f]]
1859 ) {
1860
1/2
✓ Branch 1 taken 49804 times.
✗ Branch 2 not taken.
 49804 visited[f] = 1;
1861 99608 operand_inclusion_bits[f] =
1862
3/6
✓ Branch 1 taken 49804 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49804 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 49804 times.
✗ Branch 8 not taken.
 49804 component_inclusion_bits[facet_component[f]];
1863
1/2
✓ Branch 1 taken 49804 times.
✗ Branch 2 not taken.
 49804 S.push(f);
1864 }
1865 }
1866
1867
2/2
✓ Branch 1 taken 175992 times.
✓ Branch 2 taken 30 times.
 176022 while(!S.empty()) {
1868 175992 index_t f1 = S.top();
1869 175992 S.pop();
1870 {
1871
1/2
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
 175992 index_t f2 = halfedges_.facet_alpha3(f1);
1872
6/8
✓ Branch 0 taken 175992 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 175992 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 23739 times.
✓ Branch 6 taken 152253 times.
✓ Branch 7 taken 23739 times.
✓ Branch 8 taken 152253 times.
 175992 if(f2 != NO_INDEX && !visited[f2]) {
1873
1/2
✓ Branch 1 taken 23739 times.
✗ Branch 2 not taken.
 23739 visited[f2] = true;
1874
1/2
✓ Branch 1 taken 23739 times.
✗ Branch 2 not taken.
 23739 S.push(f2);
1875 23739 operand_inclusion_bits[f2] =
1876
3/6
✓ Branch 1 taken 23739 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23739 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 23739 times.
✗ Branch 8 not taken.
 23739 operand_inclusion_bits[f1] ^ operand_bit[f1];
1877 }
1878 }
1879
2/2
✓ Branch 0 taken 527976 times.
✓ Branch 1 taken 175992 times.
 703968 for(index_t le=0; le<3; ++le) {
1880
1/2
✓ Branch 1 taken 527976 times.
✗ Branch 2 not taken.
 527976 index_t f2 = mesh_.facets.adjacent(f1,le);
1881
6/8
✓ Branch 0 taken 527976 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 527976 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 102449 times.
✓ Branch 6 taken 425527 times.
✓ Branch 7 taken 102449 times.
✓ Branch 8 taken 425527 times.
 527976 if(f2 != NO_INDEX && !visited[f2]) {
1882
1/2
✓ Branch 1 taken 102449 times.
✗ Branch 2 not taken.
 102449 visited[f2] = true;
1883
1/2
✓ Branch 1 taken 102449 times.
✗ Branch 2 not taken.
 102449 S.push(f2);
1884
2/4
✓ Branch 1 taken 102449 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 102449 times.
✗ Branch 5 not taken.
 102449 operand_inclusion_bits[f2] = operand_inclusion_bits[f1];
1885 }
1886 }
1887 }
1888 30 }
1889
1890 // Classify facets based on ther operand inclusion bits and on the
1891 // boolean expression
1892
1893
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 vector<index_t> classify_facet(mesh_.facets.nb(), 0);
1894
3/4
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
✓ Branch 4 taken 21 times.
 30 if(expr == "intersection") {
1895 // If operation is an intersection, return the neighbors of
1896 // the facets that have all their operand inclusion bit sets.
1897 9 index_t all_bits_set = (1u << nb_operands)-1u;
1898
2/2
✓ Branch 5 taken 81288 times.
✓ Branch 6 taken 9 times.
 81297 for(index_t f: mesh_.facets) {
1899 bool flipped =
1900
2/4
✓ Branch 1 taken 81288 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81288 times.
✗ Branch 5 not taken.
 81288 (max_chart_volume_in_component[facet_component[f]] < 0.0);
1901
1/2
✓ Branch 0 taken 81288 times.
✗ Branch 1 not taken.
 81288 if(flipped) {
1902 81288 classify_facet[f] =
1903
2/4
✓ Branch 1 taken 81288 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81288 times.
✗ Branch 5 not taken.
 81288 (operand_inclusion_bits[f] == all_bits_set);
1904 } else {
1905 ✗ classify_facet[f] =
1906 (
1907 ✗ operand_inclusion_bits[halfedges_.facet_alpha3(f)] ==
1908 all_bits_set
1909 );
1910 }
1911 }
1912 } else {
1913
1914
5/14
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 9 taken 21 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 21 times.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 21 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 21 BooleanExpression E(expr == "union" ? "*" : expr);
1915
1916 // For a general operation, return the facets f for which the
1917 // expression evaluates to false on f and to true on the neighbors
1918 // of f. Rember: what we want to compute is the *boundary* of the
1919 // region defined by the boolean expression, that is, the facets
1920 // for which the result of the boolean expression changes when they
1921 // are traversed by alpha3.
1922
1923 94704 auto classify = [&](index_t f)->bool {
1924 bool flipped =
1925 94704 (max_chart_volume_in_component[facet_component[f]]<0.0);
1926 94704 index_t f_in_sets = operand_inclusion_bits[f];
1927 index_t g_in_sets = operand_inclusion_bits[
1928 94704 halfedges_.facet_alpha3(f)
1929 94704 ];
1930 return flipped
1931
4/6
✓ Branch 0 taken 94704 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 23506 times.
✓ Branch 4 taken 71198 times.
✓ Branch 6 taken 23506 times.
✗ Branch 7 not taken.
 94704 ? (E(f_in_sets) && !E(g_in_sets))
1932
0/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 94704 : (E(g_in_sets) && !E(f_in_sets));
1933 21 };
1934
1935 try {
1936
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 13 times.
 21 if(mesh_.facets.nb() > 1024) {
1937
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 parallel_for(
1938 16 0, mesh_.facets.nb(), [&](index_t f) {
1939 88700 classify_facet[f] = classify(f);
1940 88700 }
1941 );
1942 } else {
1943
2/2
✓ Branch 4 taken 6004 times.
✓ Branch 5 taken 13 times.
 6017 for(index_t f: mesh_.facets) {
1944
2/4
✓ Branch 1 taken 6004 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6004 times.
✗ Branch 5 not taken.
 6004 classify_facet[f] = classify(f);
1945 }
1946 }
1947 ✗ } catch(...) {
1948 ✗ }
1949 21 }
1950
1951 // invert classification (we mark the facets that
1952 // we want to delete)
1953
2/2
✓ Branch 4 taken 175992 times.
✓ Branch 5 taken 30 times.
 176022 for(index_t f: mesh_.facets) {
1954
2/4
✓ Branch 1 taken 175992 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 175992 times.
✗ Branch 5 not taken.
 175992 classify_facet[f] = 1u - classify_facet[f];
1955 }
1956
1957
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 mesh_.facets.delete_elements(classify_facet);
1958 //remove_fins();
1959
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 mesh_.facets.connect();
1960
1961
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
 30 if(verbose_) {
1962 ✗ Logger::out("Weiler") << "Facets classified" << std::endl;
1963 }
1964 30 }
1965
1966 11 index_t MeshSurfaceIntersection::compute_component_inclusion_bits(
1967 index_t component, index_t v
1968 ) {
1969
2/4
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 11 times.
 11 if(!is_original_vertex(v)) {
1970 ✗ return compute_component_inclusion_bits_exact(component, v);
1971 }
1972
1973 Attribute<index_t> operand_bit_copy(
1974
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 33 mesh_copy_.facets.attributes(), "operand_bit"
1975
1/2
✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
 11 );
1976
1977 Attribute<index_t> facet_component_copy(
1978
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 22 mesh_copy_.facets.attributes(), "component"
1979
1/2
✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
 11 );
1980
1981
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if(verbose_) {
1982 ✗ Logger::out("Weiler") << " component" << component << std::endl;
1983 }
1984
1985 11 index_t component_inclusion_bits = 0;
1986
1987
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 11 vec3 q1 = mesh_.vertices.point(v);
1988
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 11 vec3 q2{q1.x, q1.y, q1.z + 1e6};
1989
1990
2/2
✓ Branch 5 taken 1882 times.
✓ Branch 6 taken 11 times.
 1893 for(index_t f: mesh_copy_.facets) {
1991
1992 // TODO: understand when it can happen (probably when merging
1993 // overlapping facets, one of the merged facets is no longer seen
1994 // by the mesh_ -> mesh_copy_ original_facet relation).
1995
2/4
✓ Branch 1 taken 1882 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1882 times.
 1882 if(facet_component_copy[f] == NO_INDEX) {
1996 ✗ continue;
1997 }
1998
1999
3/4
✓ Branch 1 taken 1882 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 830 times.
✓ Branch 4 taken 1052 times.
 1882 if(facet_component_copy[f] == component) {
2000 830 continue;
2001 }
2002
2003
1/2
✓ Branch 1 taken 1052 times.
✗ Branch 2 not taken.
 1052 index_t facet_operand_bits = operand_bit_copy[f];
2004
1/2
✓ Branch 1 taken 1052 times.
✗ Branch 2 not taken.
 1052 const vec3& p1 = mesh_copy_.facets.point(f,0);
2005
1/2
✓ Branch 1 taken 1052 times.
✗ Branch 2 not taken.
 1052 const vec3& p2 = mesh_copy_.facets.point(f,1);
2006
1/2
✓ Branch 1 taken 1052 times.
✗ Branch 2 not taken.
 1052 const vec3& p3 = mesh_copy_.facets.point(f,2);
2007
3/4
✓ Branch 1 taken 1052 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 1049 times.
 1052 if(segment_triangle_intersection_SOS(q1,q2,p1,p2,p3)) {
2008 3 component_inclusion_bits ^= facet_operand_bits;
2009 }
2010 }
2011 11 return component_inclusion_bits;
2012 11 }
2013
2014 ✗ index_t MeshSurfaceIntersection::compute_component_inclusion_bits_exact(
2015 index_t component, index_t v
2016 ) {
2017 Attribute<index_t> operand_bit_copy(
2018 ✗ mesh_copy_.facets.attributes(), "operand_bit"
2019 ✗ );
2020
2021 Attribute<index_t> facet_component_copy(
2022 ✗ mesh_copy_.facets.attributes(), "component"
2023 ✗ );
2024
2025 ✗ if(verbose_) {
2026 ✗ Logger::out("Weiler") << " componentE" << component << std::endl;
2027 }
2028
2029 ✗ index_t component_inclusion_bits = 0;
2030
2031 ✗ ExactPoint q1 = exact_vertex(v);
2032 ✗ ExactPoint q2 = q1;
2033 ✗ if(q2.w.is_one()) {
2034 ✗ q2.z += exact::scalar(1e6);
2035 } else {
2036 ✗ q2.z += (q2.w * exact::scalar(1e6));
2037 }
2038
2039 ✗ for(index_t f: mesh_copy_.facets) {
2040 ✗ if(facet_component_copy[f] == NO_INDEX) {
2041 ✗ continue;
2042 }
2043
2044 ✗ if(facet_component_copy[f] == component) {
2045 ✗ continue;
2046 }
2047
2048 ✗ index_t facet_operand_bits = operand_bit_copy[f];
2049 ✗ ExactPoint p1(mesh_copy_.facets.point(f,0));
2050 ✗ ExactPoint p2(mesh_copy_.facets.point(f,1));
2051 ✗ ExactPoint p3(mesh_copy_.facets.point(f,2));
2052 ✗ if(segment_triangle_intersection_SOS(q1,q2,p1,p2,p3)) {
2053 ✗ component_inclusion_bits ^= facet_operand_bits;
2054 }
2055 ✗ }
2056 ✗ return component_inclusion_bits;
2057 ✗ }
2058
2059 index_t
2060 ✗ MeshSurfaceIntersection::compute_component_inclusion_bits_exact_exact(
2061 index_t component, index_t v
2062 ) {
2063 Attribute<index_t> operand_bit(
2064 ✗ mesh_.facets.attributes(), "operand_bit"
2065 ✗ );
2066
2067 Attribute<index_t> facet_component(
2068 ✗ mesh_.facets.attributes(), "component"
2069 ✗ );
2070
2071 ✗ if(verbose_) {
2072 ✗ Logger::out("Weiler") << " componentEE" << component << std::endl;
2073 }
2074
2075 ✗ index_t component_inclusion_bits = 0;
2076
2077 ✗ ExactPoint q1 = exact_vertex(v);
2078 ✗ ExactPoint q2 = q1;
2079 ✗ if(q2.w.is_one()) {
2080 ✗ q2.z += exact::scalar(1e6);
2081 } else {
2082 ✗ q2.z += (q2.w * exact::scalar(1e6));
2083 }
2084
2085 ✗ for(index_t f: mesh_.facets) {
2086 ✗ if(facet_component[f] == NO_INDEX) {
2087 ✗ continue;
2088 }
2089
2090 ✗ if(facet_component[f] == component) {
2091 ✗ continue;
2092 }
2093
2094 ✗ index_t facet_operand_bits = operand_bit[f];
2095 ✗ ExactPoint p1 = exact_vertex(mesh_.facets.vertex(f,0));
2096 ✗ ExactPoint p2 = exact_vertex(mesh_.facets.vertex(f,1));
2097 ✗ ExactPoint p3 = exact_vertex(mesh_.facets.vertex(f,2));
2098 ✗ if(segment_triangle_intersection_SOS(q1,q2,p1,p2,p3)) {
2099 ✗ component_inclusion_bits ^= facet_operand_bits;
2100 }
2101 ✗ }
2102 ✗ return component_inclusion_bits;
2103 ✗ }
2104
2105 /*************************************************************************/
2106
2107 57 void MeshSurfaceIntersection::simplify_coplanar_facets(
2108 double angle_tolerance
2109 ) {
2110
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 57 times.
 57 if(mesh_.facets.nb() == 0) {
2111 ✗ return;
2112 }
2113
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 57 times.
 57 if(interpolate_attributes_) {
2114 ✗ Logger::warn("Intersect")
2115 ✗ << "Cannot simplify coplanar facets with interpolated attributes"
2116 ✗ << std::endl;
2117 ✗ return;
2118 }
2119
2120
2/4
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57 times.
✗ Branch 5 not taken.
 171 Stopwatch W("Coplanar",verbose_);
2121
2/4
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 57 times.
✗ Branch 6 not taken.
 57 Attribute<index_t> facet_group(mesh_.facets.attributes(), "group");
2122 57 vector<index_t> group_facet; // one facet per group
2123
2/2
✓ Branch 4 taken 123492 times.
✓ Branch 5 taken 57 times.
 123549 for(index_t f: mesh_.facets) {
2124
1/2
✓ Branch 1 taken 123492 times.
✗ Branch 2 not taken.
 123492 facet_group[f] = NO_INDEX;
2125 }
2126
2/4
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 57 times.
✗ Branch 6 not taken.
 57 Attribute<bool> keep_vertex(mesh_.vertices.attributes(), "keep");
2127
2/2
✓ Branch 4 taken 61752 times.
✓ Branch 5 taken 57 times.
 61809 for(index_t v: mesh_.vertices) {
2128
1/2
✓ Branch 2 taken 61752 times.
✗ Branch 3 not taken.
 61752 keep_vertex[v] = false;
2129 }
2130 57 index_t current_group = 0;
2131 {
2132 // clear attributes -------------v
2133
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 CoplanarFacets coplanar(*this, true, angle_tolerance);
2134
2/2
✓ Branch 5 taken 123492 times.
✓ Branch 6 taken 57 times.
 123549 for(index_t f: mesh_.facets) {
2135
3/4
✓ Branch 1 taken 123492 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 73084 times.
✓ Branch 4 taken 50408 times.
 123492 if(facet_group[f] == NO_INDEX) {
2136
1/2
✓ Branch 1 taken 73084 times.
✗ Branch 2 not taken.
 73084 coplanar.get(f, current_group); // This sets facet_group_[f]
2137
1/2
✓ Branch 1 taken 73084 times.
✗ Branch 2 not taken.
 73084 coplanar.mark_vertices_to_keep();
2138
1/2
✓ Branch 1 taken 73084 times.
✗ Branch 2 not taken.
 73084 group_facet.push_back(f);
2139 73084 ++current_group;
2140 }
2141 }
2142 57 }
2143
2144
1/2
✓ Branch 2 taken 57 times.
✗ Branch 3 not taken.
 57 vector<index_t> remove_f(mesh_.facets.nb(), 0);
2145 57 index_t nb_groups = current_group;
2146
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 57 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 57 geo_debug_assert(nb_groups == group_facet.size());
2147
2148 // Avoid to have reallocations in parallel with access by
2149 // preallocating facets (we are going to create a maximum
2150 // nb of facets that corresponds to the actual nb of facets)
2151
1/2
✓ Branch 2 taken 57 times.
✗ Branch 3 not taken.
 57 mesh_.facets.reserve(mesh_.facets.nb());
2152
2153 // Triangulate coplanar facet groups in parallel
2154 57 Process::spinlock lock = GEOGRAM_SPINLOCK_INIT;
2155
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 parallel_for_slice(
2156
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 114 0, nb_groups, [&](index_t b, index_t e) {
2157 // do not clear attributes -----v
2158
1/2
✓ Branch 1 taken 228 times.
✗ Branch 2 not taken.
 228 CoplanarFacets coplanar(*this,false,angle_tolerance);
2159
2160
2/2
✓ Branch 0 taken 73084 times.
✓ Branch 1 taken 228 times.
 73312 for(index_t group=b; group<e; ++group) {
2161
2/4
✓ Branch 1 taken 73084 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 73084 times.
✗ Branch 5 not taken.
 73084 coplanar.get(group_facet[group],group);
2162
2163
3/4
✓ Branch 1 taken 73084 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 65017 times.
✓ Branch 4 taken 8067 times.
 73084 if(coplanar.nb_facets() < 2) {
2164 65017 continue;
2165 }
2166
2167
1/2
✓ Branch 1 taken 8067 times.
✗ Branch 2 not taken.
 8067 coplanar.triangulate();
2168
2169 8067 bool OK = true;
2170
3/4
✓ Branch 1 taken 35610 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 27543 times.
✓ Branch 4 taken 8067 times.
 35610 for(index_t t=0; t<coplanar.CDT.nT(); ++t) {
2171 // If one of v1,v2,v3 is NO_INDEX,
2172 // it means that v1,v2 or v3 was one of the four
2173 // vertices of the external quad.
2174 // It means that there was probably an
2175 // inside/outside classification error.
2176
2177
2/4
✓ Branch 1 taken 27543 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27543 times.
✗ Branch 5 not taken.
 27543 index_t v1=coplanar.CDT.vertex_id(coplanar.CDT.Tv(t,0));
2178
2/4
✓ Branch 1 taken 27543 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27543 times.
✗ Branch 5 not taken.
 27543 index_t v2=coplanar.CDT.vertex_id(coplanar.CDT.Tv(t,1));
2179
2/4
✓ Branch 1 taken 27543 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27543 times.
✗ Branch 5 not taken.
 27543 index_t v3=coplanar.CDT.vertex_id(coplanar.CDT.Tv(t,2));
2180
2/4
✓ Branch 0 taken 27543 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 27543 times.
✗ Branch 3 not taken.
 27543 OK = OK && (v1 != NO_INDEX);
2181
2/4
✓ Branch 0 taken 27543 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 27543 times.
✗ Branch 3 not taken.
 27543 OK = OK && (v2 != NO_INDEX);
2182
2/4
✓ Branch 0 taken 27543 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 27543 times.
✗ Branch 3 not taken.
 27543 OK = OK && (v3 != NO_INDEX);
2183 }
2184
2185
1/2
✓ Branch 0 taken 8067 times.
✗ Branch 1 not taken.
 8067 if(OK) {
2186
1/2
✓ Branch 1 taken 8067 times.
✗ Branch 2 not taken.
 8067 coplanar.mark_facets(remove_f);
2187 8067 Process::acquire_spinlock(lock);
2188
3/4
✓ Branch 1 taken 35610 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 27543 times.
✓ Branch 4 taken 8067 times.
 35610 for(index_t t=0; t<coplanar.CDT.nT(); ++t) {
2189 index_t v1=
2190
2/4
✓ Branch 1 taken 27543 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27543 times.
✗ Branch 5 not taken.
 27543 coplanar.CDT.vertex_id(coplanar.CDT.Tv(t,0));
2191 index_t v2=
2192
2/4
✓ Branch 1 taken 27543 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27543 times.
✗ Branch 5 not taken.
 27543 coplanar.CDT.vertex_id(coplanar.CDT.Tv(t,1));
2193 index_t v3=
2194
2/4
✓ Branch 1 taken 27543 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27543 times.
✗ Branch 5 not taken.
 27543 coplanar.CDT.vertex_id(coplanar.CDT.Tv(t,2));
2195 index_t new_f =
2196
1/2
✓ Branch 1 taken 27543 times.
✗ Branch 2 not taken.
 27543 mesh_.facets.create_triangle(v1,v2,v3);
2197
1/2
✓ Branch 1 taken 27543 times.
✗ Branch 2 not taken.
 27543 facet_group[new_f] = current_group;
2198 }
2199 8067 Process::release_spinlock(lock);
2200 }
2201 }
2202 228 }
2203 );
2204
2205 // Delete temporary attributes
2206
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 facet_group.destroy();
2207
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 keep_vertex.destroy();
2208
2209
1/2
✓ Branch 2 taken 57 times.
✗ Branch 3 not taken.
 57 remove_f.resize(mesh_.facets.nb(),0);
2210
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 mesh_.facets.delete_elements(remove_f);
2211
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 mesh_.facets.connect();
2212 57 }
2213
2214
2215
2216 }
2217
2218 /******************************************************************************/
2219
2220 namespace {
2221 using namespace GEO;
2222
2223 ✗ void copy_operand(Mesh& result, const Mesh& operand, index_t operand_id) {
2224 Attribute<index_t> operand_bit(
2225 ✗ result.facets.attributes(), "operand_bit"
2226 ✗ );
2227 ✗ index_t v_ofs = result.vertices.create_vertices(operand.vertices.nb());
2228 ✗ for(index_t v: operand.vertices) {
2229 ✗ result.vertices.point(v + v_ofs) = operand.vertices.point(v);
2230 }
2231 ✗ for(index_t f1: operand.facets) {
2232 ✗ index_t N = operand.facets.nb_vertices(f1);
2233 ✗ index_t f2 = result.facets.create_polygon(N);
2234 ✗ for(index_t lv=0; lv<N; ++lv) {
2235 ✗ result.facets.set_vertex(
2236 ✗ f2,lv,operand.facets.vertex(f1,lv) + v_ofs
2237 );
2238 }
2239 ✗ operand_bit[f2] = index_t(1) << operand_id;
2240 }
2241 ✗ }
2242 }
2243
2244 namespace GEO {
2245
2246 ✗ void mesh_boolean_operation(
2247 Mesh& result, const Mesh& A, const Mesh& B,
2248 const std::string& operation,
2249 MeshBooleanOperationFlags flags
2250 ) {
2251 ✗ bool verbose = ((flags & MESH_BOOL_OPS_VERBOSE) != 0);
2252 ✗ if(&result == &A) {
2253 Attribute<index_t> operand_bit(
2254 ✗ result.facets.attributes(), "operand_bit"
2255 ✗ );
2256 ✗ for(index_t f: A.facets) {
2257 ✗ operand_bit[f] = index_t(1);
2258 }
2259 ✗ copy_operand(result,B,1);
2260 ✗ } else if(&result == &B) {
2261 ✗ mesh_boolean_operation(
2262 ✗ result, B, A, (operation=="A-B") ? "B-A" : operation, verbose
2263 );
2264 ✗ return;
2265 } else {
2266 ✗ result.clear();
2267 ✗ result.vertices.set_dimension(3);
2268 ✗ copy_operand(result,A,0);
2269 ✗ copy_operand(result,B,1);
2270 }
2271 ✗ MeshSurfaceIntersection I(result);
2272 ✗ I.set_radial_sort(true);
2273 ✗ I.set_verbose(verbose);
2274 ✗ if((flags & MESH_BOOL_OPS_ATTRIBS) != 0) {
2275 ✗ I.set_interpolate_attributes(true);
2276 }
2277 ✗ if((flags & MESH_BOOL_OPS_NO_CHECK_NEIGHBORS) != 0) {
2278 ✗ I.set_detect_intersecting_neighbors(false);
2279 }
2280 ✗ I.intersect();
2281 ✗ I.classify(operation);
2282 ✗ if(
2283 ✗ (flags & MESH_BOOL_OPS_ATTRIBS) == 0 &&
2284 ✗ (flags & MESH_BOOL_OPS_NO_SIMPLIFY) == 0
2285 ) {
2286 ✗ I.simplify_coplanar_facets();
2287 }
2288 ✗ }
2289
2290 1 void mesh_remove_intersections(Mesh& M, index_t max_iter, bool verbose) {
2291 // TODO: same as tet_meshing() (compute union) ?
2292
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
 4 for(index_t k=0; k<max_iter; ++k) {
2293
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 MeshSurfaceIntersection I(M);
2294 3 I.set_radial_sort(false);
2295 3 I.set_verbose(verbose);
2296
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 I.intersect();
2297
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 mesh_repair(M);
2298 3 }
2299 1 }
2300
2301 ✗ bool mesh_facets_have_intersection(Mesh& M, index_t f1, index_t f2) {
2302 ✗ for(auto [ p1, p2, p3] : M.facets.triangle_points(f1)) {
2303 ✗ for(auto [ q1, q2, q3] : M.facets.triangle_points(f2)) {
2304 ✗ if(triangles_intersections(p1,p2,p3,q1,q2,q3)) {
2305 ✗ return true;
2306 }
2307 }
2308 }
2309 ✗ return false;
2310 }
2311 }
2312