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/quad_cover.h>
41			#include <exploragram/hexdom/mixed_constrained_solver.h>
42			#include <geogram/mesh/mesh.h>
43			#include <geogram/mesh/mesh_geometry.h>
44			#include <geogram/mesh/mesh_frame_field.h>
45
46			namespace GEO {
47
48			namespace GlobalParam2d {
49
50			namespace Internal {
51
52			/**
53			* \brief the four 2x2 rotation matrices associated with the
54			* values of R that transform coordinates between two triangles.
55			* \details First index is Rij (in 0..3, number of 90 degrees
56			* rotations), then row and column index of the 2x2 rotation
57			* matrix, where Rij is the number of times coordinates axes are
58			* rotated by 90 degrees.
59			* \note The rotation is inversed as compared with the
60			* computation in Rij() since when the axes rotate clockwise,
61			* the coordinates rotate anticlockwise (and conversely).
62			*/
63			static double Rot[4][2][2] = {
64			{{1, 0},
65			{0, 1}},
66
67			{{ 0, 1},
68			{-1, 0}},
69
70			{{-1, 0},
71			{ 0,-1}},
72
73			{{0 ,-1},
74			{1 , 0}}
75			};
76
77		✗	void quad_cover_solve(
78			Mesh* mesh,
79			Attribute<vec3>& B, Attribute<index_t>& R,
80			Attribute<index_t>& on_border,
81			Attribute<index_t>& constraints,
82			Attribute<vec2>& U,
83			Attribute<double>& T,
84			Attribute<bool>& v_is_singular,
85			double scaling,
86			bool constrain_hard_edges,
87			bool integer_constraints
88			) {
89		✗	scaling *= surface_average_edge_length(*mesh);
90		✗	index_t nb_U = mesh->facets.nb()*3*2;
91			index_t nb_T = mesh->facets.nb()*3*2;
92		✗	MatrixMixedConstrainedSolver solver(nb_U+nb_T);
93
94			// All Tijs are even integers (entiers pairs).
95		✗	FOR(t, nb_T) {
96		✗	solver.set_multiplicity(nb_U+t,2);
97			}
98
99			// Constrained u,v coordinates
100		✗	if(constrain_hard_edges) {
101		✗	FOR(c, mesh->facet_corners.nb()) {
102		✗	if((constraints[c] & CNSTR_U) != 0) {
103		✗	solver.set_multiplicity(2*c,1);
104			}
105		✗	if((constraints[c] & CNSTR_V) != 0) {
106		✗	solver.set_multiplicity(2*c+1,1);
107			}
108			}
109			}
110
111			// Compute vertex-to-corner map (exclude vertices on the border
112			// and singular vertices)
113
114			vector<index_t> v2c(mesh->vertices.nb(), NO_CORNER);
115		✗	FOR(c, mesh->facet_corners.nb()) {
116			index_t v = mesh->facet_corners.vertex(c);
117		✗	v2c[v] = c;
118			}
119			{
120		✗	FOR(c, mesh->facet_corners.nb()) {
121			index_t v = mesh->facet_corners.vertex(c);
122			if(
123			(
124			mesh->facet_corners.adjacent_facet(c) ==
125			NO_FACET
126		✗	) || v_is_singular[v]
127			) {
128		✗	v2c[v] = NO_CORNER;
129			}
130			}
131			}
132
133
134		✗	FOR(pass, 4) {
135		✗	FOR(c, mesh->facet_corners.nb()) {
136		✗	if(mesh->facet_corners.adjacent_facet(c) == NO_FACET) {
137		✗	continue;
138			}
139
140			index_t f2 = mesh->facet_corners.adjacent_facet(c);
141		✗	index_t e2 = mesh->facets.find_adjacent(f2,c/3);
142		✗	index_t c2 = mesh->facets.corners_begin(f2)+e2;
143			index_t c3 =
144			mesh->facets.next_corner_around_facet(f2,c2);
145
146		✗	geo_assert(
147			mesh->facet_corners.vertex(c) ==
148			mesh->facet_corners.vertex(c3)
149			);
150		✗	index_t Rij = R[c];
151
152			// Chart transform for each pair of adjacent triangles
153
154			// On the border of the ball, Tij + Rij*Tji = 0
155			// (the Tij 1-form is ... a 1-form)
156		✗	if(on_border[c]) {
157			solver.begin_constraint();
158		✗	solver.add_constraint_coeff(nb_U+2*c, 1.0);
159		✗	solver.add_constraint_coeff(nb_U+2*c2, Rot[Rij][0][0]);
160		✗	solver.add_constraint_coeff(nb_U+2*c2+1, Rot[Rij][0][1]);
161			solver.end_constraint();
162			solver.begin_constraint();
163		✗	solver.add_constraint_coeff(nb_U+2*c+1, 1.0);
164		✗	solver.add_constraint_coeff(nb_U+2*c2, Rot[Rij][1][0]);
165		✗	solver.add_constraint_coeff(nb_U+2*c2+1, Rot[Rij][1][1]);
166			solver.end_constraint();
167			} else {
168			// Inside the ball, Tij = 0
169			solver.begin_constraint();
170		✗	solver.add_constraint_coeff(nb_U+2*c, 1.0);
171			solver.end_constraint();
172			solver.begin_constraint();
173		✗	solver.add_constraint_coeff(nb_U+2*c+1, 1.0);
174			solver.end_constraint();
175			}
176
177			// Setup relation between Ui - Rij*Uj - Tij = 0
178			// (dU = T)
179			solver.begin_constraint();
180		✗	solver.add_constraint_coeff(2*c , 1.0);
181		✗	solver.add_constraint_coeff(2*c3 , -Rot[Rij][0][0]);
182		✗	solver.add_constraint_coeff(2*c3+1, -Rot[Rij][0][1]);
183		✗	solver.add_constraint_coeff(nb_U+2*c,-1.0);
184			solver.end_constraint();
185
186			solver.begin_constraint();
187		✗	solver.add_constraint_coeff(2*c+1 , 1.0);
188		✗	solver.add_constraint_coeff(2*c3 , -Rot[Rij][1][0]);
189		✗	solver.add_constraint_coeff(2*c3+1, -Rot[Rij][1][1]);
190		✗	solver.add_constraint_coeff(nb_U+2*c+1,-1.0);
191			solver.end_constraint();
192			}
193
194			//Wheel compatibility constraints (the Tij 1-form is closed)
195		✗	FOR(v, mesh->vertices.nb()) {
196			// If the corner is on the border or incident
197			// to a singular vertex then it is skipped.
198		✗	if(v2c[v] == NO_CORNER) {
199		✗	continue;
200			}
201			// Enforce the constraint on the wheel
202			// neighborhood for each component of the Tijs.
203		✗	FOR(coord, 2) {
204		✗	index_t c = v2c[v];
205			index_t r = 0;
206			solver.begin_constraint();
207			do {
208		✗	solver.add_constraint_coeff(
209		✗	nb_U+2*c, Rot[r][coord][0]
210			);
211		✗	solver.add_constraint_coeff(
212			nb_U+2*c+1, Rot[r][coord][1]
213			);
214			// Accumulate the rotation.
215		✗	r = (r + R[c]) % 4;
216			// Find the next corner around the vertex.
217			index_t f =
218			mesh->facet_corners.adjacent_facet(c);
219
220		✗	geo_assert(f != NO_FACET);
221			index_t next_c = NO_CORNER;
222		✗	for(
223			next_c = mesh->facets.corners_begin(f);
224		✗	next_c<mesh->facets.corners_end(f);
225			++next_c
226			) {
227		✗	if(mesh->facet_corners.vertex(next_c)==v) {
228			break;
229			}
230			}
231		✗	geo_assert(
232			mesh->facet_corners.vertex(next_c) == v
233			);
234			c = next_c;
235		✗	} while(c != v2c[v]);
236			// On non-singular vertices, by definition,
237			// compose of all rotations = identity.
238		✗	geo_assert(r == 0);
239			solver.end_constraint();
240			}
241			}
242
243			// Constrained edges - equality between coordinates
244			// Note: sometimes, setting this constraint causes
245			// an assertion failure in Nico's mixed integer solver:
246			// Assertion failed: pass != 3 || cM0M1M2.empty().
247			// Note2: seems to be OK now that the wheel compat. cnstr.
248			// is there (to be checked).
249		✗	if(constrain_hard_edges) {
250		✗	FOR(c, mesh->facet_corners.nb()) {
251		✗	index_t f=c/3;
252			index_t c2 =
253			mesh->facets.next_corner_around_facet(f,c);
254		✗	index_t cnstr = get_edge_constraints(mesh,c,B);
255		✗	if(cnstr != 0) {
256		✗	if(cnstr == CNSTR_U) {
257			solver.begin_constraint();
258		✗	solver.add_constraint_coeff(2*c, 1.0);
259		✗	solver.add_constraint_coeff(2*c2, -1.0);
260			solver.end_constraint();
261		✗	} else if(cnstr == CNSTR_V) {
262			solver.begin_constraint();
263		✗	solver.add_constraint_coeff(2*c+1, 1.0);
264		✗	solver.add_constraint_coeff(2*c2+1, -1.0);
265			solver.end_constraint();
266			} else {
267		✗	geo_assert_not_reached;
268			}
269			}
270			}
271			}
272
273			solver.end_pass(pass);
274			}
275
276
277		✗	while (!solver.converged()) {
278		✗	plop("MIQ iter");
279		✗	solver.start_new_iter();
280
281		✗	FOR(f, mesh->facets.nb()) {
282			// setup objective function :
283			// For each edge (pi,pj):
284			// ( B * (pj-pi) - (uj-ui))^2 +
285			// ( rot90(B) * (pj-pi) - (vj-vi))^2
286		✗	vec3 N = normalize(Geom::mesh_facet_normal(*mesh,f));
287		✗	vec3 Bf = normalize(B[f]);
288			vec3 BTf = cross(N,Bf);
289
290			for(index_t c1 = mesh->facets.corners_begin(f);
291		✗	c1 < mesh->facets.corners_end(f); ++c1) {
292			index_t c2 =
293			mesh->facets.next_corner_around_facet(f,c1);
294			index_t v1 = mesh->facet_corners.vertex(c1);
295			index_t v2 = mesh->facet_corners.vertex(c2);
296			vec3 E =
297			vec3(mesh->vertices.point_ptr(v2)) -
298			vec3(mesh->vertices.point_ptr(v1));
299			solver.begin_energy();
300		✗	solver.add_energy_coeff(2*c2, scaling);
301		✗	solver.add_energy_coeff(2*c1,-scaling);
302			solver.add_energy_rhs(dot(Bf,E));
303			solver.end_energy();
304			solver.begin_energy();
305		✗	solver.add_energy_coeff(2*c2+1, scaling);
306		✗	solver.add_energy_coeff(2*c1+1,-scaling);
307			solver.add_energy_rhs(dot(BTf,E));
308			solver.end_energy();
309			}
310			}
311		✗	solver.end_iter();
312		✗	if(!integer_constraints) {
313			break;
314			}
315			}
316
317
318			// Get the result
319		✗	FOR(u, nb_U) {
320		✗	double coord = solver.value(u);
321		✗	snap_tex_coord(coord); // Required by mesh extraction
322		✗	U[u/2][u%2] = coord;
323			}
324		✗	FOR(t, nb_T) {
325		✗	T[t] = solver.value(nb_U+t);
326			}
327		✗	}
328
329			} // namespace Internal
330
331		✗	void quad_cover(
332			Mesh* mesh,
333			Attribute<vec3>& B, Attribute<vec2>& U,
334			double scaling, bool constrain_hard_edges, bool do_brush,
335			bool integer_constraints
336			) {
337			{
338		✗	Attribute<index_t> R_ff(mesh->facet_corners.attributes(),"R");
339			Attribute<index_t> c_on_border(
340		✗	mesh->facet_corners.attributes(), "on_border"
341			);
342
343		✗	if(do_brush) {
344		✗	Internal::brush(mesh,B);
345			}
346		✗	Internal::compute_R_ff(mesh,B,R_ff);
347			Attribute<bool> v_is_singular(
348		✗	mesh->vertices.attributes(), "is_singular"
349			);
350		✗	Internal::mark_singular_vertices(mesh, R_ff, v_is_singular);
351
352		✗	if(do_brush) {
353		✗	Internal::do_the_ball(mesh, R_ff, c_on_border);
354			} else {
355		✗	Internal::do_the_ball_no_brush_no_zip(mesh, c_on_border);
356			}
357
358			Attribute<double> T;
359		✗	T.bind_if_is_defined(mesh->facet_corners.attributes(),"T");
360			if(!T.is_bound()) {
361		✗	T.create_vector_attribute(
362		✗	mesh->facet_corners.attributes(), "T", 2
363			);
364			}
365			Attribute<index_t> constraint(
366		✗	mesh->facet_corners.attributes(), "cnstr"
367			);
368		✗	Internal::get_constraints(mesh, B, R_ff, constraint);
369		✗	Internal::quad_cover_solve(
370			mesh, B, R_ff, c_on_border,
371			constraint, U, T, v_is_singular,
372			scaling, constrain_hard_edges,
373			integer_constraints
374			);
375			}
376
377			// Destroy the temporary attributes
378			// I keep them for now, for debugging...
379			//
380			// mesh->facet_corners.attributes().delete_attribute_store("cnstr");
381			// mesh->facet_corners.attributes().delete_attribute_store("on_border");
382			// mesh->facet_corners.attributes().delete_attribute_store("T");
383			// mesh->facet_corners.attributes().delete_attribute_store("R");
384			// mesh->facet_corners.attributes().delete_attribute_store("UU");
385			// mesh->vertices.attributes().delete_attribute_store("is_singular");
386
387		✗	}
388
389			} // namespace GlobalParam2d
390
391			} // namespace GEO
392