GCC Code Coverage Report

Directory: ./
File: lib/geogram/points/co3ne.cpp
Date: 2026-09-07 02:37:58
Exec Total Coverage
Lines: 814 1123 72.5%
Functions: 75 90 83.3%
Branches: 712 1799 39.6%
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/points/co3ne.h>
41 #include <geogram/points/nn_search.h>
42 #include <geogram/points/principal_axes.h>
43 #include <geogram/mesh/mesh.h>
44 #include <geogram/mesh/mesh_io.h>
45 #include <geogram/mesh/index.h>
46 #include <geogram/mesh/mesh_repair.h>
47 #include <geogram/mesh/mesh_topology.h>
48 #include <geogram/mesh/mesh_reorder.h>
49 #include <geogram/basic/geometry.h>
50 #include <geogram/basic/process.h>
51 #include <geogram/basic/assert.h>
52 #include <geogram/basic/progress.h>
53 #include <geogram/basic/command_line.h>
54 #include <geogram/basic/algorithm.h>
55 #include <geogram/basic/stopwatch.h>
56 #include <stack>
57 #include <queue>
58
59 namespace {
60 using namespace GEO;
61
62 /**
63 * \brief number of elements in sine/cosine table.
64 */
65 static constexpr index_t sincos_nb = 10;
66
67 /**
68 * \brief sine/cosine table.
69 * \details We keep a small table of sines and cosines for
70 * speeding up things a little bit.
71 * Table entries are as follows:
72 * - sincos_table[i][0] = sin(2*M_PI*i/(sincos_nb-1))
73 * - sincos_table[i][1] = cos(2*M_PI*i/(sincos_nb-1))
74 */
75 static double sincos_table[10][2] = {
76 {0,1},
77 {0.642788,0.766044},
78 {0.984808,0.173648},
79 {0.866025,-0.5},
80 {0.34202,-0.939693},
81 {-0.34202,-0.939693},
82 {-0.866025,-0.5},
83 {-0.984808,0.173648},
84 {-0.642788,0.766044},
85 {-2.44929e-16,1}
86 };
87
88 /**
89 * \brief Used by the algorithm that reorients normals.
90 */
91 struct OrientNormal {
92 /**
93 * \brief OrientNormal constructor.
94 * \param[in] v_in the index of a point
95 * \param[in] dot_in the dot product between the (unit)
96 * normal vector at \p v_in and the normal vector at
97 * the point that initiated propagation to \p v_in.
98 */
99 ✗ OrientNormal(
100 index_t v_in, double dot_in
101 ✗ ) : v(v_in), dot(dot_in) {
102 ✗ }
103
104 /**
105 * \brief Compares two OrientNormal objects
106 * \retval true if \p rhs should be processed before this
107 * OrientObject.
108 * \retval false otherwise.
109 */
110 ✗ bool operator<(const OrientNormal& rhs) const {
111 ✗ return (::fabs(dot) < ::fabs(rhs.dot));
112 }
113 index_t v;
114 double dot;
115 };
116
117
118 /************************************************************/
119
120 /**
121 * \brief Extracts a manifold surface from the set of
122 * triangles reconstructed by Co3Ne.
123 */
124 class Co3NeManifoldExtraction {
125 public:
126 static constexpr index_t NO_CORNER = NO_INDEX;
127 static constexpr index_t NO_FACET = NO_INDEX;
128 static constexpr index_t NO_CNX = NO_INDEX;
129
130 /**
131 * \brief Initializes a new Co3NeManifoldExtraction with
132 * a list of triangles.
133 * \param[in,out] target the target mesh. It needs to be already
134 * initialized with the vertices.
135 * \param[in,out] good_triangles the good triangles reconstructed
136 * by Co3Ne. They are 'stealed' by the mesh (on exit, good_triangles
137 * is empty). If some non-manifold edges are detected, then all
138 * the triangles incident to any manifold edge are ignored.
139 */
140 2 Co3NeManifoldExtraction(
141 Mesh& target,
142 vector<index_t>& good_triangles
143 2 ) : M_(target) {
144
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 strict_ = CmdLine::get_arg_bool("co3ne:strict");
145
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(strict_) {
146 ✗ vector<index_t> first_triangle;
147 ✗ for(index_t i=0; i<3; ++i) {
148 ✗ first_triangle.push_back(*good_triangles.rbegin());
149 ✗ good_triangles.pop_back();
150 }
151 ✗ M_.facets.assign_triangle_mesh(first_triangle, true);
152 ✗ init_and_remove_non_manifold_edges();
153 ✗ init_connected_components();
154 ✗ add_triangles(good_triangles);
155 ✗ } else {
156
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 M_.facets.assign_triangle_mesh(good_triangles, true);
157
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 init_and_remove_non_manifold_edges();
158
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 init_connected_components();
159 }
160 2 }
161
162 /**
163 * \brief Tentatively adds triangle from the specified list.
164 * \details Some geometric and topological properties are
165 * verified by connect_and_validate_triangle() before accepting
166 * the triangle.
167 * \see connect_and_validate_triangle()
168 */
169 2 void add_triangles(const vector<index_t>& not_so_good_triangles) {
170
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 bool pretty = CmdLine::get_arg_bool("log:pretty");
171
172 2 index_t nb_triangles = not_so_good_triangles.size()/3;
173
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 Logger::out("Co3ne") << "Tentatively add "
174
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.
 2 << nb_triangles << " triangles" << std::endl;
175
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 vector<bool> t_is_classified(nb_triangles,false);
176 2 bool changed = true;
177
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 index_t max_iter = strict_ ? 5000 : 50;
178 2 index_t iter = 0;
179 2 bool first = true;
180
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 7 while(changed && iter < max_iter) {
181
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if(first) {
182
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 CmdLine::ui_clear_line();
183 } else {
184 ✗ first = false;
185 }
186
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if(pretty) {
187
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 CmdLine::ui_message(
188
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 10 "o-[Manifold Rec] Iteration:" + String::to_string(iter)
189 );
190 } else {
191 ✗ Logger::out("Manifold Rec")
192 ✗ << "Iteration:" << iter << std::endl;
193 }
194 5 changed = false;
195 5 ++iter;
196
2/2
✓ Branch 0 taken 992 times.
✓ Branch 1 taken 5 times.
 997 for(index_t t=0; t<nb_triangles; ++t) {
197
2/2
✓ Branch 2 taken 359 times.
✓ Branch 3 taken 633 times.
 992 if(!t_is_classified[t]) {
198
1/2
✓ Branch 1 taken 359 times.
✗ Branch 2 not taken.
 359 index_t i = not_so_good_triangles[3*t];
199
1/2
✓ Branch 1 taken 359 times.
✗ Branch 2 not taken.
 359 index_t j = not_so_good_triangles[3*t+1];
200
1/2
✓ Branch 1 taken 359 times.
✗ Branch 2 not taken.
 359 index_t k = not_so_good_triangles[3*t+2];
201
1/2
✓ Branch 1 taken 359 times.
✗ Branch 2 not taken.
 359 index_t new_t = add_triangle(i,j,k);
202 359 bool classified = false;
203
3/4
✓ Branch 1 taken 359 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 174 times.
✓ Branch 4 taken 185 times.
 359 if(connect_and_validate_triangle(new_t, classified)) {
204 174 changed = true;
205 } else {
206
1/2
✓ Branch 1 taken 185 times.
✗ Branch 2 not taken.
 185 rollback_triangle();
207 }
208
2/2
✓ Branch 0 taken 356 times.
✓ Branch 1 taken 3 times.
 359 if(classified) {
209 356 t_is_classified[t] = true;
210 }
211 }
212 }
213 }
214
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(pretty) {
215
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 CmdLine::ui_clear_line();
216
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 CmdLine::ui_message(
217 2 "o-[Manifold Rec] Iteration:" +
218
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.
 6 String::to_string(iter) + "\n"
219 );
220 } else {
221 ✗ Logger::out("Manifold Rec")
222 ✗ << "Iteration:" << iter << std::endl;
223 }
224 2 }
225
226 protected:
227
228 /**
229 * \brief Initializes the combinatorial data
230 * structures and deletes all facets incident
231 * to a non-manifold edge.
232 */
233 2 void init_and_remove_non_manifold_edges() {
234
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 next_c_around_v_.assign(M_.facet_corners.nb(), NO_CORNER);
235
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 v2c_.assign(M_.vertices.nb(),NO_CORNER);
236
2/2
✓ Branch 5 taken 124806 times.
✓ Branch 6 taken 2 times.
 124808 for(index_t t: M_.facets) {
237
1/2
✓ Branch 1 taken 124806 times.
✗ Branch 2 not taken.
 124806 insert(t);
238 }
239
240 2 index_t nb_non_manifold = 0;
241 2 vector<index_t> remove_t;
242
2/2
✓ Branch 5 taken 124806 times.
✓ Branch 6 taken 2 times.
 124808 for(index_t t: M_.facets) {
243
2/4
✓ Branch 1 taken 124806 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 124806 times.
 124806 if(!connect(t)) {
244 ✗ remove_t.resize(M_.facets.nb(),0);
245 ✗ remove_t[t] = 1;
246 ✗ ++nb_non_manifold;
247 }
248 }
249
250
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 mesh_reorient(M_, &remove_t);
251
252
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 if(remove_t.size() == 0) {
253
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 6 Logger::out("Co3Ne")
254
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << "All edges are manifold and well oriented"
255
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << std::endl;
256 } else {
257 ✗ index_t nb_remove_t = 0;
258 ✗ for(index_t t: M_.facets) {
259 ✗ if(remove_t[t] != 0) {
260 ✗ ++nb_remove_t;
261 }
262 }
263 ✗ index_t nb_moebius = nb_remove_t - nb_non_manifold;
264 ✗ Logger::out("Co3Ne")
265 ✗ << "Removing " << nb_remove_t
266 ✗ << " triangles ("
267 ✗ << nb_non_manifold << " non_manifold, "
268 ✗ << nb_moebius
269 ✗ << " moebius)"
270 ✗ << std::endl;
271 ✗ M_.facets.delete_elements(remove_t,false);
272 }
273
274 // We need to re-compute next_c_around_v_ and v2c_
275 // since all the indices changed in the mesh
276 // (even if remove_t is empty, because mesh_reorient() may
277 // have changed triangles orientation).
278
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 next_c_around_v_.assign(M_.facet_corners.nb(), NO_CORNER);
279
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 v2c_.assign(M_.vertices.nb(),NO_CORNER);
280
2/2
✓ Branch 5 taken 124806 times.
✓ Branch 6 taken 2 times.
 124808 for(index_t t: M_.facets) {
281
1/2
✓ Branch 1 taken 124806 times.
✗ Branch 2 not taken.
 124806 insert(t);
282 }
283 2 }
284
285 /**
286 * \brief Tentatively connects a newly added triangle
287 * to the current mesh under construction. Accepted
288 * triangles satisfy the following criteria:
289 * - each new triangle should be either incident to at least
290 * two edges of existing triangles, or to one existing triangle
291 * and one isolated point.
292 * - the normals to the new triangle and its neighbor should
293 * not point to opposite directions.
294 * - inserting the new triangle should not generate 'by-excess'
295 * non-manifold vertices. A 'by-excess' non-manifold vertex
296 * has a closed loop of triangles in its neighbors plus
297 * additional triangles.
298 * - the orientation of the surface should be coherent (no Moebius
299 * strip).
300 * \param[in] t index of the triangle
301 * \param[out] classified true if the status of the triangle
302 * (accepted/rejected) could be completely determined,
303 * false if its status may still change during subsequent iterations
304 * \retval true if all combinatorial and geometric tests succeeded
305 * \retval false otherwise
306 */
307 359 bool connect_and_validate_triangle(index_t t, bool& classified) {
308 index_t adj_c[3];
309 359 classified = false;
310
311 // Combinatorial test (I): tests whether the three
312 // candidate edges are manifold.
313
3/4
✓ Branch 1 taken 359 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 182 times.
✓ Branch 4 taken 177 times.
 359 if(!get_adjacent_corners(t,adj_c)) {
314 182 classified = true;
315 182 return false ;
316 }
317
318 // Geometric test: tests whether the angles formed with
319 // the candidate neighbors do not indicate degenerate sharp
320 // creases.
321
2/2
✓ Branch 0 taken 531 times.
✓ Branch 1 taken 177 times.
 708 for(index_t i=0; i<3; ++i) {
322
2/2
✓ Branch 0 taken 438 times.
✓ Branch 1 taken 93 times.
 531 if(adj_c[i] != NO_CORNER) {
323
1/2
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
 438 index_t t2 = c2f(adj_c[i]);
324
2/4
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 438 times.
 438 if(!triangles_normals_agree(t,t2)) {
325 ✗ classified = true;
326 ✗ return false;
327 }
328 }
329 }
330
331 177 int nb_neighbors =
332 177 (adj_c[0] != NO_CORNER) +
333 177 (adj_c[1] != NO_CORNER) +
334 177 (adj_c[2] != NO_CORNER) ;
335
336 // Combinatorial test (II)
337
2/3
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 174 times.
 177 switch(nb_neighbors) {
338 // If the candidate triangle is adjacent to no other
339 // triangle, reject it
340 ✗ case 0: {
341 ✗ return false ;
342 }
343 // If the candidate triangle is adjacent to a single
344 // triangle, reject it if the vertex opposite to
345 // the common edge is not isolated.
346 3 case 1: {
347 // If not in strict mode, we reject the triangle.
348 // Experimentally, it improves the result.
349
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if(!strict_) {
350 3 return false;
351 }
352 ✗ index_t other_vertex=NO_INDEX;
353 ✗ for(index_t i=0; i<3; ++i) {
354 ✗ if(adj_c[i] != NO_CORNER) {
355 other_vertex =
356 ✗ M_.facet_corners.vertex(
357 ✗ M_.facets.corners_begin(t) + ((i+2)%3)
358 );
359 }
360 }
361 ✗ geo_debug_assert(other_vertex != NO_INDEX);
362 // Test whether other_vertex is isolated, reject
363 // the triangle if other_vertex is NOT isolated.
364 ✗ index_t nb_incident_T = nb_incident_triangles(other_vertex);
365 ✗ geo_assert(nb_incident_T != 0); // There is at least THIS T.
366 ✗ if(nb_incident_T > 1) {
367 ✗ return false;
368 }
369 }
370 }
371
372
1/2
✓ Branch 1 taken 174 times.
✗ Branch 2 not taken.
 174 connect_adjacent_corners(t,adj_c);
373
374 // Combinatorial test (III): test non-manifold vertices
375
6/10
✓ Branch 1 taken 174 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 174 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 174 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 522 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 522 times.
✓ Branch 15 taken 174 times.
 696 for(index_t v: M_.facets.vertices(t)) {
376 522 bool moebius=false;
377
2/4
✓ Branch 1 taken 522 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 522 times.
 522 if(vertex_is_non_manifold_by_excess(v,moebius)) {
378 ✗ classified = true;
379 ✗ return false;
380 }
381 // It should not occur since we remove all Moebius configs
382 // from the T3s and forbid Moebius configs when inserting
383 // the T12s. However, some transient moebius configurations
384 // due to triangle t may appear (since the Moebius test is
385 // right after the non-manifold test).
386
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 522 times.
 522 if(moebius) {
387 ✗ Logger::warn("Co3Ne")
388 ✗ << "Encountered Moebius configuration" << std::endl;
389 ✗ classified = true;
390 ✗ return false;
391 }
392 }
393
394 // Combinatorial test (IV): orientability
395
2/4
✓ Branch 1 taken 174 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 174 times.
 174 if(!enforce_orientation_from_triangle(t)) {
396 ✗ return false;
397 }
398
399 174 classified = true;
400 174 return true;
401 }
402
403
404 /**
405 * \brief Tentatively enforces mesh orientation starting from a
406 * given triangle.
407 * \details The triangle \p t is rejected if it is incident to
408 * the same connected component with two different orientations.
409 * \param[in] t index of the triangle to start mesh orientation from
410 * \retval true if the mesh could be coherently oriented
411 * \retval false otherwise
412 */
413 174 bool enforce_orientation_from_triangle(index_t t) {
414
415 // Index of adjacent triangle
416 // (or NO_FACET if no neighbor)
417 index_t adj[3];
418
419 // Index of adjacent connected component
420 // (or NO_CNX if no neighbor)
421 index_t adj_cnx[3];
422
423 // Orientation of adjacent triangle relative to
424 // triangle t (or 0 if no neighbor)
425 signed_index_t adj_ori[3];
426
427
2/2
✓ Branch 0 taken 522 times.
✓ Branch 1 taken 174 times.
 696 for(index_t i=0; i<3; ++i) {
428
1/2
✓ Branch 1 taken 522 times.
✗ Branch 2 not taken.
 522 index_t c = M_.facets.corners_begin(t)+i;
429
1/2
✓ Branch 1 taken 522 times.
✗ Branch 2 not taken.
 522 adj[i] = index_t(M_.facet_corners.adjacent_facet(c));
430 }
431
432
433
2/2
✓ Branch 0 taken 522 times.
✓ Branch 1 taken 174 times.
 696 for(index_t i=0; i<3; ++i) {
434
2/2
✓ Branch 0 taken 87 times.
✓ Branch 1 taken 435 times.
 522 if(adj[i] == NO_FACET) {
435 87 adj_ori[i] = 0;
436 87 adj_cnx[i] = NO_CNX;
437 } else {
438 435 adj_ori[i] =
439
3/4
✓ Branch 1 taken 435 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 208 times.
✓ Branch 4 taken 227 times.
 435 (triangles_have_same_orientation(t,adj[i])) ? 1 : -1;
440
1/2
✓ Branch 1 taken 435 times.
✗ Branch 2 not taken.
 435 adj_cnx[i] = cnx_[adj[i]];
441 }
442 }
443
444 // If in the neighborhood the same connected component appears
445 // with two opposite orientations, then connecting the triangle
446 // would create a Moebius strip (the triangle is rejected)
447
2/2
✓ Branch 0 taken 522 times.
✓ Branch 1 taken 174 times.
 696 for(index_t i=0; i<3; ++i) {
448
2/2
✓ Branch 0 taken 435 times.
✓ Branch 1 taken 87 times.
 522 if(adj[i] != NO_FACET) {
449
2/2
✓ Branch 0 taken 470 times.
✓ Branch 1 taken 435 times.
 905 for(index_t j=i+1; j<3; ++j) {
450 470 if(
451
2/2
✓ Branch 0 taken 346 times.
✓ Branch 1 taken 124 times.
 470 adj_cnx[j] == adj_cnx[i] &&
452
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 346 times.
 346 adj_ori[j] != adj_ori[i]
453 ) {
454 ✗ return false;
455 }
456 }
457 }
458 }
459
460 // The triangle is accepted,
461 // now reorient all the connected components and the
462 // triangle coherently.
463
464 // Find the largest component incident to t
465 174 index_t largest_neigh_comp = NO_CNX;
466
2/2
✓ Branch 0 taken 522 times.
✓ Branch 1 taken 174 times.
 696 for(index_t i=0; i<3; ++i) {
467 522 if(
468
6/6
✓ Branch 0 taken 435 times.
✓ Branch 1 taken 87 times.
✓ Branch 2 taken 261 times.
✓ Branch 3 taken 174 times.
✓ Branch 4 taken 174 times.
✓ Branch 5 taken 348 times.
 783 adj_cnx[i] != NO_CNX && (
469 261 largest_neigh_comp == NO_CNX ||
470
1/2
✓ Branch 1 taken 261 times.
✗ Branch 2 not taken.
 261 cnx_size_[adj_cnx[i]] >
471
2/4
✓ Branch 1 taken 261 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 261 times.
 261 cnx_size_[adj_cnx[largest_neigh_comp]]
472 )
473 ) {
474
475 174 largest_neigh_comp = i;
476 }
477 }
478
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 174 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 174 geo_assert(largest_neigh_comp != NO_CNX);
479
480 // Orient t like the largest incident component
481 174 index_t comp = adj_cnx[largest_neigh_comp];
482
483
1/2
✓ Branch 3 taken 174 times.
✗ Branch 4 not taken.
 174 cnx_.resize(std::max(t+1, cnx_.size()));
484
1/2
✓ Branch 1 taken 174 times.
✗ Branch 2 not taken.
 174 cnx_[t] = comp;
485
1/2
✓ Branch 1 taken 174 times.
✗ Branch 2 not taken.
 174 ++cnx_size_[comp];
486
2/2
✓ Branch 0 taken 90 times.
✓ Branch 1 taken 84 times.
 174 if(adj_ori[largest_neigh_comp] == -1) {
487
1/2
✓ Branch 1 taken 90 times.
✗ Branch 2 not taken.
 90 flip_triangle(t);
488
2/2
✓ Branch 0 taken 270 times.
✓ Branch 1 taken 90 times.
 360 for(index_t i=0; i<3; ++i) {
489 270 adj_ori[i] = -adj_ori[i];
490 }
491 }
492
493 // Merge (and reorient if need be) all the other incident
494 // components
495
2/2
✓ Branch 0 taken 522 times.
✓ Branch 1 taken 174 times.
 696 for(index_t i=0; i<3; ++i) {
496 522 if(
497 348 i != largest_neigh_comp &&
498
9/10
✓ Branch 0 taken 348 times.
✓ Branch 1 taken 174 times.
✓ Branch 2 taken 261 times.
✓ Branch 3 taken 87 times.
✓ Branch 5 taken 261 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 260 times.
✓ Branch 9 taken 1 times.
✓ Branch 10 taken 521 times.
 522 adj[i] != NO_FACET && cnx_[adj[i]] != comp
499 ) {
500 1 merge_connected_component(
501
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 adj[i], comp, (adj_ori[i] == -1)
502 );
503 }
504 }
505
506 174 return true;
507 }
508
509
510 /**
511 * \brief Adds a new triangle to the surface and to the
512 * combinatorial data structure.
513 * \param[in] i first index of the triangle
514 * \param[in] j second index of the triangle
515 * \param[in] k third index of the triangle
516 */
517 359 index_t add_triangle(index_t i, index_t j, index_t k) {
518 359 index_t result = M_.facets.create_triangle(i,j,k);
519 359 next_c_around_v_.push_back(NO_CORNER);
520 359 next_c_around_v_.push_back(NO_CORNER);
521 359 next_c_around_v_.push_back(NO_CORNER);
522 359 insert(result);
523 359 return result;
524 }
525
526 /**
527 * \brief Removes the latest triangle from both
528 * the mesh and the combinatorial data structure.
529 */
530 185 void rollback_triangle() {
531 185 index_t t = M_.facets.nb()-1;
532 185 remove(t);
533 185 M_.facets.pop();
534 185 }
535
536
537 /**
538 * \brief Inverts the orientation of a triangle.
539 * \param[in] t the index of the triangle to be flipped.
540 */
541 90 void flip_triangle(index_t t) {
542
543 // Remove t from the additional combinatorial data structure
544 // (it is both simpler and more efficient to do that
545 // than updating it).
546 90 remove(
547 t,
548 false // disconnect is set to false because
549 // we will re-insert t right after.
550 );
551
552 90 index_t c1 = M_.facets.corners_begin(t);
553 90 index_t c2 = c1+1;
554 90 index_t c3 = c2+1;
555 90 index_t v1 = M_.facet_corners.vertex(c1);
556 90 index_t f1 = M_.facet_corners.adjacent_facet(c1);
557 90 index_t f2 = M_.facet_corners.adjacent_facet(c2);
558 90 index_t v3 = M_.facet_corners.vertex(c3);
559
560 90 M_.facet_corners.set_vertex(c1,v3);
561 90 M_.facet_corners.set_adjacent_facet(c1,f2);
562 90 M_.facet_corners.set_adjacent_facet(c2,f1);
563 90 M_.facet_corners.set_vertex(c3,v1);
564
565 // Re-insert t into the additional combinatorial data structure.
566 90 insert(t);
567 90 }
568
569 /**
570 * \brief Inserts a triangle of the mesh into the data structures
571 * used for topology checks.
572 * \param[in] t index of the triangles to be inserted
573 * \pre \p t is a valid triangle index in the mesh
574 */
575 250061 void insert(index_t t) {
576
3/4
✓ Branch 1 taken 250061 times.
✗ Branch 2 not taken.
✓ Branch 8 taken 750183 times.
✓ Branch 9 taken 250061 times.
 1000244 for(index_t c : M_.facets.corners(t)) {
577
1/2
✓ Branch 1 taken 750183 times.
✗ Branch 2 not taken.
 750183 index_t v = M_.facet_corners.vertex(c);
578
3/4
✓ Branch 1 taken 750183 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 127048 times.
✓ Branch 4 taken 623135 times.
 750183 if(v2c_[v] == NO_CORNER) {
579
1/2
✓ Branch 1 taken 127048 times.
✗ Branch 2 not taken.
 127048 v2c_[v] = c;
580
1/2
✓ Branch 1 taken 127048 times.
✗ Branch 2 not taken.
 127048 next_c_around_v_[c] = c;
581 } else {
582
3/6
✓ Branch 1 taken 623135 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 623135 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 623135 times.
✗ Branch 8 not taken.
 623135 next_c_around_v_[c] = next_c_around_v_[v2c_[v]];
583
2/4
✓ Branch 1 taken 623135 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 623135 times.
✗ Branch 5 not taken.
 623135 next_c_around_v_[v2c_[v]] = c;
584 }
585 }
586 250061 }
587
588 /**
589 * \brief Removes a triangle of the mesh from the data structures
590 * used for topology/combinatorial checks.
591 * \param[in] t index of the triangles to be removed
592 * \param[in] disconnect if true, connections from the neighbors
593 * to t are set to -1 (facet_corners.adjacent_facet).
594 * \pre \p t is a valid triangle index in the mesh
595 */
596 275 void remove(index_t t, bool disconnect=true) {
597
2/2
✓ Branch 0 taken 185 times.
✓ Branch 1 taken 90 times.
 275 if(disconnect) {
598
6/10
✓ Branch 1 taken 185 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 185 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 185 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 555 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 555 times.
✓ Branch 15 taken 185 times.
 740 for(index_t t2: M_.facets.adjacent(t)) {
599 // Disconnect facet-facet link that point to t
600
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 555 times.
 555 if(t2 != NO_FACET) {
601 ✗ for(index_t c2: M_.facets.corners(t2)) {
602 ✗ if(
603 ✗ M_.facet_corners.adjacent_facet(c2) == t
604 ) {
605 ✗ M_.facet_corners.set_adjacent_facet(
606 c2,NO_FACET
607 );
608 }
609 }
610 }
611 }
612 }
613
614
3/4
✓ Branch 1 taken 275 times.
✗ Branch 2 not taken.
✓ Branch 8 taken 825 times.
✓ Branch 9 taken 275 times.
 1100 for(index_t c : M_.facets.corners(t)) {
615 // Remove t from combinatorial data structures
616
1/2
✓ Branch 1 taken 825 times.
✗ Branch 2 not taken.
 825 index_t v = M_.facet_corners.vertex(c);
617
2/4
✓ Branch 1 taken 825 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 825 times.
 825 if(next_c_around_v_[c] == c) {
618 ✗ v2c_[v] = NO_CORNER;
619 } else {
620
1/2
✓ Branch 1 taken 825 times.
✗ Branch 2 not taken.
 825 index_t c_pred = next_c_around_v_[c];
621
3/4
✓ Branch 1 taken 4322 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3497 times.
✓ Branch 4 taken 825 times.
 4322 while(next_c_around_v_[c_pred] != c) {
622
1/2
✓ Branch 1 taken 3497 times.
✗ Branch 2 not taken.
 3497 c_pred = next_c_around_v_[c_pred];
623 }
624
2/4
✓ Branch 1 taken 825 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 825 times.
✗ Branch 5 not taken.
 825 next_c_around_v_[c_pred] = next_c_around_v_[c];
625
1/2
✓ Branch 1 taken 825 times.
✗ Branch 2 not taken.
 825 v2c_[v] = c_pred;
626 }
627 }
628 275 }
629
630 /**
631 * \brief Gets the number of triangles incident
632 * to a vertex.
633 * \param[in] v index of the vertex
634 * \return the number of triangles incident to \p v
635 */
636 522 index_t nb_incident_triangles(index_t v) const {
637 522 index_t result = 0;
638 522 index_t c = v2c_[v];
639 do {
640 3059 ++result;
641 3059 c = next_c_around_v_[c];
642
2/2
✓ Branch 1 taken 2537 times.
✓ Branch 2 taken 522 times.
 3059 } while(c != v2c_[v]);
643 522 return result;
644 }
645
646 /**
647 * \brief Tests whether a given vertex is non-manifold
648 * by excess.
649 * \details A vertex is non-manifold by-excess if its
650 * set of incident triangles contains a closed loop
651 * of triangles and additional triangles.
652 * \param[in] v index of the vertex to be tested
653 * \retval true if \p v is non-manifold by excess
654 * \retval false otherwise
655 */
656 522 bool vertex_is_non_manifold_by_excess(index_t v, bool& moebius) {
657 522 index_t nb_v_neighbors = nb_incident_triangles(v);
658 522 index_t c = v2c_[v];
659 do {
660 3059 index_t loop_size=0;
661 3059 index_t c_cur = c ;
662 do {
663 17451 ++loop_size;
664
2/2
✓ Branch 0 taken 1022 times.
✓ Branch 1 taken 16429 times.
 17451 if(c_cur == NO_CORNER) {
665 1022 break;
666 }
667
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 16429 times.
 16429 if(loop_size > 100) {
668 // Probably Moebious strip or something...
669 ✗ moebius = true;
670 ✗ break;
671 }
672 16429 c_cur = next_around_vertex_unoriented(v,c_cur);
673
2/2
✓ Branch 0 taken 14392 times.
✓ Branch 1 taken 2037 times.
 16429 } while(c_cur != c);
674
675
3/4
✓ Branch 0 taken 2037 times.
✓ Branch 1 taken 1022 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2037 times.
 3059 if(c_cur == c && loop_size < nb_v_neighbors) {
676 ✗ return true;
677 }
678 3059 c = next_c_around_v_[c];
679
2/2
✓ Branch 1 taken 2537 times.
✓ Branch 2 taken 522 times.
 3059 } while(c != v2c_[v]);
680
681 522 return false;
682 }
683
684 /**
685 * \brief Gets the next corner around a vertex from a given
686 * corner.
687 * \details This function works even for a mesh that has triangles
688 * that are not coherently oriented. In other words, for two
689 * corners c1, c2, if we have:
690 * - v1 = facet_corners.vertex(c1)
691 * - v2 = facet_corners.vertex(
692 * c1,facets.next_corner_around_facet(c2f(c1),c1)
693 * )
694 * - w1 = facet_corners.vertex(c2)
695 * - w2 = facet_corners.vertex(
696 * c2,facets.next_corner_around_facet(c2f(c2),c2)
697 * )
698 * then we can have:
699 * - v1=w2 and v2=w1 (as usual) or:
700 * - v1=v2 and w1=w2 ('inverted' configuration)
701 * \param[in] v the vertex
702 * \param[in] c1 a corner incident to \p v or pointing to \p v
703 * \return another corner incident to the \p v
704 */
705 16429 index_t next_around_vertex_unoriented(
706 index_t v, index_t c1
707 ) const {
708 16429 index_t f1 = c2f(c1);
709 16429 index_t v1 = M_.facet_corners.vertex(c1);
710 16429 index_t v2 = M_.facet_corners.vertex(
711 16429 M_.facets.next_corner_around_facet(f1,c1)
712 );
713
714
3/8
✓ Branch 0 taken 2264 times.
✓ Branch 1 taken 14165 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2264 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 16429 geo_debug_assert(v1 == v || v2 == v);
715
716 16429 index_t f2 = M_.facet_corners.adjacent_facet(c1);
717
2/2
✓ Branch 0 taken 15407 times.
✓ Branch 1 taken 1022 times.
 16429 if(f2 != NO_FACET) {
718
2/4
✓ Branch 1 taken 15407 times.
✗ Branch 2 not taken.
✓ Branch 8 taken 31156 times.
✗ Branch 9 not taken.
 31156 for(index_t c2: M_.facets.corners(f2)) {
719
1/2
✓ Branch 1 taken 31156 times.
✗ Branch 2 not taken.
 31156 index_t w1 = M_.facet_corners.vertex(c2);
720
1/2
✓ Branch 1 taken 31156 times.
✗ Branch 2 not taken.
 31156 index_t w2 = M_.facet_corners.vertex(
721
1/2
✓ Branch 1 taken 31156 times.
✗ Branch 2 not taken.
 31156 M_.facets.next_corner_around_facet(f2,c2)
722 );
723
2/2
✓ Branch 0 taken 6930 times.
✓ Branch 1 taken 24226 times.
 31156 if(
724
4/4
✓ Branch 0 taken 4473 times.
✓ Branch 1 taken 2457 times.
✓ Branch 2 taken 14615 times.
✓ Branch 3 taken 14084 times.
 31156 (v1 == w1 && v2 == w2) ||
725
2/2
✓ Branch 0 taken 12950 times.
✓ Branch 1 taken 1665 times.
 14615 (v1 == w2 && v2 == w1)
726 ) {
727
2/2
✓ Branch 0 taken 13143 times.
✓ Branch 1 taken 2264 times.
 15407 if(w2 == v) {
728
1/2
✓ Branch 1 taken 13143 times.
✗ Branch 2 not taken.
 15407 return M_.facets.next_corner_around_facet(f2,c2);
729 } else {
730
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2264 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2264 geo_debug_assert(w1 == v);
731
1/2
✓ Branch 1 taken 2264 times.
✗ Branch 2 not taken.
 2264 return M_.facets.prev_corner_around_facet(f2,c2);
732 }
733 }
734 }
735 }
736 1022 return NO_CORNER;
737 }
738
739 /**
740 * \brief Gets the three corners adjacent to a triangle.
741 * \details This function works even for a mesh that has triangles
742 * that are not coherently oriented. In other words, for two
743 * corners c1, c2, if we have:
744 * - v1 = facet_corners.vertex(c1)
745 * - v2 = facet_corners.vertex(
746 * c1,facets.next_corner_around_facet(c2f(c1),c1))
747 * - w1 = facet_corners.vertex(c2)
748 * - w2 = facet_corners.vertex(
749 * c2,facets.next_corner_around_facet(c2f(c2),c2))
750 * then c1 and c2 are adjacent if we have:
751 * - v1=w2 and v2=w1 (as usual) or:
752 * - v1=v2 and w1=w2 ('inverted' configuration)
753 * \param[in] t1 index of the triangle
754 * \param[out] adj_c index of the adjacent corners
755 * (array of 3 integers). Each entry contains a valid corner index
756 * or NO_CORNER if the corresponding edge is on the border.
757 * \retval true if the three edges are manifold
758 * \retval false otherwise (and then \p adj_c contains undefined
759 * values).
760 */
761 125165 bool get_adjacent_corners(index_t t1, index_t* adj_c) {
762
3/4
✓ Branch 1 taken 125165 times.
✗ Branch 2 not taken.
✓ Branch 7 taken 375175 times.
✓ Branch 8 taken 124983 times.
 500158 for(index_t c1: M_.facets.corners(t1)) {
763
1/2
✓ Branch 1 taken 375175 times.
✗ Branch 2 not taken.
 375175 index_t v2 = M_.facet_corners.vertex(
764
1/2
✓ Branch 1 taken 375175 times.
✗ Branch 2 not taken.
 375175 M_.facets.next_corner_around_facet(t1,c1)
765 );
766
767 375175 *adj_c = NO_CORNER;
768
769 // Traverse the circular incident edge list
770
1/2
✓ Branch 1 taken 375175 times.
✗ Branch 2 not taken.
 375175 index_t c2=next_c_around_v_[c1];
771
2/2
✓ Branch 0 taken 1898501 times.
✓ Branch 1 taken 374993 times.
 2273494 while(c2 != c1) {
772
1/2
✓ Branch 1 taken 1898501 times.
✗ Branch 2 not taken.
 1898501 index_t t2 = c2f(c2);
773
1/2
✓ Branch 1 taken 1898501 times.
✗ Branch 2 not taken.
 1898501 index_t c3 = M_.facets.prev_corner_around_facet(t2,c2);
774
1/2
✓ Branch 1 taken 1898501 times.
✗ Branch 2 not taken.
 1898501 index_t v3 = M_.facet_corners.vertex(c3);
775
2/2
✓ Branch 0 taken 159007 times.
✓ Branch 1 taken 1739494 times.
 1898501 if(v3 == v2) {
776 // Found an adjacent edge
777
2/2
✓ Branch 0 taken 158922 times.
✓ Branch 1 taken 85 times.
 159007 if(*adj_c == NO_CORNER) {
778 158922 *adj_c = c3;
779
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 158922 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 158922 geo_debug_assert(c3 != c1);
780 } else {
781 // If there was already an adjacent edge,
782 // then this is a non-manifold configuration
783 182 return false;
784 }
785 }
786
787 // Check with the other (wrong) orientation
788
1/2
✓ Branch 1 taken 1898416 times.
✗ Branch 2 not taken.
 1898416 c3 = M_.facets.next_corner_around_facet(t2,c2);
789
1/2
✓ Branch 1 taken 1898416 times.
✗ Branch 2 not taken.
 1898416 v3 = M_.facet_corners.vertex(c3);
790
2/2
✓ Branch 0 taken 213804 times.
✓ Branch 1 taken 1684612 times.
 1898416 if(v3 == v2) {
791 // Found an adjacent edge
792
2/2
✓ Branch 0 taken 213707 times.
✓ Branch 1 taken 97 times.
 213804 if(*adj_c == NO_CORNER) {
793 213707 *adj_c = c2;
794
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 213707 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 213707 geo_debug_assert(c2 != c1);
795 } else {
796 // If there was already an adjacent edge,
797 // then this is a non-manifold configuration
798 97 return false;
799 }
800 }
801
1/2
✓ Branch 1 taken 1898319 times.
✗ Branch 2 not taken.
 1898319 c2 = next_c_around_v_[c2];
802 }
803 374993 ++adj_c;
804 }
805 124983 return true;
806 }
807
808 /**
809 * \brief Tentatively connect a triangle of the mesh with its
810 * neighbors.
811 * \details This function is independent of triangles orientations,
812 * see get_adjacent_corners().
813 * \param[in] t index of the triangle to be connected
814 * \param[in] adj_c an array of three integers that indicate
815 * for each corner of the triangle the index of the adjacent
816 * corner or NO_CORNER if the corner is on the border.
817 */
818 124980 void connect_adjacent_corners(index_t t, index_t* adj_c) {
819
2/2
✓ Branch 0 taken 374940 times.
✓ Branch 1 taken 124980 times.
 499920 for(index_t i=0; i<3; ++i) {
820
2/2
✓ Branch 0 taken 372427 times.
✓ Branch 1 taken 2513 times.
 374940 if(adj_c[i] != NO_CORNER) {
821 372427 index_t c = M_.facets.corners_begin(t)+i;
822 372427 M_.facet_corners.set_adjacent_facet(c, c2f(adj_c[i]));
823 372427 M_.facet_corners.set_adjacent_facet(adj_c[i], t);
824 }
825 }
826 124980 }
827
828 /**
829 * \brief Tentatively connect a triangle of the mesh with its
830 * neighbors.
831 * \details This function is independent of triangles orientations,
832 * see get_adjacent_corners().
833 * \param[in] t index of the triangle to be connected
834 * \retval false if the connection would have created non-manifold
835 * edges
836 * \retval true otherwise
837 */
838 124806 bool connect(index_t t) {
839 124806 index_t adj_c[3] = {NO_CORNER, NO_CORNER, NO_CORNER};
840
2/4
✓ Branch 1 taken 124806 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 124806 times.
 124806 if(!get_adjacent_corners(t,adj_c)) {
841 ✗ return false;
842 }
843
1/2
✓ Branch 1 taken 124806 times.
✗ Branch 2 not taken.
 124806 connect_adjacent_corners(t, adj_c);
844 124806 return true;
845 }
846
847 /**
848 * \brief Gets a facet index by corner index.
849 * \details for a triangulated mesh, indexing is
850 * implicit, and we do not need to store a c2f array.
851 * \param[in] c corner index
852 * \return the index of the facet incident to c
853 */
854 2287795 index_t c2f(index_t c) const {
855
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2287795 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2287795 geo_debug_assert(c != NO_CORNER);
856
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2287795 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2287795 geo_debug_assert(c < M_.facet_corners.nb());
857 2287795 return c/3;
858 }
859
860
861 /**
862 * \brief Tests whether two triangles have the
863 * same orientation.
864 * \param[in] t1 first triangle
865 * \param[in] t2 second triangle
866 * \retval true if \p t1 and \p t2 have the same
867 * orientation
868 * \retval false otherwise
869 * \pre \p t1 and \p t2 share an edge
870 */
871 435 bool triangles_have_same_orientation(
872 index_t t1,
873 index_t t2
874 ) {
875 435 index_t c1 = M_.facets.corners_begin(t1);
876 435 index_t i1 = M_.facet_corners.vertex(c1);
877 435 index_t j1 = M_.facet_corners.vertex(c1+1);
878 435 index_t k1 = M_.facet_corners.vertex(c1+2);
879
880 435 index_t c2 = M_.facets.corners_begin(t2);
881 435 index_t i2 = M_.facet_corners.vertex(c2);
882 435 index_t j2 = M_.facet_corners.vertex(c2+1);
883 435 index_t k2 = M_.facet_corners.vertex(c2+2);
884
885
2/2
✓ Branch 0 taken 94 times.
✓ Branch 1 taken 341 times.
 435 if(
886
4/4
✓ Branch 0 taken 89 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 80 times.
✓ Branch 3 taken 350 times.
 435 (i1==i2 && j1==j2) ||
887
4/4
✓ Branch 0 taken 46 times.
✓ Branch 1 taken 34 times.
✓ Branch 2 taken 127 times.
✓ Branch 3 taken 269 times.
 430 (i1==k2 && j1==i2) ||
888
4/4
✓ Branch 0 taken 80 times.
✓ Branch 1 taken 47 times.
✓ Branch 2 taken 79 times.
✓ Branch 3 taken 270 times.
 396 (i1==j2 && j1==k2) ||
889
4/4
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 28 times.
✓ Branch 2 taken 120 times.
✓ Branch 3 taken 201 times.
 349 (k1==k2 && i1==i2) ||
890
4/4
✓ Branch 0 taken 94 times.
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 68 times.
✓ Branch 3 taken 227 times.
 321 (k1==j2 && i1==k2) ||
891
4/4
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 27 times.
✓ Branch 3 taken 252 times.
 295 (k1==i2 && i1==j2) ||
892
4/4
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 97 times.
✓ Branch 3 taken 168 times.
 279 (j1==j2 && k1==k2) ||
893
4/4
✓ Branch 0 taken 65 times.
✓ Branch 1 taken 32 times.
✓ Branch 2 taken 92 times.
✓ Branch 3 taken 141 times.
 265 (j1==i2 && k1==j2) ||
894
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 67 times.
 92 (j1==k2 && k1==i2)
895 ) {
896 227 return false;
897 }
898
899 208 return true;
900 }
901
902
903 /**
904 * \brief Tests whether the normals of two triangles that
905 * share an edge 'agree', i.e. whether they do not form
906 * a too sharp angle.
907 * \param[in] t1 index of the first triangle
908 * \param[in] t2 index of the second triangle
909 * \retval true if the normals of both triangles do not
910 * point in opposite directions
911 * \retval false otherwise
912 * \pre the two triangles are incident to the same edge
913 * (they have two vertices in common)
914 */
915 438 bool triangles_normals_agree(
916 index_t t1,
917 index_t t2
918 ) const {
919 const vec3* points =
920
1/2
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
 438 reinterpret_cast<const vec3*>(M_.vertices.point_ptr(0));
921
922
1/2
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
 438 index_t c1 = M_.facets.corners_begin(t1);
923
1/2
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
 438 index_t i1 = M_.facet_corners.vertex(c1);
924
1/2
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
 438 index_t j1 = M_.facet_corners.vertex(c1+1);
925
1/2
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
 438 index_t k1 = M_.facet_corners.vertex(c1+2);
926
927
1/2
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
 438 index_t c2 = M_.facets.corners_begin(t2);
928
1/2
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
 438 index_t i2 = M_.facet_corners.vertex(c2);
929
1/2
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
 438 index_t j2 = M_.facet_corners.vertex(c2+1);
930
1/2
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
 438 index_t k2 = M_.facet_corners.vertex(c2+2);
931
932
1/2
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
 438 vec3 n1 = normalize(
933 438 cross(
934 438 points[j1] - points[i1],
935 438 points[k1] - points[i1]
936 )
937 );
938
939
1/2
✓ Branch 1 taken 438 times.
✗ Branch 2 not taken.
 438 vec3 n2 = normalize(
940 438 cross(
941 438 points[j2] - points[i2],
942 438 points[k2] - points[i2]
943 )
944 );
945
946 438 double d = dot(n1,n2);
947 // Test for combinatorial orientation,
948 // if t1 and t2 have opposite orientation,
949 // then we flip one of the normals (i.e.,
950 // we simply change the sign of the dot product).
951
2/2
✓ Branch 0 taken 94 times.
✓ Branch 1 taken 344 times.
 438 if(
952
4/4
✓ Branch 0 taken 89 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 80 times.
✓ Branch 3 taken 353 times.
 438 (i1==i2 && j1==j2) ||
953
4/4
✓ Branch 0 taken 46 times.
✓ Branch 1 taken 34 times.
✓ Branch 2 taken 129 times.
✓ Branch 3 taken 270 times.
 433 (i1==k2 && j1==i2) ||
954
4/4
✓ Branch 0 taken 82 times.
✓ Branch 1 taken 47 times.
✓ Branch 2 taken 79 times.
✓ Branch 3 taken 273 times.
 399 (i1==j2 && j1==k2) ||
955
4/4
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 28 times.
✓ Branch 2 taken 121 times.
✓ Branch 3 taken 203 times.
 352 (k1==k2 && i1==i2) ||
956
4/4
✓ Branch 0 taken 95 times.
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 69 times.
✓ Branch 3 taken 229 times.
 324 (k1==j2 && i1==k2) ||
957
4/4
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 17 times.
✓ Branch 2 taken 27 times.
✓ Branch 3 taken 254 times.
 298 (k1==i2 && i1==j2) ||
958
4/4
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 99 times.
✓ Branch 3 taken 168 times.
 281 (j1==j2 && k1==k2) ||
959
4/4
✓ Branch 0 taken 66 times.
✓ Branch 1 taken 33 times.
✓ Branch 2 taken 92 times.
✓ Branch 3 taken 142 times.
 267 (j1==i2 && k1==j2) ||
960
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 67 times.
 92 (j1==k2 && k1==i2)
961 ) {
962 229 d = -d;
963 }
964 438 return (d > -0.8);
965 }
966
967 /**
968 * \brief Merges two connected components.
969 * \details The connected component incident to \p t
970 * is replaced with \p comp2.
971 * \param [in] t index of a triangle incident
972 * to the first connected component
973 * \param [in] comp2 index of the second connected
974 * component
975 * \param [in] flip if true, flip the triangles
976 * \pre At least one of the triangles adjacent to
977 * \p t (directly or not) is incident to
978 * component \p comp2
979 */
980 1 void merge_connected_component(
981 index_t t,
982 index_t comp2,
983 bool flip
984 ) {
985
2/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 1 geo_assert(comp2 != cnx_[t]);
986
987
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 std::stack<index_t> S;
988
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 index_t comp1 = cnx_[t];
989
990
991
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 cnx_[t] = comp2;
992
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 --cnx_size_[comp1];
993
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 ++cnx_size_[comp2];
994
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if(flip) {
995 ✗ flip_triangle(t);
996 }
997
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 S.push(t);
998
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 while(!S.empty()) {
999 1 index_t t1 = S.top();
1000 1 S.pop();
1001
6/10
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 3 times.
✓ Branch 15 taken 1 times.
 4 for(index_t t2: M_.facets.adjacent(t1)) {
1002
5/8
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 3 times.
 3 if(t2 != NO_FACET && cnx_[t2] == comp1) {
1003 ✗ cnx_[t2] = comp2;
1004 ✗ --cnx_size_[comp1];
1005 ✗ ++cnx_size_[comp2];
1006 ✗ if(flip) {
1007 ✗ flip_triangle(t2);
1008 }
1009 ✗ S.push(t2);
1010 }
1011 }
1012 }
1013
2/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 1 geo_assert(cnx_size_[comp1] == 0);
1014 1 }
1015
1016 /**
1017 * \brief Initializes the date structures
1018 * that represent the connected components.
1019 * \details This function computes cnx_ and
1020 * cnx_size_. The array cnx_[f] gives for each
1021 * facet f the index of the connected component
1022 * that contains f, and the array cnx_size_[comp]
1023 * gives for each connected component comp the
1024 * number of facets in comp.
1025 */
1026 2 void init_connected_components() {
1027 2 cnx_.assign(M_.facets.nb(), NO_CNX);
1028 2 cnx_size_.clear();
1029
2/2
✓ Branch 1 taken 124806 times.
✓ Branch 2 taken 2 times.
 124808 for(index_t t=0; t<M_.facets.nb(); ++t) {
1030
3/4
✓ Branch 1 taken 124806 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 124803 times.
 124806 if(cnx_[t] == NO_CNX) {
1031 3 index_t cnx_id = cnx_size_.size();
1032 3 index_t nb = 0;
1033
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 std::stack<index_t> S;
1034
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 S.push(t);
1035
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 cnx_[t] = cnx_id;
1036 3 ++nb;
1037
2/2
✓ Branch 1 taken 124806 times.
✓ Branch 2 taken 3 times.
 124809 while(!S.empty()) {
1038 124806 index_t t2 = S.top();
1039 124806 S.pop();
1040
6/10
✓ Branch 1 taken 124806 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 124806 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 124806 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 374418 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 374418 times.
✓ Branch 15 taken 124806 times.
 499224 for(index_t t3: M_.facets.adjacent(t2)) {
1041
7/8
✓ Branch 0 taken 371992 times.
✓ Branch 1 taken 2426 times.
✓ Branch 3 taken 371992 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 124803 times.
✓ Branch 6 taken 247189 times.
✓ Branch 7 taken 124803 times.
✓ Branch 8 taken 249615 times.
 374418 if(t3 != NO_FACET && cnx_[t3] != cnx_id) {
1042
2/8
✓ Branch 1 taken 124803 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 124803 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 124803 geo_assert(cnx_[t3] == NO_CNX);
1043
1/2
✓ Branch 1 taken 124803 times.
✗ Branch 2 not taken.
 124803 cnx_[t3] = cnx_id;
1044 124803 ++nb;
1045
1/2
✓ Branch 1 taken 124803 times.
✗ Branch 2 not taken.
 124803 S.push(t3);
1046 }
1047 }
1048 }
1049
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 cnx_size_.push_back(nb);
1050 3 }
1051 }
1052
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 Logger::out("Co3Ne")
1053
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
 2 << "Found " << cnx_size_.size() << " connected components"
1054
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << std::endl;
1055 2 }
1056
1057 private:
1058 Mesh& M_;
1059
1060 /**
1061 * \brief For each corner, next_c_around_v_[c]
1062 * chains the circular list of corners
1063 * incident to the same corner as c.
1064 */
1065 vector<index_t> next_c_around_v_;
1066
1067 /**
1068 * \brief For each vertex v, v2c_[v] contains a
1069 * corner incident to v, or NO_VERTEX if v is
1070 * isolated.
1071 */
1072 vector<index_t> v2c_;
1073
1074
1075 /**
1076 * \brief For each triangle t, cnx_[t] contains
1077 * the index of the connected component of the
1078 * mesh incident to t.
1079 */
1080 vector<index_t> cnx_;
1081
1082 /**
1083 * \brief For each connected component C,
1084 * cnx_size_[C] contains the number of
1085 * facets in C.
1086 */
1087 vector<index_t> cnx_size_;
1088
1089 /**
1090 * \brief In strict mode, each inserted triangle
1091 * is checked for non-manifold configuration.
1092 * In non-strict mode, only T2 and T1 triangles are
1093 * tested (those seen from only 2 or only 1 Voronoi
1094 * cell), T3 triangles are inserted without test.
1095 */
1096 bool strict_;
1097 };
1098
1099 /************************************************************/
1100
1101 /**
1102 * \brief Comparator class for sorting facets.
1103 */
1104 class CompareTriangles {
1105 public:
1106 /**
1107 * \brief Constructs a new CompareFacets.
1108 * \param[in] triangles a const reference to a vector
1109 * of indices triplets
1110 */
1111 2 explicit CompareTriangles(const vector<index_t>& triangles) :
1112 2 triangles_(triangles) {
1113 2 }
1114
1115 /**
1116 * \brief Tests the lexicographic order of two facets by their indices.
1117 * \param[in] f1 index of the first facet
1118 * \param[in] f2 index of the second facet
1119 * \return true if facet \p f1 is before facet \p f2 according to
1120 * the lexicographic order of its vertices, false otherwise.
1121 */
1122 10031635 bool is_before(index_t f1, index_t f2) const {
1123
2/2
✓ Branch 0 taken 18724043 times.
✓ Branch 1 taken 3468409 times.
 22192452 for(index_t c=0; c<3; c++) {
1124 18724043 index_t v1 = triangles_[3*f1+c];
1125 18724043 index_t v2 = triangles_[3*f2+c];
1126
2/2
✓ Branch 0 taken 3995070 times.
✓ Branch 1 taken 14728973 times.
 18724043 if(v1 > v2) {
1127 3995070 return false;
1128 }
1129
2/2
✓ Branch 0 taken 2568156 times.
✓ Branch 1 taken 12160817 times.
 14728973 if(v1 < v2) {
1130 2568156 return true;
1131 }
1132 }
1133 3468409 return false;
1134 }
1135
1136 /**
1137 * \brief Tests whether two facets are identical.
1138 * \param[in] f1 index of the first facet
1139 * \param[in] f2 index of the second facet
1140 * \return true if facets \p f1 and \p f2 have the same
1141 * vertices, false otherwise
1142 */
1143 374956 bool is_same(index_t f1, index_t f2) const {
1144
2/2
✓ Branch 0 taken 946863 times.
✓ Branch 1 taken 249796 times.
 1196659 for(index_t c=0; c<3; c++) {
1145 946863 index_t v1 = triangles_[3*f1+c];
1146 946863 index_t v2 = triangles_[3*f2+c];
1147
2/2
✓ Branch 0 taken 125160 times.
✓ Branch 1 taken 821703 times.
 946863 if(v1 != v2) {
1148 125160 return false;
1149 }
1150 }
1151 249796 return true;
1152 }
1153
1154 /**
1155 * \brief Tests the lexicographic order of two facets by their indices.
1156 * \param[in] f1 index of the first facet
1157 * \param[in] f2 index of the second facet
1158 * \return true if facet \p f1 is before facet \p f2 according to
1159 * the lexicographic order of its vertices, false otherwise.
1160 */
1161 10031635 bool operator() (index_t f1, index_t f2) const {
1162 10031635 return is_before(f1, f2);
1163 }
1164
1165 private:
1166 const vector<index_t>& triangles_;
1167 };
1168
1169
1170 /**
1171 * \brief Splits the raw list of triangles reconstructed
1172 * by the Co3Ne algorithm into two lists, good triangles
1173 * and "not so good" triangles.
1174 * \details The triangles that appear 3 times (seen from 3
1175 * different Voronoi cells) are the good ones, else they
1176 * are the "not so good" ones.
1177 * \param[in,out] triangles the input list of triangles. It
1178 * is modified by the algorithm (it is reordered).
1179 * \param[out] good_triangles the good triangles
1180 * \param[out] not_so_good_triangles the not-so-good triangles
1181 */
1182 2 static void co3ne_split_triangles_list(
1183 vector<index_t>& triangles,
1184 vector<index_t>& good_triangles,
1185 vector<index_t>& not_so_good_triangles
1186 ) {
1187 2 index_t nb_triangles = triangles.size()/3;
1188
1189 // Step 1: normalize vertices order
1190
2/2
✓ Branch 1 taken 374958 times.
✓ Branch 2 taken 2 times.
 374960 for(index_t i=0; i<triangles.size(); i+=3) {
1191
1/2
✓ Branch 1 taken 374958 times.
✗ Branch 2 not taken.
 374958 index_t* ptr = &triangles[i];
1192
1/2
✓ Branch 1 taken 374958 times.
✗ Branch 2 not taken.
 374958 std::sort(ptr, ptr+3);
1193 }
1194
1195 // Step 2: sort the triangles in lexicographic order
1196
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 vector<index_t> t_sort(nb_triangles);
1197
2/2
✓ Branch 0 taken 374958 times.
✓ Branch 1 taken 2 times.
 374960 for(index_t t=0; t<nb_triangles; ++t) {
1198
1/2
✓ Branch 1 taken 374958 times.
✗ Branch 2 not taken.
 374958 t_sort[t] = t ;
1199 }
1200 2 CompareTriangles compare_triangles(triangles);
1201
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 GEO::sort(t_sort.begin(), t_sort.end(), compare_triangles);
1202
1203
1204 // Step 3: select the triangles that appear exactly 3 times
1205 2 index_t if1 = 0;
1206
2/2
✓ Branch 0 taken 125162 times.
✓ Branch 1 taken 2 times.
 125164 while(if1 < nb_triangles) {
1207 125162 index_t if2 = if1 + 1;
1208 125162 while(
1209
4/4
✓ Branch 0 taken 374956 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 249796 times.
✓ Branch 3 taken 125162 times.
 749914 if2 < nb_triangles &&
1210
5/8
✓ Branch 1 taken 374956 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 374956 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 374956 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 249796 times.
✓ Branch 10 taken 125160 times.
 374956 compare_triangles.is_same(t_sort[if1], t_sort[if2])
1211 ) {
1212 249796 if2++;
1213 }
1214
1215
1/2
✓ Branch 1 taken 125162 times.
✗ Branch 2 not taken.
 125162 index_t t = t_sort[if1];
1216
2/2
✓ Branch 0 taken 124806 times.
✓ Branch 1 taken 356 times.
 125162 if(if2 - if1 == 3) {
1217
2/4
✓ Branch 1 taken 124806 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 124806 times.
✗ Branch 5 not taken.
 124806 good_triangles.push_back(triangles[3*t]);
1218
2/4
✓ Branch 1 taken 124806 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 124806 times.
✗ Branch 5 not taken.
 124806 good_triangles.push_back(triangles[3*t+1]);
1219
2/4
✓ Branch 1 taken 124806 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 124806 times.
✗ Branch 5 not taken.
 124806 good_triangles.push_back(triangles[3*t+2]);
1220
1/2
✓ Branch 0 taken 356 times.
✗ Branch 1 not taken.
 356 } else if(if2 - if1 <= 2) {
1221
2/4
✓ Branch 1 taken 356 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 356 times.
✗ Branch 5 not taken.
 356 not_so_good_triangles.push_back(triangles[3*t]);
1222
2/4
✓ Branch 1 taken 356 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 356 times.
✗ Branch 5 not taken.
 356 not_so_good_triangles.push_back(triangles[3*t+1]);
1223
2/4
✓ Branch 1 taken 356 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 356 times.
✗ Branch 5 not taken.
 356 not_so_good_triangles.push_back(triangles[3*t+2]);
1224 }
1225 125162 if1 = if2;
1226 }
1227 2 }
1228
1229 /************************************************************/
1230
1231 /**
1232 * \brief Used internally by the reconstruction algorithm.
1233 * Co3NeRestrictedVoronoiDiagram computes the restricted
1234 * Voronoi diagram of a set of disks.
1235 *
1236 * \details Given a point set with normals and a radius, this class
1237 * computes the intersection between the Voronoi diagram of
1238 * the points and the disks centered on the points and
1239 * orthogonal to the normals.
1240 */
1241 class Co3NeRestrictedVoronoiDiagram {
1242 public:
1243 /**
1244 * \brief Stores a 3D point and the combinatorial information
1245 * (index of the adjacent seed). The combinatorial information
1246 * is used to reconstruct the triangles at the end of the
1247 * algorithm.
1248 */
1249 class Vertex {
1250 public:
1251 /**
1252 * \brief Constructs a new uninitialized Vertex.
1253 */
1254 770936 Vertex() {
1255 770936 }
1256
1257 /**
1258 * \brief Constructs a Vertex from a 3d point.
1259 */
1260 635240 Vertex(const vec3& v) :
1261 635240 point_(v),
1262 635240 adjacent_seed_(-1) {
1263 635240 }
1264
1265 /**
1266 * \brief Gets the 3d point associated with this vertex.
1267 * \return a const reference to the 3d point
1268 */
1269 27984975 const vec3& point() const {
1270 27984975 return point_;
1271 }
1272
1273 /**
1274 * \brief Gets the 3d point associated with this vertex.
1275 * \return a const reference to the 3d point
1276 */
1277 2308008 vec3& point() {
1278 2308008 return point_;
1279 }
1280
1281 /**
1282 * \brief Gets the index of the adjacent seed associated with
1283 * this vertex.
1284 * \details Each vertex stores combinatorial information, i.e.
1285 * the index of the adjacent Voronoi seed accros the edge
1286 * starting from this vertex
1287 * \return the index of the adjacent Voronoi seed
1288 */
1289 1143556 signed_index_t adjacent_seed() const {
1290 1143556 return adjacent_seed_;
1291 }
1292
1293 /**
1294 * \brief Sets the index of the adjacent seed associated with
1295 * this vertex.
1296 * \details Each vertex stores combinatorial information, i.e.
1297 * the index of the adjacent Voronoi seed accros the edge
1298 * starting from this vertex
1299 * \param[in] x the index of the adjacent Voronoi seed
1300 */
1301 769336 void set_adjacent_seed(signed_index_t x) {
1302 769336 adjacent_seed_ = x;
1303 769336 }
1304
1305 private:
1306 vec3 point_;
1307 signed_index_t adjacent_seed_;
1308 };
1309
1310 /**
1311 * \brief Internal representation of the polygons, that represent
1312 * the intersection between the disks and the Voronoi cells.
1313 */
1314 class Polygon {
1315 public:
1316 /**
1317 * \brief Creates a new uninitialized polygon with a given
1318 * number of vertices.
1319 * \param[in] size number of vertices
1320 */
1321 16 Polygon(index_t size) :
1322 16 vertices_(size) {
1323 16 }
1324
1325 /**
1326 * \brief Gets the number of vertices.
1327 * \return the number of vertices of this Polygon
1328 */
1329 30536857 index_t nb_vertices() const {
1330 30536857 return vertices_.size();
1331 }
1332
1333 /**
1334 * \brief Adds a new vertex to this Polygon.
1335 * \param[in] v the vertex to be added.
1336 */
1337 10845872 void add_vertex(const Vertex& v) {
1338 10845872 vertices_.push_back(v);
1339 10845872 }
1340
1341 /**
1342 * \brief Gets a Vertex by its index.
1343 * \param[in] i the index of the Vertex
1344 * \return a reference to the Vertex
1345 */
1346 13288827 Vertex& vertex(index_t i) {
1347 13288827 return vertices_[i];
1348 }
1349
1350 /**
1351 * \brief Gets a Vertex by its index.
1352 * \param[in] i the index of the Vertex
1353 * \return a const reference to the Vertex
1354 */
1355 10839020 const Vertex& vertex(index_t i) const {
1356 10839020 return vertices_[i];
1357 }
1358
1359 /**
1360 * \brief Gets the index of the next vertex around
1361 * the polygon.
1362 * \param[in] i index of the vertex
1363 * \return index of the next vertex (successor of \p i)
1364 * around the Polygon.
1365 */
1366 379444 index_t next_vertex(index_t i) const {
1367
2/2
✓ Branch 1 taken 315920 times.
✓ Branch 2 taken 63524 times.
 379444 return (i == nb_vertices() - 1) ? 0 : i + 1;
1368 }
1369
1370 /**
1371 * \brief Removes all the vertices.
1372 */
1373 2127035 void clear() {
1374 2127035 vertices_.resize(0);
1375 2127035 }
1376
1377 /**
1378 * \brief Swaps the vertices of this Polygon with
1379 * the vertices of another polygon.
1380 * \param[in] P the other polygon
1381 */
1382 2063511 void swap(Polygon& P) {
1383 2063511 vertices_.swap(P.vertices_);
1384 2063511 }
1385
1386 private:
1387 vector<Vertex> vertices_;
1388 };
1389
1390 /**
1391 * \brief Constructs a new uninitialized Co3NeRestrictedVoronoiDiagram.
1392 */
1393 2 Co3NeRestrictedVoronoiDiagram() :
1394 2 nb_points_(0),
1395 2 p_(nullptr),
1396 2 p_stride_(0),
1397 2 n_(nullptr),
1398 2 n_stride_(0),
1399 2 radius_(0.0),
1400
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 NN_(NearestNeighborSearch::create(3)),
1401 2 sqROS_(0.0),
1402 2 nb_neighbors_(0)
1403 {
1404 2 }
1405
1406 /**
1407 * \brief Co3NeRestrictedVoronoiDiagram destructor.
1408 */
1409 2 ~Co3NeRestrictedVoronoiDiagram() {
1410 2 clear();
1411 2 }
1412
1413
1414 /**
1415 * \brief Sets or resets exact mode for nearest neighbor search
1416 * (default is exact).
1417 * \details Nearest neighbor search can be exact or approximate.
1418 * Note that approximate mode cannot be used for the
1419 * final reconstruction phase (that needs exact combinatorics),
1420 * but it may speedup the smoothing phase.
1421 * \param[in] x if set, nearest neighbors search are exact, else they
1422 * are approximate
1423 */
1424 3 void set_exact(bool x) {
1425 3 NN_->set_exact(x);
1426 3 }
1427
1428 /**
1429 * \brief Clears this Co3NeRestrictedVoronoiDiagram.
1430 */
1431 4 void clear() {
1432 4 NN_.reset();
1433 4 nb_points_ = 0;
1434 4 p_ = nullptr;
1435 4 p_stride_ = 0;
1436 4 n_ = nullptr;
1437 4 n_stride_ = 0;
1438 4 nb_neighbors_ = 0;
1439 4 }
1440
1441 /**
1442 * \brief Initializes this Co3NeRestrictedVoronoiDiagram from a
1443 * pointset stored in a mesh.
1444 * \details If the mesh \p M has normals, then they are used.
1445 * \param[in] M the pointset
1446 */
1447 2 void init(Mesh& M) {
1448
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2 geo_assert(M.vertices.dimension() >= 3);
1449
2/10
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 2 geo_assert(M.vertices.dimension() == 3 || NN_->stride_supported());
1450 2 double* normals_pointer = nullptr;
1451 {
1452
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Attribute<double> normal;
1453
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 2 normal.bind_if_is_defined(M.vertices.attributes(), "normal");
1454
2/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 2 if(normal.is_bound() && normal.dimension() == 3) {
1455 ✗ normals_pointer = &normal[0];
1456 }
1457 2 }
1458
1459
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(normals_pointer == nullptr) {
1460 2 init(
1461 M.vertices.nb(),
1462 M.vertices.point_ptr(0), M.vertices.dimension(),
1463 nullptr, 0
1464 );
1465 } else {
1466 ✗ init(
1467 M.vertices.nb(),
1468 M.vertices.point_ptr(0), M.vertices.dimension(),
1469 normals_pointer, 3
1470 );
1471 }
1472 2 }
1473
1474 /**
1475 * \brief Initializes this Co3NeRestrictedVoronoiDiagram from an
1476 * array of points and an array of normals.
1477 * \param[in] nb_points_in number of points
1478 * \param[in] p pointer to the coordinates of the points
1479 * \param[in] p_stride number of doubles between two consecutive points
1480 * \param[in] n pointer to the coordinates of the normals
1481 * \param[in] n_stride number of doubles between two consecutive normals
1482 */
1483 5 void init(
1484 index_t nb_points_in,
1485 double* p, index_t p_stride,
1486 double* n, index_t n_stride
1487 ) {
1488 5 nb_points_ = nb_points_in;
1489 5 p_ = p;
1490 5 p_stride_ = p_stride;
1491 5 n_ = n;
1492 5 n_stride_ = n_stride;
1493 5 NN_->set_points(nb_points(), p_, p_stride_);
1494 5 }
1495
1496 /**
1497 * \brief Reconstructs the nearest neighbors search data
1498 * structure.
1499 * \details This function needs to be called whenever
1500 * the point set changes.
1501 */
1502 3 void update() {
1503 3 init(nb_points_, p_, p_stride_, n_, n_stride_);
1504 3 }
1505
1506 /**
1507 * \brief Sets the radius of the circles used to determine
1508 * points adjacencies.
1509 * \param[in] r the radius of the circles
1510 */
1511 2 void set_circles_radius(double r) {
1512 2 radius_ = r;
1513 2 sqROS_ = 4.0 * radius_ * radius_; // squared radius of security
1514 // when a neighbor is further away than ROS, then it cannot
1515 // clip a circle of radius r
1516 2 }
1517
1518 /**
1519 * \brief Gets the number of points.
1520 * \return the number of points
1521 */
1522 11333564 index_t nb_points() const {
1523 11333564 return nb_points_;
1524 }
1525
1526 /**
1527 * \brief Gets a point by its index.
1528 * \param[in] i index of the point
1529 * \return a const reference to the point
1530 */
1531 11175752 const vec3& point(index_t i) const {
1532
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 11175752 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 11175752 geo_debug_assert(i < nb_points());
1533 11175752 return *(vec3*) (p_ + i * p_stride_);
1534 // Yes I know, this is a bit ugly...
1535 }
1536
1537 /**
1538 * \brief Gets a normal by point index.
1539 * \param[in] i index of the point
1540 * \return a const reference to the normal
1541 * associated with the point
1542 */
1543 ✗ const vec3& normal(index_t i) const {
1544 ✗ geo_debug_assert(n_ != nullptr);
1545 ✗ geo_debug_assert(i < nb_points());
1546 ✗ return *(vec3*) (n_ + i * n_stride_);
1547 // Yes I know, this is a bit ugly...
1548 }
1549
1550 /**
1551 * \brief Sets the normal associated wigth a point.
1552 * \param[in] i the index of the point
1553 * \param[in] N the normal
1554 */
1555 ✗ void set_normal(index_t i, const vec3& N) const {
1556 ✗ geo_debug_assert(n_ != nullptr);
1557 ✗ geo_debug_assert(i < nb_points());
1558 ✗ double* n = n_ + i * n_stride_;
1559 ✗ n[0] = N.x;
1560 ✗ n[1] = N.y;
1561 ✗ n[2] = N.z;
1562 ✗ }
1563
1564 /**
1565 * \brief Computes the intersection between a polygon
1566 * and the halfspace determined by the bisector
1567 * of two points.
1568 * \param[in,out] Ping polygon to be clipped
1569 * \param[in,out] Pong a temporary work variable provided
1570 * by the caller
1571 * \param[in] pi first extremity of the bisector
1572 * \param[in] pj second extremity of the bisector
1573 * \param[in] j index of the second extremity of the
1574 * bisector (used to store the combinatorial information).
1575 */
1576 2063511 static void clip_polygon_by_bisector(
1577 Polygon& Ping, Polygon& Pong,
1578 const vec3& pi, const vec3& pj, index_t j
1579 ) {
1580
2/4
✓ Branch 1 taken 2063511 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2063511 times.
 2063511 if(Ping.nb_vertices() == 0) {
1581 ✗ return;
1582 }
1583
1/2
✓ Branch 1 taken 2063511 times.
✗ Branch 2 not taken.
 2063511 Pong.clear();
1584
1585 vec3 n(
1586 4127022 pi.x - pj.x,
1587 4127022 pi.y - pj.y,
1588 4127022 pi.z - pj.z
1589 2063511 );
1590
1591 // Compute d = n . m, where n is the
1592 // normal vector of the bisector [pi,pj]
1593 // and m twice the middle point of the bisector.
1594 2063511 double d =
1595 2063511 n.x * (pi.x + pj.x) +
1596 2063511 n.y * (pi.y + pj.y) +
1597 2063511 n.z * (pi.z + pj.z);
1598
1599 // The predecessor of the first vertex is the last vertex
1600
1/2
✓ Branch 1 taken 2063511 times.
✗ Branch 2 not taken.
 2063511 index_t prev_k = Ping.nb_vertices() - 1;
1601
1/2
✓ Branch 1 taken 2063511 times.
✗ Branch 2 not taken.
 2063511 const Vertex* prev_vk = &(Ping.vertex(prev_k));
1602
1603 // We compute:
1604 // prev_l = prev_vk . n
1605 2063511 double prev_l = dot(prev_vk->point(), n);
1606
1607 // We compute:
1608 // side1(pi,pj,q) = sign(2*q.n - n.m) = sign(2*l - d)
1609
1/2
✓ Branch 1 taken 2063511 times.
✗ Branch 2 not taken.
 2063511 Sign prev_status = geo_sgn(2.0 * prev_l - d);
1610
1611
3/4
✓ Branch 1 taken 12529939 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10466428 times.
✓ Branch 4 taken 2063511 times.
 12529939 for(index_t k = 0; k < Ping.nb_vertices(); k++) {
1612
1/2
✓ Branch 1 taken 10466428 times.
✗ Branch 2 not taken.
 10466428 const Vertex* vk = &(Ping.vertex(k));
1613
1614 // We compute: l = vk . n
1615 10466428 double l = dot(vk->point(), n);
1616
1617 // We compute:
1618 // side1(pi,pj,q) = sign(2*q.n - n.m) = sign(2*l - d)
1619
1/2
✓ Branch 1 taken 10466428 times.
✗ Branch 2 not taken.
 10466428 Sign status = geo_sgn(2.0 * l - d);
1620
1621 // If status of edge extremities differ,
1622 // then there is an intersection.
1623
3/4
✓ Branch 0 taken 769336 times.
✓ Branch 1 taken 9697092 times.
✓ Branch 2 taken 769336 times.
✗ Branch 3 not taken.
 10466428 if(status != prev_status && (prev_status != 0)) {
1624
1625 // Compute lambda1 and lambda2, the
1626 // barycentric coordinates of the intersection I
1627 // in the segment [prev_vk vk]
1628 // Note that d and l (used for the predicates)
1629 // are reused here.
1630 769336 double denom = 2.0 * (prev_l - l);
1631 double lambda1, lambda2;
1632
1633 // Shit happens ! [Forrest Gump]
1634
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 769336 times.
 769336 if(::fabs(denom) < 1e-20) {
1635 ✗ lambda1 = 0.5;
1636 ✗ lambda2 = 0.5;
1637 } else {
1638 769336 lambda1 = (d - 2.0 * l) / denom;
1639 // Note: lambda2 is also given
1640 // by (2.0*l2-d)/denom
1641 // (but 1.0 - lambda1 is a bit
1642 // faster to compute...)
1643 769336 lambda2 = 1.0 - lambda1;
1644 }
1645
1/2
✓ Branch 1 taken 769336 times.
✗ Branch 2 not taken.
 769336 Vertex V;
1646 1538672 V.point().x =
1647 769336 lambda1 * prev_vk->point().x + lambda2 * vk->point().x;
1648 1538672 V.point().y =
1649 769336 lambda1 * prev_vk->point().y + lambda2 * vk->point().y;
1650 1538672 V.point().z =
1651 769336 lambda1 * prev_vk->point().z + lambda2 * vk->point().z;
1652
2/2
✓ Branch 0 taken 384668 times.
✓ Branch 1 taken 384668 times.
 769336 if(status > 0) {
1653 384668 V.set_adjacent_seed(prev_vk->adjacent_seed());
1654 } else {
1655 384668 V.set_adjacent_seed(signed_index_t(j));
1656 }
1657
1/2
✓ Branch 1 taken 769336 times.
✗ Branch 2 not taken.
 769336 Pong.add_vertex(V);
1658 }
1659
2/2
✓ Branch 0 taken 9441296 times.
✓ Branch 1 taken 1025132 times.
 10466428 if(status > 0) {
1660
1/2
✓ Branch 1 taken 9441296 times.
✗ Branch 2 not taken.
 9441296 Pong.add_vertex(*vk);
1661 }
1662 10466428 prev_vk = vk;
1663 10466428 prev_status = status;
1664 10466428 prev_k = k;
1665 10466428 prev_l = l;
1666 }
1667 2063511 Ping.swap(Pong);
1668 }
1669
1670 /**
1671 * \brief Computes the squared maximum distance between a point
1672 * and the vertices of a polygon.
1673 * \param[in] p the point
1674 * \param[in] P the polygon
1675 * \return the maximum squared distance between \p p and the vertices
1676 * of \p P
1677 */
1678 2125620 static double squared_radius(const vec3& p, const Polygon& P) {
1679 2125620 double result = 0.0;
1680
3/4
✓ Branch 1 taken 12964640 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10839020 times.
✓ Branch 4 taken 2125620 times.
 12964640 for(index_t i = 0; i < P.nb_vertices(); i++) {
1681
2/4
✓ Branch 1 taken 10839020 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 10839020 times.
✗ Branch 6 not taken.
 10839020 result = std::max(result, distance2(p, P.vertex(i).point()));
1682 }
1683 2125620 return result;
1684 }
1685
1686 /**
1687 * \brief Computes a polygon that approximates a disk centered
1688 * at a point and orthogonal to its normal vector.
1689 * \param[in] i index of the point
1690 * \param[out] P an approximation of the circle centered
1691 * at point \p i with normal vector \p N. The radius is
1692 * defined by set_circles_radius().
1693 * \param[in] N normal vector
1694 */
1695 63524 void get_circle(index_t i, Polygon& P, const vec3& N) const {
1696
1/2
✓ Branch 1 taken 63524 times.
✗ Branch 2 not taken.
 63524 P.clear();
1697
1/2
✓ Branch 1 taken 63524 times.
✗ Branch 2 not taken.
 63524 const vec3& pi = point(i);
1698
1/2
✓ Branch 1 taken 63524 times.
✗ Branch 2 not taken.
 63524 vec3 U = Geom::perpendicular(N);
1699
1/2
✓ Branch 1 taken 63524 times.
✗ Branch 2 not taken.
 63524 U = normalize(U);
1700 63524 vec3 V = cross(N, U);
1701
1/2
✓ Branch 1 taken 63524 times.
✗ Branch 2 not taken.
 63524 V = normalize(V);
1702 // We use a table for sine and cosine for speeding up things
1703 // a little bit (especially on some cell phones / handheld devices
1704 // that do not have a FPU).
1705 /*
1706 const index_t nb = 10;
1707 for(index_t k=0; k<nb; ++k) {
1708 double alpha = 2.0 * M_PI * double(k) / double(nb - 1);
1709 double s = sin(alpha);
1710 double c = cos(alpha);
1711 vec3 p = pi + c * radius_ * U + s * radius_ * V;
1712 P.add_vertex(p);
1713 }
1714 */
1715
1716
2/2
✓ Branch 0 taken 635240 times.
✓ Branch 1 taken 63524 times.
 698764 for(index_t k = 0; k < sincos_nb; ++k) {
1717 635240 double s = sincos_table[k][0];
1718 635240 double c = sincos_table[k][1];
1719 635240 vec3 p = pi + c * radius_ * U + s * radius_ * V;
1720
1/2
✓ Branch 2 taken 635240 times.
✗ Branch 3 not taken.
 635240 P.add_vertex(p);
1721 }
1722
1723 63524 }
1724
1725 /**
1726 * \brief Nearest neighbor search
1727 * \param[in] i index of the query point
1728 * \param[out] neigh array of nb signed_index_t
1729 * \param[out] sq_dist array of nb doubles
1730 * \param[in] nb number of neighbors to be searched
1731 */
1732 184788 void get_neighbors(
1733 index_t i,
1734 index_t* neigh,
1735 double* sq_dist,
1736 index_t nb
1737 ) const {
1738 184788 return NN_->get_nearest_neighbors(
1739 nb, i, neigh, sq_dist
1740 184788 );
1741 }
1742
1743 /**
1744 * \brief Nearest neighbor search
1745 * \param[in] i index of the query point
1746 * \param[out] neigh vector of signed_index_t
1747 * \param[out] sq_dist array of nb doubles
1748 * \param[in] nb number of neighbors to be searched
1749 */
1750 184788 void get_neighbors(
1751 index_t i,
1752 vector<index_t>& neigh,
1753 vector<double>& sq_dist,
1754 index_t nb
1755 ) const {
1756 184788 neigh.resize(nb);
1757 184788 sq_dist.resize(nb);
1758 184788 get_neighbors(i, neigh.data(), sq_dist.data(), nb);
1759 184788 }
1760
1761 /**
1762 * \brief Computes a Restricted Voronoi Cell (RVC), i.e.
1763 * the intersection between a disk and the Voronoi cell
1764 * of a point.
1765 * \details The temporary work variables provided by the caller
1766 * make it possible to reuse memory accros multiple calls to this
1767 * function and thus avoid multiple dynamic memory allocations.
1768 * \param[in] i index of the point that determines the Voronoi cell.
1769 * \param[out] P result
1770 * \param[in] Q work temporary variable provided by caller
1771 * \param[in] neighbor work temporary variable provided by caller
1772 * \param[in] squared_dist work temporary variable provided by caller
1773 */
1774 ✗ void get_RVC(
1775 index_t i, Polygon& P,
1776 Polygon& Q,
1777 vector<index_t>& neighbor,
1778 vector<double>& squared_dist
1779 ) const {
1780 ✗ neighbor.resize(0);
1781 ✗ squared_dist.resize(0);
1782 ✗ get_RVC(i, normal(i), P, Q, neighbor, squared_dist);
1783 ✗ }
1784
1785 /**
1786 * \brief Computes a Restricted Voronoi Cell (RVC), i.e.
1787 * the intersection between a disk and the Voronoi cell
1788 * of a point.
1789 * \details The temporary work variables provided by the caller
1790 * make it possible to reuse memory accros multiple calls to this
1791 * function and thus avoid multiple dynamic memory allocations.
1792 * \param[in] i index of the point that determines the Voronoi cell.
1793 * \param[in] N normal vector at point \p i
1794 * \param[out] P result
1795 * \param[in] Q work temporary variable, provided by caller
1796 * \param[in] neighbor initial neighbor indices
1797 * if size is not zero, contains (previously computed)
1798 * neighbor indices.
1799 * \param[in] squared_dist initial neighbor squared distances
1800 * if size is not zero, contains (previously computed)
1801 * neighbor squared distances.
1802 */
1803 63524 void get_RVC(
1804 index_t i, const vec3& N, Polygon& P,
1805 Polygon& Q,
1806 vector<index_t>& neighbor,
1807 vector<double>& squared_dist
1808 ) const {
1809
1/2
✓ Branch 1 taken 63524 times.
✗ Branch 2 not taken.
 63524 get_circle(i, P, N);
1810
1811 63524 index_t nb_neigh = std::min(index_t(nb_points() - 1), index_t(20));
1812 63524 index_t jj = 0;
1813
1814 // just in case, limit to 1000 neighbors.
1815 63524 index_t max_neigh = std::min(index_t(1000), nb_points() - 1);
1816
1817
2/2
✓ Branch 0 taken 92844 times.
✓ Branch 1 taken 1415 times.
 94259 while(nb_neigh < max_neigh) {
1818
2/4
✓ Branch 1 taken 92844 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 92844 times.
 92844 if(P.nb_vertices() < 3) {
1819 62109 return;
1820 }
1821
2/2
✓ Branch 1 taken 25978 times.
✓ Branch 2 taken 66866 times.
 92844 if(neighbor.size() < nb_neigh) {
1822
1/2
✓ Branch 1 taken 25978 times.
✗ Branch 2 not taken.
 25978 get_neighbors(i, neighbor, squared_dist, nb_neigh);
1823 }
1824
6/8
✓ Branch 0 taken 156368 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 156368 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 63524 times.
✓ Branch 6 taken 92844 times.
✓ Branch 7 taken 63524 times.
✓ Branch 8 taken 92844 times.
 156368 while(jj < nb_neigh && squared_dist[jj] < 1e-30) {
1825 63524 jj++;
1826 }
1827
2/2
✓ Branch 0 taken 2125620 times.
✓ Branch 1 taken 30735 times.
 2156355 while(jj < nb_neigh) {
1828
2/4
✓ Branch 1 taken 2125620 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2125620 times.
 2125620 if(squared_dist[jj] > sqROS_) {
1829 ✗ return;
1830 }
1831
1/2
✓ Branch 1 taken 2125620 times.
✗ Branch 2 not taken.
 2125620 index_t j = neighbor[jj];
1832
2/4
✓ Branch 1 taken 2125620 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2125620 times.
✗ Branch 5 not taken.
 2125620 double Rk = squared_radius(point(i), P);
1833
3/4
✓ Branch 1 taken 2125620 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 62109 times.
✓ Branch 4 taken 2063511 times.
 2125620 if(squared_dist[jj] > 4.0 * Rk) {
1834 62109 return;
1835 }
1836
3/6
✓ Branch 1 taken 2063511 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2063511 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2063511 times.
✗ Branch 8 not taken.
 2063511 clip_polygon_by_bisector(P, Q, point(i), point(j), j);
1837 2063511 jj++;
1838 }
1839
1/2
✓ Branch 0 taken 30735 times.
✗ Branch 1 not taken.
 30735 if(nb_neigh > 3) {
1840 30735 nb_neigh += nb_neigh / 3;
1841 } else {
1842 ✗ nb_neigh++;
1843 }
1844 30735 nb_neigh = std::min(nb_neigh, nb_points()-1);
1845 }
1846 }
1847
1848 /**
1849 * \brief Gets the number of neighbors, used for nearest neighbors
1850 * queries.
1851 * \return the number of neighbors
1852 */
1853 20 index_t nb_neighbors() const {
1854 20 return std::min(nb_neighbors_,nb_points()-1);
1855 }
1856
1857 /**
1858 * \brief Sets the number of neighbors, used for nearest neighbors
1859 * queries.
1860 * \param[in] x the number of neighbors
1861 */
1862 5 void set_nb_neighbors(index_t x) {
1863 5 nb_neighbors_ = x;
1864 5 }
1865
1866 private:
1867 friend class Co3Ne;
1868
1869 index_t nb_points_;
1870 double* p_;
1871 index_t p_stride_;
1872 double* n_;
1873 index_t n_stride_;
1874 double radius_;
1875
1876 NearestNeighborSearch_var NN_;
1877
1878 double sqROS_;
1879 index_t nb_neighbors_;
1880 };
1881
1882 /************************************************************************/
1883
1884 class Co3Ne;
1885
1886 /**
1887 * \brief Determines what a thread does in
1888 * the multithreaded implementation of the Co3Ne reconstruction algorithm.
1889 */
1890 enum Co3NeMode {
1891 CO3NE_NONE, /**< uninitialized */
1892 CO3NE_NORMALS, /**< estimate normals in pointset */
1893 CO3NE_SMOOTH, /**< smooth the pointset */
1894 CO3NE_RECONSTRUCT, /**< reconstruct the triangles */
1895 CO3NE_NORMALS_AND_RECONSTRUCT
1896 /**< combined normal estimation and reconstruction */
1897 };
1898
1899 /**
1900 * \brief Internal implementation class for Co3Ne.
1901 * Encapsulates the work of one thread.
1902 */
1903 class Co3NeThread : public Thread {
1904 public:
1905 /**
1906 * \brief Creates a new Co3NeThread
1907 * \param[in] master the Co3Ne this thread depends on
1908 * \param[in] from index of the first point to process
1909 * \param[in] to one position past the index of the last point
1910 */
1911 8 Co3NeThread(
1912 Co3Ne* master,
1913 index_t from, index_t to
1914 8 ) :
1915 8 master_(master),
1916 8 from_(from),
1917
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 to_(to) {
1918 8 mode_ = CO3NE_NONE;
1919 8 }
1920
1921 /**
1922 * \brief Sets the mode of this thread
1923 * \param[in] m the mode, that determines whether normal computation,
1924 * smoothing or reconstruction is performed
1925 */
1926 20 void set_mode(Co3NeMode m) {
1927 20 mode_ = m;
1928 20 }
1929
1930 /**
1931 * \brief Does the actual computation of this thread.
1932 * \details The actual computation is determined by set_mode().
1933 */
1934 20 void run() override {
1935
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 20 switch(mode_) {
1936 ✗ case CO3NE_NORMALS:
1937 ✗ run_normals();
1938 ✗ break;
1939 12 case CO3NE_SMOOTH:
1940 12 run_smooth();
1941 12 break;
1942 ✗ case CO3NE_RECONSTRUCT:
1943 ✗ run_reconstruct();
1944 ✗ break;
1945 8 case CO3NE_NORMALS_AND_RECONSTRUCT:
1946 8 run_normals_and_reconstruct();
1947 8 break;
1948 ✗ case CO3NE_NONE:
1949 ✗ break;
1950 }
1951 20 }
1952
1953 /**
1954 * \brief Gets the reconstructed triangles.
1955 * \return a reference to a vector of indices
1956 */
1957 24 vector<index_t>& triangles() {
1958 24 return triangles_;
1959 }
1960
1961
1962 /**
1963 * \brief Gets the number of reconstructed triangles.
1964 * \return the number of reconstructed triangles
1965 */
1966 8 index_t nb_triangles() const {
1967 8 return triangles_.size()/3;
1968 }
1969
1970 protected:
1971 /**
1972 * \brief Estimates the normals in the pointset.
1973 */
1974 void run_normals();
1975
1976 /**
1977 * \brief Smoothes the pointset.
1978 */
1979 void run_smooth();
1980
1981 /**
1982 * \brief Reconstructs the triangles.
1983 */
1984 void run_reconstruct();
1985
1986 /**
1987 * \brief Estimates the normals and reconstructs the triangles.
1988 */
1989 void run_normals_and_reconstruct();
1990
1991 private:
1992 Co3Ne* master_;
1993 index_t from_;
1994 index_t to_;
1995 Co3NeMode mode_;
1996 PrincipalAxes3d least_squares_normal_;
1997 vector<index_t> triangles_;
1998 };
1999
2000 /************************************************************************/
2001
2002 /**
2003 * \brief Reconstructs a mesh from a set of point with
2004 * the Co3Ne algorithm (concurrent co-cones).
2005 * This class also implements point set smoothing and point
2006 * set normal estimation.
2007 */
2008 class Co3Ne {
2009 public:
2010 /**
2011 * \brief Constructs a new Co3Ne.
2012 * \param[in] M the pointset
2013 */
2014 2 Co3Ne(Mesh& M) :
2015
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 mesh_(M) {
2016 // TODO: interlace threads (more cache friendly)
2017
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 RVD_.init(mesh_);
2018
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 index_t nb = Process::maximum_concurrent_threads();
2019 2 thread_.clear();
2020 2 index_t batch_size = RVD_.nb_points() / nb;
2021 2 index_t cur = 0;
2022 2 index_t remaining = RVD_.nb_points();
2023
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
 10 for(index_t i = 0; i < nb; i++) {
2024 8 index_t this_batch_size = batch_size;
2025
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
 8 if(i == nb - 1) {
2026 2 this_batch_size = remaining;
2027 }
2028
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 thread_.push_back(
2029
1/4
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 16 new Co3NeThread(
2030 this, cur, cur + this_batch_size
2031
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
 8 )
2032 );
2033 8 cur += this_batch_size;
2034 8 remaining -= this_batch_size;
2035 }
2036
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2 geo_assert(remaining == 0);
2037
2038 // TODO: pass it as an argument and let Vorpaline's main.cpp
2039 // communicate with CmdLine.
2040
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 double alpha = CmdLine::get_arg_double("co3ne:max_N_angle");
2041 2 alpha = alpha * M_PI / 180.0;
2042 2 set_max_angle(alpha);
2043 2 }
2044
2045 /**
2046 * \brief Runs the threads.
2047 */
2048 5 void run_threads() {
2049 5 Process::run_threads(thread_);
2050 5 }
2051
2052 /**
2053 * \brief Estimates the normals of the point set.
2054 * \details They are stored in the "normal" vertex attribute.
2055 * \param[in] nb_neighbors number of neighbors to be
2056 * used for normal estimation
2057 */
2058 ✗ void compute_normals(index_t nb_neighbors) {
2059 ✗ Attribute<double> normals;
2060 ✗ normals.bind_if_is_defined(mesh_.vertices.attributes(), "normal");
2061 ✗ if(!normals.is_bound()) {
2062 ✗ normals.create_vector_attribute(
2063 ✗ mesh_.vertices.attributes(), "normal", 3
2064 );
2065 }
2066 ✗ RVD_.init(mesh_);
2067 ✗ RVD_.set_nb_neighbors(nb_neighbors);
2068 ✗ for(index_t t = 0; t < thread_.size(); t++) {
2069 ✗ thread_[t]->set_mode(CO3NE_NORMALS);
2070 }
2071 ✗ run_threads();
2072 ✗ }
2073
2074 ✗ static inline double cos_angle(
2075 Attribute<double>& normal, index_t v1, index_t v2
2076 ) {
2077 ✗ vec3 V1(normal[3*v1],normal[3*v1+1],normal[3*v1+2]);
2078 ✗ vec3 V2(normal[3*v2],normal[3*v2+1],normal[3*v2+2]);
2079 ✗ return Geom::cos_angle(V1,V2);
2080 }
2081
2082 ✗ static inline void flip(Attribute<double>& normal, index_t v) {
2083 ✗ normal[3*v] = -normal[3*v];
2084 ✗ normal[3*v+1] = -normal[3*v+1];
2085 ✗ normal[3*v+2] = -normal[3*v+2];
2086 ✗ }
2087
2088 /**
2089 * \brief Tentatively enforces a coherent orientation of normals
2090 * using a breadth-first traveral of the K-nearest-neighbor graph.
2091 * \retval true if normals where computed
2092 * \retval false otherwise (when the user pushes the cancel button).
2093 */
2094 ✗ bool reorient_normals() {
2095 ✗ Attribute<double> normal;
2096 ✗ normal.bind_if_is_defined(mesh_.vertices.attributes(), "normal");
2097 ✗ geo_assert(normal.is_bound());
2098
2099 // To resist noisy inputs, propagation is prioritized to the points
2100 // that have smallest normal deviations.
2101
2102 ✗ std::priority_queue<OrientNormal> S;
2103 ✗ vector<index_t> neighbors(RVD_.nb_neighbors());
2104 ✗ vector<double> dist(RVD_.nb_neighbors());
2105
2106 ✗ index_t nb=0;
2107 ✗ ProgressTask progress("Reorient");
2108
2109 try {
2110 ✗ std::vector<bool> visited(mesh_.vertices.nb(), false);
2111 ✗ for(index_t v=0; v<mesh_.vertices.nb(); ++v) {
2112 ✗ if(!visited[v]) {
2113 ✗ S.push(OrientNormal(v,0.0));
2114 ✗ visited[v] = true;
2115 ✗ while(!S.empty()) {
2116 ✗ OrientNormal top = S.top();
2117 ✗ ++nb;
2118 ✗ progress.progress(nb*100/mesh_.vertices.nb());
2119 ✗ S.pop();
2120 ✗ if(top.dot < 0.0) {
2121 ✗ flip(normal,top.v);
2122 }
2123 ✗ RVD_.get_neighbors(
2124 top.v,
2125 neighbors.data(),dist.data(),RVD_.nb_neighbors()
2126 );
2127 ✗ for(index_t i=0; i<RVD_.nb_neighbors(); ++i) {
2128 ✗ index_t neigh = neighbors[i];
2129 ✗ if(!visited[neigh]) {
2130 ✗ visited[neigh] = true;
2131 double dot =
2132 ✗ cos_angle(normal, top.v, neigh);
2133 ✗ S.push(OrientNormal(neigh,dot));
2134 }
2135 }
2136 }
2137 }
2138 }
2139 ✗ } catch(const TaskCanceled&) {
2140 ✗ return false;
2141 ✗ }
2142 ✗ return true;
2143 ✗ }
2144
2145 /**
2146 * \brief Smoothes a point set by projection
2147 * onto the nearest neighbors best
2148 * approximating planes.
2149 * \param[in] nb_neighbors number of neighbors to be
2150 * used for best approximating plane estimation
2151 */
2152 3 void smooth(index_t nb_neighbors) {
2153 3 new_vertices_.resize(mesh_.vertices.nb() * 3);
2154 3 RVD_.set_nb_neighbors(nb_neighbors);
2155
2/2
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 3 times.
 15 for(index_t t = 0; t < thread_.size(); t++) {
2156 12 thread_[t]->set_mode(CO3NE_SMOOTH);
2157 }
2158 3 run_threads();
2159 /*
2160 // TODO: once 'steal-arg' mode works for vertices,
2161 // we can use this one.
2162 if(RVD_.p_stride_ == 3) {
2163 MeshMutator::vertices(mesh_).swap(new_vertices_);
2164 } else */ {
2165 3 index_t idx = 0;
2166
2/2
✓ Branch 1 taken 95286 times.
✓ Branch 2 taken 3 times.
 95289 for(index_t i = 0; i < mesh_.vertices.nb(); i++) {
2167 95286 double* p = mesh_.vertices.point_ptr(i);
2168
2/2
✓ Branch 0 taken 285858 times.
✓ Branch 1 taken 95286 times.
 381144 for(coord_index_t c = 0; c < 3; c++) {
2169 285858 p[c] = new_vertices_[idx];
2170 285858 idx++;
2171 }
2172 }
2173 }
2174 3 }
2175
2176 /**
2177 * \brief This function needs to be called after the
2178 * last iteration of smoothing.
2179 * \details Deallocates the temporary
2180 * variables used for smoothing.
2181 */
2182 1 void end_smooth() {
2183 1 new_vertices_.clear();
2184 1 }
2185
2186 /**
2187 * \brief Reconstructs a mesh from a point set.
2188 * \details If the mesh has a "normal" vertex attribute,
2189 * then the existing normals are used, else normals are estimated.
2190 * \param[in] r maximum distance used to determine
2191 * points adjacencies.
2192 */
2193 2 void reconstruct(double r) {
2194 2 bool has_normals = false;
2195 {
2196
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Attribute<double> normal;
2197
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 2 normal.bind_if_is_defined(mesh_.vertices.attributes(),"normal");
2198 2 has_normals = (
2199
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2 normal.is_bound() && normal.dimension() == 3
2200 );
2201 2 }
2202
2203
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 ProgressTask progress("reconstruct",100);
2204
2205
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(has_normals) {
2206 ✗ Stopwatch W("Co3Ne recons");
2207 ✗ RVD_.set_circles_radius(r);
2208 ✗ for(index_t t = 0; t < thread_.size(); t++) {
2209 ✗ thread_[t]->set_mode(CO3NE_RECONSTRUCT);
2210 ✗ thread_[t]->triangles().clear();
2211 }
2212 ✗ progress.progress(1);
2213 ✗ run_threads();
2214 ✗ progress.progress(50);
2215 ✗ } else {
2216
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 Stopwatch W("Co3Ne recons");
2217
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 Logger::out("Co3Ne")
2218
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << "using combined \'normals and reconstruct\'"
2219
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << std::endl;
2220
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 RVD_.set_nb_neighbors(
2221
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 CmdLine::get_arg_uint("co3ne:nb_neighbors")
2222 );
2223 2 RVD_.set_circles_radius(r);
2224
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
 10 for(index_t t = 0; t < thread_.size(); t++) {
2225
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 thread_[t]->set_mode(CO3NE_NORMALS_AND_RECONSTRUCT);
2226
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 thread_[t]->triangles().clear();
2227 }
2228
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 progress.progress(1);
2229
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 run_threads();
2230
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 progress.progress(50);
2231 2 }
2232
2233 {
2234
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 Stopwatch W("Co3Ne manif.");
2235
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 RVD_.clear(); // reclaim memory used by ANN
2236
2237 2 index_t nb_triangles = 0;
2238
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
 10 for(index_t t = 0; t < thread_.size(); t++) {
2239
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
 8 nb_triangles += thread_[t]->nb_triangles();
2240 }
2241
2242
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 Logger::out("Co3Ne") << "Raw triangles: "
2243
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << nb_triangles
2244
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << std::endl;
2245
2246 2 vector<index_t> raw_triangles;
2247
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 raw_triangles.reserve(nb_triangles * 3);
2248
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
 10 for(index_t th = 0; th < thread_.size(); th++) {
2249
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 vector<index_t>& triangles = thread_[th]->triangles();
2250
1/2
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
 16 raw_triangles.insert(
2251 8 raw_triangles.end(),
2252 triangles.begin(), triangles.end()
2253 );
2254
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 thread_[th]->triangles().clear();
2255 }
2256
2257
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
 4 if(CmdLine::get_arg_bool("dbg:co3ne")) {
2258 ✗ Logger::out("Co3Ne") << ">> co3ne_raw.geogram"
2259 ✗ << std::endl;
2260 ✗ Mesh M;
2261 ✗ M.vertices.assign_points(
2262 ✗ mesh_.vertices.point_ptr(0),
2263 ✗ mesh_.vertices.dimension(),
2264 ✗ mesh_.vertices.nb()
2265 );
2266 ✗ M.facets.assign_triangle_mesh(raw_triangles, false);
2267 ✗ M.vertices.set_dimension(3);
2268 ✗ mesh_save(M, "co3ne_raw.geogram");
2269 ✗ }
2270
2271 2 vector<index_t> good_triangles;
2272 2 vector<index_t> not_so_good_triangles;
2273
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 co3ne_split_triangles_list(
2274 raw_triangles, good_triangles, not_so_good_triangles
2275 );
2276
2277
2278
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
 4 if(CmdLine::get_arg_bool("dbg:co3ne")) {
2279 ✗ Logger::out("Co3Ne") << ">> co3ne_T3.geogram"
2280 ✗ << std::endl;
2281 ✗ Mesh M;
2282 ✗ M.vertices.assign_points(
2283 ✗ mesh_.vertices.point_ptr(0),
2284 ✗ mesh_.vertices.dimension(),
2285 ✗ mesh_.vertices.nb()
2286 );
2287 ✗ M.facets.assign_triangle_mesh(good_triangles, false);
2288 ✗ M.vertices.set_dimension(3);
2289 ✗ mesh_save(M, "co3ne_T3.geogram");
2290 ✗ }
2291
2292
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
 4 if(CmdLine::get_arg_bool("dbg:co3ne")) {
2293 ✗ Logger::out("Co3Ne") << ">> co3ne_T12.geogram"
2294 ✗ << std::endl;
2295 ✗ Mesh M;
2296 ✗ M.vertices.assign_points(
2297 ✗ mesh_.vertices.point_ptr(0),
2298 ✗ mesh_.vertices.dimension(),
2299 ✗ mesh_.vertices.nb()
2300 );
2301 ✗ M.facets.assign_triangle_mesh(not_so_good_triangles, false);
2302 ✗ M.vertices.set_dimension(3);
2303 ✗ mesh_save(M, "co3ne_T12.geogram");
2304 ✗ }
2305
2306
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 progress.progress(53);
2307
2308 Co3NeManifoldExtraction manifold_extraction(
2309 mesh_, good_triangles
2310
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 );
2311
2312
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 progress.progress(55);
2313
2314
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 if(CmdLine::get_arg_bool("co3ne:T12")) {
2315
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 manifold_extraction.add_triangles(not_so_good_triangles);
2316 }
2317
2318
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 progress.progress(57);
2319
2320
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 mesh_reorient(mesh_);
2321
2322
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 progress.progress(60);
2323
2324
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
 4 if(CmdLine::get_arg_bool("dbg:co3ne")) {
2325 ✗ Logger::out("Co3Ne") << ">> co3ne_manif.geogram"
2326 ✗ << std::endl;
2327 ✗ mesh_save(mesh_, "co3ne_manif.geogram");
2328 }
2329 2 }
2330
2331
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 if(CmdLine::get_arg_bool("co3ne:repair")) {
2332
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 Stopwatch W("Co3Ne post.");
2333
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 mesh_repair(mesh_,
2334 MeshRepairMode(
2335 MESH_REPAIR_DEFAULT | MESH_REPAIR_RECONSTRUCT
2336 )
2337 );
2338
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
 4 if(CmdLine::get_arg_bool("dbg:co3ne")) {
2339 ✗ Logger::out("Co3Ne") << ">> co3ne_post.geogram"
2340 ✗ << std::endl;
2341 ✗ mesh_save(mesh_, "co3ne_post.geogram");
2342 }
2343 2 }
2344
2345
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 progress.progress(100);
2346
2347
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 Logger::out("Topology")
2348
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.
 2 << "nb components=" << mesh_nb_connected_components(mesh_)
2349
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.
 2 << " nb borders=" << mesh_nb_borders(mesh_)
2350
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << std::endl;
2351
2352 2 }
2353
2354 /**
2355 * \brief Gets the Co3NeRestrictedVoronoiDiagram associated
2356 * with this Co3Ne.
2357 * \return a reference to the Co3NeRestrictedVoronoiDiagram
2358 */
2359 26 Co3NeRestrictedVoronoiDiagram& RVD() {
2360 26 return RVD_;
2361 }
2362
2363 /**
2364 * \brief Sets a point
2365 * \param[in] i the index of the point
2366 * \param[in] P the new geometric location at the point
2367 */
2368 95286 void set_point(index_t i, const vec3& P) {
2369
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 95286 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 95286 geo_debug_assert(new_vertices_.size() > 3 * i + 2);
2370 95286 new_vertices_[3 * i] = P.x;
2371 95286 new_vertices_[3 * i + 1] = P.y;
2372 95286 new_vertices_[3 * i + 2] = P.z;
2373 95286 }
2374
2375 /**
2376 * \brief Sets a normal vector
2377 * \param[in] i the index of the point
2378 * \param[in] N the new normal vector associated with the point
2379 */
2380 ✗ void set_normal(index_t i, const vec3& N) {
2381 ✗ RVD_.set_normal(i, N);
2382 ✗ }
2383
2384 /**
2385 * \brief Sets the maximum angle for determining admissible triangles.
2386 * \details Admissible triangles have a deviation between their normals
2387 * and the normals estimated in the pointset smaller than a given
2388 * threshold \p alpha.
2389 * \param[in] alpha the maximum normal angle deviation
2390 */
2391 2 void set_max_angle(double alpha) {
2392 2 min_cos_angle_ = ::cos(alpha);
2393 2 }
2394
2395 8 Mesh& mesh() {
2396 8 return mesh_;
2397 }
2398
2399 private:
2400 Mesh& mesh_;
2401 vector<double> new_vertices_;
2402 Co3NeRestrictedVoronoiDiagram RVD_;
2403 TypedThreadGroup<Co3NeThread> thread_;
2404 double min_cos_angle_;
2405 };
2406
2407 /************************************************************************/
2408
2409 ✗ void Co3NeThread::run_normals() {
2410 ✗ Co3NeRestrictedVoronoiDiagram& RVD = master_->RVD();
2411 ✗ index_t nb_neigh = RVD.nb_neighbors();
2412 ✗ vector<index_t> neigh(nb_neigh);
2413 ✗ vector<double> sq_dist(nb_neigh);
2414
2415 ✗ for(index_t i = from_; i < to_; i++) {
2416 ✗ RVD.get_neighbors(
2417 i, neigh, sq_dist, nb_neigh
2418 );
2419 ✗ least_squares_normal_.begin();
2420 ✗ for(index_t jj = 0; jj < neigh.size(); jj++) {
2421 ✗ least_squares_normal_.add_point(RVD.point(neigh[jj]));
2422 }
2423 ✗ least_squares_normal_.end();
2424 ✗ master_->set_normal(i, least_squares_normal_.normal());
2425 }
2426 ✗ }
2427
2428 12 void Co3NeThread::run_smooth() {
2429 12 Co3NeRestrictedVoronoiDiagram& RVD = master_->RVD();
2430
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 index_t nb_neigh = RVD.nb_neighbors();
2431
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 vector<index_t> neigh(nb_neigh);
2432
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 vector<double> sq_dist(nb_neigh);
2433
2434
2/2
✓ Branch 0 taken 95286 times.
✓ Branch 1 taken 12 times.
 95298 for(index_t i = from_; i < to_; i++) {
2435
1/2
✓ Branch 1 taken 95286 times.
✗ Branch 2 not taken.
 95286 RVD.get_neighbors(
2436 i, neigh, sq_dist, nb_neigh
2437 );
2438
1/2
✓ Branch 1 taken 95286 times.
✗ Branch 2 not taken.
 95286 least_squares_normal_.begin();
2439
2/2
✓ Branch 1 taken 2858580 times.
✓ Branch 2 taken 95286 times.
 2953866 for(index_t jj = 0; jj < neigh.size(); jj++) {
2440
3/6
✓ Branch 1 taken 2858580 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2858580 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2858580 times.
✗ Branch 8 not taken.
 2858580 least_squares_normal_.add_point(RVD.point(neigh[jj]));
2441 }
2442
1/2
✓ Branch 1 taken 95286 times.
✗ Branch 2 not taken.
 95286 least_squares_normal_.end();
2443
1/2
✓ Branch 2 taken 95286 times.
✗ Branch 3 not taken.
 95286 vec3 N = normalize(least_squares_normal_.normal());
2444
1/2
✓ Branch 1 taken 95286 times.
✗ Branch 2 not taken.
 95286 vec3 g = least_squares_normal_.center();
2445
1/2
✓ Branch 1 taken 95286 times.
✗ Branch 2 not taken.
 95286 vec3 d = RVD.point(i) - g;
2446 95286 d -= dot(d, N) * N;
2447
1/2
✓ Branch 2 taken 95286 times.
✗ Branch 3 not taken.
 95286 master_->set_point(i, g + d);
2448 }
2449 12 }
2450
2451 ✗ void Co3NeThread::run_reconstruct() {
2452 ✗ Co3NeRestrictedVoronoiDiagram& RVD = master_->RVD();
2453 ✗ vector<index_t> neigh(100);
2454 ✗ vector<double> sq_dist(100);
2455 ✗ Co3NeRestrictedVoronoiDiagram::Polygon P(100);
2456 ✗ Co3NeRestrictedVoronoiDiagram::Polygon Q(100);
2457
2458 ✗ for(index_t i = from_; i < to_; i++) {
2459 ✗ RVD.get_RVC(i, P, Q, neigh, sq_dist);
2460 ✗ for(index_t v1 = 0; v1 < P.nb_vertices(); v1++) {
2461 ✗ index_t v2 = P.next_vertex(v1);
2462 ✗ signed_index_t j = P.vertex(v1).adjacent_seed();
2463 ✗ signed_index_t k = P.vertex(v2).adjacent_seed();
2464 ✗ if(
2465 ✗ j >= 0 && k >= 0 && j != k
2466 ) {
2467 ✗ triangles_.push_back(i);
2468 ✗ triangles_.push_back(index_t(j));
2469 ✗ triangles_.push_back(index_t(k));
2470 }
2471 }
2472 }
2473 ✗ }
2474
2475 8 void Co3NeThread::run_normals_and_reconstruct() {
2476
2477
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 Attribute<double> normal;
2478
3/6
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
 16 if(CmdLine::get_arg_bool("co3ne:use_normals")) {
2479
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 normal.bind_if_is_defined(
2480
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 24 master_->mesh().vertices.attributes(), "normal"
2481 );
2482 }
2483
2484
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 std::ofstream RVD_file;
2485 8 bool debug_RVD = false;
2486 8 if(
2487
3/6
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 8 times.
 16 CmdLine::get_arg_bool("dbg:co3neRVD")
2488 ) {
2489 ✗ if(CmdLine::get_arg_bool("sys:multithread")) {
2490 ✗ Logger::warn("Co3Ne")
2491 ✗ << "dbg:Co3NeRVD cannot work in multithread mode"
2492 ✗ << std::endl;
2493 ✗ Logger::warn("Co3Ne")
2494 ✗ << "use sys:multithread=false"
2495 ✗ << std::endl;
2496 } else {
2497 ✗ Logger::out("Co3Ne") << "Saving RVD in co3neRVD.obj"
2498 ✗ << std::endl;
2499 ✗ RVD_file.open("co3neRVD.obj");
2500 ✗ debug_RVD=true;
2501 }
2502 }
2503 8 index_t cur_v = 0;
2504
2505 8 Co3NeRestrictedVoronoiDiagram& RVD = master_->RVD();
2506
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 index_t nb_neigh = RVD.nb_neighbors();
2507
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 vector<index_t> neigh(100);
2508
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 vector<double> sq_dist(100);
2509
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 Co3NeRestrictedVoronoiDiagram::Polygon P(100);
2510
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 Co3NeRestrictedVoronoiDiagram::Polygon Q(100);
2511
2512
2/2
✓ Branch 0 taken 63524 times.
✓ Branch 1 taken 8 times.
 63532 for(index_t i = from_; i < to_; i++) {
2513
2514 63524 vec3 N;
2515
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 63524 times.
 63524 if(normal.is_bound()) {
2516 ✗ RVD.get_neighbors(
2517 i, neigh, sq_dist, nb_neigh
2518 );
2519 ✗ N = vec3(normal[3*i], normal[3*i+1], normal[3*i+2]);
2520 } else {
2521
1/2
✓ Branch 1 taken 63524 times.
✗ Branch 2 not taken.
 63524 RVD.get_neighbors(
2522 i, neigh, sq_dist, nb_neigh
2523 );
2524
1/2
✓ Branch 1 taken 63524 times.
✗ Branch 2 not taken.
 63524 least_squares_normal_.begin();
2525
2/2
✓ Branch 1 taken 1905720 times.
✓ Branch 2 taken 63524 times.
 1969244 for(index_t jj = 0; jj < neigh.size(); jj++) {
2526
3/6
✓ Branch 1 taken 1905720 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1905720 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1905720 times.
✗ Branch 8 not taken.
 1905720 least_squares_normal_.add_point(RVD.point(neigh[jj]));
2527 }
2528
1/2
✓ Branch 1 taken 63524 times.
✗ Branch 2 not taken.
 63524 least_squares_normal_.end();
2529 63524 N = least_squares_normal_.normal();
2530 }
2531
2532
1/2
✓ Branch 1 taken 63524 times.
✗ Branch 2 not taken.
 63524 RVD.get_RVC(i, N, P, Q, neigh, sq_dist);
2533
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 63524 times.
 63524 if(debug_RVD) {
2534 ✗ for(index_t v = 0; v < P.nb_vertices(); ++v) {
2535 ✗ RVD_file << "v "
2536 ✗ << P.vertex(v).point().x
2537 ✗ << " "
2538 ✗ << P.vertex(v).point().y
2539 ✗ << " "
2540 ✗ << P.vertex(v).point().z
2541 ✗ << std::endl;
2542 }
2543 ✗ RVD_file << "f ";
2544 ✗ for(index_t v = 0; v < P.nb_vertices(); ++v) {
2545 ✗ ++cur_v;
2546 ✗ RVD_file << cur_v << " ";
2547 }
2548 ✗ RVD_file << std::endl;
2549 ✗ RVD_file << "#" << i << " ";
2550 ✗ for(index_t v1 = 0; v1 < P.nb_vertices(); ++v1) {
2551 ✗ RVD_file << P.vertex(v1).adjacent_seed() << " ";
2552 }
2553 ✗ RVD_file << std::endl;
2554 }
2555
3/4
✓ Branch 1 taken 442968 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 379444 times.
✓ Branch 4 taken 63524 times.
 442968 for(index_t v1 = 0; v1 < P.nb_vertices(); v1++) {
2556
1/2
✓ Branch 1 taken 379444 times.
✗ Branch 2 not taken.
 379444 index_t v2 = P.next_vertex(v1);
2557
1/2
✓ Branch 1 taken 379444 times.
✗ Branch 2 not taken.
 379444 signed_index_t j = P.vertex(v1).adjacent_seed();
2558
1/2
✓ Branch 1 taken 379444 times.
✗ Branch 2 not taken.
 379444 signed_index_t k = P.vertex(v2).adjacent_seed();
2559
2/2
✓ Branch 0 taken 377016 times.
✓ Branch 1 taken 2428 times.
 379444 if(
2560
3/4
✓ Branch 0 taken 374958 times.
✓ Branch 1 taken 2058 times.
✓ Branch 2 taken 374958 times.
✗ Branch 3 not taken.
 377016 j >= 0 && k >= 0 && j != k
2561 ) {
2562
1/2
✓ Branch 1 taken 374958 times.
✗ Branch 2 not taken.
 374958 triangles_.push_back(i);
2563
1/2
✓ Branch 1 taken 374958 times.
✗ Branch 2 not taken.
 374958 triangles_.push_back(index_t(j));
2564
1/2
✓ Branch 1 taken 374958 times.
✗ Branch 2 not taken.
 374958 triangles_.push_back(index_t(k));
2565 }
2566 }
2567 }
2568
2569
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
 8 if(normal.is_bound()) {
2570 ✗ normal.unbind();
2571 }
2572 8 }
2573 }
2574
2575 /****************************************************************************/
2576
2577 namespace GEO {
2578
2579 ✗ void Co3Ne_smooth(Mesh& M, index_t nb_neighbors, index_t nb_iterations) {
2580 ✗ Co3Ne co3ne(M);
2581 try {
2582 ✗ ProgressTask progress("Smoothing", nb_iterations);
2583 ✗ for(index_t i = 0; i < nb_iterations; i++) {
2584 ✗ co3ne.smooth(nb_neighbors);
2585 ✗ if(i != nb_iterations - 1) {
2586 ✗ co3ne.RVD().update();
2587 }
2588 ✗ progress.next();
2589 }
2590 ✗ co3ne.end_smooth();
2591 ✗ }
2592 ✗ catch(const TaskCanceled&) {
2593 ✗ }
2594 ✗ }
2595
2596 ✗ bool Co3Ne_compute_normals(Mesh& M, index_t nb_neighbors, bool reorient) {
2597 {
2598 ✗ Attribute<double> normal;
2599 ✗ normal.bind_if_is_defined(M.vertices.attributes(), "normal");
2600 ✗ if(!normal.is_bound()) {
2601 ✗ normal.create_vector_attribute(
2602 ✗ M.vertices.attributes(), "normal", 3
2603 );
2604 }
2605 ✗ }
2606 ✗ Co3Ne co3ne(M);
2607 ✗ Logger::out("Co3Ne") << "Computing normals" << std::endl;
2608 ✗ co3ne.compute_normals(nb_neighbors);
2609 ✗ if(reorient) {
2610 ✗ Logger::out("Co3Ne") << "Orienting normals" << std::endl;
2611 ✗ if(!co3ne.reorient_normals()) {
2612 ✗ return false;
2613 }
2614 }
2615 ✗ return true;
2616 ✗ }
2617
2618 ✗ void Co3Ne_reconstruct(Mesh& M, double radius) {
2619 ✗ Co3Ne co3ne(M);
2620 ✗ co3ne.reconstruct(radius);
2621 ✗ }
2622
2623 2 void Co3Ne_smooth_and_reconstruct(
2624 Mesh& M, index_t nb_neighbors, index_t nb_iterations, double radius
2625 ) {
2626
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 Stopwatch W("Co3Ne total");
2627
2628
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 if(CmdLine::get_arg_bool("co3ne:use_normals")) {
2629
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Attribute<double> normal;
2630
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 2 normal.bind_if_is_defined(M.vertices.attributes(), "normal");
2631
2/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 2 if(normal.is_bound() && normal.dimension() == 3) {
2632 ✗ Logger::out("Co3Ne") << "Using existing normal attribute"
2633 ✗ << std::endl;
2634 } else {
2635
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.
 6 Logger::out("Co3Ne") << "No \'normal\' vertex attribute found"
2636
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << std::endl;
2637
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.
 6 Logger::out("Co3Ne") << "(estimating normals)"
2638
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << std::endl;
2639 }
2640 2 }
2641
2642
2643
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Co3Ne co3ne(M);
2644
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if(nb_iterations != 0) {
2645 try {
2646
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 co3ne.RVD().set_exact(false);
2647
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 ProgressTask progress("Co3Ne smooth", nb_iterations);
2648
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
 4 for(index_t i = 0; i < nb_iterations; i++) {
2649
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 co3ne.smooth(nb_neighbors);
2650
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 co3ne.RVD().update();
2651
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 progress.next();
2652 }
2653 1 co3ne.end_smooth();
2654 1 }
2655 ✗ catch(const TaskCanceled&) {
2656 // TODO_CANCEL
2657 ✗ }
2658 }
2659
4/8
✓ 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.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
 2 Logger::out("Co3Ne") << "Reconstruction..." << std::endl;
2660
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 co3ne.RVD().set_exact(true);
2661
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 co3ne.reconstruct(radius);
2662 2 }
2663 }
2664