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
41			#include <exploragram/basic/common.h>
42			#include <exploragram/hexdom/quadmesher.h>
43			#include <exploragram/hexdom/geometry.h>
44			#include <exploragram/hexdom/mesh_utils.h>
45			#include <geogram/mesh/mesh_io.h>
46			#include <exploragram/hexdom/extra_connectivity.h>
47			#include <geogram/NL/nl.h>
48			#include <geogram/mesh/mesh_geometry.h>
49
50			#include <deque>
51			#include <assert.h>
52			#include <cmath>
53			namespace GEO {
54
55
56			template<class T>
57			T lower_abs_modulo(T value, T module = 1.) {
58			return value - module * std::floor(value / module + .5);
59			}
60			template<class T>
61			T lower_positive_modulo(T value, T module = 1.) {
62			return value - module * std::floor(value / module);
63			}
64
65			typedef vecng<6, double> vec6;
66
67		✗	void add_edge(Mesh* m, vec3 P0, vec3 P1, double value ) {
68		✗	Attribute<double> val(m->edges.attributes(), "val");
69			index_t off_v = m->vertices.create_vertices(2);
70		✗	X(m)[off_v] = P0;
71		✗	X(m)[off_v + 1] = P1;
72		✗	val[m->edges.create_edge(off_v, off_v + 1)] = value;
73		✗	}
74
75		✗	void add_triangle(Mesh* m, vec3 P0, vec3 P1, vec3 P2, double value ) {
76		✗	Attribute<double> val(m->facets.attributes(), "val");
77			index_t off_v = m->vertices.create_vertices(3);
78		✗	X(m)[off_v] = P0; X(m)[off_v + 1] = P1; X(m)[off_v + 2] = P2;
79		✗	val[m->facets.create_triangle(off_v, off_v + 1, off_v + 2)] = value;
80		✗	}
81
82		✗	static mat2 angle_to_mat(double p_alpha, double p_beta) {
83		✗	double rot_angle = (p_beta + p_alpha )/ 2.;
84		✗	double rot_values[4] = { cos(rot_angle), sin(rot_angle), -sin(rot_angle) ,cos(rot_angle)};
85			mat2 rot(rot_values);
86		✗	mat2 inv_rot = rot.inverse();
87
88		✗	double diag_angle = (p_beta - p_alpha) / 2.;
89		✗	double diag_values[4] = { cos(diag_angle), 0,0, sin(diag_angle) };
90			mat2 diag(diag_values);
91		✗	return inv_rot*diag*rot;
92			}
93
94
95			static vec2 operator*(mat2& M, vec2& v) {
96		✗	return vec2(M(0, 0)*v[0] + M(0, 1)*v[1], M(1, 0)*v[0] + M(1, 1)*v[1] );
97			}
98
99
100			struct FF2D {
101		✗	FF2D(Mesh* p_m) : fec(p_m){
102		✗	m = p_m;
103		✗	feature_edge.bind(m->facet_corners.attributes(), "feature_edge");
104		✗	alpha.bind(m->facets.attributes(), "alpha");
105		✗	beta.bind(m->facets.attributes(), "delta");
106		✗	aniso.bind(m->facets.attributes(), "aniso");
107		✗	}
108
109		✗	void init_feature_edge() {
110		✗	FOR(h, m->facet_corners.nb()) {
111		✗	feature_edge[h] = false;
112		✗	if (fec.opposite(h) == NOT_AN_ID) {
113			feature_edge[h] = true;
114		✗	continue;
115			}
116			vec3 n[2] = {
117		✗	Geom::triangle_normal(X(m)[fec.org(fec.prev(h))], X(m)[fec.org(h)], X(m)[fec.dest(h)]),
118		✗	Geom::triangle_normal(X(m)[fec.dest(fec.next(fec.opposite(h)))], X(m)[fec.dest(h)], X(m)[fec.org(h)])
119			};
120		✗	FOR(f, 2) n[f] = normalize(n[f]);
121		✗	if (acos(dot(n[0], n[1])) > M_PI / 3.) {
122			feature_edge[h] = true;
123		✗	feature_edge[fec.opposite(h)] = true;
124			}
125			}
126		✗	}
127
128
129		✗	void local_basis(index_t h, vec3& x, vec3& y, vec3& z) {
130		✗	z = normalize(Geom::mesh_facet_normal(*m, fec.facet(h)));
131		✗	x = normalize(X(m)[fec.dest(h)] - X(m)[fec.org(h)]);
132		✗	y = normalize(cross(z, x));
133		✗	}
134
135		✗	void param_per_triangle(Mesh* debug_mesh) {
136		✗	FOR(f, m->facets.nb()) {
137		✗	vec3 G = Geom::mesh_facet_center(*m, f);
138			vec3 x, y, z;
139		✗	local_basis(m->facets.corner(f, 0), x, y, z);
140		✗	double scale = std::sqrt(Geom::mesh_facet_area(*m, f));
141		✗	FOR(i, 2) {
142		✗	double rot = i?beta[f] : alpha[f];
143		✗	if (rot> 100) continue;
144		✗	vec3 vect = .2*scale*(x*cos(rot) + y*sin(rot));
145		✗	add_edge(debug_mesh, G - vect, G + vect, i);
146			}
147		✗	FOR(i, 32) {
148		✗	vec2 ref(cos(double(i)/5.) , sin(double(i) / 5.));
149			ref = aniso[f]*ref;
150		✗	add_edge(debug_mesh, G , G +0.1*scale*(ref[0]*x+ref[1]*y), i);
151			}
152
153			}
154		✗	}
155
156		✗	vec3 edge_geom(index_t h) {
157		✗	return X(m)[fec.dest(h)] - X(m)[fec.org(h)];
158			}
159		✗	double vector_angle(vec3 v0, vec3 v1) { return atan2(cross(v0, v1).length(), dot(v0, v1)); }
160		✗	double angle_w_r_t_ref(index_t h) {
161		✗	index_t h_ref = m->facets.corner(fec.facet(h),0);
162		✗	if (h == h_ref) return 0;
163		✗	double angle = vector_angle(edge_geom(h), edge_geom(h_ref));
164		✗	if (h == m->facets.corner(fec.facet(h), 2)) return angle;
165		✗	if (h == m->facets.corner(fec.facet(h), 1)) return 2.*M_PI - angle;
166			return angle; // [BL seems to be missing !!]
167			}
168
169		✗	double corner_angle(index_t h) { return vector_angle(X(m)[fec.dest(h)] - X(m)[fec.org(h)], X(m)[fec.org(fec.prev(h))] - X(m)[fec.org(h)]);}
170
171			// if // transport: alpha[fec.facet(fec.opposite(h)) ] = basis_change(h) + alpha[fec.facet(h) ];
172		✗	double basis_change(index_t h) {
173		✗	return M_PI+ angle_w_r_t_ref(h) - angle_w_r_t_ref(fec.opposite(h));
174			}
175
176
177
178		✗	void naive_LS_blur_delta(Mesh* debug_mesh) {
179			geo_argused(debug_mesh);
180		✗	nlNewContext();
181		✗	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
182		✗	nlSolverParameteri(NL_NB_VARIABLES, NLint(3*m->facets.nb()));
183		✗	nlBegin(NL_SYSTEM);
184
185		✗	FOR(f, m->facets.nb()) if (beta[f] < 100) {
186		✗	Matrix<2, double> M = angle_to_mat(alpha[f], beta[f]);
187		✗	nlSetVariable(3 * f, M(0, 0)); nlLockVariable(3 * f);
188		✗	nlSetVariable(3 * f+1, M(0,1)); nlLockVariable(3 * f+1);
189		✗	nlSetVariable(3 * f+2, M(1, 1)); nlLockVariable(3 * f+2);
190			}
191
192		✗	nlBegin(NL_MATRIX);
193		✗	FOR(h, 3 * m->facets.nb()) {
194		✗	if (feature_edge[h]) continue;
195		✗	geo_assert(fec.opposite(h) != NOT_AN_ID);
196		✗	double angle = -basis_change(h);
197		✗	double c = cos(angle);
198		✗	double s = sin(angle);
199
200		✗	index_t f = fec.facet(h);
201		✗	index_t opp = fec.facet(fec.opposite(h));
202		✗	index_t a0 = 3 * f;
203		✗	index_t b0 = 3 * f + 1;
204		✗	index_t c0 = 3 * f + 2;
205		✗	index_t a1 = 3 * opp;
206		✗	index_t b1 = 3 * opp + 1;
207		✗	index_t c1 = 3 * opp + 2;
208		✗	nlBegin(NL_ROW);
209		✗	nlCoefficient(a0, -1);
210		✗	nlCoefficient(a1, c*c);
211		✗	nlCoefficient(b1, -2 * s*c);
212		✗	nlCoefficient(c1, s*s);
213		✗	nlEnd(NL_ROW);
214		✗	nlBegin(NL_ROW);
215		✗	nlCoefficient(b0, -1);
216		✗	nlCoefficient(a1, s*c);
217		✗	nlCoefficient(b1, c*c -s*s);
218		✗	nlCoefficient(c1, -s*c);
219		✗	nlEnd(NL_ROW);
220
221		✗	nlBegin(NL_ROW);
222		✗	nlCoefficient(c0, -1);
223		✗	nlCoefficient(a1, s*s);
224		✗	nlCoefficient(b1, 2.*s*c);
225		✗	nlCoefficient(c1, c*c);
226		✗	nlEnd(NL_ROW);
227
228
229
230			// data fitting term
231
232			// double scale = .1; [BL unused]
233			//nlBegin(NL_ROW); nlCoefficient(a0, -scale); nlRightHandSide(scale); nlEnd(NL_ROW);
234			//nlBegin(NL_ROW); nlCoefficient(b0, -scale); nlRightHandSide(0); nlEnd(NL_ROW);
235			//nlBegin(NL_ROW); nlCoefficient(c0, -scale); nlRightHandSide(scale); nlEnd(NL_ROW);
236
237			}
238		✗	nlEnd(NL_MATRIX);
239		✗	nlEnd(NL_SYSTEM);
240		✗	nlSolve();
241		✗	Attribute<double> A(m->facets.attributes(), "a");
242		✗	Attribute<double> B(m->facets.attributes(), "b");
243		✗	Attribute<double> C(m->facets.attributes(), "c");
244		✗	FOR(f, m->facets.nb()) {
245			mat2 M;
246		✗	M(0, 0) = nlGetVariable(3 * f);
247		✗	M(0, 1) = nlGetVariable(3 * f + 1);
248		✗	M(1, 0) = nlGetVariable(3 * f + 1);
249		✗	M(1, 1) = nlGetVariable(3 * f + 2);
250		✗	aniso[f] = M;
251
252			}
253
254		✗	}
255
256			void naive_LS_smooth() {
257			static const double N = 4.;// N sym dir field... N=4, just change it for debug
258			nlNewContext();
259			nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
260			nlSolverParameteri(NL_NB_VARIABLES, NLint(2 * m->facets.nb()));
261			nlBegin(NL_SYSTEM);
262
263			FOR(h, 3 * m->facets.nb()) {
264			if (!feature_edge[h]) continue;
265			index_t f = fec.facet(h);
266			double angle = -angle_w_r_t_ref(h);//alpha[f] + M_PI / 4. + .5*delta[f];
267			double avt[2] = { cos(angle),sin(angle) };
268			double ap[2];
269			mat2 inv = aniso[f].inverse();
270			mult(inv, avt,ap);
271			angle = atan2(ap[1], ap[0]);
272			nlSetVariable(f * 2, cos(N*angle));
273			nlLockVariable(f * 2);
274			nlSetVariable(f * 2 + 1, sin(N*angle));
275			nlLockVariable(f * 2 + 1);
276			}
277
278			nlBegin(NL_MATRIX);
279			FOR(h, 3 * m->facets.nb()) {
280			if (feature_edge[h]) continue;
281			geo_assert(fec.opposite(h) != NOT_AN_ID);
282
283			double angle = basis_change(h);
284			angle *= N;
285			double rot[2][2] = { { cos(angle),sin(angle) },{ -sin(angle),cos(angle) } };
286			FOR(d, 2) {
287			nlBegin(NL_ROW);
288			nlCoefficient(fec.facet(h) * 2 + d, -1);
289			FOR(dd, 2) nlCoefficient(fec.facet(fec.opposite(h)) * 2 + dd, rot[d][dd]);
290			nlEnd(NL_ROW);
291			}
292			}
293			nlEnd(NL_MATRIX);
294			nlEnd(NL_SYSTEM);
295			nlSolve();
296			FOR(f, m->facets.nb()) {
297			double angle = (1. / N)*atan2(nlGetVariable(f * 2 + 1), nlGetVariable(f * 2));
298			double avt[2] = { cos(angle),sin(angle) };
299			double ap[2];
300			mult(aniso[f], avt, ap);
301			angle = atan2(ap[1], ap[0]);
302
303			alpha[f] = (1. / N)*atan2(nlGetVariable(f * 2 + 1), nlGetVariable(f * 2));
304			}
305			}
306
307
308		✗	void apply(Mesh* debug_mesh) {
309		✗	init_feature_edge();
310		✗	FOR(f, m->facets.nb()) alpha[f] = beta[f] = 1000;
311
312			// for each fature edge:
313			// -> compute the #angu geom to the next feature edge around vertex
314			// -> define #angu in map
315			// -> if #angu in map est impair:
316			// => alpha = #angu in map - #angu geom
317			// => paralell transport first constraint
318			// ===> RESULT : fix an angle + possible delta a some triangles
319
320		✗	FOR(h, 3*m->facets.nb()) {
321		✗	if (!feature_edge[h]) continue;
322		✗	if (beta[fec.facet(h)]<20) continue; // already constrained by previous halfedge
323		✗	alpha[fec.facet(h)] = -angle_w_r_t_ref(h);
324
325			// find all corners sharing org(h)
326			vector<index_t> edges;
327			edges.push_back(h);
328		✗	while (!feature_edge[fec.prev(edges.back())])
329		✗	edges.push_back(fec.opposite(fec.prev(edges.back())));
330
331			// compute their angles
332			double sum = 0;
333		✗	FOR(i, edges.size()) sum+= corner_angle(edges[i]);
334		✗	double wanted_sum = 0.5*M_PI * floor(sum / (0.5*M_PI) + .5);
335		✗	if (wanted_sum < .1) wanted_sum = 0.5*M_PI;
336		✗	int nb_angu = int(std::floor(wanted_sum / (0.5*M_PI)));
337
338
339			// parallel transport
340		✗	FOR(i, edges.size()-1) alpha[fec.facet(edges[i+1])] = basis_change(fec.prev(edges[i])) + alpha[fec.facet(edges[i ])];
341
342		✗	if (nb_angu % 2 == 0)
343		✗	FOR(i, edges.size() - 1) alpha[fec.facet(edges[i + 1])] += (sum - wanted_sum)*double(i + 1) / double(edges.size() - 1);
344		✗	else FOR(i, edges.size())
345		✗	beta[fec.facet(edges[i])] = alpha[fec.facet(edges[i])] + M_PI / 2. +(sum - wanted_sum);
346
347			//plop(edges.size());
348			//plop(wanted_sum);
349			}
350		✗	naive_LS_blur_delta(debug_mesh);
351
352			//naive_LS_smooth();
353
354			//FOR(f, m->facets.nb()) alpha[f] = delta[f] = 1000;
355			//FOR(f, m->facets.nb()) FOR(lv, 3) {
356			// index_t h = m->facets.corner(f, lv);
357			// if (feature_edge[h]) alpha[f] = -angle_w_r_t_ref(h);
358			//}
359		✗	}
360
361
362
363
364
365
366			FacetsExtraConnectivity fec;
367			Attribute<bool> feature_edge; // halfedge boolean
368			Attribute<double> alpha;
369			Attribute<double> beta;
370			Attribute<mat2> aniso;
371			Mesh* m;
372			};
373
374
375		✗	void current_test(Mesh* m, Mesh* debug_mesh) {
376		✗	FF2D ff2d(m);
377		✗	ff2d.apply(debug_mesh);
378		✗	ff2d.param_per_triangle(debug_mesh);
379		✗	return;
380
381		✗	}
382			}
383