GCC Code Coverage Report


Directory: ./
File: lib/geogram/mesh/mesh_degree3_vertices.cpp
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 61 98 62.2%
Functions: 3 4 75.0%
Branches: 43 114 37.7%

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_degree3_vertices.h>
41 #include <geogram/mesh/mesh_geometry.h>
42 #include <geogram/mesh/index.h>
43 #include <geogram/basic/geometry_nd.h>
44 #include <geogram/basic/logger.h>
45
46 namespace {
47
48 using namespace GEO;
49
50 // TODO: include in Doxygen documentation.
51 /*
52 * There are two different numerotation for triangles:
53 *
54 * o GEO::Mesh numerotation:
55 * - Used everywhere in Vorpaline
56 * - Edge ei is right after vertex vi when turning
57 * around the facet.
58 *
59 * o Standard triangle numerotation:
60 * - Used in this file, by t_xxx() functions
61 * and Degree3Vertex class.
62 * - Edge ei is opposite to vertex vi.
63 *
64 * GEO::Mesh numerotation:
65 * v2
66 * e2 / \e1
67 * / \
68 * v0----v1
69 * e0
70 *
71 * Standard triangle numerotation:
72 * v2
73 * e1 / \e0
74 * / \
75 * v0----v1
76 * e2
77 */
78
79 /**
80 * \brief Gets the vertex of a triangle by its local index, using
81 * standard triangle numerotation (edge i is opposite to vertex i).
82 * \param[in] M the mesh
83 * \param[in] t the index of the triangle
84 * \param[in] li the local index (0,1 or 2) of the vertex in \p t
85 * \return the global index of the vertex in \p M
86 * \note Everywhere else in Vorpaline, triangle edge numerotation is
87 * different from standard triangle numerotation used here.
88 */
89 inline index_t t_vertex(
90 const Mesh& M, index_t t, index_t li
91 ) {
92 geo_debug_assert(M.facets.nb_vertices(t) == 3);
93 geo_debug_assert(li < 3);
94 1 return M.facet_corners.vertex(M.facets.corners_begin(t) + li);
95 }
96
97 /**
98 * \brief Gets the index of an adjacent triangle, using
99 * \param[in] M the mesh
100 * \param[in] t the index of the triangle
101 * \param[in] le the local index (0,1 or 2) of the edge in \p t, using
102 * standard triangle numerotation (edge i is opposite to vertex i).
103 * \return the global index of the triangle adjacent to \p t accros
104 * edge \p le in \p M
105 * \note Everywhere else in Vorpaline, triangle edge numerotation is
106 * different from standard triangle numerotation used here.
107 */
108
109 inline index_t t_adjacent(const Mesh& M, index_t t, index_t le) {
110 geo_debug_assert(M.facets.nb_vertices(t) == 3);
111 geo_debug_assert(le < 3);
112 2 le = (le + 1) % 3; // convert corner numerotation
113 // to standard triangle numerotation
114 2 return M.facet_corners.adjacent_facet(M.facets.corners_begin(t) + le);
115 }
116
117 /**
118 * \brief Gets the local vertex index in a triangle by its global index,
119 * using standard triangle numerotation (edge i is opposite to vertex i).
120 * \param[in] M the mesh
121 * \param[in] t the index of the triangle
122 * \param[in] v the global index of the vertex in \p M
123 * \return the local index (0,1 or 2) of the vertex \p v in \p t
124 * \note Everywhere else in Vorpaline, triangle edge numerotation is
125 * different from standard triangle numerotation used here.
126 */
127 3 inline index_t t_index(const Mesh& M, index_t t, index_t v) {
128 geo_debug_assert(M.facets.nb_vertices(t) == 3);
129
1/2
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
7 for(index_t li = 0; li < 3; li++) {
130
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 4 times.
7 if(M.facet_corners.vertex(M.facets.corners_begin(t) + li) == v) {
131 3 return li;
132 }
133 }
134 geo_assert_not_reached;
135 }
136
137 /**
138 * \brief Sets a triangle vertices and neighbors
139 * using standard triangle numerotation (edge i is opposite to vertex i).
140 * \param[in] M the mesh
141 * \param[in] t the index of the triangle
142 * \param[in] v1 the global index of the first vertex
143 * \param[in] v2 the global index of the second vertex
144 * \param[in] v3 the global index of the third vertex
145 * \param[in] adj1 the index of the adjacent triangle opposite to \p v1
146 * \param[in] adj2 the index of the adjacent triangle opposite to \p v2
147 * \param[in] adj3 the index of the adjacent triangle opposite to \p v3
148 * \note Everywhere else in Vorpaline, triangle edge numerotation is
149 * different from standard triangle numerotation used here.
150 */
151 inline void t_set(
152 Mesh& M,
153 index_t t,
154 index_t v1, index_t v2, index_t v3,
155 index_t adj1, index_t adj2, index_t adj3
156 ) {
157 geo_debug_assert(M.facets.nb_vertices(t) == 3);
158 geo_debug_assert(adj1 != NO_FACET);
159 geo_debug_assert(adj2 != NO_FACET);
160 geo_debug_assert(adj3 != NO_FACET);
161 index_t c1 = M.facets.corners_begin(t);
162 index_t c2 = c1 + 1;
163 index_t c3 = c1 + 2;
164 M.facet_corners.set_vertex(c1, v1);
165 M.facet_corners.set_vertex(c2, v2);
166 M.facet_corners.set_vertex(c3, v3);
167 M.facet_corners.set_adjacent_facet(c1, adj3);
168 M.facet_corners.set_adjacent_facet(c2, adj1);
169 M.facet_corners.set_adjacent_facet(c3, adj2);
170 }
171
172 // TODO: ascii-art needed here !!
173
174 /**
175 * \brief Stores the three facets incident to a degree 3 vertex
176 * and the three adjacent facets.
177 * \details Implements operator< (according to distance between degree 3
178 * vertex and neighbor vertices supporting plane).
179 * \note Uses standard triangle numerotation (edge i is opposite to
180 * vertex i). Everywhere else in Vorpaline, triangle edge numerotation is
181 * different from standard triangle numerotation used here.
182 */
183 struct Degree3Vertex {
184
185 /**
186 * \brief Constructs a new Degree3Vertex
187 * \param[in] M the mesh
188 * \param[in] v_in the index of the degree 3 vertex in \p M
189 * \param[in] t_in the index of a triangle incident to \p v_in
190 * \pre there are exactly three triangles incident to \p v_in
191 */
192 1 Degree3Vertex(const Mesh& M, index_t v_in, index_t t_in) {
193 1 v[3] = v_in;
194 1 t[0] = t_in;
195 {
196 1 index_t i = t_index(M, t[0], v[3]);
197 1 index_t j = (i + 1) % 3;
198 1 index_t k = (j + 1) % 3;
199
200
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 t[1] = index_t(t_adjacent(M, t[0], j));
201 geo_debug_assert(t[1] != NO_INDEX);
202 1 t[2] = index_t(t_adjacent(M, t[0], k));
203 geo_debug_assert(t[2] != NO_INDEX);
204 1 v[1] = t_vertex(M, t[0], j);
205 1 v[2] = t_vertex(M, t[0], k);
206 1 adj[0] = t_adjacent(M, t[0], i);
207 }
208
209 {
210 1 index_t i = t_index(M, t[1], v[3]);
211 1 index_t j = (i + 1) % 3;
212 1 index_t k = (j + 1) % 3;
213
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
2 v[0] = t_vertex(M, t[1], k);
214 geo_debug_assert(t_vertex(M, t[1], j) == v[2]);
215 geo_debug_assert(t_adjacent(M, t[1], j) == t[2]);
216 geo_debug_assert(t_adjacent(M, t[1], k) == t[0]);
217 1 adj[1] = t_adjacent(M, t[1], i);
218 }
219
220 {
221 1 index_t i = t_index(M, t[2], v[3]);
222 #ifdef GEO_DEBUG
223 index_t j = (i + 1) % 3;
224 index_t k = (j + 1) % 3;
225 geo_debug_assert(t_vertex(M, t[2], j) == v[0]);
226 geo_debug_assert(t_vertex(M, t[2], k) == v[1]);
227 geo_debug_assert(t_adjacent(M, t[2], j) == t[0]);
228 geo_debug_assert(t_adjacent(M, t[2], k) == t[1]);
229 #endif
230
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 adj[2] = t_adjacent(M, t[2], i);
231 }
232
233 1 const vec3& p0 = M.vertices.point(v[0]);
234 1 const vec3& p1 = M.vertices.point(v[1]);
235 1 const vec3& p2 = M.vertices.point(v[2]);
236 1 const vec3& p3 = M.vertices.point(v[3]);
237
238 1 dist = ::sqrt(
239 Geom::point_triangle_squared_distance(p3, p0, p1, p2)
240 );
241 1 }
242
243 /**
244 * \brief Compares two Degree3Vertex by the distance to the
245 * supporting planes of the neighbors.
246 * \param[in] rhs the comparand
247 * \return true if this Degree3Vertex can be removed before
248 * \p rhs, false otherwise
249 */
250 bool operator< (const Degree3Vertex& rhs) const {
251 return dist < rhs.dist;
252 }
253
254 double dist;
255 index_t v[4];
256 index_t t[3];
257 index_t adj[3];
258 };
259 }
260
261 /****************************************************************************/
262
263 namespace GEO {
264
265 8 index_t remove_degree3_vertices(Mesh& M, double max_dist) {
266
267
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
8 if(max_dist == 0.0) {
268 return 0;
269 }
270
271 // Step 1: detect degree3 vertices
272
273 8 vector<signed_index_t> vertex_degree(M.vertices.nb(), 0);
274 // or -1 if v is on border or if v has an incident facet that
275 // is not a triangle.
276
277
2/2
✓ Branch 0 taken 61104 times.
✓ Branch 1 taken 8 times.
61112 for(index_t f: M.facets) {
278 bool f_is_triangle = (M.facets.nb_vertices(f) == 3);
279
2/2
✓ Branch 0 taken 183312 times.
✓ Branch 1 taken 61104 times.
244416 for(index_t c: M.facets.corners(f)) {
280 index_t v = M.facet_corners.vertex(c);
281 if(
282
3/4
✓ Branch 0 taken 183312 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 194 times.
✓ Branch 3 taken 183118 times.
183312 !f_is_triangle ||
283 (M.facet_corners.adjacent_facet(c) == NO_FACET)
284 ) {
285 194 vertex_degree[v] = -1;
286 } else {
287
2/2
✓ Branch 0 taken 182981 times.
✓ Branch 1 taken 137 times.
183118 if(vertex_degree[v] != -1) {
288 182981 vertex_degree[v]++;
289 }
290 }
291 }
292 }
293
294 // Step 2: count degree3 vertices
295
296 index_t nb_degree3_vertices = 0;
297
2/2
✓ Branch 0 taken 30654 times.
✓ Branch 1 taken 8 times.
30662 for(index_t v: M.vertices) {
298
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 30653 times.
30654 if(vertex_degree[v] == 3) {
299 1 nb_degree3_vertices++;
300 }
301 }
302
303
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 1 times.
8 if(nb_degree3_vertices == 0) {
304
2/6
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
7 Logger::out("Degree3")
305 << "Does not have any degree 3 vertex (good)" << std::endl;
306 7 return 0;
307 }
308
309 // Step 3: v2f[v] is one of the facets adjacent to v
310
311 vector<index_t> v2f(M.vertices.nb());
312
2/2
✓ Branch 0 taken 9996 times.
✓ Branch 1 taken 1 times.
9997 for(index_t f: M.facets) {
313
2/2
✓ Branch 0 taken 29988 times.
✓ Branch 1 taken 9996 times.
39984 for(index_t c: M.facets.corners(f)) {
314 index_t v = M.facet_corners.vertex(c);
315 29988 v2f[v] = f;
316 }
317 }
318
319 // Step 4: compute and sort degree3 vertices configurations
320
321 vector<Degree3Vertex> degree3vertices;
322
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 degree3vertices.reserve(nb_degree3_vertices);
323
2/2
✓ Branch 0 taken 5000 times.
✓ Branch 1 taken 1 times.
5001 for(index_t v: M.vertices) {
324
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4999 times.
5000 if(vertex_degree[v] == 3) {
325
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 Degree3Vertex V(M, v, v2f[v]);
326
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(V.dist < max_dist) {
327 degree3vertices.push_back(V);
328 }
329 }
330 }
331
332
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
2 Logger::out("Degree3")
333 << "Removing " << degree3vertices.size()
334 << "/" << nb_degree3_vertices
335 << " degree 3 vertices (within max_deg3_dist)"
336 << std::endl;
337
338 1 std::sort(degree3vertices.begin(), degree3vertices.end());
339
340 // Step 5: remove degree3 vertices in order
341
342 enum FacetStatus {
343 F_UNTOUCHED = 0,
344 F_TO_REMOVE = 1,
345 F_MODIFIED = 2
346 };
347
348
1/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
1 vector<index_t> facet_status(M.facets.nb(), F_UNTOUCHED);
349 index_t nb_removed = 0;
350
351
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 for(index_t i = 0; i < degree3vertices.size(); i++) {
352 const Degree3Vertex& V = degree3vertices[i];
353 index_t t1 = V.t[0];
354 index_t t2 = V.t[1];
355 index_t t3 = V.t[2];
356
357 // Skip vertex if one of its incident facet was already modified
358 // when removing a neighboring degree 3 vertex.
359 if(
360 facet_status[t1] != F_UNTOUCHED ||
361 facet_status[t2] != F_UNTOUCHED ||
362 facet_status[t3] != F_UNTOUCHED
363 ) {
364 continue;
365 }
366
367 // Skip vertex if one of its neighboring facet
368 // was already modified...
369 index_t adj1 = V.adj[0];
370 index_t adj2 = V.adj[1];
371 index_t adj3 = V.adj[2];
372 if(adj1 != NO_FACET && facet_status[adj1] != F_UNTOUCHED) {
373 continue;
374 }
375 if(adj2 != NO_FACET && facet_status[adj2] != F_UNTOUCHED) {
376 continue;
377 }
378 if(adj3 != NO_FACET && facet_status[adj3] != F_UNTOUCHED) {
379 continue;
380 }
381
382 facet_status[t2] = F_TO_REMOVE; // t2 will be deleted.
383 facet_status[t3] = F_TO_REMOVE; // t3 will be deleted.
384 facet_status[t1] = F_MODIFIED; // t1 is recycled and
385 // used by the new facet.
386 t_set(
387 M, t1, V.v[0], V.v[1], V.v[2],
388 index_t(V.adj[0]), index_t(V.adj[1]), index_t(V.adj[2])
389 );
390
391 // connect the neighbors of t2 and t3 to t1
392 for(index_t j = 1; j <= 2; j++) {
393 if(V.adj[j] != NO_FACET) {
394 index_t f = index_t(V.adj[j]);
395 for(index_t c: M.facets.corners(f)) {
396 if(M.facet_corners.adjacent_facet(c) == V.t[j]) {
397 M.facet_corners.set_adjacent_facet(c, V.t[0]);
398 }
399 }
400 }
401 }
402 nb_removed++;
403 }
404
405 // Step 6: delete dangling triangles
406
407 // Do not delete the facets that were recycled.
408
2/2
✓ Branch 0 taken 9996 times.
✓ Branch 1 taken 1 times.
19993 for(index_t i = 0; i < facet_status.size(); i++) {
409
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9996 times.
9996 if(facet_status[i] == F_MODIFIED) {
410 facet_status[i] = 0;
411 }
412 }
413
414 // Note: remove_facets() calls remove_isolated_vertices()
415
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 M.facets.delete_elements(facet_status);
416
417 return nb_removed;
418 }
419 }
420