GCC Code Coverage Report

Directory: ./
File: lib/geogram/mesh/mesh_repair.cpp
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 298 393 75.8%
Functions: 20 26 76.9%
Branches: 358 653 54.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_repair.h>
41 #include <geogram/mesh/mesh_geometry.h>
42 #include <geogram/mesh/index.h>
43 #include <geogram/mesh/mesh_halfedges.h>
44 #include <geogram/mesh/mesh_io.h>
45 #include <geogram/mesh/mesh_preprocessing.h>
46 #include <geogram/points/colocate.h>
47 #include <geogram/basic/geometry_nd.h>
48 #include <geogram/basic/stopwatch.h>
49 #include <geogram/basic/command_line.h>
50 #include <geogram/basic/argused.h>
51 #include <geogram/basic/algorithm.h>
52
53 #include <stack>
54 #include <queue>
55
56 namespace {
57
58 using namespace GEO;
59
60 /**
61 * \brief Tests whether a facet is degenerate.
62 * \param[in] M the mesh that the facet belongs to
63 * \param[in] f the index of the facet in \p M
64 * \return true if facet \p f has duplicated vertices,
65 * false otherwise
66 */
67
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 942622 times.
 942622 bool facet_is_degenerate(const Mesh& M, index_t f) {
68 index_t nb_vertices = M.facets.nb_vertices(f);
69 ✗ if(nb_vertices != 3) {
70 ✗ index_t* vertices = (index_t*)alloca(nb_vertices*sizeof(index_t));
71 ✗ for(index_t lv=0; lv<nb_vertices; ++lv) {
72 ✗ vertices[lv] = M.facets.vertex(f,lv);
73 }
74 ✗ std::sort(vertices, vertices + nb_vertices);
75 return std::unique(
76 vertices, vertices + nb_vertices
77 ✗ ) != vertices + nb_vertices;
78 }
79 index_t c1 = M.facets.corners_begin(f);
80 942622 index_t c2 = c1 + 1;
81
2/2
✓ Branch 0 taken 942521 times.
✓ Branch 1 taken 101 times.
 942622 index_t c3 = c2 + 1;
82 index_t v1 = M.facet_corners.vertex(c1);
83 index_t v2 = M.facet_corners.vertex(c2);
84 index_t v3 = M.facet_corners.vertex(c3);
85
4/4
✓ Branch 0 taken 942521 times.
✓ Branch 1 taken 101 times.
✓ Branch 2 taken 79 times.
✓ Branch 3 taken 942442 times.
 942622 return v1 == v2 || v2 == v3 || v3 == v1;
86 }
87
88 /**
89 * \brief Generates a unique ordering of the vertices of
90 * a facet.
91 * \details Shifts and inverts the order of f's vertices in
92 * such a way that f's first vertex has the smallest index
93 * and it's predecessor->successor have increasing vertex
94 * indices. This ensures that the same facet has a unique
95 * representation (used to detect duplicated facets).
96 * \param[in] M the mesh that the facet belongs to
97 * \param[in] f the index of the facet in \p M
98 * \return true if the facet was flipped, false otherwise
99 */
100
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 848937 times.
 848937 bool normalize_facet_vertices_order(Mesh& M, index_t f) {
101 index_t d = M.facets.nb_vertices(f);
102
103 // Step 1: corner-to-vertex connections
104 // ------------------------------------
105
106 // Determine index c_min of corner with smallest vertex id
107 index_t c0 = M.facets.corners_begin(f);
108 index_t c_min = c0;
109
2/2
✓ Branch 0 taken 1697874 times.
✓ Branch 1 taken 848937 times.
 5093622 for(index_t c = c0 + 1; c < M.facets.corners_end(f); ++c) {
110
2/2
✓ Branch 0 taken 685435 times.
✓ Branch 1 taken 1012439 times.
 1697874 if(M.facet_corners.vertex(c) < M.facet_corners.vertex(c_min)) {
111 c_min = c;
112 }
113 }
114
115 // Determine whether facet should be flipped
116
2/2
✓ Branch 0 taken 378988 times.
✓ Branch 1 taken 469949 times.
 848937 index_t c_prev = M.facets.prev_corner_around_facet(f, c_min);
117 index_t c_next = M.facets.next_corner_around_facet(f, c_min);
118 bool direct = (
119 M.facet_corners.vertex(c_next) >= M.facet_corners.vertex(c_prev)
120 );
121
122 // Assign corner-to-vertex links
123 848937 index_t* f_vertex = (index_t*) alloca(sizeof(index_t) * d);
124 {
125 index_t c = c_min;
126
2/2
✓ Branch 0 taken 2546811 times.
✓ Branch 1 taken 848937 times.
 3395748 for(index_t i = 0; i < d; i++) {
127
2/2
✓ Branch 0 taken 1278246 times.
✓ Branch 1 taken 1268565 times.
 2546811 f_vertex[i] = M.facet_corners.vertex(c);
128
2/2
✓ Branch 0 taken 1278246 times.
✓ Branch 1 taken 1268565 times.
 2546811 c = direct ? M.facets.next_corner_around_facet(f, c)
129 1268565 : M.facets.prev_corner_around_facet(f, c);
130 }
131 }
132
2/2
✓ Branch 0 taken 2546811 times.
✓ Branch 1 taken 848937 times.
 3395748 for(index_t i = 0; i < d; i++) {
133 2546811 M.facet_corners.set_vertex(c0 + i, f_vertex[i]);
134 }
135
136 // Step 2: permute corner attributes, using the function
137 // that swaps attributes between two elements.
138 // -----------------------------------------------------
139
140 // Compute permutation P and inverse permutation Pinv
141 // P[i]: from where we fetch the attributes of the i-th corner
142 // P_inv[i]: where we want to put the attributes of the i-th corner
143
144 848937 index_t* P = (index_t*) alloca(sizeof(index_t) * d);
145 848937 index_t* P_inv = (index_t*) alloca(sizeof(index_t) * d);
146 {
147 848937 index_t cur = c_min - c0;
148
2/2
✓ Branch 0 taken 2546811 times.
✓ Branch 1 taken 848937 times.
 3395748 for(index_t i=0; i<d; ++i) {
149 2546811 P[i] = cur;
150
2/2
✓ Branch 0 taken 1278246 times.
✓ Branch 1 taken 1268565 times.
 2546811 if(direct) {
151
2/2
✓ Branch 0 taken 852164 times.
✓ Branch 1 taken 426082 times.
 1278246 cur = (cur == d-1) ? 0 : cur+1;
152 } else {
153
2/2
✓ Branch 0 taken 422855 times.
✓ Branch 1 taken 845710 times.
 1268565 cur = (cur == 0) ? d-1 : cur-1;
154 }
155 }
156
2/2
✓ Branch 0 taken 2546811 times.
✓ Branch 1 taken 848937 times.
 3395748 for(index_t i=0; i<d; ++i) {
157 2546811 P_inv[P[i]] = i;
158 }
159 }
160
161 // Permute attributes (and update P and P_inv accordingly)
162
2/2
✓ Branch 0 taken 1697874 times.
✓ Branch 1 taken 848937 times.
 2546811 for(index_t i=0; i<d-1; ++i) {
163 1697874 index_t j = P[i];
164 1697874 index_t j_inv = P_inv[i];
165
2/2
✓ Branch 0 taken 860821 times.
✓ Branch 1 taken 837053 times.
 1697874 if(i != j) {
166 860821 M.facet_corners.attributes().swap_items(c0+i,c0+j);
167 }
168
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1697874 times.
 1697874 std::swap(P[i],P[j_inv]);
169 1697874 std::swap(P_inv[i], P_inv[j]);
170
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1697874 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1697874 geo_assert(P[i] == i);
171
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1697874 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1697874 geo_assert(P_inv[i] == i);
172 }
173
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 848937 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 848937 geo_assert(P[d-1] == d-1);
174
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 848937 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 848937 geo_assert(P_inv[d-1] == d-1);
175 848937 return !direct;
176 }
177
178 /**
179 * \brief Comparator class for sorting facets.
180 */
181 class CompareFacets {
182 public:
183 /**
184 * \brief Constructs a new CompareFacets.
185 * \param[in] M the mesh
186 */
187 250 explicit CompareFacets(const Mesh& M) :
188
1/2
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
 250 mesh_(M) {
189 }
190
191 /**
192 * \brief Tests the lexicographic order of two facets by their indices.
193 * \param[in] f1 index of the first facet
194 * \param[in] f2 index of the second facet
195 * \return true if facet \p f1 is before facet \p f2 according to
196 * the lexicographic order of its vertices, false otherwise.
197 */
198 16089472 bool is_before(index_t f1, index_t f2) const {
199
1/2
✓ Branch 0 taken 16089472 times.
✗ Branch 1 not taken.
 16089472 index_t c1 = mesh_.facets.corners_begin(f1);
200 index_t c2 = mesh_.facets.corners_begin(f2);
201 while(
202
3/4
✓ Branch 0 taken 17650196 times.
✓ Branch 1 taken 40245 times.
✓ Branch 2 taken 17650196 times.
✗ Branch 3 not taken.
 35340637 c1 != mesh_.facets.corners_end(f1) &&
203 c2 != mesh_.facets.corners_end(f2)
204 ) {
205 index_t v1 = mesh_.facet_corners.vertex(c1);
206 index_t v2 = mesh_.facet_corners.vertex(c2);
207
2/2
✓ Branch 0 taken 3823127 times.
✓ Branch 1 taken 13827069 times.
 17650196 if(v1 > v2) {
208 return false;
209 }
210
2/2
✓ Branch 0 taken 1600969 times.
✓ Branch 1 taken 12226100 times.
 13827069 if(v1 < v2) {
211 return true;
212 }
213 1600969 c1++;
214 1600969 c2++;
215 }
216 return (
217
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 40245 times.
✓ Branch 2 taken 40245 times.
✗ Branch 3 not taken.
 80490 c1 == mesh_.facets.corners_end(f1) &&
218 c2 != mesh_.facets.corners_end(f2)
219 ) ;
220 }
221
222 /**
223 * \brief Tests whether two facets are identical.
224 * \param[in] f1 index of the first facet
225 * \param[in] f2 index of the second facet
226 * \return true if facets \p f1 and \p f2 have the same
227 * vertices, false otherwise
228 */
229 848687 bool is_same(index_t f1, index_t f2) const {
230
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 848687 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 848687 times.
 1697374 if(mesh_.facets.nb_vertices(f1) != mesh_.facets.nb_vertices(f2)) {
231 return false;
232 }
233 index_t c1 = mesh_.facets.corners_begin(f1);
234 index_t c2 = mesh_.facets.corners_begin(f2);
235
2/2
✓ Branch 0 taken 1620252 times.
✓ Branch 1 taken 29881 times.
 1650133 while(c1 != mesh_.facets.corners_end(f1)) {
236 geo_debug_assert(c2 != mesh_.facets.corners_end(f2));
237 index_t v1 = mesh_.facet_corners.vertex(c1);
238 index_t v2 = mesh_.facet_corners.vertex(c2);
239
2/2
✓ Branch 0 taken 818806 times.
✓ Branch 1 taken 801446 times.
 1620252 if(v1 != v2) {
240 return false;
241 }
242 801446 c1++;
243 801446 c2++;
244 }
245 return true;
246 }
247
248 /**
249 * \brief Tests the lexicographic order of two facets by their indices.
250 * \param[in] f1 index of the first facet
251 * \param[in] f2 index of the second facet
252 * \return true if facet \p f1 is before facet \p f2 according to
253 * the lexicographic order of its vertices, false otherwise.
254 */
255 bool operator() (index_t f1, index_t f2) const {
256
22/36
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 12165 times.
✓ Branch 13 taken 61133 times.
✓ Branch 14 taken 192295 times.
✓ Branch 15 taken 184318 times.
✓ Branch 16 taken 1902777 times.
✓ Branch 17 taken 848124 times.
✓ Branch 18 taken 4380783 times.
✓ Branch 19 taken 1323193 times.
✓ Branch 20 taken 5607362 times.
✓ Branch 21 taken 1323193 times.
✓ Branch 22 taken 52197 times.
✓ Branch 23 taken 47273 times.
✓ Branch 24 taken 26279 times.
✓ Branch 25 taken 25918 times.
✓ Branch 26 taken 13115 times.
✓ Branch 27 taken 12803 times.
✓ Branch 28 taken 21518 times.
✓ Branch 29 taken 25755 times.
✓ Branch 30 taken 17046 times.
✓ Branch 31 taken 8709 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✓ Branch 34 taken 563 times.
✓ Branch 35 taken 2953 times.
 15938329 return is_before(f1, f2);
257 }
258
259 private:
260 const Mesh& mesh_;
261 };
262
263
264 /**
265 * \brief Finds the non-duplicated vertices of a facet
266 * \param[in] M a const reference to a mesh
267 * \param[in] f a facet index in \p M
268 * \param[out] new_polygon on exit, where to append the
269 * non-duplicated vertices of facet \p f
270 * and a terminal NO_INDEX
271 * \retval true if the facet has three non-duplicated
272 * vertices and more
273 * \retval false otherwise
274 */
275 ✗ bool find_facet_non_duplicated_vertices(
276 const Mesh& M, index_t f, vector<index_t>& new_polygon
277 ) {
278 index_t first_corner = NO_INDEX;
279
280 // Find the first vertex that is different from
281 // its predecessor around the facet.
282 ✗ for(index_t c1: M.facets.corners(f)) {
283 index_t c2 = M.facets.next_corner_around_facet(f,c1);
284 ✗ if(M.facet_corners.vertex(c1) != M.facet_corners.vertex(c2)) {
285 first_corner = c2;
286 break;
287 }
288 }
289
290 // All the vertices may be identical (if the facet
291 // is completely degenerate).
292 ✗ if(first_corner == NO_INDEX) {
293 return false;
294 }
295
296 index_t c = first_corner;
297 index_t cur_v = NO_INDEX;
298 index_t nb = 0;
299
300 do {
301 ✗ index_t v = M.facet_corners.vertex(c);
302 ✗ if(v != cur_v) {
303 ✗ new_polygon.push_back(v);
304 cur_v = v;
305 ✗ ++nb;
306 }
307 c = M.facets.next_corner_around_facet(f,c);
308 ✗ } while(c != first_corner);
309 ✗ new_polygon.push_back(NO_INDEX);
310
311 // If there were only 2 non-duplicated vertices, then
312 // the facet is degenerate and is "rolled back" (we do
313 // not want to generate facets with two vertices only).
314 ✗ if(nb == 2) {
315 ✗ new_polygon.resize(new_polygon.size()-3);
316 ✗ return false;
317 }
318
319 return true;
320 }
321
322
323 /**
324 * \brief Detects degenerate facets in a mesh.
325 * \param[in] M the mesh
326 * \param[in] check_duplicates if true, duplicated facets are
327 * detected and all but one instance of each is marked as to be
328 * removed.
329 * \param[out] remove_f indicates for each facet whether it should be
330 * removed. If remove_f[f] != 0 if f should be removed, else f
331 * should be kept. If remove_f.size() == 0, then there is
332 * no facet to remove, else remove_f.size() == M.facets.nb().
333 * \param[out] old_polygons if non zero, on exit contains the
334 * indices of the polygonal facets that had duplicated vertices
335 * \param[out] new_polygons if non zero, on exit contains the
336 * polygonal facets to be created to replace the input polygonal
337 * facets that have duplicated vertices. Each individual facet
338 * is terminated by NO_INDEX
339 */
340 319 void detect_bad_facets(
341 Mesh& M, bool check_duplicates, vector<index_t>& remove_f,
342 vector<index_t>* old_polygons = nullptr,
343 vector<index_t>* new_polygons = nullptr,
344 bool verbose = false
345 ) {
346 index_t nb_duplicates = 0;
347 index_t nb_degenerate = 0;
348
2/2
✓ Branch 0 taken 250 times.
✓ Branch 1 taken 69 times.
 319 if(check_duplicates) {
349 vector<char> flipped(M.facets.nb());
350
351 // Used by boolean operations
352 Attribute<index_t> operand_bit;
353
1/2
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
 250 operand_bit.bind_if_is_defined(
354
2/4
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 250 times.
 500 M.facets.attributes(),"operand_bit"
355 );
356
357 // Reorder vertices around each facet to make
358 // it easier to compare two facets, and memorize
359 // initial facet orientation.
360
361
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 250 times.
 250 if(M.facets.nb() > 65535) { // Do that in parallel if mesh is large
362 ✗ parallel_for(
363 0, M.facets.nb(),
364 ✗ [&](index_t f) {
365 ✗ flipped[f] = normalize_facet_vertices_order(M, f);
366 }
367 );
368 } else {
369
2/2
✓ Branch 0 taken 848937 times.
✓ Branch 1 taken 250 times.
 849187 for(index_t f: M.facets) {
370
1/2
✓ Branch 1 taken 848937 times.
✗ Branch 2 not taken.
 848937 flipped[f] = normalize_facet_vertices_order(M, f);
371 }
372 }
373
374 // Indirect-sort the facets in lexicographic
375 // order.
376 vector<index_t> f_sort(M.facets.nb());
377
2/2
✓ Branch 0 taken 848937 times.
✓ Branch 1 taken 250 times.
 849187 for(index_t f: M.facets) {
378 848937 f_sort[f] = f;
379 }
380 CompareFacets compare_facets(M);
381
1/2
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
 250 GEO::sort(f_sort.begin(), f_sort.end(), compare_facets);
382
383 // Now f_sort[0] ... fsort[nb_facets-1] contains the indices
384 // of the sorted facets. This ensures that the indices of the
385 // facets with the same vertices (i.e. duplicated facets)
386 // appear at contiguous sequences in fsort.
387
388 // Traverse in fsort the sequences of duplicate facets.
389 // The algorithm detects the sequence of indices
390 // f_sort[if1] ... f_sort[if2-1] that contain facets
391 // with the same indices.
392 index_t if1 = 0;
393
2/2
✓ Branch 0 taken 819056 times.
✓ Branch 1 taken 250 times.
 819306 while(if1 < M.facets.nb()) {
394 819056 index_t if2 = if1 + 1;
395 819056 while(
396
2/2
✓ Branch 0 taken 848687 times.
✓ Branch 1 taken 250 times.
 848937 if2 < M.facets.nb() &&
397
2/2
✓ Branch 0 taken 29881 times.
✓ Branch 1 taken 818806 times.
 848687 compare_facets.is_same(f_sort[if1], f_sort[if2])
398 ) {
399
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 29881 times.
 29881 nb_duplicates++;
400 // Tag all facets in f_sort[if1+1] ... f_sort[if2-1] as
401 // 'to be removed' (because they all have the same vertices
402 // as f_sort[if1]).
403
2/2
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 29855 times.
 29881 if(remove_f.size() == 0) {
404
1/4
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 26 remove_f.resize(M.facets.nb(), 0);
405 }
406
1/2
✓ Branch 0 taken 29881 times.
✗ Branch 1 not taken.
 29881 remove_f[f_sort[if2]] = 1;
407 // Used by boolean operations
408 // ^= instead of |= because there can be "fins" in the
409 // input of boolean operators, and ^= discards duplicated
410 // facets in fins.
411 if(operand_bit.is_bound()) {
412 29881 operand_bit[f_sort[if1]] ^= operand_bit[f_sort[if2]];
413 }
414 29881 if2++;
415 }
416 if1 = if2;
417 }
418
419 // Restore initial facets orientation
420
2/2
✓ Branch 0 taken 848937 times.
✓ Branch 1 taken 250 times.
 849187 for(index_t f: M.facets) {
421
2/2
✓ Branch 0 taken 422855 times.
✓ Branch 1 taken 426082 times.
 848937 if(flipped[f]) {
422
1/2
✓ Branch 1 taken 422855 times.
✗ Branch 2 not taken.
 422855 M.facets.flip(f);
423 }
424 }
425 }
426
427 // Now, we tag the degenerate facets as 'to be removed'. A
428 // facet is degenerate if it is incident to the same vertex several
429 // times.
430
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 972503 times.
✓ Branch 2 taken 972503 times.
✓ Branch 3 taken 319 times.
 972822 for(index_t f: M.facets) {
431 if(
432
6/6
✓ Branch 0 taken 178915 times.
✓ Branch 1 taken 793588 times.
✓ Branch 2 taken 149034 times.
✓ Branch 3 taken 29881 times.
✓ Branch 4 taken 942442 times.
✓ Branch 5 taken 180 times.
 1915125 (remove_f.size() == 0 || remove_f[f] == 0) &&
433 942622 facet_is_degenerate(M, f)
434 ) {
435
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 180 times.
 180 nb_degenerate++;
436
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 177 times.
 180 if(remove_f.size() == 0) {
437 3 remove_f.resize(M.facets.nb(), 0);
438 }
439 180 remove_f[f] = 1;
440
441 // If we found a degenerate polygonal facet and
442 // we want to regenerate a valid one:
443 if(
444 180 old_polygons != nullptr &&
445
1/4
✓ Branch 0 taken 180 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 180 new_polygons != nullptr &&
446 M.facets.nb_vertices(f) > 3
447 ) {
448 ✗ if(find_facet_non_duplicated_vertices(
449 M,f,*new_polygons
450 )) {
451 ✗ old_polygons->push_back(f);
452 }
453 }
454 }
455 }
456
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 319 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 319 if(verbose && (nb_duplicates != 0 || nb_degenerate != 0)) {
457 ✗ Logger::out("Validate")
458 << "Detected " << nb_duplicates << " duplicate and "
459 << nb_degenerate << " degenerate facets"
460 << std::endl;
461 }
462 319 }
463
464 /************************************************************************/
465
466 /**
467 * \brief Connects the facets in a mesh.
468 * \details Reconstructs the corners.adjacent_facet links.
469 * Note that the Moebius law is not respected by this
470 * function (adjacent facets may have incoherent orientations).
471 * This function outputs a mesh with possibly not coherently
472 * oriented triangles. In other words, for two
473 * corners c1, c2, if we have:
474 * - v1 = facet_corners.vertex(c1)
475 * - v2 = facet_corners.vertex(
476 * c1,facets.next_corner_around_facet(c2f(c1),c1)
477 * )
478 * - w1 = facet_corners.vertex(c2)
479 * - w2 = facet_corners.vertex(
480 * c2,facets.next_corner_around_facet(c2f(c2),c2)
481 * )
482 * then c1 and c2 are adjacent if we have:
483 * - v1=w2 and v2=w1 (as usual) or:
484 * - v1=v2 and w1=w2 ('inverted' configuration)
485 * The output of this function can be then post-processed by
486 * repair_reorient_facets_anti_moebius() to recover coherent
487 * orientations.
488 * \param[in] M the mesh to repair
489 */
490 199 void repair_connect_facets(Mesh& M) {
491 const index_t NON_MANIFOLD=index_t(-2);
492
493 // Reset all facet-facet adjacencies.
494
2/2
✓ Branch 0 taken 2233176 times.
✓ Branch 1 taken 199 times.
 2233375 for(index_t c: M.facet_corners) {
495 M.facet_corners.set_adjacent_facet(c,NO_FACET);
496 }
497
498 // For each vertex v, v2c[v] gives the index of a
499 // corner incident to vertex v.
500 vector<index_t> v2c(M.vertices.nb(),NO_CORNER);
501
502 // For each corner c, next_c_around_v[c] is the
503 // linked list of all the corners incident to
504 // vertex v.
505 vector<index_t> next_c_around_v(M.facet_corners.nb(),NO_CORNER);
506
507 // For each corner c, c2f[c] is the index of
508 // the facet incident to c (or use c/3 if
509 // M is triangulated).
510 vector<index_t> c2f;
511
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 199 times.
 199 if(!M.facets.are_simplices()) {
512 ✗ c2f.assign(M.facet_corners.nb(), NO_FACET);
513 }
514
515 // Compute v2c and next_c_around_v
516
2/2
✓ Branch 0 taken 2233176 times.
✓ Branch 1 taken 199 times.
 2233375 for(index_t c: M.facet_corners) {
517 index_t v = M.facet_corners.vertex(c);
518 2233176 next_c_around_v[c] = v2c[v];
519 2233176 v2c[v] = c;
520 }
521
522 // Compute f2c (only if M is not triangulated,
523 // because if M is triangulated, we have f2c(c) = c/3).
524
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 199 times.
 199 if(!M.facets.are_simplices()) {
525 ✗ for(index_t f: M.facets) {
526 ✗ for(index_t c: M.facets.corners(f)) {
527 ✗ c2f[c]=f;
528 }
529 }
530 }
531
532
2/2
✓ Branch 0 taken 744392 times.
✓ Branch 1 taken 199 times.
 744591 for(index_t f1: M.facets) {
533
2/2
✓ Branch 0 taken 2233176 times.
✓ Branch 1 taken 744392 times.
 2977568 for(index_t c1: M.facets.corners(f1)) {
534
2/2
✓ Branch 0 taken 1121680 times.
✓ Branch 1 taken 1111496 times.
 2233176 if(M.facet_corners.adjacent_facet(c1) == NO_FACET) {
535 index_t adj_corner = NO_CORNER;
536 index_t v1=M.facet_corners.vertex(c1);
537 index_t v2=M.facet_corners.vertex(
538 M.facets.next_corner_around_facet(f1,c1)
539 );
540
541 1121680 index_t c2 = v2c[v1];
542
543 // Lookup candidate adjacent edges from incident
544 // edges list.
545
2/2
✓ Branch 0 taken 7003983 times.
✓ Branch 1 taken 1121680 times.
 8125663 while(c2 != NO_CORNER) {
546
2/2
✓ Branch 0 taken 5882303 times.
✓ Branch 1 taken 1121680 times.
 7003983 if(c2 != c1) {
547 index_t f2 =
548
1/2
✓ Branch 0 taken 5882303 times.
✗ Branch 1 not taken.
 5882303 M.facets.are_simplices() ? c2/3 : c2f[c2];
549 index_t c3 =
550
2/2
✓ Branch 0 taken 1074136 times.
✓ Branch 1 taken 4808167 times.
 5882303 M.facets.prev_corner_around_facet(f2,c2);
551 index_t v3 = M.facet_corners.vertex(c3);
552 // Check with standard orientation.
553
2/2
✓ Branch 0 taken 1074136 times.
✓ Branch 1 taken 4808167 times.
 5882303 if(v3 == v2) {
554
2/2
✓ Branch 0 taken 245 times.
✓ Branch 1 taken 1073891 times.
 1074136 if(adj_corner == NO_CORNER) {
555 adj_corner = c3;
556 } else {
557 adj_corner = NON_MANIFOLD;
558 }
559 } else {
560 // Check with the other ("wrong") orientation
561 c3 = M.facets.next_corner_around_facet(f2,c2);
562 v3 = M.facet_corners.vertex(c3);
563
2/2
✓ Branch 0 taken 37996 times.
✓ Branch 1 taken 4770171 times.
 4808167 if(v3 == v2) {
564
2/2
✓ Branch 0 taken 163 times.
✓ Branch 1 taken 37833 times.
 37996 if(adj_corner == NO_CORNER) {
565 adj_corner = c2;
566 } else {
567 adj_corner = NON_MANIFOLD;
568 }
569 }
570 }
571 }
572 7003983 c2 = next_c_around_v[c2];
573 }
574
2/2
✓ Branch 0 taken 1111496 times.
✓ Branch 1 taken 10184 times.
 1121680 if(
575 adj_corner != NO_CORNER &&
576 adj_corner != NON_MANIFOLD
577 ) {
578 M.facet_corners.set_adjacent_facet(adj_corner,f1);
579
1/2
✓ Branch 0 taken 1111496 times.
✗ Branch 1 not taken.
 1111496 index_t f2 = M.facets.are_simplices() ?
580 adj_corner/3 :
581 1111496 c2f[adj_corner] ;
582 M.facet_corners.set_adjacent_facet(c1,f2);
583 }
584 }
585 }
586 }
587 199 }
588
589 /************************************************************************/
590
591 /**
592 * \brief Tests the relative orientation of two adjacent facets
593 * \param[in] M the mesh
594 * \param[in] f1 index of the first facet
595 * \param[in] c11 index of a corner in facet \p f1
596 * \param[in] f2 index of the second facet
597 * \return 1 if \p f1 and \p f2 have compatible orientations, -1 if
598 * they have incompatible orientations, 0 if they are not adjacent
599 */
600
1/2
✓ Branch 0 taken 1483923 times.
✗ Branch 1 not taken.
 1483923 inline signed_index_t repair_relative_orientation(
601 Mesh& M, index_t f1, index_t c11, index_t f2
602 ) {
603 index_t c12 = M.facets.next_corner_around_facet(f1, c11);
604 index_t v11 = M.facet_corners.vertex(c11);
605 index_t v12 = M.facet_corners.vertex(c12);
606
1/2
✓ Branch 0 taken 2941899 times.
✗ Branch 1 not taken.
 2941899 for(index_t c21: M.facets.corners(f2)) {
607 index_t c22 = M.facets.next_corner_around_facet(f2, c21);
608 index_t v21 = M.facet_corners.vertex(c21);
609 index_t v22 = M.facet_corners.vertex(c22);
610
2/2
✓ Branch 0 taken 2802673 times.
✓ Branch 1 taken 139226 times.
 2941899 if(v11 == v21 && v12 == v22) {
611 return -1;
612 }
613
2/2
✓ Branch 0 taken 1457976 times.
✓ Branch 1 taken 1344697 times.
 2802673 if(v11 == v22 && v12 == v21) {
614 return 1;
615 }
616 }
617 return 0;
618 }
619
620 /**
621 * \brief Removes an adjacency connections between two facets in a mesh
622 * \param[in] M the mesh
623 * \param[in] f1 index of the first facet
624 * \param[in] f2 index of the second facet
625 */
626 ✗ void repair_dissociate(
627 Mesh& M, index_t f1, index_t f2
628 ) {
629 ✗ for(index_t c: M.facets.corners(f1)) {
630 ✗ if(M.facet_corners.adjacent_facet(c) == f2) {
631 M.facet_corners.set_adjacent_facet(c, NO_FACET);
632 }
633 }
634 ✗ for(index_t c: M.facets.corners(f2)) {
635 ✗ if(M.facet_corners.adjacent_facet(c) == f1) {
636 M.facet_corners.set_adjacent_facet(c, NO_FACET);
637 }
638 }
639 ✗ }
640
641 /**
642 * \brief Greedily propagates facet reorientation in a mesh
643 * \details Whenever a Moebius loop is encountered, the involved
644 * facets are disconnected from their neighbors.
645 * \param[in] M the mesh
646 * \param[in] f index of the current facet
647 * \param[in,out] visited a vector used to mark facets that were
648 * already traversed
649 * \param[out] moebius_count number of Moebius loops encountered
650 * \param[out] moebius_facets a pointer to a vector. On exit,
651 * *moebius_facets[f] has a non-zero value if facet f is
652 * incident to an edge that could not be consistently oriented.
653 * If nullptr, then this information is not returned.
654 */
655
1/2
✓ Branch 0 taken 993991 times.
✗ Branch 1 not taken.
 993991 void repair_propagate_orientation(
656 Mesh& M, index_t f, const std::vector<bool>& visited,
657 index_t& moebius_count,
658 vector<index_t>* moebius_facets = nullptr
659 ) {
660 index_t nb_plus = 0;
661 index_t nb_minus = 0;
662
2/2
✓ Branch 0 taken 2981973 times.
✓ Branch 1 taken 993991 times.
 3975964 for(index_t c: M.facets.corners(f)) {
663 index_t f2 = M.facet_corners.adjacent_facet(c);
664
4/4
✓ Branch 0 taken 2967323 times.
✓ Branch 1 taken 14650 times.
✓ Branch 2 taken 1483923 times.
✓ Branch 3 taken 1483400 times.
 2981973 if(f2 != NO_FACET && visited[index_t(f2)]) {
665 signed_index_t ori =
666 1483923 repair_relative_orientation(M, f, c, f2);
667
2/3
✓ Branch 0 taken 1344697 times.
✓ Branch 1 taken 139226 times.
✗ Branch 2 not taken.
 1483923 switch(ori) {
668 1344697 case 1:
669 1344697 nb_plus++;
670 1344697 break;
671 139226 case -1:
672 139226 nb_minus++;
673 139226 break;
674 case 0:
675 ✗ geo_assert_not_reached;
676 }
677 }
678 }
679
3/4
✓ Branch 0 taken 900759 times.
✓ Branch 1 taken 93232 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 900759 times.
 993991 if(nb_plus != 0 && nb_minus != 0) {
680 ✗ moebius_count++;
681 ✗ if(moebius_facets != nullptr) {
682 ✗ moebius_facets->resize(M.facets.nb(), 0);
683 ✗ (*moebius_facets)[f] = 1;
684 ✗ for(index_t c: M.facets.corners(f)) {
685 index_t f2 = M.facet_corners.adjacent_facet(c);
686 ✗ if(f2 != NO_FACET) {
687 ✗ (*moebius_facets)[f2] = 1;
688 }
689 }
690 }
691 ✗ if(nb_plus > nb_minus) {
692 nb_minus = 0;
693 ✗ for(index_t c: M.facets.corners(f)) {
694 index_t f2 = M.facet_corners.adjacent_facet(c);
695 if(
696 ✗ f2 != NO_FACET && visited[f2] &&
697 ✗ repair_relative_orientation(M, f, c, f2) < 0
698 ) {
699 ✗ repair_dissociate(M, f, f2);
700 }
701 }
702 } else {
703 nb_plus = 0;
704 ✗ for(index_t c: M.facets.corners(f)) {
705 index_t f2 = M.facet_corners.adjacent_facet(c);
706 if(
707 ✗ f2 != NO_FACET && visited[index_t(f2)] &&
708 ✗ repair_relative_orientation(M, f, c, f2) > 0
709 ) {
710 ✗ repair_dissociate(M, f, f2);
711 }
712 }
713 }
714 }
715 geo_argused(nb_plus);
716
1/2
✓ Branch 0 taken 93232 times.
✗ Branch 1 not taken.
 93232 if(nb_minus != 0) {
717 93232 M.facets.flip(f);
718 }
719 993991 }
720
721 /**
722 * \brief Tests whether a facet of a mesh is on the border.
723 * \param[in] M the mesh
724 * \param[in] f index of the facet
725 * \return true if \p f is on the border of \p M, false otherwise
726 */
727
1/2
✓ Branch 0 taken 994178 times.
✗ Branch 1 not taken.
 994178 bool facet_is_on_border(Mesh& M, index_t f) {
728
2/2
✓ Branch 0 taken 2965806 times.
✓ Branch 1 taken 979726 times.
 3945532 for(index_t c: M.facets.corners(f)) {
729
2/2
✓ Branch 0 taken 2951354 times.
✓ Branch 1 taken 14452 times.
 2965806 if(M.facet_corners.adjacent_facet(c) == NO_FACET) {
730 return true;
731 }
732 }
733 return false;
734 }
735
736 /**
737 * \brief Used to represent graph distance to border.
738 * Since values are clamped to a small number (typically 5),
739 * this fits in a single byte.
740 */
741 typedef Numeric::uint8 facet_distance_t;
742
743 /**
744 * \brief Computes for each facet its facet-graph distance to
745 * the border of the mesh, clamped to max_iter.
746 * \param[in] M the mesh
747 * \param[out] D for each facet, its graph distance to the border
748 * \param[in] max_iter maximumm number of iterations (determines the
749 * largest possible computed graph distance).
750 */
751 203 void compute_border_distance(
752 Mesh& M, vector<facet_distance_t>& D, index_t max_iter
753 ) {
754
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 203 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 203 geo_assert(max_iter < 256);
755 203 D.assign(M.facets.nb(), facet_distance_t(max_iter));
756
2/2
✓ Branch 0 taken 994178 times.
✓ Branch 1 taken 203 times.
 994381 for(index_t f: M.facets) {
757
2/2
✓ Branch 0 taken 14452 times.
✓ Branch 1 taken 979726 times.
 994178 if(facet_is_on_border(M, f)) {
758 14452 D[f] = facet_distance_t(0);
759 }
760 }
761
2/2
✓ Branch 0 taken 812 times.
✓ Branch 1 taken 203 times.
 1015 for(signed_index_t i = 1; i < signed_index_t(max_iter); i++) {
762
2/2
✓ Branch 0 taken 3976712 times.
✓ Branch 1 taken 812 times.
 3977524 for(index_t f: M.facets) {
763
2/2
✓ Branch 0 taken 3866806 times.
✓ Branch 1 taken 109906 times.
 3976712 if(D[f] == signed_index_t(max_iter)) {
764
2/2
✓ Branch 0 taken 11561115 times.
✓ Branch 1 taken 3832869 times.
 15393984 for(index_t c: M.facets.corners(f)) {
765 index_t g = M.facet_corners.adjacent_facet(c);
766
3/4
✓ Branch 0 taken 11561115 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 33937 times.
✓ Branch 3 taken 11527178 times.
 11561115 if(g != NO_FACET && D[g] == facet_distance_t(i - 1)) {
767 33937 D[f] = facet_distance_t(i);
768 33937 break;
769 }
770 }
771 }
772 }
773 }
774 203 }
775
776 /**
777 * \brief A priority queue specialized to
778 * the specific case where priorities can
779 * take a small number of distinct values.
780 * \details
781 * It is implemented as an array of stacks.
782 */
783 203 class SimplePriorityQueue {
784 public:
785 /**
786 * \param[in] D priorities
787 * \param[in] max_distance max value in D
788 */
789 SimplePriorityQueue(
790 const vector<facet_distance_t>& D,
791 facet_distance_t max_distance
792
1/2
✓ Branch 1 taken 203 times.
✗ Branch 2 not taken.
 203 ) :
793 stacks_(index_t(max_distance + 1)),
794 203 D_(D) {
795 }
796
797 /**
798 * \brief Pushes a facet onto the priority queue
799 * \param[in] f index of the facet to push
800 */
801 994178 void push(index_t f) {
802 geo_debug_assert(D_[f] < stacks_.size());
803
2/2
✓ Branch 0 taken 986459 times.
✓ Branch 1 taken 7719 times.
 994178 stacks_[D_[f]].push(f);
804 994178 }
805
806 /**
807 * \brief Pops a facet from the priority queue
808 * \return the index of the popped facet
809 */
810
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 994178 times.
 994178 index_t pop() {
811 158295 for(
812 994178 signed_index_t i = signed_index_t(stacks_.size()) - 1;
813
1/2
✓ Branch 0 taken 1152473 times.
✗ Branch 1 not taken.
 1152473 i >= 0; i--
814 ) {
815
2/2
✓ Branch 0 taken 994178 times.
✓ Branch 1 taken 158295 times.
 1152473 if(!stacks_[i].empty()) {
816
2/2
✓ Branch 0 taken 986459 times.
✓ Branch 1 taken 7719 times.
 994178 index_t result = stacks_[i].top();
817 stacks_[i].pop();
818 994178 return result;
819 }
820 }
821 ✗ geo_assert_not_reached;
822 }
823
824 /**
825 * \brief Tests whether this SimplePriorityQueue is empty
826 * \return true if this SimplePriorityQueue is empty, false
827 * otherwise
828 */
829 bool empty() {
830
2/2
✓ Branch 0 taken 5141935 times.
✓ Branch 1 taken 187 times.
 10284244 for(index_t i = 0; i < stacks_.size(); i++) {
831
2/2
✓ Branch 0 taken 4147757 times.
✓ Branch 1 taken 994178 times.
 5141935 if(!stacks_[i].empty()) {
832 return false;
833 }
834 }
835 return true;
836 }
837
838 private:
839 vector<std::stack<index_t> > stacks_;
840 const vector<facet_distance_t>& D_;
841 };
842
843 /**
844 * \brief Reorients the facets with a heuristic that reduces
845 * the impact of Moebius loops.
846 * \param[in,out] M the mesh to repair
847 * \param[out] moebius_facets a pointer to a vector. On exit,
848 * *moebius_facets[f] has a non-zero value if facet f is
849 * incident to an edge that could not be consistently oriented.
850 * If nullptr, then this information is not returned.
851 */
852
1/2
✓ Branch 1 taken 203 times.
✗ Branch 2 not taken.
 203 void repair_reorient_facets_anti_moebius(
853 Mesh& M, vector<index_t>* moebius_facets=nullptr
854 ) {
855 const int max_iter = 5;
856 vector<facet_distance_t> D;
857
1/4
✓ Branch 1 taken 203 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 203 std::vector<bool> visited(M.facets.nb(), false);
858
1/2
✓ Branch 1 taken 203 times.
✗ Branch 2 not taken.
 203 compute_border_distance(M, D, max_iter);
859 SimplePriorityQueue Q(D, max_iter);
860
861 203 index_t moebius_count = 0;
862 index_t nb_visited = 0;
863
2/2
✓ Branch 0 taken 1218 times.
✓ Branch 1 taken 203 times.
 1421 for(signed_index_t i = max_iter; i >= 0; i--) {
864
2/2
✓ Branch 0 taken 380959 times.
✓ Branch 1 taken 453 times.
 381412 for(index_t f: M.facets) {
865
4/4
✓ Branch 0 taken 8187 times.
✓ Branch 1 taken 372772 times.
✓ Branch 2 taken 187 times.
✓ Branch 3 taken 8000 times.
 380959 if(!visited[f] && D[f] == i) {
866
1/2
✓ Branch 1 taken 187 times.
✗ Branch 2 not taken.
 187 Q.push(f);
867 visited[f] = true;
868 187 nb_visited++;
869
2/2
✓ Branch 0 taken 994178 times.
✓ Branch 1 taken 187 times.
 994552 while(!Q.empty()) {
870
1/2
✓ Branch 1 taken 994178 times.
✗ Branch 2 not taken.
 994178 index_t f1 = Q.pop();
871
2/2
✓ Branch 0 taken 2982534 times.
✓ Branch 1 taken 994178 times.
 3976712 for(index_t c: M.facets.corners(f1)) {
872 index_t f2 = M.facet_corners.adjacent_facet(c);
873
4/4
✓ Branch 0 taken 2967846 times.
✓ Branch 1 taken 14688 times.
✓ Branch 2 taken 993991 times.
✓ Branch 3 taken 1973855 times.
 2982534 if(f2 != NO_FACET && !visited[f2]) {
874 visited[f2] = true;
875 993991 nb_visited++;
876
1/2
✓ Branch 1 taken 993991 times.
✗ Branch 2 not taken.
 993991 repair_propagate_orientation(
877 M, f2, visited,
878 moebius_count, moebius_facets
879 );
880
1/2
✓ Branch 1 taken 993991 times.
✗ Branch 2 not taken.
 993991 Q.push(f2);
881 }
882 }
883 }
884 }
885
2/2
✓ Branch 0 taken 380194 times.
✓ Branch 1 taken 765 times.
 380959 if(nb_visited == M.facets.nb()) {
886 break;
887 }
888 }
889 }
890
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 203 times.
 203 if(moebius_count != 0) {
891 ✗ Logger::out("Validate")
892 << "Encountered " << moebius_count
893 << " ambiguous facet orientation (Moebius)"
894 << std::endl;
895 }
896 203 }
897
898 /************************************************************************/
899
900 /**
901 * \brief Finds the corner by facet and vertex index
902 * \param[in] M the mesh
903 * \param[in] f the facet index
904 * \param[in] v the vertex index
905 * \return the index of the corner that corresponds to \p f and \p v
906 * \pre such a corner does not exist in \p M
907 */
908
1/2
✓ Branch 0 taken 1855840 times.
✗ Branch 1 not taken.
 1855840 inline index_t find_corner(
909 const Mesh& M, index_t f, index_t v
910 ) {
911
1/2
✓ Branch 0 taken 3713537 times.
✗ Branch 1 not taken.
 3713537 for(index_t c: M.facets.corners(f)) {
912
2/2
✓ Branch 0 taken 1855840 times.
✓ Branch 1 taken 1857697 times.
 3713537 if(M.facet_corners.vertex(c) == v) {
913 1855840 return c;
914 }
915 }
916 ✗ geo_assert_not_reached;
917 }
918
919 /**
920 * \brief Splits the non-manifold vertices
921 * \param[in] M the mesh to repair
922 */
923 199 void repair_split_non_manifold_vertices(Mesh& M, bool verbose=false) {
924 199 std::vector<bool> c_is_visited(M.facet_corners.nb(), false);
925
1/2
✓ Branch 1 taken 199 times.
✗ Branch 2 not taken.
 199 std::vector<bool> v_is_used(M.vertices.nb(), false);
926 // new vertices are stored separately to avoid
927 // too large vector growth that would occur if
928 // pushed back to M.vertices_.
929 vector<double> new_vertices;
930 index_t nb_vertices = M.vertices.nb();
931
2/2
✓ Branch 0 taken 744392 times.
✓ Branch 1 taken 199 times.
 744591 for(index_t f: M.facets) {
932
2/2
✓ Branch 0 taken 2233176 times.
✓ Branch 1 taken 744392 times.
 2977568 for(index_t c: M.facets.corners(f)) {
933
2/2
✓ Branch 0 taken 377336 times.
✓ Branch 1 taken 1855840 times.
 2233176 if(!c_is_visited[c]) {
934 index_t cur_f = f;
935 index_t cur_c = c;
936 index_t old_v = M.facet_corners.vertex(c);
937 index_t new_v = old_v;
938
2/2
✓ Branch 0 taken 168 times.
✓ Branch 1 taken 377168 times.
 377336 if(v_is_used[old_v]) {
939 new_v = nb_vertices;
940 168 nb_vertices++;
941 672 for(
942
2/2
✓ Branch 0 taken 504 times.
✓ Branch 1 taken 168 times.
 672 index_t coord = 0; coord < M.vertices.dimension();
943 coord++
944 ) {
945 new_vertices.push_back(
946
2/2
✓ Branch 0 taken 464 times.
✓ Branch 1 taken 40 times.
 504 M.vertices.point_ptr(old_v)[coord]
947 );
948 }
949 } else {
950 v_is_used[old_v] = true;
951 }
952
953 index_t count = 0;
954 for(;;) {
955
2/2
✓ Branch 0 taken 1845823 times.
✓ Branch 1 taken 377336 times.
 2223159 c_is_visited[cur_c] = true;
956 // cannot use corners.set_vertex
957 // since vertices are not created yet
958 // (would generate an assertion fail).
959 M.facet_corners.set_vertex_no_check(cur_c,new_v);
960 cur_f = M.facet_corners.adjacent_facet(cur_c);
961
2/2
✓ Branch 0 taken 1845823 times.
✓ Branch 1 taken 377336 times.
 2223159 if(cur_f == NO_FACET || cur_f == f) {
962 break;
963 }
964
1/2
✓ Branch 1 taken 1845823 times.
✗ Branch 2 not taken.
 1845823 cur_c = find_corner(M, index_t(cur_f), old_v);
965 1845823 count++;
966
1/8
✓ Branch 0 taken 1845823 times.
✗ Branch 1 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 1845823 geo_assert(count < 10000);
967 }
968
969
2/2
✓ Branch 0 taken 10184 times.
✓ Branch 1 taken 367152 times.
 377336 if(cur_f == NO_FACET) {
970 cur_f = f;
971 cur_c = c;
972 count = 0;
973 for(;;) {
974
2/2
✓ Branch 0 taken 10017 times.
✓ Branch 1 taken 10184 times.
 20201 cur_c = M.facets.prev_corner_around_facet(
975 index_t(cur_f), cur_c
976 );
977 cur_f = M.facet_corners.adjacent_facet(cur_c);
978
2/2
✓ Branch 0 taken 10017 times.
✓ Branch 1 taken 10184 times.
 20201 if(cur_f == NO_FACET) {
979 break;
980 }
981
1/2
✓ Branch 1 taken 10017 times.
✗ Branch 2 not taken.
 10017 cur_c = find_corner(M, index_t(cur_f), old_v);
982
1/2
✓ Branch 0 taken 10017 times.
✗ Branch 1 not taken.
 10017 c_is_visited[cur_c] = true;
983 // cannot use corners.set_vertex
984 // since size is not updated yet
985 // (would generate an assertion fail).
986 M.facet_corners.set_vertex_no_check(cur_c,new_v);
987 10017 count++;
988
1/8
✓ Branch 0 taken 10017 times.
✗ Branch 1 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 10017 geo_assert(count < 10000);
989 }
990 }
991 }
992 }
993 }
994
995
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 194 times.
 199 if(new_vertices.size() != 0) {
996
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if(verbose) {
997
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 5 Logger::out("Validate")
998 << "Detected non-manifold vertices" << std::endl;
999
2/6
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 10 Logger::out("Validate") << " (fixed by generating "
1000
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 << nb_vertices - M.vertices.nb()
1001 << " new vertices)"
1002 << std::endl;
1003 }
1004
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 index_t first_v = M.vertices.create_vertices(
1005 new_vertices.size() / M.vertices.dimension()
1006 );
1007
1008
2/2
✓ Branch 0 taken 504 times.
✓ Branch 1 taken 5 times.
 509 for(index_t i=0; i<new_vertices.size(); ++i) {
1009 504 M.vertices.point_ptr(first_v)[i] = new_vertices[i];
1010 }
1011 }
1012 199 }
1013 }
1014
1015 /****************************************************************************/
1016
1017 namespace GEO {
1018
1019 ✗ void mesh_connect_and_reorient_facets_no_check(
1020 Mesh& M
1021 ) {
1022 ✗ repair_connect_facets(M);
1023 ✗ repair_reorient_facets_anti_moebius(M);
1024 ✗ }
1025
1026 191 void mesh_repair(
1027 Mesh& M, MeshRepairMode mode, double colocate_epsilon
1028 ) {
1029
2/2
✓ Branch 0 taken 160 times.
✓ Branch 1 taken 31 times.
 191 bool verbose = ((mode & MESH_REPAIR_QUIET) == 0);
1030
1031 index_t nb_vertices_in = M.vertices.nb();
1032 index_t nb_facets_in = M.facets.nb();
1033
1034
2/2
✓ Branch 0 taken 160 times.
✓ Branch 1 taken 31 times.
 191 if(mode & MESH_REPAIR_COLOCATE) {
1035 160 mesh_colocate_vertices_no_check(M, colocate_epsilon, verbose);
1036 }
1037
2/2
✓ Branch 0 taken 114 times.
✓ Branch 1 taken 77 times.
 191 if(mode & MESH_REPAIR_TRIANGULATE) {
1038 114 M.facets.triangulate();
1039 }
1040 191 mesh_remove_bad_facets_no_check(
1041 M, (mode & MESH_REPAIR_DUP_F) != 0
1042 );
1043
1044 191 repair_connect_facets(M);
1045 191 repair_reorient_facets_anti_moebius(M);
1046 191 repair_split_non_manifold_vertices(M,verbose);
1047
1048
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 189 times.
 191 if(
1049 (mode & MESH_REPAIR_RECONSTRUCT) != 0
1050 ) {
1051 2 double Marea = Geom::mesh_area(M,3);
1052
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 4 remove_small_connected_components(
1053 M,
1054
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 4 CmdLine::get_arg_percent("co3ne:min_comp_area",Marea),
1055 2 CmdLine::get_arg_uint("co3ne:min_comp_facets")
1056 );
1057
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 4 fill_holes(
1058 M,
1059
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 4 CmdLine::get_arg_percent("co3ne:max_hole_area",Marea),
1060 2 CmdLine::get_arg_uint("co3ne:max_hole_edges")
1061 );
1062 // We do that one more time, to remove the small
1063 // connected components
1064 // yielded by the detected non-manifold edges.
1065
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 4 remove_small_connected_components(
1066 M,
1067
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 4 CmdLine::get_arg_percent("co3ne:min_comp_area",Marea),
1068 2 CmdLine::get_arg_uint("co3ne:min_comp_facets")
1069 );
1070
1071 // We need to do that one more time after removing the
1072 // small component, to ensure that everything is correct.
1073 2 repair_connect_facets(M);
1074 2 repair_reorient_facets_anti_moebius(M);
1075 2 repair_split_non_manifold_vertices(M,verbose);
1076
1077 }
1078
1079
2/2
✓ Branch 0 taken 162 times.
✓ Branch 1 taken 29 times.
 191 if((mode & MESH_REPAIR_QUIET) == 0) {
1080 if(
1081
4/4
✓ Branch 0 taken 129 times.
✓ Branch 1 taken 33 times.
✓ Branch 2 taken 20 times.
✓ Branch 3 taken 109 times.
 162 M.vertices.nb() != nb_vertices_in ||
1082 M.facets.nb() != nb_facets_in
1083 ) {
1084
1/2
✓ Branch 2 taken 53 times.
✗ Branch 3 not taken.
 106 M.show_stats("Validate");
1085 }
1086 }
1087
1088
1/2
✓ Branch 0 taken 191 times.
✗ Branch 1 not taken.
 191 if(M.vertices.dimension() >= 3) {
1089 191 orient_normals(M);
1090 }
1091 191 }
1092
1093 6 void mesh_postprocess_RDT(
1094 Mesh& M, bool verbose
1095 ) {
1096 6 vector<index_t> f_is_bad(M.facets.nb(), 0);
1097
1/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 6 vector<signed_index_t> v_nb_incident(M.vertices.nb(), 0);
1098
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 detect_bad_facets(M, true, f_is_bad, nullptr, nullptr, verbose);
1099 bool changed = false;
1100
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 6 do {
1101 changed = false;
1102
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 v_nb_incident.assign(M.vertices.nb(), 0);
1103
2/2
✓ Branch 0 taken 59804 times.
✓ Branch 1 taken 6 times.
 59810 for(index_t f: M.facets) {
1104
1/2
✓ Branch 0 taken 59804 times.
✗ Branch 1 not taken.
 59804 if(f_is_bad[f] == 0) {
1105
2/2
✓ Branch 0 taken 179412 times.
✓ Branch 1 taken 59804 times.
 239216 for(index_t c: M.facets.corners(f)) {
1106 179412 ++v_nb_incident[M.facet_corners.vertex(c)];
1107 }
1108 }
1109 }
1110
2/2
✓ Branch 0 taken 59804 times.
✓ Branch 1 taken 6 times.
 59810 for(index_t f: M.facets) {
1111
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 59804 times.
 59804 if(f_is_bad[f] == 0) {
1112
2/2
✓ Branch 0 taken 59804 times.
✓ Branch 1 taken 179412 times.
 239216 for(index_t c: M.facets.corners(f)) {
1113
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 179412 times.
 179412 if(v_nb_incident[M.facet_corners.vertex(c)] == 1) {
1114 ✗ f_is_bad[f] = 1;
1115 changed = true;
1116 ✗ break;
1117 }
1118 }
1119 }
1120 }
1121 } while(changed);
1122
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 M.facets.delete_elements(f_is_bad);
1123
1124
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 repair_connect_facets(M);
1125
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 repair_reorient_facets_anti_moebius(M);
1126
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 repair_split_non_manifold_vertices(M,verbose);
1127
1128
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 6 if(verbose) {
1129 ✗ M.show_stats("Validate");
1130 }
1131 6 }
1132
1133 4 void mesh_reorient(Mesh& M, vector<index_t>* moebius_facets) {
1134 4 repair_reorient_facets_anti_moebius(M, moebius_facets);
1135 4 }
1136
1137 ✗ void mesh_detect_colocated_vertices(
1138 const Mesh& M, vector<index_t>& v_colocated_index,
1139 double colocate_epsilon
1140 ) {
1141 ✗ Geom::colocate(
1142 M.vertices.point_ptr(0),
1143 coord_index_t(M.vertices.dimension()),
1144 M.vertices.nb(),
1145 v_colocated_index,
1146 colocate_epsilon
1147 );
1148 ✗ }
1149
1150 ✗ void mesh_detect_isolated_vertices(
1151 const Mesh& M, vector<index_t>& v_is_isolated
1152 ) {
1153 ✗ v_is_isolated.assign(M.vertices.nb(),1);
1154 ✗ for(index_t e: M.edges) {
1155 ✗ v_is_isolated[M.edges.vertex(e,0)] = 0;
1156 ✗ v_is_isolated[M.edges.vertex(e,1)] = 0;
1157 }
1158 ✗ for(index_t f: M.facets) {
1159 ✗ for(index_t lv=0; lv<M.facets.nb_vertices(f); ++lv) {
1160 ✗ v_is_isolated[M.facets.vertex(f,lv)] = 0;
1161 }
1162 }
1163 ✗ for(index_t c: M.cells) {
1164 ✗ for(index_t lv=0; lv<M.cells.nb_vertices(c); ++lv) {
1165 ✗ v_is_isolated[M.cells.vertex(c,lv)] = 0;
1166 }
1167 }
1168 ✗ }
1169
1170 ✗ void mesh_detect_degenerate_facets(
1171 const Mesh& M, vector<index_t>& f_is_degenerate
1172 ) {
1173 ✗ f_is_degenerate.resize(M.facets.nb());
1174 ✗ for(index_t f: M.facets) {
1175 ✗ f_is_degenerate[f] = facet_is_degenerate(M,f);
1176 }
1177 ✗ }
1178
1179
1/2
✓ Branch 0 taken 221 times.
✗ Branch 1 not taken.
 221 void mesh_colocate_vertices_no_check(
1180 Mesh& M, double colocate_epsilon, bool verbose
1181 ) {
1182 vector<index_t> old2new;
1183
1184
1/2
✓ Branch 0 taken 221 times.
✗ Branch 1 not taken.
 221 if(M.vertices.nb() == 0) {
1185 return;
1186 }
1187
1188 index_t nb_new_vertices = 0;
1189
2/2
✓ Branch 0 taken 178 times.
✓ Branch 1 taken 43 times.
 221 if(colocate_epsilon == 0.0) {
1190
1/2
✓ Branch 1 taken 178 times.
✗ Branch 2 not taken.
 178 nb_new_vertices = Geom::colocate_by_lexico_sort(
1191 M.vertices.point_ptr(0), 3, M.vertices.nb(),
1192 old2new, M.vertices.dimension()
1193 );
1194 } else {
1195
2/4
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
 129 nb_new_vertices = Geom::colocate(
1196 M.vertices.point_ptr(0), 3, M.vertices.nb(),
1197 old2new, colocate_epsilon, M.vertices.dimension()
1198 );
1199 }
1200
1201
2/2
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 179 times.
 221 if(nb_new_vertices == M.vertices.nb()) {
1202 return;
1203 }
1204
1205
2/2
✓ Branch 0 taken 29 times.
✓ Branch 1 taken 13 times.
 42 if(verbose) {
1206
2/6
✓ Branch 1 taken 29 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 29 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 58 Logger::out("Validate") << "Removed "
1207
1/2
✓ Branch 1 taken 29 times.
✗ Branch 2 not taken.
 29 << M.vertices.nb() - nb_new_vertices
1208 << " duplicated vertices" << std::endl;
1209 }
1210
1211 // Replace vertex indices for edges
1212
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
 42 for(index_t e: M.edges) {
1213 ✗ M.edges.set_vertex(e, 0, old2new[M.edges.vertex(e,0)]);
1214 ✗ M.edges.set_vertex(e, 1, old2new[M.edges.vertex(e,1)]);
1215 }
1216
1217 // Replace vertex indices for facets
1218
2/2
✓ Branch 0 taken 262272 times.
✓ Branch 1 taken 42 times.
 262314 for(index_t c: M.facet_corners) {
1219 262272 M.facet_corners.set_vertex(c, old2new[M.facet_corners.vertex(c)]);
1220 }
1221
1222 // Replace vertex indices for cells
1223
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
 42 for(index_t ce: M.cells) {
1224 ✗ for(index_t c: M.cells.corners(ce)) {
1225 ✗ M.cell_corners.set_vertex(c, old2new[M.cell_corners.vertex(c)]);
1226 }
1227 }
1228
1229 // Now old2new is "recycled" for marking vertices that
1230 // need to be removed.
1231
2/2
✓ Branch 0 taken 96398 times.
✓ Branch 1 taken 42 times.
 192838 for(index_t i = 0; i < old2new.size(); i++) {
1232
2/2
✓ Branch 0 taken 47156 times.
✓ Branch 1 taken 49242 times.
 96398 if(old2new[i] == i) {
1233 47156 old2new[i] = 0;
1234 } else {
1235 49242 old2new[i] = 1;
1236 }
1237 }
1238
1/2
✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
 42 M.vertices.delete_elements(old2new);
1239
2/2
✓ Branch 0 taken 262272 times.
✓ Branch 1 taken 42 times.
 262314 for(index_t c: M.facet_corners) {
1240 M.facet_corners.set_adjacent_facet(c, NO_INDEX);
1241 }
1242 }
1243
1244 /*************************************************************************/
1245
1246
1/2
✓ Branch 1 taken 313 times.
✗ Branch 2 not taken.
 313 void mesh_remove_bad_facets_no_check(Mesh& M, bool check_duplicates) {
1247 vector<index_t> remove_f;
1248 vector<index_t> old_polygons;
1249 vector<index_t> new_polygons;
1250
1/2
✓ Branch 1 taken 313 times.
✗ Branch 2 not taken.
 313 detect_bad_facets(
1251 M, check_duplicates, remove_f, &old_polygons, &new_polygons
1252 );
1253 index_t current_old_polygon=0;
1254
2/2
✓ Branch 0 taken 29 times.
✓ Branch 1 taken 284 times.
 313 if(remove_f.size() != 0) {
1255 // Create the new facets that correspond to input polygonal
1256 // facets that had duplicated vertices.
1257 // This needs to be done before deleting the bad facets,
1258 // else some vertices will become isolated and will be
1259 // discarded.
1260 index_t b=0;
1261 index_t e=0;
1262
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 29 times.
 29 while(b < new_polygons.size()) {
1263 ✗ while(new_polygons[e] != NO_INDEX) {
1264 ✗ ++e;
1265 }
1266 ✗ index_t new_f = M.facets.create_polygon(e-b);
1267 ✗ M.facets.attributes().copy_item(
1268 new_f, old_polygons[current_old_polygon]
1269 );
1270 ✗ ++current_old_polygon;
1271 // We created a new facet that we want to keep !!
1272 ✗ remove_f.push_back(0);
1273 ✗ for(index_t lv=0; lv<e-b; ++lv) {
1274 ✗ M.facets.set_vertex(new_f,lv,new_polygons[b+lv]);
1275 }
1276 ✗ ++e;
1277 b=e;
1278 }
1279
1/2
✓ Branch 1 taken 29 times.
✗ Branch 2 not taken.
 29 M.facets.delete_elements(remove_f);
1280 }
1281
2/2
✓ Branch 0 taken 2647914 times.
✓ Branch 1 taken 313 times.
 2648227 for(index_t c: M.facet_corners) {
1282 M.facet_corners.set_adjacent_facet(c, NO_INDEX);
1283 }
1284 313 }
1285
1286 /*************************************************************************/
1287
1288 }
1289