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/frame.h>
41
42			namespace GEO {
43
44		✗	mat3 normalize_columns(const mat3& B) {
45		✗	mat3 res;
46		✗	vec3 n;
47		✗	FOR(j, 3) n[j] = col(B, j).length(); //= std::sqrt(pow(B(0, j), 2) + pow(B(1, j), 2) + pow(B(2, j), 2));
48		✗	FOR(j, 3) geo_assert(n[j] > 1e-20);
49		✗	FOR(i, 3)FOR(j, 3) res(i, j) = B(i, j) / n[j];
50		✗	return res;
51			}
52
53		✗	mat3 invert_columns_norm(const mat3& B) {
54		✗	mat3 res;
55		✗	vec3 n;
56		✗	FOR(j, 3) n[j] = col(B, j).length2();// pow(B(0, j), 2) + pow(B(1, j), 2) + pow(B(2, j), 2);
57		✗	FOR(j, 3) geo_assert(n[j] > 1e-20);
58		✗	FOR(i, 3)FOR(j, 3) res(i, j) = B(i, j) / n[j];
59		✗	return res;
60			}
61
62			/**************************************************************************************************/
63
64		✗	mat3 rotx(double angle) {
65		✗	double c = cos(angle);
66		✗	double s = sin(angle);
67		✗	mat3 res = mat3_from_coeffs(
68			1,0,0,
69			0,c,-s,
70			0,s,c
71			);
72		✗	return res;
73			}
74
75		✗	mat3 roty(double angle) {
76		✗	double c = cos(angle);
77		✗	double s = sin(angle);
78		✗	mat3 res = mat3_from_coeffs(
79			c, 0, s,
80			0, 1, 0,
81			-s, 0, c
82			);
83		✗	return res;
84			}
85
86		✗	mat3 rotz(double angle) {
87		✗	double c = cos(angle);
88		✗	double s = sin(angle);
89		✗	mat3 res = mat3_from_coeffs(
90			c, -s, 0,
91			s, c, 0,
92			0, 0, 1
93			);
94		✗	return res;
95			}
96
97			// non optimized version is "return rotz(xyz[2]) *roty(xyz[1]) *rotx(xyz[0]);"
98		✗	mat3 euler_to_mat3(vec3 xyz) {
99		✗	double ca = cos(xyz[0]), sa = sin(xyz[0]);
100		✗	double cb = cos(xyz[1]), sb = sin(xyz[1]);
101		✗	double cg = cos(xyz[2]), sg = sin(xyz[2]);
102		✗	mat3 res = mat3_from_coeffs(
103		✗	cb*cg, cg*sa*sb - ca*sg, ca*cg*sb + sa*sg,
104		✗	cb*sg, sa*sb*sg + ca*cg, ca*sb*sg - cg*sa,
105			-sb, cb*sa, ca*cb
106			);
107		✗	return res;
108			}
109
110		✗	vec3 mat3_to_euler(const mat3& r) {//http://www.staff.city.ac.uk/~sbbh653/publications/euler.pdf
111		✗	vec3 res;
112		✗	if (std::abs(std::abs(r(2, 0)) - 1) > 1e-5) {
113		✗	res[1] = -asin(r(2, 0));
114		✗	res[0] = atan2(r(2, 1), r(2, 2));
115		✗	res[2] = atan2(r(1, 0), r(0, 0));
116			} else {
117		✗	res[2] = 0;
118		✗	if (std::abs(r(2, 0) + 1) < 1e-5) {
119		✗	res[1] = M_PI / 2.;
120		✗	res[0] = atan2(r(0, 1), r(0, 2));
121			} else {
122		✗	res[1] = -M_PI / 2.;
123		✗	res[0] = atan2(-r(0, 1), -r(0, 2));
124			}
125			}
126		✗	return res;
127			}
128
129
130			/***********************************************************************************************************************************/
131
132			//void generate_AxisPermutations_init_code() {
133			// mat3 r90z = mat3_from_coeffs(0, -1, 0, 1, 0, 0, 0, 0, 1);
134			// FOR(inv, 2) FOR(perm, 3) FOR(rotz, 4) {
135			// mat3 m;
136			// m.load_identity();
137			// if (inv == 1) m = mat3_from_coeffs(0, 0, -1, 0, -1, 0, -1, 0, 0);
138
139			// FOR(p, perm) FOR(d, 3) {
140			// double tmp = m(d, 0);
141			// m(d, 0) = m(d, 1);
142			// m(d, 1) = m(d, 2);
143			// m(d, 2) = tmp;
144			// }
145			// FOR(r, rotz) m = r90z*m;
146			// FOR(i, 3)FOR(j, 3) m(i, j) = std::floor(m(i, j) + 0.5);
147			// AxisPermutations[12 * inv + 4 * perm + rotz] = m;
148			// }
149			// FOR(i, 24)FOR(j, 24) if ((AxisPermutations[i] * AxisPermutations[j]).is_identity()) AxisPermutations_inv[i] = j;
150
151
152			// std::cerr << "static mat3 AxisPermutations[24] = { ";
153			// FOR(i, 24) {
154			// double* ptr = AxisPermutations[i].data();
155			// std::cerr << "mat3_from_coeffs(";
156			// std::cerr << ptr[0];
157			// FOR(d, 8) std::cerr << "," << ptr[d + 1];
158			// if (i<23)std::cerr << ")," << ((i % 4) ? "" : "\n");
159			// else std::cerr << ")};\n";
160			// }
161
162			// std::cerr << "static int AxisPermutations_inv[24] = { ";
163			// FOR(i, 24) {
164			// if (i<23)std::cerr << AxisPermutations_inv[i] << ",";
165			// else std::cerr << AxisPermutations_inv[i] << "};\n";
166			// }
167			//};
168
169			static mat3 AxisPermutations[24] = {
170			mat3_from_coeffs(1,0,0,0,1,0,0,0,1),
171			mat3_from_coeffs(0,-1,0,1,0,0,0,0,1),mat3_from_coeffs(-1,0,0,0,-1,0,0,0,1),mat3_from_coeffs(0,1,0,-1,0,0,0,0,1),mat3_from_coeffs(0,0,1,1,0,0,0,1,0),
172			mat3_from_coeffs(-1,0,0,0,0,1,0,1,0),mat3_from_coeffs(0,0,-1,-1,0,0,0,1,0),mat3_from_coeffs(1,0,0,0,0,-1,0,1,0),mat3_from_coeffs(0,1,0,0,0,1,1,0,0),
173			mat3_from_coeffs(0,0,-1,0,1,0,1,0,0),mat3_from_coeffs(0,-1,0,0,0,-1,1,0,0),mat3_from_coeffs(0,0,1,0,-1,0,1,0,0),mat3_from_coeffs(0,0,-1,0,-1,0,-1,0,0),
174			mat3_from_coeffs(0,1,0,0,0,-1,-1,0,0),mat3_from_coeffs(0,0,1,0,1,0,-1,0,0),mat3_from_coeffs(0,-1,0,0,0,1,-1,0,0),mat3_from_coeffs(0,-1,0,-1,0,0,0,0,-1),
175			mat3_from_coeffs(1,0,0,0,-1,0,0,0,-1),mat3_from_coeffs(0,1,0,1,0,0,0,0,-1),mat3_from_coeffs(-1,0,0,0,1,0,0,0,-1),mat3_from_coeffs(-1,0,0,0,0,-1,0,-1,0),
176			mat3_from_coeffs(0,0,1,-1,0,0,0,-1,0),mat3_from_coeffs(1,0,0,0,0,1,0,-1,0),mat3_from_coeffs(0,0,-1,1,0,0,0,-1,0)
177			};
178
179			static index_t AxisPermutations_inv[24] = { 0,3,2,1,8,5,15,22,4,14,21,11,12,23,9,6,16,17,18,19,20,10,7,13 };
180
181		✗	void AxisPermutation::aligns_B_wrt_ref(mat3 ref, mat3 B) {
182		✗	ref= normalize_columns(ref);
183		✗	B = normalize_columns(B);
184		✗	index_t best_i = 0;
185		✗	double best_score = 1e20;
186		✗	FOR(i, 24) {
187		✗	double score = Frobenius_norm(B*AxisPermutations[i]-ref);
188		✗	if (score < best_score) {
189		✗	best_i = i;
190		✗	best_score = score;
191			}
192			}
193		✗	mid = best_i;
194		✗	}
195
196		✗	void AxisPermutation::make_col2_equal_to_z(mat3 B, vec3 z) {
197		✗	B = normalize_columns(B);
198		✗	index_t best_i = 0;
199		✗	double best_score = 1e20;
200		✗	FOR(i, 24) {
201		✗	double score = (col(B*AxisPermutations[i],2) - z).length2();
202		✗	if (score < best_score) {
203		✗	best_i = i;
204		✗	best_score = score;
205			}
206			}
207		✗	mid = best_i;
208		✗	}
209
210		✗	const mat3& AxisPermutation::get_mat() const {
211		✗	return AxisPermutations[mid];
212			}
213
214		✗	AxisPermutation AxisPermutation::inverse() {
215		✗	return AxisPermutation(AxisPermutations_inv[mid]);
216			}
217
218			/***********************************************************************************************************************************/
219
220		✗	void Frame::make_z_equal_to(vec3 z) {
221		✗	z = normalize(z);
222		✗	vec3 x(1, 0, 0);
223		✗	vec3 y = cross(z, x);
224		✗	if (y.length2() < .1) {
225		✗	x = vec3(0, 1, 0);
226		✗	y = cross(z, x);
227			}
228		✗	y = normalize(y);
229		✗	x = cross(y, z);
230		✗	FOR(d, 3) r(d, 0) = x[d];
231		✗	FOR(d, 3) r(d, 1) = y[d];
232		✗	FOR(d, 3) r(d, 2) = z[d];
233		✗	}
234
235		✗	mat3 Frame::average_frame(vector<mat3>& data) {
236		✗	SphericalHarmonicL4 sum;
237		✗	FOR(i, data.size()) {
238		✗	SphericalHarmonicL4 nv;
239		✗	nv[4] = std::sqrt(7. / 12.);
240		✗	nv[8] = std::sqrt(5. / 12.);
241		✗	nv.euler_rot(mat3_to_euler(data[i]));
242		✗	sum = sum + nv;
243			}
244		✗	sum = sum*(1. / sum.norm());
245		✗	return sum.project_mat3();
246			}
247
248		✗	mat3 Frame::representative_frame(vector<vec3>& bunch_of_vectors, vector<double>& w) {
249		✗	SphericalHarmonicL4 sum;
250		✗	FOR(i,bunch_of_vectors.size()) {
251		✗	vec3 n = normalize(bunch_of_vectors[i]);
252		✗	mat3 r;
253		✗	Frame(r).make_z_equal_to(n);
254		✗	SphericalHarmonicL4 nv;
255		✗	nv[4] = std::sqrt(7. / 12.);
256		✗	nv.euler_rot(mat3_to_euler(r));
257		✗	nv = nv * w[i];
258		✗	sum = sum + nv;
259			}
260		✗	sum = sum*(1. / sum.norm());
261		✗	return sum.project_mat3();
262			}
263
264		✗	mat3 Frame::representative_frame(vector<vec3>& bunch_of_vectors) {
265		✗	vector<double> w(bunch_of_vectors.size(), 1);
266		✗	return representative_frame(bunch_of_vectors, w);
267		✗	}
268
269			/***********************************************************************************************************************************/
270
271		✗	AxisPermutation Rij(Mesh* m, Attribute<mat3>& B, index_t i, index_t j) {
272		✗	if (i > j) return Rij(m, B, j, i).inverse(); // façon obscure de rendre le bazar symmetrique ?
273		✗	AxisPermutation permut;
274		✗	permut.aligns_B_wrt_ref(B[i], B[j]);
275		✗	return permut;
276			}
277
278		✗	bool triangle_is_frame_singular(Mesh* m, Attribute<mat3>& B, index_t c, index_t cf) {
279		✗	mat3 ac_rot;
280		✗	ac_rot.load_identity();
281		✗	FOR(e, 3) ac_rot = Rij(m, B, m->cells.facet_vertex(c, cf, e), m->cells.facet_vertex(c, cf, next_mod(e, 3))).get_mat()*ac_rot;
282		✗	return !is_identity_auvp(ac_rot);
283			}
284
285		✗	bool triangle_is_frame_singular___give_stable_direction(Mesh* m, int& stable_dir_index, Attribute<mat3>& B, index_t c, index_t cf, index_t cfv) {
286		✗	mat3 ac_rot,id3D;
287		✗	stable_dir_index = -1;
288		✗	ac_rot.load_identity();
289		✗	FOR(e, 3) ac_rot = Rij(m, B, m->cells.facet_vertex(c, cf, (e+ cfv)%3), m->cells.facet_vertex(c, cf, (e + cfv+1) % 3)).get_mat()*ac_rot;
290		✗	if (is_identity_auvp(ac_rot)) return false;
291		✗	FOR(ax, 3){
292		✗	double sum = 0;
293		✗	FOR(d, 3) sum += std::abs(ac_rot(d, ax) - id3D(d, ax));
294		✗	if (sum<1e-10) stable_dir_index = int(ax);
295			}
296		✗	return true;
297			}
298		✗	bool triangle_is_frame_singular___give_stable_direction(Mesh* m, vec3& stable_dir_geom, Attribute<mat3>& B, index_t c, index_t cf, index_t cfv) {
299			int stable_dir_index;
300		✗	if (!triangle_is_frame_singular___give_stable_direction(m, stable_dir_index, B, c, cf,cfv)) return false;
301		✗	stable_dir_geom = vec3(0, 0, 0);
302		✗	if (stable_dir_index == -1) return true;
303		✗	FOR(d, 3) stable_dir_geom[d] = B[m->cells.facet_vertex(c, cf, cfv)](d, index_t(stable_dir_index));
304		✗	stable_dir_geom = normalize(stable_dir_geom);
305		✗	return true;
306			}
307
308			}
309