GCC Code Coverage Report

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 <exploragram/hexdom/PGP.h>
41			#include <exploragram/hexdom/frame.h>
42			#include <exploragram/hexdom/extra_connectivity.h>
43			#include <exploragram/hexdom/geometry.h>
44			#include <exploragram/hexdom/mesh_utils.h>
45			#include <exploragram/hexdom/time_log.h>
46
47
48
49
50			#include <exploragram/hexdom/quadmesher.h>
51			#include <geogram/NL/nl.h>
52
53			#include <algorithm>
54			#include <cmath>
55			#include <queue>
56
57			namespace GEO {
58
59		✗	PGPopt::PGPopt(Mesh* p_m) : m(p_m) {
60			// init fast acces to edges from vertices
61		✗	if (m->edges.nb() == 0) compute_tet_edge_graph(m, v2e, false);
62		✗	else restore_v2e(m, v2e);
63		✗	v2eopp = vector<vector<index_t> >(m->edges.nb(), vector<index_t>());
64		✗	FOR(e, m->edges.nb()) v2eopp[m->edges.vertex(e, 1)].push_back(e);
65
66			//bind attributes
67		✗	U.bind(m->vertices.attributes(), "U");
68		✗	B.bind(m->vertices.attributes(), "B");
69		✗	corr.bind(m->edges.attributes(), "corr");
70		✗	tij.bind(m->edges.attributes(), "tij");
71		✗	}
72
73		✗	bool PGPopt::is_PGP_singular(index_t c, index_t lf) {
74			index_t v[3];
75		✗	FOR(i, 3) v[i] = m->cells.facet_vertex(c, lf, i);
76			vec3i t(0, 0, 0);
77			mat3 R;
78			R.load_identity();
79		✗	FOR(i, 3) {
80		✗	AxisPermutation r = Rij(m, B, v[i], v[(i + 1) % 3]);
81			bool inv;
82		✗	index_t e = edge_from_vertices(v[i], v[(i + 1) % 3], inv);
83		✗	int i1 = int(tij[e][0]);
84		✗	int i2 = int(tij[e][1]);
85		✗	int i3 = int(tij[e][2]);
86			vec3i ltij(i1,i2,i3);
87		✗	t += R*(inv ? -(r.inverse()*ltij) : ltij);
88		✗	R = R* r.inverse().get_mat();
89			}
90
91		✗	geo_assert(R.is_identity());
92		✗	return t[0] || t[1] || t[2];
93			}
94
95
96		✗	bool PGPopt::tet_is_PGP_singular_fct(index_t t) {
97			bool is_sing = false;
98		✗	FOR(f, 4) is_sing = is_sing || is_PGP_singular(t, f);
99		✗	return is_sing;
100			}
101
102
103		✗	bool PGPopt::face_is_resp(index_t c, index_t lf) {
104		✗	if (m->cells.adjacent(c, lf) == NO_CELL) return true;
105		✗	return m->cells.adjacent(c, lf) < c;
106			}
107
108			struct TriangleEdgesInSameBasis {
109			bool try_new_triangle(PGPopt* pgp, index_t c, index_t lf) {
110			if (!pgp->face_is_resp(c, lf)) { return false; } // avoid doing twice the same work
111			if (triangle_is_frame_singular(pgp->m, pgp->B, c, lf)) { return false; } // curl correction on singular face is meaningless
112			// without chain basis change
113			index_t vid[3]; // vertices index
114			// init vertex indices
115			FOR(e, 3) vid[e] = pgp->m->cells.facet_vertex(c, lf, e);
116			// init edges
117			FOR(e, 3) {
118			edge[e] = pgp->edge_from_vertices(vid[e], vid[next_mod(e, 3)], inv[e]);
119			geo_assert(edge[e] != NOT_AN_ID);
120			}
121			// express corr of all edges in a common basis with edge_ap[e]
122			FOR(e, 3) {
123			if (inv[e]) edge_ap[e] = Rij(pgp->m, pgp->B, vid[0], vid[next_mod(e, 3)]);
124			else edge_ap[e] = Rij(pgp->m, pgp->B, vid[0], vid[e]);
125			}
126			return true;
127			}
128
129
130			index_t edge[3];
131			bool inv[3];
132			AxisPermutation edge_ap[3];
133			};
134
135
136
137
138			// ___ __
139			// / _ \ _ _ ___ / _|___ _ _ _ __ ___
140			// | (_) | ' \/ -_) | _/ _ \ '_| ' \(_-<
141			// \___/|_||_\___| |_| \___/_| |_|_|_/__/
142			//
143
144
145
146		✗	vec3 PGPopt::wish_angle_corr(index_t e, bool inv) {
147			vec3 c(0, 0, 0);
148			if (corr.is_bound()) { // CubeCover ne l'utilise pas forcement
149		✗	c = corr[e];
150		✗	if (inv) {
151		✗	AxisPermutation chg = Rij(m, B, m->edges.vertex(e, 0), m->edges.vertex(e, 1));
152		✗	c = -(chg.inverse() *c);
153			}
154			}
155		✗	return c;
156			}
157
158		✗	vec3 PGPopt::wish_angle_edge_geom(index_t e, bool inv) {
159		✗	index_t org = m->edges.vertex(e, 0);
160			index_t dest = m->edges.vertex(e, 1);
161		✗	if (inv) std::swap(org, dest);
162
163		✗	AxisPermutation ap = Rij(m, B, org, dest);
164
165		✗	mat3 frame = B[org] + Frame(B[dest]).apply_permutation(ap);
166			frame *= 0.5;
167		✗	frame = invert_columns_norm(frame);
168		✗	vec3 angle = frame.transpose() * (m->vertices.point(dest) - m->vertices.point(org));
169		✗	return 2.*M_PI *angle; // one cycle length is edgelength_
170			}
171
172		✗	vec3 PGPopt::wish_angle(index_t e, bool inv) {
173		✗	return wish_angle_edge_geom(e, inv) + PGPopt::wish_angle_corr(e, inv);
174			}
175
176
177
178
179
180			// _ _ _
181			// ___ _ __| |_(_)_ __ (_)______ __ ___ _ _ _ _
182			// / _ \ '_ \ _| | ' \| |_ / -_) / _/ _ \ '_| '_|
183			// \___/ .__/\__|_|_|_|_|_/__\___|_\__\___/_| |_|
184			// |_|
185
186
187			///////////////////////////////////////////////////
188			//optimize_corr
189			///////////////////////////////////////////////////
190			// problem: the objective one form derived from the frame field is not close, leading to PGP singularities (T_junctions)
191			// solution: this function computes a one form such that :
192			// *its norm is minimal
193			// *adding it to the original objective one form, make it close... everywhere but on frame field singularities
194			// param : max_corr_prop typically in [0,.5] (0 => no correction, .5 => take initial field into account for only 50%
195			//
196			// note: it is the same energy as in cubecover, exept that there is no integer constraints, and the frame singularities are ignored
197
198
199		✗	void PGPopt::optimize_corr(double max_corr_prop) {
200		✗	FOR(e, m->edges.nb()) corr[e] = vec3(0, 0, 0);
201		✗	if (max_corr_prop == 0) {
202			return;
203			}
204		✗	cubcover(true);
205			// clamp the result: very naive way to avoid too large corrections. May be improved.
206		✗	FOR(e, m->edges.nb()) {
207			// BasisChg chg = edge_basis_change(e, false);
208		✗	vec3 geom_form = wish_angle_edge_geom(e, false);
209		✗	vec3 corr_form = wish_angle_corr(e, false);
210			double scale = 1.;
211			double cl = corr_form.length();
212			double gl = geom_form.length();
213		✗	if (cl > max_corr_prop* gl) scale = max_corr_prop* gl / cl;
214		✗	FOR(d, 3) corr[e][d] = scale * corr[e][d];
215			}
216			}
217
218
219
220		✗	void PGPopt::move_U_to_corner() {
221		✗	Attribute<vec3> UC(m->cell_corners.attributes(), "U");
222		✗	Attribute<bool> has_param(m->cell_facets.attributes(), "has_param");
223
224		✗	FOR(c, m->cells.nb()) {
225
226			bool all_faces_have_param = true;
227			// flag has param
228		✗	FOR(lf, 4) {
229		✗	index_t f = m->cells.facet(c, lf);
230			has_param[f] = true;
231		✗	if (triangle_is_frame_singular(m, B, c, lf))has_param[f] = false;
232		✗	else if (is_PGP_singular(c, lf)) has_param[f] = false;
233		✗	all_faces_have_param = all_faces_have_param && has_param[f];
234			}
235
236
237			index_t org = 0;
238
239			// find a seed that will properly reconstruct all valid triangle param
240		✗	if (!all_faces_have_param) {
241		✗	bool can_be_ref[4] = { true, true, true, true };
242		✗	FOR(lf, 4) {
243		✗	index_t f = m->cells.facet(c, lf);
244		✗	if (has_param[f]) {
245		✗	bool local_can_be_ref[4] = { false, false, false, false };
246		✗	FOR(lv, 3) local_can_be_ref[m->cells.descriptor(c).facet_vertex[lf][lv]] = true;
247		✗	FOR(i, 4) can_be_ref[i] = can_be_ref[i] && local_can_be_ref[i];
248			}
249			}
250			org = NOT_AN_ID;
251		✗	FOR(i, 4) if (can_be_ref[i]) org = i;
252			}
253
254
255		✗	if (org != NOT_AN_ID) FOR(i, 4) {
256		✗	index_t corner = m->cells.corner(c, i);
257		✗	UC[corner] = U[m->cells.vertex(c, i)];
258		✗	if (i != org) {
259		✗	AxisPermutation change = Rij(m, B, m->cells.vertex(c, org), m->cells.vertex(c, i));
260		✗	UC[corner] = change.inverse() * U[m->cells.vertex(c, i)];
261			bool inv;
262		✗	index_t e = edge_from_vertices(m->cells.vertex(c, org), m->cells.vertex(c, i), inv);
263			geo_assert(e != NOT_AN_ID);
264		✗	vec3 t2 = !inv ? tij[e] : -(change.inverse()*tij[e]);
265			UC[corner] += vec3(t2[0], t2[1], t2[2]);
266
267			}
268			}
269		✗	else geo_assert(!has_param[m->cells.facet(c, 0)] && !has_param[m->cells.facet(c, 1)] && !has_param[m->cells.facet(c, 2)] && !has_param[m->cells.facet(c, 3)]);
270			}
271
272			// TODO REMOVE THAT !!!!! just to TEST
273			//FOR(c, m->cell_corners.nb()) UC[c] = 2*UC[c];
274		✗	}
275
276			/*
277			* _ _ _ ___ ___ ___
278			* ___ _ __| |_(_)_ __ (_)______ | _ \/ __| _ \
279			* / _ \ '_ \ _| | ' \| |_ / -_) | _/ (_ | _/
280			* \___/ .__/\__|_|_|_|_|_/__\___|_|_| \___|_|
281			* |_|
282			*/
283
284		✗	void PGPopt::optimize_PGP() {
285		✗	Attribute<vec3> lockU(m->vertices.attributes(), "lockU");
286			// Create and initialize OpenNL context
287		✗	nlNewContext();
288		✗	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
289		✗	nlSolverParameteri(NL_NB_VARIABLES, NLint(6 * m->vertices.nb()));
290
291		✗	nlBegin(NL_SYSTEM);
292		✗	FOR(v, m->vertices.nb())FOR(d, 3) if (std::fabs(lockU[v][d]) > 0)FOR(c, 2) {
293		✗	nlSetVariable(6 * v + 2 * d + c, 1 - c);
294		✗	nlLockVariable(6 * v + 2 * d + c);
295			}
296
297		✗	nlBegin(NL_MATRIX);
298		✗	FOR(e, m->edges.nb()) {
299		✗	AxisPermutation ap = Rij(m, B, m->edges.vertex(e, 0), m->edges.vertex(e, 1));
300
301		✗	vec3 theta = wish_angle(e, false);
302		✗	FOR(d, 3) {
303		✗	double c = cos(theta[d]);
304		✗	double s = sin(theta[d]);
305		✗	index_t off0 = 6 * m->edges.vertex(e, 0) + 2 * d;
306
307		✗	nlBegin(NL_ROW);
308		✗	FOR(dd, 3) if (ap.get_mat()(dd, d) != 0)
309		✗	nlCoefficient(6 * m->edges.vertex(e, 1) + 2 * dd, -1.);
310		✗	nlCoefficient(off0, c);
311		✗	nlCoefficient(off0 + 1, s);
312		✗	nlEnd(NL_ROW);
313		✗	nlBegin(NL_ROW);
314		✗	FOR(dd, 3)
315		✗	nlCoefficient(6 * m->edges.vertex(e, 1) + 2 * dd + 1, -ap.get_mat()(dd, d));
316		✗	nlCoefficient(off0, -s);
317		✗	nlCoefficient(off0 + 1, c);
318		✗	nlEnd(NL_ROW);
319			}
320			}
321
322		✗	nlEnd(NL_MATRIX);
323		✗	nlEnd(NL_SYSTEM);
324			// Solve and get solution
325		✗	nlSolve();
326
327		✗	FOR(v, m->vertices.nb()) FOR(d, 3)
328		✗	U[v][d] = (.5 / M_PI) * atan2(nlGetVariable(6 * v + 2 * d + 1), nlGetVariable(6 * v + 2 * d));
329
330		✗	nlDeleteContext(nlGetCurrent());
331
332
333		✗	FOR(e, m->edges.nb()) {
334			index_t i = m->edges.vertex(e, 0);
335			index_t j = m->edges.vertex(e, 1);
336		✗	AxisPermutation rij = Rij(m, B, i, j);
337
338		✗	vec3 gij = wish_angle(e, false) / (2.*M_PI);
339
340		✗	FOR(d, 3) {
341		✗	tij[e][d] = int(round(-(rij.inverse().get_mat()*U[j])[d] - gij[d] + U[i][d]));
342			}
343
344			}
345		✗	move_U_to_corner();
346		✗	}
347
348
349
350
351
352
353			// _____ _ _____
354			// / ____| | | / ____|
355			// | | _ _| |__ ___| | _____ _____ _ __
356			// | | | | | | '_ \ / _ \ | / _ \ \ / / _ \ '__|
357			// | |___| |_| | |_) | __/ |___| (_) \ V / __/ |
358			// \_____\__,_|_.__/ \___|\_____\___/ \_/ \___|_|
359			//
360			//
361			struct TetHalfedge {
362		✗	TetHalfedge(index_t p_cell, index_t p_org, index_t p_dest) {
363		✗	cell = p_cell; org = p_org; dest = p_dest;
364			}
365			index_t cell;
366			index_t org;
367			index_t dest;
368			};
369		✗	index_t PGPopt::get_non_nulledge(index_t c, index_t cf, Attribute<bool>& nulledge){
370			int nbzero = 0;
371		✗	FOR(cfv, 3) {
372			bool inv;
373		✗	index_t e = edge_from_vertices(m->cells.facet_vertex(c, cf, cfv), m->cells.facet_vertex(c, cf, next_mod(cfv, 3)), inv);
374			geo_assert(e != NOT_AN_ID);
375		✗	if (nulledge[e]) nbzero++;
376			}
377		✗	if (nbzero == 2 && !triangle_is_frame_singular(m, B, c, cf)) {
378		✗	FOR(cfv, 3) {
379			bool inv;
380		✗	index_t e = edge_from_vertices(m->cells.facet_vertex(c, cf, cfv), m->cells.facet_vertex(c, cf, next_mod(cfv, 3)), inv);
381			geo_assert(e != NOT_AN_ID);
382		✗	if (!nulledge[e]) {
383			nulledge[e] = true;
384		✗	return cfv;
385			}
386			}
387			}
388			return NOT_AN_ID;
389			}
390
391
392
393
394		✗	void PGPopt::mark_null_edges(Attribute<bool>& nulledge) {
395
396		✗	Attribute<vec3> lockU(m->vertices.attributes(), "lockU");
397		✗	index_t alot = std::numeric_limits<index_t>::max() - 10;
398		✗	FOR(e, m->edges.nb()) nulledge[e] = false;
399
400			vector<index_t> offset_from_org;
401			vector<index_t> dest;
402
403			// compute covering tree
404		✗	cell_edges_in_RCS(m, offset_from_org, dest);
405		✗	std::vector<index_t > dist(m->vertices.nb(), alot);
406
407			{
408		✗	FOR(seed, m->vertices.nb()) {
409		✗	if (dist[seed] != alot) continue;
410		✗	dist[seed] = 0;
411			std::queue<index_t> queue;
412			queue.push(seed);
413		✗	while (!queue.empty()) {
414		✗	index_t v = queue.front();
415			queue.pop();
416		✗	for (index_t ed = offset_from_org[v]; ed < offset_from_org[v + 1]; ed++) {
417		✗	index_t vopp = dest[ed];
418		✗	if (dist[vopp] == alot) {
419			bool inv;
420		✗	index_t e = edge_from_vertices(v, vopp, inv);
421			geo_assert(e != NOT_AN_ID);
422			nulledge[e] = true;
423		✗	dist[vopp] = dist[v] + 1;
424			queue.push(vopp);
425			}
426			}
427			}
428			}
429			}
430
431			std::queue<TetHalfedge> queue;
432			// init
433		✗	FOR(c, m->cells.nb()) FOR(cf, 4) {
434		✗	index_t cfe = get_non_nulledge(c, cf, nulledge);
435		✗	if (cfe != NOT_AN_ID)
436		✗	queue.push(TetHalfedge(c, m->cells.facet_vertex(c, cf, cfe), m->cells.facet_vertex(c, cf, next_mod(cfe, 3))));
437			}
438		✗	while (!queue.empty()) {
439		✗	TetHalfedge cur = queue.front();
440			queue.pop();
441			index_t cir = cur.cell;
442			do {
443		✗	FOR(cf, 4) {
444		✗	index_t cfe = get_non_nulledge(cir, cf, nulledge);
445		✗	if (cfe != NOT_AN_ID)
446		✗	queue.push(TetHalfedge(cir, m->cells.facet_vertex(cir, cf, cfe), m->cells.facet_vertex(cir, cf, next_mod(cfe, 3))));
447			}
448		✗	cir = next_cell_around_oriented_edge(m, cir, cur.org, cur.dest);
449
450		✗	if (cir == NOT_AN_ID) {
451			cir = cur.cell;
452		✗	while (next_cell_around_oriented_edge(m, cir, cur.dest, cur.org) != NOT_AN_ID)
453		✗	cir = next_cell_around_oriented_edge(m, cir, cur.dest, cur.org);
454			}
455		✗	} while (cir != cur.cell);
456			}
457
458
459		✗	}
460
461			// // _____ _ _____
462			// // / ____| | | / ____|
463			// // | | _ _| |__ ___| | _____ _____ _ __
464			// // | | | | | | '_ \ / _ \ | / _ \ \ / / _ \ '__|
465			// // | |___| |_| | |_) | __/ |___| (_) \ V / __/ |
466			// // \_____\__,_|_.__/ \___|\_____\___/ \_/ \___|_|
467			// //
468			// //
469
470		✗	void PGPopt::grow_ball(Attribute<bool>& tet_in_ball) {
471
472		✗	vector<bool> touched(m->vertices.nb(), false);
473		✗	vector<bool> tet_is_singular(m->cells.nb(), false);
474		✗	FOR(c, m->cells.nb()) FOR(lf, 4) tet_is_singular[c] = tet_is_singular[c] || triangle_is_frame_singular(m, B, c, lf);
475
476		✗	FOR(c, m->cells.nb()) tet_in_ball[c] = false;
477		✗	FOR(seed, m->cells.nb()) {
478		✗	if (tet_in_ball[seed]) continue;
479		✗	if (tet_is_singular[seed]) continue;
480		✗	if (touched[m->cells.vertex(seed, 0)]) continue;
481		✗	if (touched[m->cells.vertex(seed, 1)]) continue;
482		✗	if (touched[m->cells.vertex(seed, 2)]) continue;
483		✗	if (touched[m->cells.vertex(seed, 3)]) continue;
484
485			tet_in_ball[seed] = true;
486		✗	FOR(lv, 4) touched[m->cells.vertex(seed, lv)] = true;
487
488			std::queue<index_t> queue;
489			queue.push(seed);
490		✗	while (!queue.empty()) {
491		✗	index_t c = queue.front();
492			queue.pop();
493		✗	FOR(lf, 4) {
494		✗	index_t copp = m->cells.adjacent(c, lf);
495		✗	if (copp == NOT_AN_ID) continue;
496		✗	if (tet_in_ball[copp]) continue;
497		✗	if (tet_is_singular[copp]) continue;
498
499			// try to reach a new vertex
500			bool all_are_fine = true;
501		✗	FOR(lvopp, 4) {
502			index_t vopp = m->cells.vertex(copp, lvopp);
503			bool is_in_c = false;
504		✗	FOR(lv, 4) is_in_c = is_in_c || (vopp == m->cells.vertex(c, lv));
505		✗	all_are_fine = all_are_fine && (is_in_c || !touched[vopp]);
506			}
507			// try to glue to two faces
508			int nb_done_neigs = 0;
509		✗	FOR(lfopp, 4)
510		✗	if (m->cells.adjacent(copp, lfopp) != NOT_AN_ID)
511		✗	if (tet_in_ball[m->cells.adjacent(copp, lfopp)])
512		✗	nb_done_neigs++;
513		✗	if (!all_are_fine && nb_done_neigs ==1) continue;
514
515		✗	FOR(lvopp, 4) touched[m->cells.vertex(copp, lvopp)] = true;
516			tet_in_ball[copp] = true;
517			queue.push(copp);
518			}
519			}
520			}
521		✗	Attribute<bool> null_edge(m->edges.attributes(), "null");
522		✗	FOR(e, m->edges.nb()) null_edge[e] = false;
523		✗	FOR(c, m->cells.nb()) if (tet_in_ball[c]) {
524			bool inv;
525		✗	FOR(lv, 3)for (index_t lv2 = lv + 1; lv2 < 4; lv2++) FOR(d, 3)
526		✗	null_edge[edge_from_vertices(m->cells.vertex(c, lv), m->cells.vertex(c, lv2), inv)] = true;
527			}
528		✗	}
529
530
531		✗	void PGPopt::cubcover(bool compute_only_corr) {
532			geo_argused(compute_only_corr);
533		✗	Logger::warn("PGP") << "Cubecover not implemented in the public version"
534			<< std::endl;
535		✗	}
536
537			}
538