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/FF.h>
41			#include <exploragram/hexdom/frame.h>
42			#include <exploragram/hexdom/basic.h>
43			#include <exploragram/hexdom/extra_connectivity.h>
44			#include <geogram/NL/nl.h>
45			#include <geogram/numerics/optimizer.h>
46
47			#ifdef GEO_OPENMP
48			#include <omp.h>
49			#endif
50			#include <queue>
51
52			namespace GEO {
53
54		✗	FFopt::FFopt(Mesh* p_m) {
55		✗	m = p_m;
56		✗	compute_tet_edge_graph(m,v2e, true); // here need a bidirectionl edge graph to speed up the LBFGS part
57		✗	Attribute<vec3> lockB(m->vertices.attributes(), "lockB");
58		✗	num_l_v = m->vertices.nb();
59		✗	num_ln_v = m->vertices.nb();
60		✗	FOR(inv_v, m->vertices.nb()) {
61		✗	index_t v = m->vertices.nb()-1 - inv_v;
62		✗	if (lockB[v][0] <.5) num_l_v = v;
63		✗	if (lockB[v][2] <.5) num_ln_v = v;
64			}
65		✗	if (num_ln_v == 0) num_ln_v = m->vertices.nb();
66		✗	}
67
68		✗	FFopt::~FFopt() {
69		✗	m->edges.clear();
70		✗	}
71
72
73
74
75		✗	void FFopt::FF_init(bool generate_sh) {
76		✗	Attribute<mat3> B(m->vertices.attributes(), "B");
77		✗	Attribute<vec3> lockB(m->vertices.attributes(), "lockB");
78		✗	Attribute<SphericalHarmonicL4> sh;
79		✗	if (generate_sh) sh.bind(m->vertices.attributes(), "sh");
80
81		✗	double smooth_coeff = 1.;
82		✗	double normal_coeff = 100.;
83
84		✗	plop(num_l_v);
85		✗	plop(num_ln_v);
86		✗	plop("construct system");
87		✗	nlNewContext();
88		✗	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
89		✗	nlSolverParameteri(NL_NB_VARIABLES, NLint(2 * (num_ln_v - num_l_v) + 9 * m->vertices.nb()));
90		✗	nlBegin(NL_SYSTEM);
91
92			// lock frames
93		✗	FOR(v, num_l_v) {
94		✗	SphericalHarmonicL4 sh48;
95		✗	sh48[4] = std::sqrt(7. / 12.);
96		✗	sh48[8] = std::sqrt(5. / 12.);
97		✗	sh48.euler_rot(mat3_to_euler(normalize_columns(B[v])));
98		✗	FOR(i, 9) {
99		✗	nlSetVariable(v * 9 + i, sh48[i]);
100		✗	nlLockVariable(v * 9 + i);
101			}
102			}
103		✗	nlBegin(NL_MATRIX);
104			// smoothing equations
105		✗	FOR(e, m->edges.nb()) {
106		✗	FOR(i, 9) {
107		✗	nlRowScaling(smooth_coeff);
108		✗	nlBegin(NL_ROW);
109		✗	nlCoefficient(m->edges.vertex(e, 0) * 9 + i, -1.);
110		✗	nlCoefficient(m->edges.vertex(e, 1) * 9 + i, 1.);
111		✗	nlEnd(NL_ROW);
112			}
113			}
114			// boundary condition enforced by barrier equations
115		✗	for (index_t v = num_l_v; v < num_ln_v; v++) {
116		✗	SphericalHarmonicL4 sh0, sh4, sh8;
117		✗	sh4[4] = std::sqrt(7. / 12.);
118		✗	sh0[0] = std::sqrt(5. / 12.);
119		✗	sh8[8] = std::sqrt(5. / 12.);
120		✗	vec3 xyz = mat3_to_euler(normalize_columns(B[v]));
121		✗	sh4.euler_rot(xyz);
122		✗	sh0.euler_rot(xyz);
123		✗	sh8.euler_rot(xyz);
124		✗	FOR(i, 9) {
125		✗	nlRowScaling(normal_coeff);
126		✗	nlBegin(NL_ROW);
127		✗	nlCoefficient(v * 9 + i, 1.);
128		✗	nlCoefficient(m->vertices.nb() * 9 + v - num_l_v, sh0[i]);
129		✗	nlCoefficient(m->vertices.nb() * 9 + (num_ln_v - num_l_v) + v - num_l_v, sh8[i]);
130		✗	nlRightHandSide(sh4[i]);
131		✗	nlEnd(NL_ROW);
132			}
133			}
134
135
136
137		✗	nlEnd(NL_MATRIX);
138		✗	nlEnd(NL_SYSTEM);
139		✗	nlSolve();
140
141		✗	plop("project SH");
142			// convert spherical harmonic coefficients to a rotation
143			#ifdef GEO_OPENMP
144		✗	#pragma omp parallel
145			#endif
146			{
147			get_thread_range(m->vertices.nb(), start, end);
148			for (index_t v = start; v < end; v++) {
149			SphericalHarmonicL4 fv;
150			FOR(i, 9) fv[i] = nlGetVariable(v * 9 + i);
151			if (generate_sh) sh[v] = fv;
152			if (v >= num_l_v) {
153			vec3 oldz = col(B[v], 2);
154			if (v > start && v > num_l_v) {
155			vec3 prev = mat3_to_euler(normalize_columns(B[v - 1]));
156			B[v] = fv.project_mat3(1e-3, 1e-5, &prev);
157			} else
158			B[v] = fv.project_mat3(1e-3, 1e-5, nullptr);
159			if (v <= num_ln_v) {
160			AxisPermutation ap;
161			ap.make_col2_equal_to_z(B[v], normalize(oldz));
162			B[v] = Frame(B[v]).apply_permutation(ap);
163			FOR(d, 3) B[v](d, 2) = oldz[d];// restore size as well
164			}
165			}
166			}
167			}
168		✗	nlDeleteContext(nlGetCurrent());
169
170		✗	}
171
172
173			// place older with constant size
174		✗	void FFopt::compute_Bid_norm() {
175		✗	Attribute<mat3> B(m->vertices.attributes(), "B");
176		✗	double scale = col(B[0], 2).length();
177		✗	FOR(v,m->vertices.nb()) {
178		✗	FOR(a, 3) {
179		✗	vec3 co = scale * normalize(col(B[v], a));
180		✗	FOR(d, 3) B[v](d, a) = co[d];// restore size as well
181			}
182			}
183		✗	}
184			}
185
186
187			namespace {
188
189			using namespace GEO;
190
191			namespace FF_LBFGS {
192			FFopt* ffopt_ptr;
193			index_t Num_ln_v;
194			index_t Num_l_v;
195			double lastf;
196			double NRJ_threshold = 1e-5;
197			int nb_iters;
198			GEO::Optimizer *solver;
199			}
200
201		✗	void new_iteration_cb(index_t N, const double* x, double f, const double* g, double gnorm) {
202		✗	FF_LBFGS::nb_iters++;
203		✗	double stop_crit = std::abs(FF_LBFGS::lastf - f) / std::abs(f);
204		✗	FF_LBFGS::lastf = f;
205		✗	std::cerr << ".";
206		✗	if (stop_crit < FF_LBFGS::NRJ_threshold) {
207		✗	GEO::Logger::out("HexDom") << " LBFGS iter " << N << " f " << f << " gnorm " << gnorm << " trash " << x[0] * g[0] << std::endl;
208		✗	GEO::Logger::out("HexDom") << "stop_crit < NRJ_threshold " << std::endl; throw 1;
209			}
210		✗	}
211
212		✗	void compute_gradient_cb2(unsigned int N, double* x, double& f, double* g) {
213		✗	mat3 mEx = mat3_from_coeffs( 0, 0, 0, 0, 0, -1, 0, 1, 0 );
214		✗	mat3 mEy = mat3_from_coeffs(0, 0, 1, 0, 0, 0, -1, 0, 0 );
215		✗	mat3 mEz = mat3_from_coeffs(0, -1, 0, 1, 0, 0, 0, 0, 0 );
216
217		✗	Attribute<mat3> B(FF_LBFGS::ffopt_ptr->m->vertices.attributes(), "B");
218		✗	index_t nverts = FF_LBFGS::ffopt_ptr->m->vertices.nb();
219		✗	geo_assert(N == 3 * (nverts - FF_LBFGS::Num_ln_v) + FF_LBFGS::Num_ln_v);
220
221
222		✗	Attribute<bool> border_vertex(FF_LBFGS::ffopt_ptr->m->vertices.attributes(), "border_vertex");
223		✗	FOR(v, FF_LBFGS::ffopt_ptr->m->vertices.nb()) border_vertex[v] = false;
224		✗	FOR(c, FF_LBFGS::ffopt_ptr->m->cells.nb()) FOR(cf, 4) if (FF_LBFGS::ffopt_ptr->m->cells.adjacent(c, cf) == NOT_AN_ID)
225		✗	FOR(cfv, 3) border_vertex[FF_LBFGS::ffopt_ptr->m->cells.facet_vertex(c, cf, cfv)] = true;
226
227			#ifdef GEO_OPENMP
228		✗	int max_threads = omp_get_max_threads();
229			#else
230			int max_threads = 1;
231			#endif
232			// f_chunks is initialized to be zero
233		✗	double *f_chunks = new double[size_t(max_threads)]();
234
235			#ifdef GEO_OPENMP
236		✗	#pragma omp parallel
237			#endif
238			{
239			#ifdef GEO_OPENMP
240			int thread_id = omp_get_thread_num();
241			#else
242			int thread_id = 0;
243			#endif
244			get_thread_range(nverts, istart, iend);
245
246
247			mat3 mJR[3], mR, mSinv, mPst, mJPst[3];
248			for (index_t v1 = istart; v1 < iend; v1++) {
249			if (v1 >= FF_LBFGS::Num_ln_v) {
250			index_t idx = FF_LBFGS::Num_ln_v + (v1 - FF_LBFGS::Num_ln_v) * 3;
251			FOR(i,3) g[idx + i] = 0.;
252			mR = euler_to_mat3(vec3(x[idx], x[idx+1], x[idx+2]));
253			mat3 mRx = rotx(x[idx]);
254			mat3 mRy = roty(x[idx+1]);
255			mat3 mRz = rotz(x[idx+2]);
256
257			mJR[0] = mR*mEx;
258			mJR[2] = mEz*mR;
259			mJR[1] = mRz*mRy*mEy*mRx;
260			}
261			else {
262			g[v1] = 0.;
263			mR = normalize_columns(B[v1]) * rotz(x[v1]);
264			mJR[0] = mR* mEz; // init JR[0] = R * Ez; JR[1] and JR[2] are not initialized
265			}
266
267
268			FOR(iv2, int(FF_LBFGS::ffopt_ptr->nb_neigs(v1))) {
269			index_t v2 = FF_LBFGS::ffopt_ptr->neig(v1, iv2);
270			if (v2 >= FF_LBFGS::Num_ln_v) { // init S = Rz Ry Rx
271			index_t idx = FF_LBFGS::Num_ln_v + (v2 - FF_LBFGS::Num_ln_v) * 3;
272			mSinv = euler_to_mat3(vec3(x[idx], x[idx + 1], x[idx + 2]));
273			} else // init S = constraint * Rz
274			mSinv = normalize_columns(B[v2]) * rotz(x[v2]);
275
276
277			mSinv = mSinv.transpose();
278			mPst = mSinv* mR; // Pst = S^{-1} * R
279
280			double scale = 1.;
281			if (HexdomParam::FF.rigid_border) {
282			if (border_vertex[v1])scale += 100.;
283			if (border_vertex[v2])scale += 100.;
284			}
285			if (v1 > v2) FOR(i, 3)
286			f_chunks[thread_id] += scale *(10. / 3.*(pow(mPst(0,i) * mPst(1,i), 2) + pow(mPst(0,i) * mPst(2,i), 2) + pow(mPst(1,i )* mPst(2,i), 2)));
287
288			if (v1 >= FF_LBFGS::Num_ln_v) {
289			index_t idx = FF_LBFGS::Num_ln_v + (v1 - FF_LBFGS::Num_ln_v) * 3;
290			FOR(d,3) {
291			mJPst[d] = mSinv* mJR[d]; // JPst[d] = S^{-1} * JR[d]
292			FOR(i,3)FOR(j,3)
293			g[idx + d] += scale *(20. / 3.*mPst(i,j) * (pow(mPst(i, (j + 1)%3), 2) + pow(mPst(i , (j + 2) % 3), 2))*mJPst[d](i , j));
294			}
295			} else if (v1 >= FF_LBFGS::Num_l_v) {
296			mJPst[0] = mSinv* mJR[0]; // JPst[0] = S^{-1} * JR[0] ; JPst[1] and JPst[2] are not initialized
297			FOR(i, 3)FOR(j, 3)
298			g[v1] += scale *(20. / 3.*mPst(i ,j) * (pow(mPst(i ,(j + 1) % 3), 2) + pow(mPst(i , (j + 2) % 3), 2))*mJPst[0](i , j));
299			}
300			} // v2
301			} // v1
302			} // omp parallel
303		✗	f = 0.;
304		✗	for (int i = max_threads; i--; f += f_chunks[i]);
305		✗	delete[] f_chunks;
306		✗	}
307			}
308
309			namespace GEO {
310
311		✗	void FFopt::FF_smooth() {
312
313		✗	Attribute<mat3> B(m->vertices.attributes(), "B");
314
315		✗	Attribute<SphericalHarmonicL4> sh(m->vertices.attributes(), "sh");
316
317			// init global variables that must be visible in callbacks
318		✗	FF_LBFGS::ffopt_ptr = this;
319		✗	FF_LBFGS::Num_ln_v = num_ln_v;
320		✗	FF_LBFGS::Num_l_v = num_l_v;
321		✗	FF_LBFGS::lastf = 1e20;
322		✗	FF_LBFGS::NRJ_threshold = 1e-5;
323		✗	FF_LBFGS::nb_iters = 0;
324
325			// create LBFGS solver and unknown vector
326		✗	index_t nverts = m->vertices.nb();
327
328			// unknown vetor is packed as follows:
329			// FF_LBFGS::Num_ln_v coordinates: 1 rotation angle around the constrained axis
330			// nverts - FF_LBFGS::Num_ln_v coordinates: 3 euler angles
331			//
332			// WARNING: note that locked frames (v<num_l_v) are associated to a useless variables
333		✗	double *x = new double[nverts * 3 - FF_LBFGS::Num_ln_v * 2];
334		✗	GEO::Optimizer *solver = GEO::Optimizer::create();
335		✗	FF_LBFGS::solver = solver;
336		✗	solver->set_N((nverts - FF_LBFGS::Num_ln_v) * 3 + FF_LBFGS::Num_ln_v);
337		✗	solver->set_M(3);
338		✗	solver->set_epsf(1e-5);
339		✗	solver->set_epsx(1e-5);
340		✗	solver->set_epsg(1e-5);
341		✗	solver->set_funcgrad_callback(compute_gradient_cb2);
342		✗	solver->set_newiteration_callback(new_iteration_cb);
343
344			// init variables
345		✗	for (index_t i = FF_LBFGS::Num_ln_v; i < nverts; i++) {
346		✗	index_t idx = (i - FF_LBFGS::Num_ln_v) * 3 + FF_LBFGS::Num_ln_v;
347		✗	vec3 xyz = mat3_to_euler(normalize_columns(B[i]));
348		✗	FOR(d, 3) x[idx+d] = xyz[d];
349			}
350		✗	FOR(i, FF_LBFGS::Num_ln_v) x[i] = 0.;
351
352
353
354			// solve until we run out of time
355		✗	solver->set_max_iter(1000000);
356		✗	try { solver->optimize(x); }
357		✗	catch (...) {
358
359		✗	}
360
361			// apply a euler rotation to rot...
362		✗	for (index_t i = nverts; i--;) {
363		✗	if (i >= FF_LBFGS::Num_ln_v) {
364		✗	index_t idx = (i - FF_LBFGS::Num_ln_v) * 3 + FF_LBFGS::Num_ln_v;
365		✗	B[i] = euler_to_mat3(vec3 (x[idx], x[idx + 1], x[idx + 2]));
366			}
367		✗	else B[i] = B[i] * rotz(x[i]);
368
369		✗	FOR(d, 9)sh[i][d] = 0;
370		✗	sh[i][4] = std::sqrt(7. / 12.);
371		✗	sh[i][8] = std::sqrt(5. / 12.);
372		✗	sh[i].euler_rot(mat3_to_euler(normalize_columns(B[i])));
373			}
374		✗	delete[] x;
375		✗	}
376
377
378
379			// ___ _ ___ _ _ _
380			// | _ )_ _ _ _ __| |_ |_ ) ___ _ __| |_(_)_ __ (_)______
381			// | _ \ '_| || (_-< ' \ / / / _ \ '_ \ _| | ' \| |_ / -_)
382			// |___/_| \_,_/__/_||_| /___| \___/ .__/\__|_|_|_|_|_/__\___|
383			// |_|
384
385
386
387		✗	void FFopt::brush_frame() {
388		✗	plop("brushing");
389		✗	Attribute<vec3> lockU(m->vertices.attributes(), "lockU");// how many dimensions are locked
390		✗	Attribute<mat3> B(m->vertices.attributes(), "B");
391		✗	vector<bool> seen(m->vertices.nb(), false);
392		✗	FOR(seed, m->vertices.nb()) { // multiple components?
393		✗	if (seen[seed]) continue;
394		✗	seen[seed] = true;
395		✗	std::deque<index_t> Q;
396		✗	Q.push_back(seed);
397		✗	while (Q.size()) { // start a breadth-first brushing
398		✗	index_t cur = Q.front();
399		✗	Q.pop_front();
400		✗	FOR(lv, nb_neigs(cur)) {
401		✗	index_t v = neig(cur, lv);
402		✗	seen[v] = true;
403		✗	AxisPermutation M=Rij(m,B,v,cur);
404		✗	if (M.mid != 0) {
405		✗	lockU[v] = M.inverse().get_mat()* lockU[v];
406		✗	B[v] = B[v]* M.get_mat();
407		✗	FOR(d, 3) if (std::abs(lockU[v][d]) < .1) lockU[v][d] = 0;
408			}
409		✗	Q.push_back(v);
410			}
411			}
412		✗	}
413		✗	}
414
415			}
416