GCC Code Coverage Report

Directory: ./
File: lib/geogram/mesh/mesh_remesh.cpp
Date: 2026-09-07 02:25:23
Exec Total Coverage
Lines: 160 192 83.3%
Functions: 5 5 100.0%
Branches: 161 312 51.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/mesh/mesh_remesh.h>
41 #include <geogram/mesh/mesh.h>
42 #include <geogram/mesh/mesh_repair.h>
43 #include <geogram/mesh/mesh_geometry.h>
44 #include <geogram/mesh/mesh_halfedges.h>
45 #include <geogram/mesh/mesh_reorder.h>
46 #include <geogram/mesh/mesh_preprocessing.h>
47 #include <geogram/mesh/mesh_io.h>
48 #include <geogram/mesh/mesh_AABB.h>
49 #include <geogram/voronoi/CVT.h>
50 #include <geogram/basic/command_line.h>
51 #include <geogram/basic/stopwatch.h>
52 #include <geogram/basic/progress.h>
53 #include <geogram/bibliography/bibliography.h>
54 #include <geogram/mesh/mesh_io.h>
55 #include <geogram/NL/nl.h>
56
57
58 /****************************************************************************/
59
60 namespace GEO {
61
62 6 void remesh_smooth(
63 Mesh& M_in, Mesh& M_out,
64 index_t nb_points,
65 coord_index_t dim,
66 index_t nb_Lloyd_iter,
67 index_t nb_Newton_iter,
68 index_t Newton_m,
69 bool adjust,
70 double adjust_max_edge_distance,
71 double adjust_border_importance
72
73 ) {
74
75 6 geo_cite("DBLP:journals/cgf/YanLLSW09");
76 6 geo_cite("DBLP:conf/imr/LevyB12");
77
78 if(dim == 0) {
79 dim = coord_index_t(M_in.vertices.dimension());
80 }
81
82 geo_argused(dim);
83
84
2/4
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
 12 Stopwatch W("Remesh");
85
86
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 CentroidalVoronoiTesselation CVT(&M_in);
87
88
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 6 if(nb_points == 0) {
89 nb_points = M_in.vertices.nb();
90 }
91
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 CVT.compute_initial_sampling(nb_points, true); // true: for verbose
92
93 try {
94
3/6
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
 12 ProgressTask progress("Lloyd", 100);
95 CVT.set_progress_logger(&progress);
96
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 CVT.Lloyd_iterations(nb_Lloyd_iter);
97 6 }
98 ✗ catch(const TaskCanceled&) {
99 // TODO_CANCEL
100 ✗ }
101
102
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 6 if(nb_Newton_iter != 0) {
103 try {
104 ✗ ProgressTask progress("Newton", 100);
105 CVT.set_progress_logger(&progress);
106 ✗ CVT.Newton_iterations(nb_Newton_iter, Newton_m);
107 ✗ }
108 ✗ catch(const TaskCanceled&) {
109 // TODO_CANCEL
110 ✗ }
111 }
112
113
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 6 if(M_in.vertices.dimension() == 6 &&
114
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 6 times.
 10 CmdLine::get_arg_bool("dbg:save_6d")
115 ) {
116 ✗ Logger::out("Remesh")
117 << "Saving source mesh into mesh6.obj6" << std::endl;
118 ✗ mesh_save(M_in, "mesh6.obj6");
119 ✗ Logger::out("Remesh")
120 << "Saving sampling into points6.txt" << std::endl;
121 ✗ std::ofstream out("points6.txt");
122 out << CVT.delaunay()->nb_vertices() << std::endl;
123 ✗ for(index_t i = 0; i < CVT.delaunay()->nb_vertices(); i++) {
124 ✗ for(coord_index_t c = 0; c < 6; c++) {
125 ✗ out << CVT.delaunay()->vertex_ptr(i)[c] << " ";
126 }
127 out << std::endl;
128 }
129 ✗ }
130
131 // Delete auxiliary storage used for each threads (it uses a lot of RAM,
132 // we need this RAM to create the surface now...)
133
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 CVT.RVD()->delete_threads();
134
135 CVT.set_use_RVC_centroids(
136
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 CmdLine::get_arg_bool("remesh:RVC_centroids")
137 );
138
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 bool multi_nerve = CmdLine::get_arg_bool("remesh:multi_nerve");
139
140
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 Logger::out("Remesh") << "Computing RVD..." << std::endl;
141
142
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 CVT.compute_surface(&M_out, multi_nerve);
143
3/6
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 6 times.
 12 if(CmdLine::get_arg_bool("dbg:save_ANN_histo")) {
144 ✗ Logger::out("ANN")
145 << "Saving histogram to ANN_histo.dat" << std::endl;
146 ✗ std::ofstream out("ANN_histo.dat");
147 ✗ CVT.delaunay()->save_histogram(out);
148 ✗ }
149
150
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 6 if(adjust) {
151
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 mesh_adjust_surface(
152 M_out, M_in, adjust_max_edge_distance,
153 false, adjust_border_importance
154 );
155 }
156 6 }
157
158 /************************************************************************/
159
160 /**
161 * \brief Gets the middle segment of a quad
162 * \param[in] AABB a reference to a MeshFacetsAABB
163 * \param[in] f a facet
164 * \param[out] q1 , q2 the extremities of the middle segment
165 * \details The quad has vertices p1, p2, p3, p4, on exit
166 * q1 = 1/2(p1+p4) and q2 = 1/2(p2+p3)
167 */
168
1/2
✓ Branch 0 taken 388 times.
✗ Branch 1 not taken.
 388 inline void get_quad_middle_segment(
169 const MeshFacetsAABB& AABB, index_t f, vec3& q1, vec3& q2
170 ) {
171
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 388 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 388 geo_assert(AABB.mesh()->facets.nb_vertices(f) == 4);
172 388 vec3 p1 = AABB.mesh()->facets.point(f,0);
173 388 vec3 p2 = AABB.mesh()->facets.point(f,1);
174 388 vec3 p3 = AABB.mesh()->facets.point(f,2);
175 388 vec3 p4 = AABB.mesh()->facets.point(f,3);
176 388 q1 = 0.5*(p1+p4);
177 388 q2 = 0.5*(p2+p3);
178 388 }
179
180 /**
181 * \brief Gets the nearest point on a surface along a ray
182 * \param[in] AABB the facets of the surface
183 * as a MeshFacetsAABB
184 * \param[in] R1 the ray, both directions are tested to find
185 * the nearest point
186 * \param[in] max_dist if the nearest point is further away
187 * than \p max_dist, then the origin of the ray is returned
188 * \param[in] ribbon_mode if set, project on middle segment
189 * of quad
190 */
191 89998 inline vec3 nearest_along_bidirectional_ray(
192 const MeshFacetsAABB& AABB, const Ray& R1,
193 double max_dist, bool ribbon_mode=false
194 ) {
195 89998 vec3 p = R1.origin;
196 89998 vec3 result = p;
197 Ray R2(R1.origin, -R1.direction);
198 MeshFacetsAABB::Intersection I1;
199 MeshFacetsAABB::Intersection I2;
200 89998 bool has_I1 = AABB.ray_nearest_intersection(R1,I1);
201 89998 bool has_I2 = AABB.ray_nearest_intersection(R2,I2);
202
203
2/2
✓ Branch 0 taken 388 times.
✓ Branch 1 taken 89610 times.
 89998 if(has_I1 && ribbon_mode) {
204 vec3 q,q1,q2;
205 double l1,l2;
206 388 get_quad_middle_segment(AABB, I1.f, q1, q2);
207 388 Geom::point_segment_squared_distance(I1.p, q1, q2, q, l1, l2);
208 388 I1.p=q;
209 }
210
211
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 89998 times.
 89998 if(has_I2 && ribbon_mode) {
212 vec3 q,q1,q2;
213 double l1,l2;
214 ✗ get_quad_middle_segment(AABB, I2.f, q1, q2);
215 ✗ Geom::point_segment_squared_distance(I2.p, q1, q2, q, l1, l2);
216 ✗ I2.p=q;
217 }
218
219
2/2
✓ Branch 0 taken 18245 times.
✓ Branch 1 taken 71753 times.
 89998 if(has_I1 && !has_I2) {
220 18245 result = I1.p;
221 }
222
2/2
✓ Branch 0 taken 5636 times.
✓ Branch 1 taken 84362 times.
 89998 if(!has_I1 && has_I2) {
223 5636 result = I2.p;
224 }
225
2/2
✓ Branch 0 taken 66117 times.
✓ Branch 1 taken 23881 times.
 89998 if(has_I1 && has_I2) {
226
2/2
✓ Branch 0 taken 24368 times.
✓ Branch 1 taken 41749 times.
 66117 if(Geom::distance2(p,I1.p) < Geom::distance2(p,I2.p)) {
227 24368 result = I1.p;
228 } else {
229 41749 result = I2.p;
230 }
231 }
232
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 89998 times.
 89998 if(Geom::distance2(result,p) > max_dist*max_dist) {
233 ✗ result = p;
234 }
235 89998 return result;
236 }
237
238
239 /**
240 * \brief Creates a ribbon orthogonal to the border of a
241 * surface mesh.
242 * \param[in] M the input surface mesh
243 * \param[out] ribbon the generated ribbon
244 * \param[in] height the height of the ribbon
245 */
246 1 static void create_ribbon_on_border(
247 const Mesh& M,
248 Mesh& ribbon,
249 double height
250 ) {
251 index_t nb_border_edges=0;
252 1 vector<vec3> Nv(M.vertices.nb(), vec3(0.0, 0.0, 0.0));
253
2/2
✓ Branch 0 taken 3000 times.
✓ Branch 1 taken 1 times.
 3001 for(index_t f: M.facets) {
254
2/4
✓ Branch 1 taken 3000 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3000 times.
✗ Branch 4 not taken.
 3000 vec3 N = Geom::mesh_facet_normal(M, f);
255
2/2
✓ Branch 0 taken 9000 times.
✓ Branch 1 taken 3000 times.
 12000 for(index_t c1: M.facets.corners(f)) {
256
2/2
✓ Branch 0 taken 178 times.
✓ Branch 1 taken 8822 times.
 9000 if(M.facet_corners.adjacent_facet(c1) == NO_FACET) {
257 index_t c2 = M.facets.next_corner_around_facet(f, c1);
258 index_t v1 = M.facet_corners.vertex(c1);
259 index_t v2 = M.facet_corners.vertex(c2);
260 Nv[v1] += N;
261 Nv[v2] += N;
262 178 ++nb_border_edges;
263 }
264 }
265 }
266
267
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 ribbon.vertices.set_dimension(3);
268
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 ribbon.vertices.create_vertices(4*nb_border_edges);
269
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 ribbon.facets.create_quads(nb_border_edges);
270
271 index_t cur_border_e = 0;
272
2/2
✓ Branch 0 taken 3000 times.
✓ Branch 1 taken 1 times.
 3001 for(index_t f: M.facets) {
273
2/2
✓ Branch 0 taken 9000 times.
✓ Branch 1 taken 3000 times.
 12000 for(index_t c1: M.facets.corners(f)) {
274
2/2
✓ Branch 0 taken 178 times.
✓ Branch 1 taken 8822 times.
 9000 if(M.facet_corners.adjacent_facet(c1) == NO_FACET) {
275 index_t c2 = M.facets.next_corner_around_facet(f, c1);
276 index_t v1 = M.facet_corners.vertex(c1);
277 index_t v2 = M.facet_corners.vertex(c2);
278 const vec3& p1 = M.vertices.point(v1);
279 const vec3& p2 = M.vertices.point(v2);
280
281 178 vec3 U1 = 0.5*height * normalize(Nv[v1]);
282
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 178 times.
 178 vec3 U2 = 0.5*height * normalize(Nv[v2]);
283
284 vec3 q1 = p1 + U1;
285 vec3 q2 = p2 + U2;
286 vec3 q3 = p2 - U2;
287 vec3 q4 = p1 - U1;
288
289 178 ribbon.vertices.point(4*cur_border_e ) = q1;
290
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 178 times.
 178 ribbon.vertices.point(4*cur_border_e+1) = q2;
291 178 ribbon.vertices.point(4*cur_border_e+2) = q3;
292
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 178 times.
 178 ribbon.vertices.point(4*cur_border_e+3) = q4;
293
294 ribbon.facets.set_vertex(cur_border_e, 0, 4*cur_border_e );
295 ribbon.facets.set_vertex(cur_border_e, 1, 4*cur_border_e+1);
296 ribbon.facets.set_vertex(cur_border_e, 2, 4*cur_border_e+2);
297 ribbon.facets.set_vertex(cur_border_e, 3, 4*cur_border_e+3);
298 178 ++cur_border_e;
299 }
300 }
301 }
302 1 }
303
304 /************************************************************************/
305
306 6 void GEOGRAM_API mesh_adjust_surface(
307 Mesh& surface,
308 Mesh& reference,
309 double max_edge_distance,
310 bool project_borders,
311 double border_importance
312 ) {
313 // The algorithm:
314 // 1) For each surface vertex v (located at Pv) with normal Nv,
315 // we determine a "target point" Qv as the intersection between
316 // the ray R(Pv,Nv) and the reference surface.
317 // 2) For each facet f with center point Pf, we determine a "target
318 // point" Qf that corresponds to the intersection between the ray
319 // R(Pf, Nf) and the reference surface, where Nf is the sum of
320 // the directions Nv associated with the vertices of the facet.
321 // 3) We optimize for the lambda's in the relations that correspond to
322 // 1) and 2), that is, the sum of the squared distances
323 // || P + lambda N - Q ||^2 for each v, for each f
324 // 4) For each v, Pv = Pv + lambda_v Nv
325 //
326 // For surfaces with borders, there is a subtlety: instead of using
327 // the normal to the surface, we use for Nv a direction tangent to
328 // the surface and normal to the border of the surface, and to find
329 // the reference point Qv, we construct a "ribbon", obtained by
330 // sweeping a segment normal to the reference surface along the
331 // border of the reference surface. The point Qv is obtained by
332 // - compute the intersection I between the ray R(Pv, Nv) and the
333 // ribbon. This intersection is in a quad supported by a segment
334 // [q1,q2] on the border of the reference surface
335 // - Qv is determined as the nearest point to I on [q1,q2]
336 //
337 // For each border edge e=(i,j), the center point Pe=0.5(Pi+Pj), the
338 // direction Ne = 0.5(Ni+Nj), a point Qe is determined (using the same
339 // algorithm as in the previous point) and a least-squares term
340 // || Pe - Qe ||^2 is added to the quantity to be minimized.
341 //
342 // An optional final step (brutally) assigns Pv = Qv for each border
343 // vertex (but in general it gives a worse result on the facets
344 // adjacent to the border, so the option project_borders is deactivated
345 // by default).
346
347 // Geometric search uses larger maximum
348 // distance for points on the border.
349 const double border_distance_factor = 10.0;
350
351 6 MeshFacetsAABB AABB(reference);
352
353 // For each vertex, direction along which the neighbor
354 // on the surface is searched (or the neighbor on the
355 // ribbon if it is a vertex on the border)
356 ✗ vector<vec3> Nv(surface.vertices.nb(), vec3(0.0, 0.0, 0.0));
357
358 // average edge length incident to vertex
359
1/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 6 vector<double> Lv(surface.vertices.nb(), 0.0);
360
361 // number of edges incident to a vertex
362
1/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 6 vector<index_t> Cv(surface.vertices.nb(), 0);
363
364 // Compute directions Nv for inner vertices
365 // (use surface normal)
366
2/2
✓ Branch 0 taken 59804 times.
✓ Branch 1 taken 6 times.
 59810 for(index_t f: surface.facets) {
367
1/2
✓ Branch 1 taken 59804 times.
✗ Branch 2 not taken.
 59804 vec3 n = Geom::mesh_facet_normal(surface, f);
368 index_t d = surface.facets.nb_vertices(f);
369
2/2
✓ Branch 0 taken 179412 times.
✓ Branch 1 taken 59804 times.
 239216 for(index_t lv=0;lv<d;++lv) {
370 index_t v1 = surface.facets.vertex(f,lv);
371 119608 index_t v2 = surface.facets.vertex(
372
2/2
✓ Branch 0 taken 119608 times.
✓ Branch 1 taken 59804 times.
 179412 f,(lv==d-1)?0:lv+1
373 );
374 Nv[v1] += n;
375 double l = Geom::distance2(
376 surface.vertices.point(v1),
377 surface.vertices.point(v2)
378 );
379 179412 l = ::sqrt(l);
380 179412 Lv[v1] += l;
381 179412 Lv[v2] += l;
382 179412 Cv[v1]++;
383 179412 Cv[v2]++;
384 }
385 }
386
387 // Normalize average incident edge length
388
2/2
✓ Branch 0 taken 30000 times.
✓ Branch 1 taken 6 times.
 30006 for(index_t v: surface.vertices) {
389
1/2
✓ Branch 0 taken 30000 times.
✗ Branch 1 not taken.
 30000 if(Cv[v] != 0) {
390 30000 Lv[v] /= double(Cv[v]);
391 }
392 }
393
394
395 // Compute Nv for vertices on the border: the used vector
396 // is normal to the border (tangent to the surface), because
397 // we will find nearest neighbor from the "ribbon" mesh.
398 index_t nb_v_on_border = 0;
399 6 vector<bool> v_on_border(surface.vertices.nb(), false);
400
2/2
✓ Branch 0 taken 59804 times.
✓ Branch 1 taken 6 times.
 59810 for(index_t f: surface.facets) {
401
2/2
✓ Branch 0 taken 179412 times.
✓ Branch 1 taken 59804 times.
 239216 for(index_t c1: surface.facets.corners(f)) {
402
2/2
✓ Branch 0 taken 194 times.
✓ Branch 1 taken 179218 times.
 179412 if(surface.facet_corners.adjacent_facet(c1) == NO_FACET) {
403 index_t v = surface.facet_corners.vertex(c1);
404 v_on_border[v] = true;
405 194 ++nb_v_on_border;
406 194 Nv[v] = vec3(0.0, 0.0, 0.0);
407 }
408 }
409 }
410
411
2/2
✓ Branch 0 taken 59804 times.
✓ Branch 1 taken 6 times.
 59810 for(index_t f: surface.facets) {
412
1/2
✓ Branch 1 taken 59804 times.
✗ Branch 2 not taken.
 59804 vec3 N = Geom::mesh_facet_normal(surface, f);
413
2/2
✓ Branch 0 taken 179412 times.
✓ Branch 1 taken 59804 times.
 239216 for(index_t c1: surface.facets.corners(f)) {
414
2/2
✓ Branch 0 taken 194 times.
✓ Branch 1 taken 179218 times.
 179412 if(surface.facet_corners.adjacent_facet(c1) == NO_FACET) {
415 index_t c2 = surface.facets.next_corner_around_facet(f, c1);
416 index_t v1 = surface.facet_corners.vertex(c1);
417 index_t v2 = surface.facet_corners.vertex(c2);
418 const vec3& p1 = surface.vertices.point(v1);
419 const vec3& p2 = surface.vertices.point(v2);
420 vec3 Ne = cross(p2 - p1, N);
421 Nv[v1] += Ne;
422 Nv[v2] += Ne;
423 }
424 }
425 }
426
427 bool reference_has_borders = false;
428
2/2
✓ Branch 0 taken 15870 times.
✓ Branch 1 taken 5 times.
 15875 for(index_t c: reference.facet_corners) {
429
2/2
✓ Branch 0 taken 15869 times.
✓ Branch 1 taken 1 times.
 15870 if(reference.facet_corners.adjacent_facet(c) == NO_FACET) {
430 reference_has_borders = true;
431 break;
432 }
433 }
434
435 // Create ribbon and AABB for the ribbon if the reference surface
436 // has vertices on the border
437
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 Mesh border_ribbon;
438 MeshFacetsAABB border_ribbon_AABB;
439
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 5 times.
 6 if(nb_v_on_border != 0 && reference_has_borders) {
440
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 create_ribbon_on_border(
441 reference, border_ribbon,
442
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 max_edge_distance*4.0*surface_average_edge_length(surface)
443 );
444
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 border_ribbon_AABB.initialize(border_ribbon);
445 }
446
447 // nearest point along Nv
448 vector<vec3> Qv(surface.vertices.nb());
449
2/2
✓ Branch 0 taken 30000 times.
✓ Branch 1 taken 6 times.
 30006 for(index_t v: surface.vertices) {
450
2/2
✓ Branch 0 taken 194 times.
✓ Branch 1 taken 29806 times.
 30000 vec3 p = surface.vertices.point(v);
451
3/4
✓ Branch 0 taken 194 times.
✓ Branch 1 taken 29806 times.
✓ Branch 2 taken 194 times.
✗ Branch 3 not taken.
 30000 if(v_on_border[v] && reference_has_borders) {
452
1/2
✓ Branch 1 taken 194 times.
✗ Branch 2 not taken.
 194 Qv[v] = nearest_along_bidirectional_ray(
453 194 border_ribbon_AABB, Ray(p, Nv[v]),
454
1/2
✓ Branch 1 taken 194 times.
✗ Branch 2 not taken.
 194 border_distance_factor*max_edge_distance*Lv[v],
455 true
456 );
457 } else {
458
1/2
✓ Branch 1 taken 29806 times.
✗ Branch 2 not taken.
 29806 Qv[v] = nearest_along_bidirectional_ray(
459 29806 AABB, Ray(p, Nv[v]),max_edge_distance*Lv[v]
460 );
461 }
462 }
463
464
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 nlNewContext();
465
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
466
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 nlSolverParameteri(
467 NL_NB_VARIABLES, NLint(surface.vertices.nb())
468 );
469
470
471
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 nlBegin(NL_SYSTEM);
472
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 nlBegin(NL_MATRIX);
473
474 // For each vertex v, least squares constraint that
475 // makes the vertex "attracted" by Qv[v]
476
2/2
✓ Branch 0 taken 30000 times.
✓ Branch 1 taken 6 times.
 30006 for(index_t v: surface.vertices) {
477
2/2
✓ Branch 0 taken 194 times.
✓ Branch 1 taken 29806 times.
 30000 if(v_on_border[v]) {
478
1/2
✓ Branch 1 taken 194 times.
✗ Branch 2 not taken.
 194 nlRowScaling(border_importance);
479 }
480
481 // p + lambda_v * Nv = q ---> lambda_v * Nv = q - v
482
2/2
✓ Branch 0 taken 90000 times.
✓ Branch 1 taken 30000 times.
 120000 for(index_t c=0; c<3; ++c) {
483
1/2
✓ Branch 1 taken 90000 times.
✗ Branch 2 not taken.
 90000 nlBegin(NL_ROW);
484
1/2
✓ Branch 1 taken 90000 times.
✗ Branch 2 not taken.
 90000 nlCoefficient(v,Nv[v][c]);
485
1/2
✓ Branch 1 taken 90000 times.
✗ Branch 2 not taken.
 90000 nlRightHandSide(Qv[v][c] - surface.vertices.point_ptr(v)[c]);
486
1/2
✓ Branch 1 taken 90000 times.
✗ Branch 2 not taken.
 90000 nlEnd(NL_ROW);
487 }
488 }
489
490 // For each facet f, least squares constraint that
491 // makes the center of the facet "attracted" by the
492 // nearest point along the averaged directions of
493 // all vertices of the facet
494
2/2
✓ Branch 0 taken 59804 times.
✓ Branch 1 taken 6 times.
 59810 for(index_t f: surface.facets) {
495 index_t d = surface.facets.nb_vertices(f);
496
497 vec3 Nf(0.0, 0.0, 0.0);
498 vec3 Pf;
499 double Lf=0.0;
500
501
2/2
✓ Branch 0 taken 179412 times.
✓ Branch 1 taken 59804 times.
 239216 for(index_t lv=0; lv<d; ++lv) {
502 index_t v= surface.facets.vertex(f,lv);
503 Nf += Nv[v];
504 Pf += surface.vertices.point(v);
505 179412 Lf += Lv[v];
506 }
507
1/2
✓ Branch 1 taken 59804 times.
✗ Branch 2 not taken.
 59804 Pf = (1.0 / double(d))*Pf;
508 59804 Lf = (1.0 / double(d))*Lf;
509
1/2
✓ Branch 1 taken 59804 times.
✗ Branch 2 not taken.
 59804 vec3 Qf = nearest_along_bidirectional_ray(
510 59804 AABB, Ray(Pf, Nf), max_edge_distance*Lf
511 );
512
513 // 1/d(Sum Pv + lambda_v Nv) = Qf
514 // --> Pf + 1/d(Sum lambda_v Nv) = Qf
515 // --> Sum (1/d lambda_v Nv) = Qf - Pf
516
2/2
✓ Branch 0 taken 179412 times.
✓ Branch 1 taken 59804 times.
 239216 for(index_t c=0; c<3; ++c) {
517
1/2
✓ Branch 1 taken 179412 times.
✗ Branch 2 not taken.
 179412 nlBegin(NL_ROW);
518
2/2
✓ Branch 0 taken 538236 times.
✓ Branch 1 taken 179412 times.
 717648 for(index_t lv=0; lv<d; ++lv) {
519 index_t v = surface.facets.vertex(f,lv);
520
1/2
✓ Branch 1 taken 538236 times.
✗ Branch 2 not taken.
 538236 nlCoefficient(v,Nv[v][c]/double(d));
521 }
522
1/2
✓ Branch 1 taken 179412 times.
✗ Branch 2 not taken.
 179412 nlRightHandSide(Qf[c]-Pf[c]);
523
1/2
✓ Branch 1 taken 179412 times.
✗ Branch 2 not taken.
 179412 nlEnd(NL_ROW);
524 }
525 }
526
527 // For each edge on the border,
528 // makes the center of the edge attracted by the
529 // nearest point on the ribbon along the averated directions
530 // of the two vertices of the edge
531
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 5 times.
 6 if(nb_v_on_border != 0 && reference_has_borders) {
532
2/2
✓ Branch 0 taken 9804 times.
✓ Branch 1 taken 1 times.
 9805 for(index_t f: surface.facets) {
533
2/2
✓ Branch 0 taken 29412 times.
✓ Branch 1 taken 9804 times.
 39216 for(index_t c1: surface.facets.corners(f)) {
534
2/2
✓ Branch 0 taken 194 times.
✓ Branch 1 taken 29218 times.
 29412 if(surface.facet_corners.adjacent_facet(c1) == NO_FACET) {
535 index_t c2 = surface.facets.next_corner_around_facet(
536 f, c1
537 );
538 index_t v1 = surface.facet_corners.vertex(c1);
539 index_t v2 = surface.facet_corners.vertex(c2);
540 const vec3& p1 = surface.vertices.point(v1);
541 const vec3& p2 = surface.vertices.point(v2);
542 vec3 p = 0.5*(p1+p2);
543 vec3 N = 0.5*(Nv[v1] + Nv[v2]);
544 194 vec3 q = nearest_along_bidirectional_ray(
545 ✗ border_ribbon_AABB, Ray(p, N),
546 194 border_distance_factor*max_edge_distance*
547
1/2
✓ Branch 1 taken 194 times.
✗ Branch 2 not taken.
 194 0.5*(Lv[v1]+Lv[v2]),
548 true
549 );
550
551 // 1/2(p1 + p2 + lambda_1 N1 + lambda_2 N2) = q
552 // --> p + 1/2(lambda_1 N1 + lambda_2 N2) = q
553 // --> (1/2 lambda_1 N1 + 1/2 lambda_2 N2) = q - p
554
2/2
✓ Branch 0 taken 582 times.
✓ Branch 1 taken 194 times.
 776 for(index_t c=0; c<3; ++c) {
555
1/2
✓ Branch 1 taken 582 times.
✗ Branch 2 not taken.
 582 nlRowScaling(0.5*border_importance);
556
1/2
✓ Branch 1 taken 582 times.
✗ Branch 2 not taken.
 582 nlBegin(NL_ROW);
557
1/2
✓ Branch 1 taken 582 times.
✗ Branch 2 not taken.
 582 nlCoefficient(v1,0.5*N[c]);
558 nlCoefficient(v2,0.5*N[c]);
559
1/2
✓ Branch 1 taken 582 times.
✗ Branch 2 not taken.
 582 nlRightHandSide(q[c]-p[c]);
560
1/2
✓ Branch 1 taken 582 times.
✗ Branch 2 not taken.
 582 nlEnd(NL_ROW);
561 }
562
563 }
564 }
565 }
566 }
567
568
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 nlEnd(NL_MATRIX);
569
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 nlEnd(NL_SYSTEM);
570
571
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 nlSolve();
572
573 // Displace each vertex v along Nv[v] by
574 // the solution of the least squares problem
575 // at v
576
2/2
✓ Branch 0 taken 30000 times.
✓ Branch 1 taken 6 times.
 30006 for(index_t v: surface.vertices) {
577 vec3& p = surface.vertices.point(v);
578 p += nlGetVariable(v)*Nv[v];
579 }
580
581 // (Brutally) project border vertices
582
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 if(project_borders && nb_v_on_border != 0 && reference_has_borders) {
583 ✗ for(index_t v: surface.vertices) {
584 ✗ if(v_on_border[v]) {
585 ✗ vec3 p = surface.vertices.point(v);
586 ✗ vec3 q = nearest_along_bidirectional_ray(
587 ✗ border_ribbon_AABB, Ray(p, Nv[v]),
588 ✗ border_distance_factor*max_edge_distance*0.5*(Lv[v]),
589 true
590 );
591 ✗ surface.vertices.point(v) = q;
592 }
593 }
594 }
595
596
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 nlDeleteContext(nlGetCurrent());
597 12 }
598 }
599