GCC Code Coverage Report

Directory: ./
File: lib/geogram/mesh/mesh_fill_holes.cpp
Date: 2026-09-07 02:36:43
Exec Total Coverage
Lines: 192 310 61.9%
Functions: 8 11 72.7%
Branches: 195 542 36.0%
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_fill_holes.h>
41 #include <geogram/mesh/mesh.h>
42 #include <geogram/mesh/mesh_repair.h>
43 #include <geogram/mesh/mesh_halfedges.h>
44 #include <geogram/mesh/mesh_io.h>
45 #include <geogram/mesh/index.h>
46 #include <geogram/basic/command_line.h>
47 #include <geogram/basic/logger.h>
48
49 #undef geo_debug_assert
50 #define geo_debug_assert(x) geo_assert(x)
51
52 namespace {
53
54 using namespace GEO;
55
56 /**
57 * \brief Checks whether a halfedge is
58 * incident to a vertex.
59 * \details Checks whether the origin of \p H
60 * is adjacent to \p v. In other words,
61 * returns true in one of the following
62 * configurations (the origin of \p H is
63 * denoted by \p x).
64 * \code
65 *
66 * H
67 * v--->x--->*
68 *
69 * H
70 * *--->x--->v
71 *
72 * \endcode
73 * \param[in] MH the mesh, wrapped with halfedge accessors
74 * \param[in] H the halfedge
75 * \param[in] v the index of the vertex
76 */
77 120834 bool halfedge_has_neighbor(
78 const MeshHalfedges& MH,
79 const MeshHalfedges::Halfedge& H,
80 index_t v
81 ) {
82
83
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 120834 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 120834 geo_debug_assert(MH.halfedge_is_valid(H));
84
85 120834 const Mesh& M = MH.mesh();
86 120834 index_t f = H.facet;
87 120834 index_t c = H.corner;
88
89 {
90 120834 index_t cnext = M.facets.next_corner_around_facet(f, c);
91
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 120834 times.
 120834 if(M.facet_corners.vertex(cnext) == v) {
92 ✗ return true;
93 }
94 }
95
96 {
97 120834 index_t cprev = M.facets.prev_corner_around_facet(f, c);
98
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 120834 times.
 120834 if(M.facet_corners.vertex(cprev) == v) {
99 ✗ return true;
100 }
101 }
102
103 120834 return false;
104 }
105
106 /**
107 * \brief Checks whether an halfedge exists between
108 * the two origins of two halfedges.
109 * \param[in] MH the mesh, wrapped with halfedge accessors
110 * \param[in] h1 first halfedge
111 * \param[in] h2 second halfedge
112 * \return true if an halfedge exists between the origins
113 * of \p h1 and \p h2, false otherwise
114 */
115 48019 bool halfedge_exists_between_vertices(
116 const MeshHalfedges& MH,
117 const MeshHalfedges::Halfedge& h1,
118 const MeshHalfedges::Halfedge& h2
119 ) {
120
1/2
✓ Branch 2 taken 48019 times.
✗ Branch 3 not taken.
 48019 index_t v2 = MH.mesh().facet_corners.vertex(h2.corner);
121 48019 MeshHalfedges::Halfedge H = h1;
122 do {
123
2/4
✓ Branch 1 taken 120834 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 120834 times.
 120834 if(halfedge_has_neighbor(MH, H, v2)) {
124 ✗ return true;
125 }
126
3/4
✓ Branch 1 taken 120834 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 47743 times.
✓ Branch 4 taken 73091 times.
 120834 if(!MH.move_to_prev_around_vertex(H)) {
127 47743 break;
128 }
129
2/2
✓ Branch 1 taken 72815 times.
✓ Branch 2 taken 276 times.
 73091 } while(H != h1);
130 48019 return false;
131 }
132
133 /**
134 * \brief Internal representation of a Hole.
135 * \details A Hole is an ordered sequence of Halfedge.
136 */
137 typedef vector<MeshHalfedges::Halfedge> Hole;
138
139 /**
140 * \brief Computes a vector orthogonal to the border of a surface and in the
141 * tangent plane of the surface.
142 * \param[in] MH the mesh, wrapped with halfedge accessors
143 * \param[in] H the halfedge
144 * \return a 3d vector orthogonal to \p H, in the plange of the
145 * surface triangle incident to \p H and pointing towards the
146 * exterior of the surface.
147 */
148 ✗ vec3 border_normal(
149 const MeshHalfedges& MH,
150 const MeshHalfedges::Halfedge& H
151 ) {
152 ✗ const Mesh& M = MH.mesh();
153 ✗ index_t c1 = H.corner;
154 ✗ index_t f = H.facet;
155 ✗ index_t c2 = M.facets.next_corner_around_facet(f, c1);
156 ✗ vec3 E = M.facet_corners.point(c2) - M.facet_corners.point(c1);
157 ✗ vec3 N = Geom::mesh_facet_normal(M, f);
158 ✗ return cross(E, N);
159 }
160
161 /**
162 * \brief Splits a hole into two.
163 * \details This function is used recursively to triangulate the holes.
164 * \param[in] MH the mesh, wrapped with halfedge accessors
165 * \param[in] hole the hole to be split
166 * \param[out] hole1 one of the computed halves
167 * \param[out] hole2 the other half
168 * \param[in] use_normals if set, then couples of vertices
169 * (v1,v2) that have their border normals that match
170 * the new segment [v1,v2], are connected in priority.
171 * This improves the result on 788_raw_hand.off (but
172 * may cause some triangles overlap)
173 * \retval true on success
174 * \retval false otherwise (for instance, if no valid edge
175 * could be found to split the hole)
176 */
177 314 bool split_hole(
178 const MeshHalfedges& MH, const Hole& hole,
179 Hole& hole1, Hole& hole2, bool use_normals
180 ) {
181
182 // Step 0: compute normals to border
183 314 vector<vec3> N;
184
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 314 times.
 314 if(use_normals) {
185 ✗ N.assign(hole.size(), vec3(0.0, 0.0, 0.0));
186 ✗ for(index_t i = 0; i < hole.size(); i++) {
187 ✗ index_t j = i + 1;
188 ✗ if(j == hole.size()) {
189 ✗ j = 0;
190 }
191 ✗ vec3 n = border_normal(MH, hole[i]);
192 ✗ N[i] += n;
193 ✗ N[j] += n;
194 }
195 ✗ for(index_t i = 0; i < N.size(); i++) {
196 ✗ N[i] = normalize(N[i]);
197 }
198 }
199
200 // Step 1: compute total curve length and curvilinear abscissa
201
1/2
✓ Branch 2 taken 314 times.
✗ Branch 3 not taken.
 314 vector<double> s(hole.size());
202 314 double cur_s = 0.0;
203
1/2
✓ Branch 1 taken 314 times.
✗ Branch 2 not taken.
 314 s[0] = cur_s;
204
2/2
✓ Branch 1 taken 2670 times.
✓ Branch 2 taken 314 times.
 2984 for(index_t i = 1; i < hole.size(); i++) {
205
2/4
✓ Branch 1 taken 2670 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 2670 times.
✗ Branch 6 not taken.
 2670 const vec3& p1 = Geom::halfedge_vertex_from(MH.mesh(), hole[i - 1]);
206
2/4
✓ Branch 1 taken 2670 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 2670 times.
✗ Branch 6 not taken.
 2670 const vec3& p2 = Geom::halfedge_vertex_from(MH.mesh(), hole[i]);
207
1/2
✓ Branch 2 taken 2670 times.
✗ Branch 3 not taken.
 2670 cur_s += length(p2 - p1);
208
1/2
✓ Branch 1 taken 2670 times.
✗ Branch 2 not taken.
 2670 s[i] = cur_s;
209 }
210
2/4
✓ Branch 1 taken 314 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 314 times.
✗ Branch 6 not taken.
 314 const vec3& p1 = Geom::halfedge_vertex_from(
211 314 MH.mesh(), hole[hole.size() - 1]
212 );
213
2/4
✓ Branch 1 taken 314 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 314 times.
✗ Branch 6 not taken.
 314 const vec3& p2 = Geom::halfedge_vertex_from(MH.mesh(), hole[0]);
214
1/2
✓ Branch 2 taken 314 times.
✗ Branch 3 not taken.
 314 double total_length = cur_s + length(p2 - p1);
215
216 // Step 2: find best pair to connect
217 314 double best_rij = Numeric::max_float64();
218 314 signed_index_t v1 = -1;
219 314 signed_index_t v2 = -1;
220
2/2
✓ Branch 1 taken 2984 times.
✓ Branch 2 taken 314 times.
 3298 for(index_t i = 0; i < hole.size(); i++) {
221
2/2
✓ Branch 1 taken 48333 times.
✓ Branch 2 taken 2984 times.
 51317 for(index_t j = i + 2; j < hole.size(); j++) {
222
223 // Do not split using vertices
224 // already connected by an edge.
225 48647 if(
226
6/6
✓ Branch 0 taken 2356 times.
✓ Branch 1 taken 45977 times.
✓ Branch 3 taken 2042 times.
✓ Branch 4 taken 314 times.
✓ Branch 5 taken 314 times.
✓ Branch 6 taken 48019 times.
 96352 (i == 0 && j == hole.size() - 1) ||
227
4/8
✓ Branch 1 taken 48019 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48019 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 48019 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 48019 times.
 48019 halfedge_exists_between_vertices(MH, hole[i], hole[j])
228 ) {
229 314 continue;
230 }
231
232 96038 double dsij = std::min(
233
4/8
✓ Branch 1 taken 48019 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48019 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 48019 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 48019 times.
✗ Branch 11 not taken.
 48019 s[j] - s[i], total_length - (s[j] - s[i])
234 48019 );
235
2/4
✓ Branch 1 taken 48019 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 48019 times.
✗ Branch 6 not taken.
 48019 const vec3& pi = Geom::halfedge_vertex_from(MH.mesh(), hole[i]);
236
2/4
✓ Branch 1 taken 48019 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 48019 times.
✗ Branch 6 not taken.
 48019 const vec3& pj = Geom::halfedge_vertex_from(MH.mesh(), hole[j]);
237
1/2
✓ Branch 2 taken 48019 times.
✗ Branch 3 not taken.
 48019 double dxij = length(pj - pi);
238
239 48019 dsij = std::max(dsij, 1e-6);
240 48019 dxij = std::max(dxij, 1e-6);
241 48019 double rij = dxij / dsij;
242
243
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 48019 times.
 48019 if(use_normals) {
244 ✗ const vec3& Pi = Geom::halfedge_vertex_from(
245 MH.mesh(), hole[i]
246 );
247 ✗ const vec3& Pj = Geom::halfedge_vertex_from(
248 MH.mesh(), hole[j]
249 );
250 ✗ vec3 Dij = normalize(Pj - Pi);
251
252 // between -1 (worse) and 1 (best)
253 double angle_factor =
254 ✗ 0.5 * (dot(Dij, N[i]) - dot(Dij, N[j]));
255
256 // between 0 (best) and 1 (worse)
257 ✗ angle_factor = 0.5 * (1.0 - angle_factor);
258
259 ✗ rij *= angle_factor;
260 }
261
262
2/2
✓ Branch 0 taken 1869 times.
✓ Branch 1 taken 46150 times.
 48019 if(rij < best_rij) {
263 1869 best_rij = rij;
264 1869 v1 = signed_index_t(i);
265 1869 v2 = signed_index_t(j);
266 }
267 }
268 }
269
270
2/4
✓ Branch 0 taken 314 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 314 times.
 314 if(v1 == -1 || v2 == -1) {
271 ✗ return false;
272 }
273
274 // Now I do not think this can happen
275 // (to be checked)
276
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 314 times.
 314 if(v2 < v1) {
277 ✗ std::swap(v1, v2);
278 }
279
280 // Step 3: copy the two "sub-holes"
281 314 hole1.clear();
282 314 hole2.clear();
283
2/2
✓ Branch 1 taken 2984 times.
✓ Branch 2 taken 314 times.
 3298 for(signed_index_t i = 0; i < signed_index_t(hole.size()); i++) {
284
4/4
✓ Branch 0 taken 2313 times.
✓ Branch 1 taken 671 times.
✓ Branch 2 taken 897 times.
✓ Branch 3 taken 1416 times.
 2984 if(i <= v1 || i >= v2) {
285
2/4
✓ Branch 1 taken 1568 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1568 times.
✗ Branch 5 not taken.
 1568 hole1.push_back(hole[i]);
286 }
287
4/4
✓ Branch 0 taken 2627 times.
✓ Branch 1 taken 357 times.
✓ Branch 2 taken 2044 times.
✓ Branch 3 taken 583 times.
 2984 if(i >= v1 && i <= v2) {
288
2/4
✓ Branch 1 taken 2044 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2044 times.
✗ Branch 5 not taken.
 2044 hole2.push_back(hole[i]);
289 }
290 }
291
292 314 return true;
293 314 }
294
295 /**
296 * \brief Triangulates a hole.
297 * \param[in] MH the mesh, wrapped with halfedge accessors
298 * \param[in] hole the hole, represented by a vector of halfedges
299 * \param[out] triangles the generated triangles
300 * \param[in] use_normals if set, then couples of vertices
301 * (v1,v2) that have their border normals that match
302 * the new segment [v1,v2], are connected in priority.
303 * This improves the result on 788_raw_hand.off (but
304 * may cause some triangles overlap)
305 * \param[in] clear if set, \p triangles is cleared
306 * \retval true on success
307 * \retval false otherwise (for instance, if no valid edge
308 * could be found to split the hole)
309 */
310 358 bool triangulate_hole_loop_splitting(
311 const MeshHalfedges& MH, const Hole& hole,
312 vector<trindex>& triangles, bool use_normals,
313 bool clear = true
314 ) {
315 358 bool ok = true;
316
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 350 times.
 358 if(clear) {
317 8 triangles.clear();
318 }
319
2/2
✓ Branch 1 taken 183 times.
✓ Branch 2 taken 175 times.
 358 if(hole.size() <= 3) {
320
1/2
✓ Branch 1 taken 183 times.
✗ Branch 2 not taken.
 183 if(hole.size() == 3) {
321 trindex T(
322 549 MH.mesh().facet_corners.vertex(hole[0].corner),
323 549 MH.mesh().facet_corners.vertex(hole[1].corner),
324 183 MH.mesh().facet_corners.vertex(hole[2].corner),
325 trindex::KEEP_ORDER
326
7/14
✓ Branch 1 taken 183 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 183 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 183 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 183 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 183 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 183 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 183 times.
✗ Branch 22 not taken.
 183 );
327
1/2
✓ Branch 1 taken 183 times.
✗ Branch 2 not taken.
 183 triangles.push_back(T);
328 }
329 } else {
330 175 Hole hole1, hole2;
331
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 175 ok = split_hole(MH, hole, hole1, hole2, use_normals);
332
3/6
✓ Branch 0 taken 175 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 175 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 175 times.
✗ Branch 6 not taken.
 175 ok = ok && triangulate_hole_loop_splitting(
333 MH, hole1, triangles, use_normals, false
334 );
335
3/6
✓ Branch 0 taken 175 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 175 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 175 times.
✗ Branch 6 not taken.
 175 ok = ok && triangulate_hole_loop_splitting(
336 MH, hole2, triangles, use_normals, false
337 );
338 175 }
339 358 return ok;
340 }
341
342 /************************************************************************/
343
344 /**
345 * \brief Computes the score obtained when generating a triangle
346 * in the ear cutting algorithm triangulate_hole_ear_cutting().
347 * \param[in] M the mesh
348 * \param[in] T1 first triangle, encoded as a vertex indices triplet
349 * \param[in] T2 second triangle, encoded as a vertex indices triplet
350 * \return the score obtained when generating an ear from \p T1 to \p T2
351 */
352 ✗ double ear_score(
353 const Mesh& M,
354 const trindex& T1,
355 const trindex& T2
356 ) {
357 ✗ geo_debug_assert(T1.indices[1] == T2.indices[0]);
358 ✗ const vec3& p10 = M.vertices.point(T1.indices[0]);
359 ✗ const vec3& p11 = M.vertices.point(T1.indices[1]);
360 ✗ const vec3& p12 = M.vertices.point(T1.indices[2]);
361 ✗ const vec3& p20 = M.vertices.point(T2.indices[0]);
362 ✗ const vec3& p21 = M.vertices.point(T2.indices[1]);
363 ✗ const vec3& p22 = M.vertices.point(T2.indices[2]);
364 ✗ vec3 n = normalize(
365 ✗ Geom::triangle_normal(p10, p11, p12) +
366 ✗ Geom::triangle_normal(p20, p21, p22)
367 );
368 ✗ vec3 a = normalize(p11 - p10);
369 ✗ vec3 b = normalize(p21 - p20);
370 ✗ return -::atan2(dot(n, cross(a, b)), dot(a, b));
371 }
372
373 /**
374 * \brief Triangulates a hole using the ear cutting algorithm.
375 * \param[in] MH the mesh, wrapped with halfedge accessors
376 * \param[in] hole_in the hole, represented by a vector of halfedges
377 * \param[out] triangles the generated triangles
378 * \param[out] clear if set, \p triangles is cleared
379 */
380 ✗ void triangulate_hole_ear_cutting(
381 const MeshHalfedges& MH, const Hole& hole_in,
382 vector<trindex>& triangles, bool clear = true
383 ) {
384 ✗ if(clear) {
385 ✗ triangles.clear();
386 }
387 ✗ if(hole_in.size() <= 3) {
388 ✗ if(hole_in.size() == 3) {
389 trindex T(
390 ✗ MH.mesh().facet_corners.vertex(hole_in[0].corner),
391 ✗ MH.mesh().facet_corners.vertex(hole_in[1].corner),
392 ✗ MH.mesh().facet_corners.vertex(hole_in[2].corner),
393 trindex::KEEP_ORDER
394 ✗ );
395 ✗ triangles.push_back(T);
396 }
397 } else {
398 ✗ const Mesh& M = MH.mesh();
399
400 // Step 1: convert hole into easier-to-manipulate representation.
401 ✗ vector<trindex> hole;
402 ✗ hole.reserve(hole_in.size());
403 ✗ for(index_t i = 0; i < hole_in.size(); i++) {
404 ✗ const MeshHalfedges::Halfedge& H = hole_in[i];
405 ✗ geo_debug_assert(H.facet != MeshHalfedges::Halfedge::NO_FACET);
406 ✗ index_t c = H.corner;
407 ✗ index_t v1 = M.facet_corners.vertex(c);
408 ✗ c = M.facets.next_corner_around_facet(H.facet, c);
409 ✗ index_t v2 = M.facet_corners.vertex(c);
410 ✗ c = M.facets.next_corner_around_facet(H.facet, c);
411 ✗ index_t v3 = M.facet_corners.vertex(c);
412 ✗ hole.push_back(trindex(v1, v2, v3, trindex::KEEP_ORDER));
413 }
414
415 // Step 2: ear cutting
416 ✗ while(hole.size() > 3) {
417 ✗ signed_index_t best_i1 = -1;
418 ✗ double best_score = Numeric::min_float64();
419 // TODO: take existing edges into account.
420 ✗ for(index_t i1 = 0; i1 < hole.size(); i1++) {
421 ✗ index_t i2 = i1 + 1;
422 ✗ if(i2 == hole.size()) {
423 ✗ i2 = 0;
424 }
425 ✗ double score = ear_score(M, hole[i1], hole[i2]);
426 ✗ if(score > best_score) {
427 ✗ best_i1 = signed_index_t(i1);
428 ✗ best_score = score;
429 }
430 }
431 ✗ geo_assert(best_i1 != -1);
432 ✗ index_t best_i2 = index_t(best_i1) + 1;
433 ✗ if(best_i2 == hole.size()) {
434 ✗ best_i2 = 0;
435 }
436 ✗ const trindex& T1 = hole[best_i1];
437 ✗ const trindex& T2 = hole[best_i2];
438 ✗ geo_debug_assert(T1.indices[1] == T2.indices[0]);
439 trindex T(
440 ✗ T1.indices[0], T2.indices[1], T1.indices[1],
441 trindex::KEEP_ORDER
442 ✗ );
443 ✗ hole[best_i1] = T;
444 ✗ hole.erase(hole.begin() + std::ptrdiff_t(best_i2));
445 ✗ triangles.push_back(T);
446 }
447
448 // Step 3: last triangle
449 ✗ geo_assert(hole.size() == 3);
450 trindex T(
451 ✗ hole[0].indices[0],
452 ✗ hole[1].indices[0],
453 ✗ hole[2].indices[0],
454 trindex::KEEP_ORDER
455 ✗ );
456 ✗ triangles.push_back(T);
457 ✗ }
458 ✗ }
459
460 /************************************************************************/
461
462 /**
463 * \brief Computes the area of a hole, i.e. the area of the generated
464 * triangles that will fill the hole.
465 * \param[in] M the mesh
466 * \param[in] triangles the triangles that will fill the hole
467 * \return the summed ares of the triangles in \p triangles
468 */
469 8 double hole_area(
470 const Mesh& M, const vector<trindex>& triangles
471 ) {
472 8 double result = 0.0;
473
2/2
✓ Branch 1 taken 183 times.
✓ Branch 2 taken 8 times.
 191 for(index_t t = 0; t < triangles.size(); t++) {
474 183 index_t i = triangles[t].indices[0];
475 183 index_t j = triangles[t].indices[1];
476 183 index_t k = triangles[t].indices[2];
477 183 const vec3& p1 = M.vertices.point(i);
478 183 const vec3& p2 = M.vertices.point(j);
479 183 const vec3& p3 = M.vertices.point(k);
480 183 result += Geom::triangle_area(p1, p2, p3);
481 }
482 8 return result;
483 }
484
485 /**
486 * \brief Strategy used to fill the holes.
487 */
488 enum HoleFilling {
489 LOOP_SPLIT, /**< Splits loops by generating small segments */
490 NLOOP_SPLIT, /**< Takes normals into account */
491 EAR_CUT /**< Uses the "ear cutting" strategy */
492 };
493
494 /************************************************************************/
495
496 /* // commented-out for now, see issue #72
497 * \brief Removes all the facets of a mesh that are
498 * on a bridge.
499 * \details A facet is said to be on a bridge if it is
500 * incident to a border and if when turning around the
501 * border it is encountered more than once.
502 *
503 *
504 *
505 void remove_bridges(Mesh& M) {
506 MeshHalfedges MH(M);
507 vector<bool> corner_is_visited(M.facet_corners.nb(),false);
508 vector<index_t> f_status(M.facets.nb(),0);
509 index_t f_stamp=1;
510 const index_t BRIDGE = NO_INDEX;
511
512 for(index_t f: M.facets) {
513 for(index_t c: M.facets.corners(f)) {
514 if(
515 M.facet_corners.adjacent_facet(c) == NO_FACET &&
516 !corner_is_visited[c]
517 ) {
518 MeshHalfedges::Halfedge first(f, c);
519 MeshHalfedges::Halfedge H(f, c);
520 do {
521 corner_is_visited[H.corner] = true;
522 MH.move_to_next_around_facet(H);
523 while(MH.move_to_next_around_vertex(H)) {
524 if(f_status[H.facet] == f_stamp) {
525 f_status[H.facet] = BRIDGE;
526 } else if(
527 f_status[H.facet] != BRIDGE &&
528 f_status[H.facet] != f_stamp) {
529 f_status[H.facet] = f_stamp;
530 }
531 }
532 } while(H != first);
533 ++f_stamp;
534 }
535 }
536 }
537 index_t nb_bridges = 0;
538 for(index_t f: M.facets) {
539 if(f_status[f] == BRIDGE) {
540 ++nb_bridges;
541 } else {
542 f_status[f] = 0;
543 }
544 }
545 if(nb_bridges != 0) {
546 M.facets.delete_elements(f_status);
547 Logger::out("Bridges")
548 << "Removed " << nb_bridges << " bridge(s)"
549 << std::endl;
550 }
551 }
552 */
553 }
554
555 /****************************************************************************/
556
557 namespace GEO {
558
559 7 void fill_holes(
560 Mesh& M, double max_area, index_t max_edges, bool repair
561 ) {
562
563
2/4
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 7 times.
 7 if(max_area == 0.0 || max_edges == 0) {
564 4 return;
565 }
566
567 // remove_bridges(M); // commented out for now, see issue #72
568
569
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 MeshHalfedges MH(M);
570
571 7 vector<Hole> holes;
572
573 7 index_t nb_filled_holes = 0;
574 7 index_t nb_skipped_by_edges = 0;
575 7 index_t nb_skipped_by_area = 0;
576 7 index_t nb_could_not_fill = 0;
577
578 // Step 1: detect holes
579 {
580
1/2
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
 7 vector<bool> corner_is_visited(M.facet_corners.nb(), false);
581
2/2
✓ Branch 5 taken 132340 times.
✓ Branch 6 taken 7 times.
 132347 for(index_t f: M.facets) {
582
3/4
✓ Branch 1 taken 132340 times.
✗ Branch 2 not taken.
✓ Branch 8 taken 397020 times.
✓ Branch 9 taken 132340 times.
 529360 for(index_t c: M.facets.corners(f)) {
583 397020 if(
584
3/4
✓ Branch 1 taken 397020 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2256 times.
✓ Branch 4 taken 394764 times.
 399276 M.facet_corners.adjacent_facet(c) == NO_FACET &&
585
4/4
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 2246 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 397010 times.
 399276 !corner_is_visited[c]
586 ) {
587
1/2
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
 10 holes.push_back(Hole());
588 10 MeshHalfedges::Halfedge first(f, c);
589 10 MeshHalfedges::Halfedge H(f, c);
590 do {
591
2/4
✓ Branch 2 taken 2256 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2256 times.
✗ Branch 6 not taken.
 2256 holes.rbegin()->push_back(H);
592 2256 corner_is_visited[H.corner] = true;
593
1/2
✓ Branch 1 taken 2256 times.
✗ Branch 2 not taken.
 2256 MH.move_to_next_around_border(H);
594
2/2
✓ Branch 1 taken 2246 times.
✓ Branch 2 taken 10 times.
 2256 } while(H != first);
595
596
3/4
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 8 times.
 10 if(holes.rbegin()->size() > max_edges) {
597 2 ++nb_skipped_by_edges;
598 2 holes.pop_back();
599 }
600 }
601 }
602 }
603 7 }
604
605
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 3 times.
 7 if(holes.size() == 0) {
606 4 return;
607 }
608
609
4/8
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 3 times.
✗ Branch 12 not taken.
 6 Logger::out("FillHoles") << "Found " << holes.size()
610
2/4
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
 3 << " holes" << std::endl;
611
612 3 HoleFilling algo = LOOP_SPLIT;
613
2/4
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
 3 std::string algo_name = CmdLine::get_arg("algo:hole_filling");
614
2/4
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
 3 if(algo_name == "loop_split") {
615 3 algo = LOOP_SPLIT;
616 ✗ } else if(algo_name == "Nloop_split") {
617 ✗ algo = NLOOP_SPLIT;
618 ✗ } else if(algo_name == "ear_cut") {
619 ✗ algo = EAR_CUT;
620 } else {
621 ✗ Logger::warn("FillHoles")
622 ✗ << algo_name << ": no such hole filling method"
623 ✗ << std::endl;
624 ✗ Logger::warn("FillHoles")
625 ✗ << "falling back to \'loop_split\'"
626 ✗ << std::endl;
627 }
628
629
630
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 3 times.
 11 for(index_t i = 0; i < holes.size(); i++) {
631 8 vector<trindex> triangles;
632 8 bool ok = true;
633
1/4
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 8 switch(algo) {
634 8 case LOOP_SPLIT:
635 16 ok = triangulate_hole_loop_splitting(
636
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
 8 MH, holes[i], triangles, false
637 );
638 8 break;
639 ✗ case NLOOP_SPLIT:
640 ✗ ok = triangulate_hole_loop_splitting(
641 ✗ MH, holes[i], triangles, true
642 );
643 ✗ break;
644 ✗ case EAR_CUT:
645 ✗ triangulate_hole_ear_cutting(MH, holes[i], triangles);
646 ✗ break;
647 }
648
649
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
 8 if(ok) {
650
3/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7 times.
✓ Branch 4 taken 1 times.
 8 if(hole_area(M, triangles) < max_area) {
651
2/2
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 7 times.
 14 for(index_t j = 0; j < triangles.size(); j++) {
652 14 M.facets.create_triangle(
653
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 triangles[j].indices[2],
654
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 triangles[j].indices[1],
655
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
 7 triangles[j].indices[0]
656 );
657 }
658 7 ++nb_filled_holes;
659 } else {
660 1 ++nb_skipped_by_area;
661 }
662 } else {
663 ✗ ++nb_could_not_fill;
664 }
665
666 8 }
667
668
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
 3 if(nb_skipped_by_area != 0) {
669
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 3 Logger::out("FillHoles")
670
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 << "Skipped " << nb_skipped_by_area
671
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 << " holes (area too large)" << std::endl;
672 }
673
674
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 if(nb_skipped_by_edges != 0) {
675
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 6 Logger::out("FillHoles")
676
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 << "Skipped " << nb_skipped_by_edges
677
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 << " holes (too many edges)" << std::endl;
678 }
679
680
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 if(nb_could_not_fill != 0) {
681 ✗ Logger::out("FillHoles")
682 ✗ << "Skipped " << nb_could_not_fill
683 ✗ << " holes (could not fill)" << std::endl;
684 }
685
686
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 3 if(nb_filled_holes != 0 && repair) {
687 // Needed because we may generate zero-length edges
688 // and zero-area facets that need to be eliminated.
689 // Note: this also reconstructs the connections between the facets.
690 2 MeshRepairMode mode = MESH_REPAIR_DEFAULT;
691
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 mesh_repair(M, mode);
692 }
693
4/4
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 4 times.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 4 times.
 11 }
694
695 /**
696 * \brief Tessellates a hole.
697 * \param[in] MH the MeshHalfHeddes the hole belongs to.
698 * \param[in] H the hole.
699 * \param[in] max_nb_vertices maximum number of vertices in
700 * the new facets to create.
701 * \param[in] copy_facet_attrib an optional facet index. If
702 * specified, all the attributes of this facet will be copied
703 * to the created facets.
704 */
705 417 static void tessellate_hole(
706 MeshHalfedges& MH, Hole& H, index_t max_nb_vertices,
707 index_t copy_facet_attrib = NO_INDEX
708 ) {
709 417 Mesh& M = MH.mesh();
710
2/2
✓ Branch 1 taken 278 times.
✓ Branch 2 taken 139 times.
 417 if(H.size() <= max_nb_vertices) {
711 278 index_t f = M.facets.create_polygon(H.size());
712
2/2
✓ Branch 1 taken 834 times.
✓ Branch 2 taken 278 times.
 1112 FOR(i,H.size()) {
713 834 index_t v = M.facet_corners.vertex(H[i].corner);
714 834 M.facets.set_vertex(f,i,v);
715 }
716
1/2
✓ Branch 0 taken 278 times.
✗ Branch 1 not taken.
 278 if(copy_facet_attrib != NO_INDEX) {
717 278 M.facets.attributes().copy_item(f, copy_facet_attrib);
718 }
719 } else {
720 139 Hole H1,H2;
721
1/2
✓ Branch 1 taken 139 times.
✗ Branch 2 not taken.
 139 split_hole(MH,H,H1,H2,false);
722
1/2
✓ Branch 1 taken 139 times.
✗ Branch 2 not taken.
 139 tessellate_hole(MH,H1,max_nb_vertices,copy_facet_attrib);
723
1/2
✓ Branch 1 taken 139 times.
✗ Branch 2 not taken.
 139 tessellate_hole(MH,H2,max_nb_vertices,copy_facet_attrib);
724 139 }
725 417 }
726
727 6 void tessellate_facets(
728 Mesh& M, index_t max_nb_vertices
729 ) {
730
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 MeshHalfedges MH(M);
731
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 6 vector<index_t> delete_f(M.facets.nb(),0);
732
2/2
✓ Branch 5 taken 143 times.
✓ Branch 6 taken 6 times.
 149 for(index_t f: M.facets) {
733
3/4
✓ Branch 1 taken 143 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 139 times.
✓ Branch 4 taken 4 times.
 143 if(M.facets.nb_vertices(f) > max_nb_vertices) {
734
1/2
✓ Branch 1 taken 139 times.
✗ Branch 2 not taken.
 139 delete_f[f] = 1;
735 139 Hole h;
736
3/4
✓ Branch 1 taken 139 times.
✗ Branch 2 not taken.
✓ Branch 7 taken 556 times.
✓ Branch 8 taken 139 times.
 695 for(index_t c: M.facets.corners(f)) {
737
1/2
✓ Branch 2 taken 556 times.
✗ Branch 3 not taken.
 556 h.push_back(MeshHalfedges::Halfedge(f,c));
738 }
739
1/2
✓ Branch 1 taken 139 times.
✗ Branch 2 not taken.
 139 tessellate_hole(MH, h, max_nb_vertices, f);
740 139 }
741 }
742
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 6 delete_f.resize(M.facets.nb());
743
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 M.facets.delete_elements(delete_f);
744
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 M.facets.connect();
745
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 6 if(max_nb_vertices == 3) {
746
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 M.facets.is_simplicial();
747 }
748 6 }
749
750 }
751