GCC Code Coverage Report

Directory: ./
File: mesh/mesh_remesh.cpp
Date: 2026-09-27 03:12:47
Exec Total Coverage
Lines: 168 197 85.3%
Functions: 5 5 100.0%
Branches: 167 316 52.8%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #include <geogram/mesh/mesh_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 7 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 7 geo_cite("DBLP:journals/cgf/YanLLSW09");
76 7 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 7 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 7 times.
✗ Branch 6 not taken.
 14 Stopwatch W("Remesh");
85
86
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 CentroidalVoronoiTesselation CVT(&M_in);
87
88
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 7 if(nb_points == 0) {
89 nb_points = M_in.vertices.nb();
90 }
91
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 CVT.compute_initial_sampling(nb_points, true); // true: for verbose
92
93 try {
94
3/6
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
 14 ProgressTask progress("Lloyd", 100);
95 CVT.set_progress_logger(&progress);
96
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 CVT.Lloyd_iterations(nb_Lloyd_iter);
97 7 }
98 ✗ catch(const TaskCanceled&) {
99 // TODO_CANCEL
100 ✗ }
101
102
1/2
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
 7 if(nb_Newton_iter != 0) {
103 try {
104
3/6
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
 14 ProgressTask progress("Newton", 100);
105 CVT.set_progress_logger(&progress);
106
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 CVT.Newton_iterations(nb_Newton_iter, Newton_m);
107 7 }
108 ✗ catch(const TaskCanceled&) {
109 // TODO_CANCEL
110 ✗ }
111 }
112
113
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 2 times.
 7 if(M_in.vertices.dimension() == 6 &&
114
3/6
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 7 times.
 12 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 7 times.
✗ Branch 2 not taken.
 7 CVT.RVD()->delete_threads();
134
135 CVT.set_use_RVC_centroids(
136
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
 7 CmdLine::get_arg_bool("remesh:RVC_centroids")
137 );
138
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
 7 bool multi_nerve = CmdLine::get_arg_bool("remesh:multi_nerve");
139
140
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
 7 Logger::out("Remesh") << "Computing RVD..." << std::endl;
141
142
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 CVT.compute_surface(&M_out, multi_nerve);
143
3/6
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 7 times.
 14 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 7 times.
✗ Branch 1 not taken.
 7 if(adjust) {
151
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 mesh_adjust_surface(
152 M_out, M_in, adjust_max_edge_distance,
153 false, adjust_border_importance
154 );
155 }
156
157
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 7 times.
 7 if(Geom::mesh_enclosed_volume(M_out) < 0.0) {
158 ✗ M_out.facets.flip();
159 }
160 7 }
161
162 /************************************************************************/
163
164 /**
165 * \brief Gets the middle segment of a quad
166 * \param[in] AABB a reference to a MeshFacetsAABB
167 * \param[in] f a facet
168 * \param[out] q1 , q2 the extremities of the middle segment
169 * \details The quad has vertices p1, p2, p3, p4, on exit
170 * q1 = 1/2(p1+p4) and q2 = 1/2(p2+p3)
171 */
172
1/2
✓ Branch 0 taken 526 times.
✗ Branch 1 not taken.
 526 inline void get_quad_middle_segment(
173 const MeshFacetsAABB& AABB, index_t f, vec3& q1, vec3& q2
174 ) {
175
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 526 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 526 geo_assert(AABB.mesh()->facets.nb_vertices(f) == 4);
176 526 vec3 p1 = AABB.mesh()->facets.point(f,0);
177 526 vec3 p2 = AABB.mesh()->facets.point(f,1);
178 526 vec3 p3 = AABB.mesh()->facets.point(f,2);
179 526 vec3 p4 = AABB.mesh()->facets.point(f,3);
180 526 q1 = 0.5*(p1+p4);
181 526 q2 = 0.5*(p2+p3);
182 526 }
183
184 /**
185 * \brief Gets the nearest point on a surface along a ray
186 * \param[in] AABB the facets of the surface
187 * as a MeshFacetsAABB
188 * \param[in] R1 the ray, both directions are tested to find
189 * the nearest point
190 * \param[in] max_dist if the nearest point is further away
191 * than \p max_dist, then the origin of the ray is returned
192 * \param[in] ribbon_mode if set, project on middle segment
193 * of quad
194 */
195 104996 inline vec3 nearest_along_bidirectional_ray(
196 const MeshFacetsAABB& AABB, const Ray& R1,
197 double max_dist, bool ribbon_mode=false
198 ) {
199 104996 vec3 p = R1.origin;
200 104996 vec3 result = p;
201 Ray R2(R1.origin, -R1.direction);
202 MeshFacetsAABB::Intersection I1;
203 MeshFacetsAABB::Intersection I2;
204 104996 bool has_I1 = AABB.ray_nearest_intersection(R1,I1);
205 104996 bool has_I2 = AABB.ray_nearest_intersection(R2,I2);
206
207
2/2
✓ Branch 0 taken 526 times.
✓ Branch 1 taken 104470 times.
 104996 if(has_I1 && ribbon_mode) {
208 vec3 q,q1,q2;
209 double l1,l2;
210 526 get_quad_middle_segment(AABB, I1.f, q1, q2);
211 526 Geom::point_segment_squared_distance(I1.p, q1, q2, q, l1, l2);
212 526 I1.p=q;
213 }
214
215
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 104996 times.
 104996 if(has_I2 && ribbon_mode) {
216 vec3 q,q1,q2;
217 double l1,l2;
218 ✗ get_quad_middle_segment(AABB, I2.f, q1, q2);
219 ✗ Geom::point_segment_squared_distance(I2.p, q1, q2, q, l1, l2);
220 ✗ I2.p=q;
221 }
222
223
2/2
✓ Branch 0 taken 26556 times.
✓ Branch 1 taken 78440 times.
 104996 if(has_I1 && !has_I2) {
224 26556 result = I1.p;
225 }
226
2/2
✓ Branch 0 taken 12922 times.
✓ Branch 1 taken 92074 times.
 104996 if(!has_I1 && has_I2) {
227 12922 result = I2.p;
228 }
229
2/2
✓ Branch 0 taken 65518 times.
✓ Branch 1 taken 39478 times.
 104996 if(has_I1 && has_I2) {
230
2/2
✓ Branch 0 taken 33773 times.
✓ Branch 1 taken 31745 times.
 65518 if(Geom::distance2(p,I1.p) < Geom::distance2(p,I2.p)) {
231 33773 result = I1.p;
232 } else {
233 31745 result = I2.p;
234 }
235 }
236
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 104995 times.
 104996 if(Geom::distance2(result,p) > max_dist*max_dist) {
237 1 result = p;
238 }
239 104996 return result;
240 }
241
242
243 /**
244 * \brief Creates a ribbon orthogonal to the border of a
245 * surface mesh.
246 * \param[in] M the input surface mesh
247 * \param[out] ribbon the generated ribbon
248 * \param[in] height the height of the ribbon
249 */
250 2 static void create_ribbon_on_border(
251 const Mesh& M,
252 Mesh& ribbon,
253 double height
254 ) {
255 index_t nb_border_edges=0;
256 2 vector<vec3> Nv(M.vertices.nb(), vec3(0.0, 0.0, 0.0));
257
2/2
✓ Branch 0 taken 65467 times.
✓ Branch 1 taken 2 times.
 65469 for(index_t f: M.facets) {
258
2/4
✓ Branch 1 taken 65467 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 65467 times.
✗ Branch 4 not taken.
 65467 vec3 N = Geom::mesh_facet_normal(M, f);
259
2/2
✓ Branch 0 taken 196401 times.
✓ Branch 1 taken 65467 times.
 261868 for(index_t c1: M.facets.corners(f)) {
260
2/2
✓ Branch 0 taken 1233 times.
✓ Branch 1 taken 195168 times.
 196401 if(M.facet_corners.adjacent_facet(c1) == NO_FACET) {
261 index_t c2 = M.facets.next_corner_around_facet(f, c1);
262 index_t v1 = M.facet_corners.vertex(c1);
263 index_t v2 = M.facet_corners.vertex(c2);
264 Nv[v1] += N;
265 Nv[v2] += N;
266 1233 ++nb_border_edges;
267 }
268 }
269 }
270
271
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 ribbon.vertices.set_dimension(3);
272
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 ribbon.vertices.create_vertices(4*nb_border_edges);
273
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 ribbon.facets.create_quads(nb_border_edges);
274
275 index_t cur_border_e = 0;
276
2/2
✓ Branch 0 taken 65467 times.
✓ Branch 1 taken 2 times.
 65469 for(index_t f: M.facets) {
277
2/2
✓ Branch 0 taken 196401 times.
✓ Branch 1 taken 65467 times.
 261868 for(index_t c1: M.facets.corners(f)) {
278
2/2
✓ Branch 0 taken 1233 times.
✓ Branch 1 taken 195168 times.
 196401 if(M.facet_corners.adjacent_facet(c1) == NO_FACET) {
279 index_t c2 = M.facets.next_corner_around_facet(f, c1);
280 index_t v1 = M.facet_corners.vertex(c1);
281 index_t v2 = M.facet_corners.vertex(c2);
282 const vec3& p1 = M.vertices.point(v1);
283 const vec3& p2 = M.vertices.point(v2);
284
285 1233 vec3 U1 = 0.5*height * normalize(Nv[v1]);
286
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1233 times.
 1233 vec3 U2 = 0.5*height * normalize(Nv[v2]);
287
288 vec3 q1 = p1 + U1;
289 vec3 q2 = p2 + U2;
290 vec3 q3 = p2 - U2;
291 vec3 q4 = p1 - U1;
292
293 1233 ribbon.vertices.point(4*cur_border_e ) = q1;
294 1233 ribbon.vertices.point(4*cur_border_e+1) = q2;
295 1233 ribbon.vertices.point(4*cur_border_e+2) = q3;
296
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1233 times.
 1233 ribbon.vertices.point(4*cur_border_e+3) = q4;
297
298 ribbon.facets.set_vertex(cur_border_e, 0, 4*cur_border_e );
299 ribbon.facets.set_vertex(cur_border_e, 1, 4*cur_border_e+1);
300 ribbon.facets.set_vertex(cur_border_e, 2, 4*cur_border_e+2);
301 ribbon.facets.set_vertex(cur_border_e, 3, 4*cur_border_e+3);
302 1233 ++cur_border_e;
303 }
304 }
305 }
306 2 }
307
308 /************************************************************************/
309
310 7 void GEOGRAM_API mesh_adjust_surface(
311 Mesh& surface,
312 Mesh& reference,
313 double max_edge_distance,
314 bool project_borders,
315 double border_importance
316 ) {
317 // The algorithm:
318 // 1) For each surface vertex v (located at Pv) with normal Nv,
319 // we determine a "target point" Qv as the intersection between
320 // the ray R(Pv,Nv) and the reference surface.
321 // 2) For each facet f with center point Pf, we determine a "target
322 // point" Qf that corresponds to the intersection between the ray
323 // R(Pf, Nf) and the reference surface, where Nf is the sum of
324 // the directions Nv associated with the vertices of the facet.
325 // 3) We optimize for the lambda's in the relations that correspond to
326 // 1) and 2), that is, the sum of the squared distances
327 // || P + lambda N - Q ||^2 for each v, for each f
328 // 4) For each v, Pv = Pv + lambda_v Nv
329 //
330 // For surfaces with borders, there is a subtlety: instead of using
331 // the normal to the surface, we use for Nv a direction tangent to
332 // the surface and normal to the border of the surface, and to find
333 // the reference point Qv, we construct a "ribbon", obtained by
334 // sweeping a segment normal to the reference surface along the
335 // border of the reference surface. The point Qv is obtained by
336 // - compute the intersection I between the ray R(Pv, Nv) and the
337 // ribbon. This intersection is in a quad supported by a segment
338 // [q1,q2] on the border of the reference surface
339 // - Qv is determined as the nearest point to I on [q1,q2]
340 //
341 // For each border edge e=(i,j), the center point Pe=0.5(Pi+Pj), the
342 // direction Ne = 0.5(Ni+Nj), a point Qe is determined (using the same
343 // algorithm as in the previous point) and a least-squares term
344 // || Pe - Qe ||^2 is added to the quantity to be minimized.
345 //
346 // An optional final step (brutally) assigns Pv = Qv for each border
347 // vertex (but in general it gives a worse result on the facets
348 // adjacent to the border, so the option project_borders is deactivated
349 // by default).
350
351 // Geometric search uses larger maximum
352 // distance for points on the border.
353 const double border_distance_factor = 10.0;
354
355 7 MeshFacetsAABB AABB(reference);
356
357 // For each vertex, direction along which the neighbor
358 // on the surface is searched (or the neighbor on the
359 // ribbon if it is a vertex on the border)
360 ✗ vector<vec3> Nv(surface.vertices.nb(), vec3(0.0, 0.0, 0.0));
361
362 // average edge length incident to vertex
363
1/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 7 vector<double> Lv(surface.vertices.nb(), 0.0);
364
365 // number of edges incident to a vertex
366
1/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 7 vector<index_t> Cv(surface.vertices.nb(), 0);
367
368 // Compute directions Nv for inner vertices
369 // (use surface normal)
370
2/2
✓ Branch 0 taken 69733 times.
✓ Branch 1 taken 7 times.
 69740 for(index_t f: surface.facets) {
371
1/2
✓ Branch 1 taken 69733 times.
✗ Branch 2 not taken.
 69733 vec3 n = Geom::mesh_facet_normal(surface, f);
372 index_t d = surface.facets.nb_vertices(f);
373
2/2
✓ Branch 0 taken 209199 times.
✓ Branch 1 taken 69733 times.
 278932 for(index_t lv=0;lv<d;++lv) {
374 index_t v1 = surface.facets.vertex(f,lv);
375 139466 index_t v2 = surface.facets.vertex(
376
2/2
✓ Branch 0 taken 139466 times.
✓ Branch 1 taken 69733 times.
 209199 f,(lv==d-1)?0:lv+1
377 );
378 Nv[v1] += n;
379 double l = Geom::distance2(
380 surface.vertices.point(v1),
381 surface.vertices.point(v2)
382 );
383 209199 l = ::sqrt(l);
384 209199 Lv[v1] += l;
385 209199 Lv[v2] += l;
386 209199 Cv[v1]++;
387 209199 Cv[v2]++;
388 }
389 }
390
391 // Normalize average incident edge length
392
2/2
✓ Branch 0 taken 35000 times.
✓ Branch 1 taken 7 times.
 35007 for(index_t v: surface.vertices) {
393
1/2
✓ Branch 0 taken 35000 times.
✗ Branch 1 not taken.
 35000 if(Cv[v] != 0) {
394 35000 Lv[v] /= double(Cv[v]);
395 }
396 }
397
398
399 // Compute Nv for vertices on the border: the used vector
400 // is normal to the border (tangent to the surface), because
401 // we will find nearest neighbor from the "ribbon" mesh.
402 index_t nb_v_on_border = 0;
403 7 vector<bool> v_on_border(surface.vertices.nb(), false);
404
2/2
✓ Branch 0 taken 69733 times.
✓ Branch 1 taken 7 times.
 69740 for(index_t f: surface.facets) {
405
2/2
✓ Branch 0 taken 209199 times.
✓ Branch 1 taken 69733 times.
 278932 for(index_t c1: surface.facets.corners(f)) {
406
2/2
✓ Branch 0 taken 263 times.
✓ Branch 1 taken 208936 times.
 209199 if(surface.facet_corners.adjacent_facet(c1) == NO_FACET) {
407 index_t v = surface.facet_corners.vertex(c1);
408 263 v_on_border[v] = true;
409 263 ++nb_v_on_border;
410 263 Nv[v] = vec3(0.0, 0.0, 0.0);
411 }
412 }
413 }
414
415
2/2
✓ Branch 0 taken 69733 times.
✓ Branch 1 taken 7 times.
 69740 for(index_t f: surface.facets) {
416
1/2
✓ Branch 1 taken 69733 times.
✗ Branch 2 not taken.
 69733 vec3 N = Geom::mesh_facet_normal(surface, f);
417
2/2
✓ Branch 0 taken 209199 times.
✓ Branch 1 taken 69733 times.
 278932 for(index_t c1: surface.facets.corners(f)) {
418
2/2
✓ Branch 0 taken 263 times.
✓ Branch 1 taken 208936 times.
 209199 if(surface.facet_corners.adjacent_facet(c1) == NO_FACET) {
419 index_t c2 = surface.facets.next_corner_around_facet(f, c1);
420 index_t v1 = surface.facet_corners.vertex(c1);
421 index_t v2 = surface.facet_corners.vertex(c2);
422 const vec3& p1 = surface.vertices.point(v1);
423 const vec3& p2 = surface.vertices.point(v2);
424 vec3 Ne = cross(p2 - p1, N);
425 Nv[v1] += Ne;
426 Nv[v2] += Ne;
427 }
428 }
429 }
430
431 bool reference_has_borders = false;
432
2/2
✓ Branch 0 taken 46802 times.
✓ Branch 1 taken 5 times.
 46807 for(index_t c: reference.facet_corners) {
433
2/2
✓ Branch 0 taken 46800 times.
✓ Branch 1 taken 2 times.
 46802 if(reference.facet_corners.adjacent_facet(c) == NO_FACET) {
434 reference_has_borders = true;
435 break;
436 }
437 }
438
439 // Create ribbon and AABB for the ribbon if the reference surface
440 // has vertices on the border
441
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 Mesh border_ribbon;
442 MeshFacetsAABB border_ribbon_AABB;
443
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 5 times.
 7 if(nb_v_on_border != 0 && reference_has_borders) {
444
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 create_ribbon_on_border(
445 reference, border_ribbon,
446
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 max_edge_distance*4.0*surface_average_edge_length(surface)
447 );
448
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 border_ribbon_AABB.initialize(border_ribbon);
449 }
450
451 // nearest point along Nv
452 vector<vec3> Qv(surface.vertices.nb());
453
2/2
✓ Branch 0 taken 35000 times.
✓ Branch 1 taken 7 times.
 35007 for(index_t v: surface.vertices) {
454 35000 vec3 p = surface.vertices.point(v);
455
3/4
✓ Branch 0 taken 263 times.
✓ Branch 1 taken 34737 times.
✓ Branch 2 taken 263 times.
✗ Branch 3 not taken.
 35000 if(v_on_border[v] && reference_has_borders) {
456
1/2
✓ Branch 1 taken 263 times.
✗ Branch 2 not taken.
 263 Qv[v] = nearest_along_bidirectional_ray(
457 263 border_ribbon_AABB, Ray(p, Nv[v]),
458
1/2
✓ Branch 1 taken 263 times.
✗ Branch 2 not taken.
 263 border_distance_factor*max_edge_distance*Lv[v],
459 true
460 );
461 } else {
462
1/2
✓ Branch 1 taken 34737 times.
✗ Branch 2 not taken.
 34737 Qv[v] = nearest_along_bidirectional_ray(
463 34737 AABB, Ray(p, Nv[v]),max_edge_distance*Lv[v]
464 );
465 }
466 }
467
468
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 nlNewContext();
469
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
470
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 nlSolverParameteri(
471 NL_NB_VARIABLES, NLint(surface.vertices.nb())
472 );
473
474
475
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 nlBegin(NL_SYSTEM);
476
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 nlBegin(NL_MATRIX);
477
478 // For each vertex v, least squares constraint that
479 // makes the vertex "attracted" by Qv[v]
480
2/2
✓ Branch 0 taken 35000 times.
✓ Branch 1 taken 7 times.
 35007 for(index_t v: surface.vertices) {
481
2/2
✓ Branch 0 taken 263 times.
✓ Branch 1 taken 34737 times.
 35000 if(v_on_border[v]) {
482
1/2
✓ Branch 1 taken 263 times.
✗ Branch 2 not taken.
 263 nlRowScaling(border_importance);
483 }
484
485 // p + lambda_v * Nv = q ---> lambda_v * Nv = q - v
486
2/2
✓ Branch 0 taken 105000 times.
✓ Branch 1 taken 35000 times.
 140000 for(index_t c=0; c<3; ++c) {
487
1/2
✓ Branch 1 taken 105000 times.
✗ Branch 2 not taken.
 105000 nlBegin(NL_ROW);
488
1/2
✓ Branch 1 taken 105000 times.
✗ Branch 2 not taken.
 105000 nlCoefficient(v,Nv[v][c]);
489
1/2
✓ Branch 1 taken 105000 times.
✗ Branch 2 not taken.
 105000 nlRightHandSide(Qv[v][c] - surface.vertices.point_ptr(v)[c]);
490
1/2
✓ Branch 1 taken 105000 times.
✗ Branch 2 not taken.
 105000 nlEnd(NL_ROW);
491 }
492 }
493
494 // For each facet f, least squares constraint that
495 // makes the center of the facet "attracted" by the
496 // nearest point along the averaged directions of
497 // all vertices of the facet
498
2/2
✓ Branch 0 taken 69733 times.
✓ Branch 1 taken 7 times.
 69740 for(index_t f: surface.facets) {
499 index_t d = surface.facets.nb_vertices(f);
500
501 vec3 Nf(0.0, 0.0, 0.0);
502 vec3 Pf;
503 double Lf=0.0;
504
505
2/2
✓ Branch 0 taken 209199 times.
✓ Branch 1 taken 69733 times.
 278932 for(index_t lv=0; lv<d; ++lv) {
506 index_t v= surface.facets.vertex(f,lv);
507 Nf += Nv[v];
508 Pf += surface.vertices.point(v);
509 209199 Lf += Lv[v];
510 }
511
1/2
✓ Branch 1 taken 69733 times.
✗ Branch 2 not taken.
 69733 Pf = (1.0 / double(d))*Pf;
512 69733 Lf = (1.0 / double(d))*Lf;
513
1/2
✓ Branch 1 taken 69733 times.
✗ Branch 2 not taken.
 69733 vec3 Qf = nearest_along_bidirectional_ray(
514 69733 AABB, Ray(Pf, Nf), max_edge_distance*Lf
515 );
516
517 // 1/d(Sum Pv + lambda_v Nv) = Qf
518 // --> Pf + 1/d(Sum lambda_v Nv) = Qf
519 // --> Sum (1/d lambda_v Nv) = Qf - Pf
520
2/2
✓ Branch 0 taken 209199 times.
✓ Branch 1 taken 69733 times.
 278932 for(index_t c=0; c<3; ++c) {
521
1/2
✓ Branch 1 taken 209199 times.
✗ Branch 2 not taken.
 209199 nlBegin(NL_ROW);
522
2/2
✓ Branch 0 taken 627597 times.
✓ Branch 1 taken 209199 times.
 836796 for(index_t lv=0; lv<d; ++lv) {
523 index_t v = surface.facets.vertex(f,lv);
524
1/2
✓ Branch 1 taken 627597 times.
✗ Branch 2 not taken.
 627597 nlCoefficient(v,Nv[v][c]/double(d));
525 }
526
1/2
✓ Branch 1 taken 209199 times.
✗ Branch 2 not taken.
 209199 nlRightHandSide(Qf[c]-Pf[c]);
527
1/2
✓ Branch 1 taken 209199 times.
✗ Branch 2 not taken.
 209199 nlEnd(NL_ROW);
528 }
529 }
530
531 // For each edge on the border,
532 // makes the center of the edge attracted by the
533 // nearest point on the ribbon along the averated directions
534 // of the two vertices of the edge
535
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 5 times.
 7 if(nb_v_on_border != 0 && reference_has_borders) {
536
2/2
✓ Branch 0 taken 19733 times.
✓ Branch 1 taken 2 times.
 19735 for(index_t f: surface.facets) {
537
2/2
✓ Branch 0 taken 59199 times.
✓ Branch 1 taken 19733 times.
 78932 for(index_t c1: surface.facets.corners(f)) {
538
2/2
✓ Branch 0 taken 263 times.
✓ Branch 1 taken 58936 times.
 59199 if(surface.facet_corners.adjacent_facet(c1) == NO_FACET) {
539 index_t c2 = surface.facets.next_corner_around_facet(
540 f, c1
541 );
542 index_t v1 = surface.facet_corners.vertex(c1);
543 index_t v2 = surface.facet_corners.vertex(c2);
544 const vec3& p1 = surface.vertices.point(v1);
545 const vec3& p2 = surface.vertices.point(v2);
546 vec3 p = 0.5*(p1+p2);
547 vec3 N = 0.5*(Nv[v1] + Nv[v2]);
548 263 vec3 q = nearest_along_bidirectional_ray(
549 ✗ border_ribbon_AABB, Ray(p, N),
550 263 border_distance_factor*max_edge_distance*
551
1/2
✓ Branch 1 taken 263 times.
✗ Branch 2 not taken.
 263 0.5*(Lv[v1]+Lv[v2]),
552 true
553 );
554
555 // 1/2(p1 + p2 + lambda_1 N1 + lambda_2 N2) = q
556 // --> p + 1/2(lambda_1 N1 + lambda_2 N2) = q
557 // --> (1/2 lambda_1 N1 + 1/2 lambda_2 N2) = q - p
558
2/2
✓ Branch 0 taken 789 times.
✓ Branch 1 taken 263 times.
 1052 for(index_t c=0; c<3; ++c) {
559
1/2
✓ Branch 1 taken 789 times.
✗ Branch 2 not taken.
 789 nlRowScaling(0.5*border_importance);
560
1/2
✓ Branch 1 taken 789 times.
✗ Branch 2 not taken.
 789 nlBegin(NL_ROW);
561
1/2
✓ Branch 1 taken 789 times.
✗ Branch 2 not taken.
 789 nlCoefficient(v1,0.5*N[c]);
562
1/2
✓ Branch 1 taken 789 times.
✗ Branch 2 not taken.
 789 nlCoefficient(v2,0.5*N[c]);
563
1/2
✓ Branch 1 taken 789 times.
✗ Branch 2 not taken.
 789 nlRightHandSide(q[c]-p[c]);
564
1/2
✓ Branch 1 taken 789 times.
✗ Branch 2 not taken.
 789 nlEnd(NL_ROW);
565 }
566
567 }
568 }
569 }
570 }
571
572
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 nlEnd(NL_MATRIX);
573
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 nlEnd(NL_SYSTEM);
574
575
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 nlSolve();
576
577 // Displace each vertex v along Nv[v] by
578 // the solution of the least squares problem
579 // at v
580
2/2
✓ Branch 0 taken 35000 times.
✓ Branch 1 taken 7 times.
 35007 for(index_t v: surface.vertices) {
581 vec3& p = surface.vertices.point(v);
582
1/2
✓ Branch 1 taken 35000 times.
✗ Branch 2 not taken.
 35000 p += nlGetVariable(v)*Nv[v];
583 }
584
585 // (Brutally) project border vertices
586
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 7 if(project_borders && nb_v_on_border != 0 && reference_has_borders) {
587 ✗ for(index_t v: surface.vertices) {
588 ✗ if(v_on_border[v]) {
589 ✗ vec3 p = surface.vertices.point(v);
590 ✗ vec3 q = nearest_along_bidirectional_ray(
591 ✗ border_ribbon_AABB, Ray(p, Nv[v]),
592 ✗ border_distance_factor*max_edge_distance*0.5*(Lv[v]),
593 true
594 );
595 ✗ surface.vertices.point(v) = q;
596 }
597 }
598 }
599
600
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
 7 nlDeleteContext(nlGetCurrent());
601 14 }
602 }
603