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 <geogram/parameterization/mesh_global_param.h>
41			#include <geogram/mesh/mesh.h>
42			#include <geogram/mesh/mesh_geometry.h>
43			#include <geogram/mesh/mesh_frame_field.h>
44
45			#include <deque>
46			#include <stack>
47
48
49			namespace {
50			using namespace GEO;
51
52			/**
53			* \brief Computes the angle between two vectors associated to two adjacent
54			* facets.
55			* \param[in] mesh a pointer to a surface mesh
56			* \param[in] f1 the first facet
57			* \param[in] B1 the first 3d vector, in the plane of f1
58			* \param[in] f2 the second facet
59			* \param[in] B2 the second 3d vector, in the plane of f2
60			* \return the angle in degrees between -180 and 180 required to
61			* transform vector \p B2 to vector \p B1 expressed in the local frame
62			* of facet f1 with the common edge of the two facets f1 and f2 as
63			* the X axis
64			*/
65		✗	double angle(
66			Mesh* mesh, index_t f1, const vec3& B1, index_t f2, const vec3& B2
67			) {
68		✗	index_t lf2 = mesh->facets.find_adjacent(f1,f2);
69		✗	geo_assert(lf2 != NO_FACET);
70		✗	index_t c1 = mesh->facets.corners_begin(f1) + lf2;
71			index_t c2 = mesh->facets.next_corner_around_facet(f1,c1);
72			index_t v1 = mesh->facet_corners.vertex(c1);
73			index_t v2 = mesh->facet_corners.vertex(c2);
74		✗	if(v2 < v1) {
75			std::swap(v1,v2);
76			}
77			const vec3& p1 = mesh->vertices.point(v1);
78			const vec3& p2 = mesh->vertices.point(v2);
79		✗	vec3 E = normalize(p2-p1);
80		✗	vec3 N1 = normalize(Geom::mesh_facet_normal(*mesh,f1));
81		✗	vec3 N2 = normalize(Geom::mesh_facet_normal(*mesh,f2));
82			vec3 Y1 = cross(N1,E);
83			vec3 Y2 = cross(N2,E);
84			double x1 = dot(B1,E);
85			double y1 = dot(B1,Y1);
86			double x2 = dot(B2,E);
87			double y2 = dot(B2,Y2);
88		✗	double a1 = atan2(y1,x1) * 180.0 / M_PI;
89		✗	double a2 = atan2(y2,x2) * 180.0 / M_PI;
90		✗	double result = a1-a2;
91
92		✗	while(result < -180.0) {
93		✗	result += 360.0;
94			}
95		✗	while(result > 180.0) {
96		✗	result -= 360.0;
97			}
98
99		✗	return result;
100			}
101
102
103			/**
104			* \brief Gets the number of 90 degrees rotations required to minimize the
105			* angle between the vectors attached to two adjacent facets.
106			* \param[in] mesh a pointer to a surface mesh
107			* \param[in] f1 the first facet
108			* \param[in] B1 the first 3D vector in the plane of f1
109			* \param[in] f2 the second facet
110			* \param[in] B2 the second 3D vector in the plane of f2
111			* \return the number of times \p B2 should be rotated around the
112			* normal vector of \p f2 to minimize its angle with \p B1 in angus
113			* (in 0,1,2,3).
114			*/
115		✗	index_t Rij(
116			Mesh* mesh, index_t f1, const vec3& B1, index_t f2, const vec3& B2
117			) {
118		✗	if(f1 > f2) {
119		✗	index_t result = Rij(mesh, f2, B2, f1, B1);
120		✗	return GlobalParam2d::Internal::inverse_R(result);
121			}
122		✗	vec3 N2 = normalize(Geom::mesh_facet_normal(*mesh,f2));
123		✗	vec3 cur_B2 = B2;
124		✗	double best_angle = ::fabs(angle(mesh, f1, B1, f2, cur_B2));
125			index_t best_i = 0;
126		✗	for(index_t i=1; i<4; ++i) {
127		✗	cur_B2 = cross(N2, cur_B2);
128		✗	double cur_angle = ::fabs(angle(mesh, f1, B1, f2, cur_B2));
129
130		✗	if(cur_angle < best_angle) {
131			best_angle = cur_angle;
132			best_i = i;
133			}
134			}
135			return best_i;
136			}
137
138			/**
139			* \brief Sets an attribute on both corners adjacent to the same edge.
140			* \param[in] mesh a pointer to a surface mesh.
141			* \param[out] attr a facet corner attribute
142			* \param[in] f1 , f2 the two facets that share the edge
143			* \param[in] val the new value of the attribute
144			*/
145		✗	void set_edge_attr(
146			Mesh* mesh, Attribute<index_t>& attr,
147			index_t f1, index_t f2, index_t val
148			) {
149		✗	index_t e1 = mesh->facets.find_adjacent(f1,f2);
150		✗	index_t e2 = mesh->facets.find_adjacent(f2,f1);
151		✗	attr[mesh->facets.corners_begin(f1)+e1] = val;
152		✗	attr[mesh->facets.corners_begin(f2)+e2] = val;
153		✗	}
154
155			}
156
157			namespace GEO {
158			namespace GlobalParam2d {
159			namespace Internal {
160
161		✗	void compute_R_ff(
162			Mesh* mesh, Attribute<vec3>& B, Attribute<index_t>& R_ff
163			) {
164		✗	for(index_t f1: mesh->facets) {
165		✗	FOR(e1, mesh->facets.nb_vertices(f1)) {
166			index_t f2 = mesh->facets.adjacent(f1,e1);
167			index_t c = mesh->facets.corners_begin(f1) + e1;
168		✗	if(f2 != NO_FACET) {
169		✗	index_t rij = Rij(mesh, f1, B[f1], f2, B[f2]);
170		✗	index_t rji = Rij(mesh, f2, B[f2], f1, B[f1]);
171		✗	geo_assert(rij == inverse_R(rji));
172		✗	R_ff[c] = rij;
173			}
174			}
175			}
176		✗	}
177
178
179		✗	void compute_R_fv(
180			Mesh* mesh,
181			Attribute<index_t>& R_ff, Attribute<index_t>& R_fv
182			) {
183			// - Each vertex has a reference corner (v2c[v])
184			// - Each corner c knows the number of rotations Rc[c] required
185			// to make the B of its triangle match the B of the triangle
186			// of the reference corner attached to its vertex (clear
187			// enough ?)
188
189			// Step 1: Compute v2c
190			// Note: if there exists a corner that has his previous corner
191			// around the facet that is on the surface border for a
192			// given v, then use this one.
193			// Later, when we turn around the vertices,
194			// it will be easier to start from such an halfedge for all
195			// vertices that are on the border.
196
197			vector<index_t> v2c(mesh->vertices.nb(), NO_CORNER);
198			{
199		✗	for(index_t c: mesh->facet_corners) {
200		✗	index_t f = c/3;
201			index_t c_prev =
202		✗	mesh->facets.prev_corner_around_facet(f,c);
203		✗	if(mesh->facet_corners.adjacent_facet(c_prev) ==
204			NO_FACET
205			) {
206			index_t v = mesh->facet_corners.vertex(c);
207		✗	v2c[v] = c;
208			}
209			}
210		✗	for(index_t c: mesh->facet_corners) {
211			index_t v = mesh->facet_corners.vertex(c);
212		✗	if(v2c[v] == NO_CORNER) {
213		✗	v2c[v] = c;
214			}
215			}
216			}
217
218			// Step 2: Compute Rv by turning around the facets that
219			// share a vertex.
220			// Yes, it is painful, I hate doing that, but R only works
221			// for pairs of *adjacent* facets (and I cannot think about
222			// a way of making work it *reliably* for any pair of facets).
223			{
224		✗	for(index_t v: mesh->vertices) {
225			index_t prev_c = NO_CORNER;
226		✗	index_t c = v2c[v];
227
228			// Isolated vertex, ignore
229		✗	if(c == NO_CORNER) {
230		✗	continue;
231			}
232
233			do {
234			index_t
235			next_f = mesh->facet_corners.adjacent_facet(c);
236			index_t next_c = NO_CORNER;
237		✗	if(next_f != NO_FACET) {
238		✗	for(index_t c2: mesh->facets.corners(next_f)) {
239		✗	if(mesh->facet_corners.vertex(c2) == v) {
240			next_c = c2;
241			break;
242			}
243			}
244		✗	geo_assert(next_c != NO_CORNER);
245			}
246		✗	if(prev_c != NO_CORNER) {
247		✗	R_fv[c] = (R_fv[prev_c] + R_ff[prev_c]) % 4;
248			}
249			prev_c = c;
250			c = next_c;
251		✗	} while(c != NO_CORNER && c != v2c[v]);
252			}
253			}
254		✗	}
255
256		✗	void mark_singular_vertices(
257			Mesh* mesh,
258			Attribute<index_t>& R_ff, Attribute<bool>& v_is_singular
259			) {
260		✗	vector<index_t> Rsum(mesh->vertices.nb(),0);
261		✗	for(index_t f: mesh->facets) {
262		✗	for(index_t c: mesh->facets.corners(f)) {
263		✗	if(
264			mesh->facet_corners.adjacent_facet(c) !=
265			NO_INDEX
266			) {
267			index_t v = mesh->facet_corners.vertex(c);
268		✗	Rsum[v] += R_ff[c];
269			}
270			}
271			}
272		✗	for(index_t v: mesh->vertices) {
273		✗	v_is_singular[v] = ((Rsum[v] % 4) != 0);
274			}
275			// Vertices on border can have non-zero Rsum without being
276			// singular.
277		✗	for(index_t c: mesh->facet_corners) {
278		✗	if(mesh->facet_corners.adjacent_facet(c) == NO_FACET) {
279			v_is_singular[mesh->facet_corners.vertex(c)] = false;
280			}
281			}
282		✗	}
283
284		✗	void brush(Mesh* mesh, Attribute<vec3>& B) {
285		✗	std::vector<bool> visited(mesh->facets.nb(),false);
286			std::deque<index_t> S;
287		✗	S.push_back(0);
288			visited[0] = true;
289		✗	while(!S.empty()) {
290		✗	index_t f1 = S.front();
291		✗	S.pop_front();
292		✗	FOR(e1, mesh->facets.nb_vertices(f1)) {
293		✗	index_t f2 = mesh->facets.adjacent(f1,e1);
294		✗	if(f2 != NO_FACET && !visited[f2]) {
295		✗	vec3 N2 = normalize(
296		✗	Geom::mesh_facet_normal(*mesh,f2)
297			);
298		✗	index_t Rc1 = Rij(mesh,f1,B[f1],f2,B[f2]);
299		✗	FOR(i, Rc1) {
300		✗	B[f2] = cross(N2,B[f2]);
301			}
302			visited[f2] = true;
303			S.push_back(f2);
304			}
305			}
306			}
307		✗	}
308
309
310		✗	void do_the_ball(
311			Mesh* mesh,
312			Attribute<index_t>& R_ff, Attribute<index_t>& c_on_border
313			) {
314		✗	for(index_t c: mesh->facet_corners) {
315		✗	c_on_border[c] = 1;
316			}
317
318			// Covering tree
319
320		✗	std::vector<bool> visited(mesh->facets.nb(),false);
321			std::deque<index_t> S;
322		✗	S.push_back(0);
323			visited[0] = true;
324		✗	while(!S.empty()) {
325		✗	index_t f1 = S.front();
326		✗	S.pop_front();
327		✗	for(index_t c1: mesh->facets.corners(f1)) {
328		✗	index_t f2 = mesh->facet_corners.adjacent_facet(c1);
329		✗	if(f2 != NO_FACET && !visited[f2] && R_ff[c1] == 0) {
330		✗	set_edge_attr(mesh, c_on_border, f1, f2, 0);
331			visited[f2] = true;
332			S.push_back(f2);
333			}
334			}
335			}
336
337			// Zipping
338
339		✗	vector<index_t> v_nb_borders(mesh->vertices.nb(), 0);
340		✗	for(index_t c: mesh->facet_corners) {
341		✗	if(c_on_border[c]) {
342		✗	++v_nb_borders[mesh->facet_corners.vertex(c)];
343			}
344			}
345
346			bool there_are_degree1_vertices = true;
347		✗	while(there_are_degree1_vertices) {
348			there_are_degree1_vertices = false;
349		✗	for(index_t c1: mesh->facet_corners) {
350			index_t v1 = mesh->facet_corners.vertex(c1);
351		✗	if(v_nb_borders[v1] == 1 &&
352		✗	(c_on_border[c1] != 0) && R_ff[c1] == 0
353			) {
354			index_t f2 = mesh->facet_corners.adjacent_facet(c1);
355		✗	if(f2 != NO_FACET) {
356			there_are_degree1_vertices = true;
357		✗	index_t f1 = c1/3;
358			index_t c2 = mesh->facets.corners_begin(f2) +
359		✗	mesh->facets.find_adjacent(f2,f1);
360			index_t v2 = mesh->facet_corners.vertex(c2);
361		✗	c_on_border[c1] = 0;
362		✗	c_on_border[c2] = 0;
363		✗	--v_nb_borders[v1];
364		✗	--v_nb_borders[v2];
365			}
366			}
367			}
368			}
369		✗	}
370
371
372		✗	void do_the_ball_no_brush_no_zip(
373			Mesh* mesh, Attribute<index_t>& c_on_border
374			) {
375		✗	for(index_t c: mesh->facet_corners) {
376		✗	c_on_border[c] = 1;
377			}
378
379			// Covering tree
380
381		✗	std::vector<bool> visited(mesh->facets.nb(),false);
382			std::deque<index_t> S;
383		✗	S.push_back(0);
384			visited[0] = true;
385		✗	while(!S.empty()) {
386		✗	index_t f1 = S.front();
387		✗	S.pop_front();
388		✗	for(index_t c1: mesh->facets.corners(f1)) {
389		✗	index_t f2 = mesh->facet_corners.adjacent_facet(c1);
390		✗	if(f2 != NO_FACET && !visited[f2]) {
391		✗	set_edge_attr(mesh, c_on_border, f1, f2, 0);
392			visited[f2] = true;
393			S.push_back(f2);
394			}
395			}
396			}
397		✗	}
398
399		✗	void get_B_on_edge(
400			Mesh* mesh, Attribute<vec3>& B, Attribute<index_t>& R_ff,
401			index_t f, index_t c,
402			vec3& Bc, vec3& BTc
403			) {
404		✗	vec3 Nf = normalize(Geom::mesh_facet_normal(*mesh,f));
405		✗	Bc = B[f];
406		✗	BTc = cross(Nf,Bc);
407			index_t f2 = mesh->facet_corners.adjacent_facet(c);
408		✗	if(f2 != NO_FACET) {
409		✗	vec3 N2 = normalize(Geom::mesh_facet_normal(*mesh,f2));
410		✗	vec3 B2 = B[f2];
411			vec3 BT2 = cross(N2,B2);
412		✗	FOR(i,R_ff[c]) {
413			B2 = cross(N2,B2);
414			BT2 = cross(N2,BT2);
415			}
416			Bc += B2;
417			BTc += BT2;
418			}
419		✗	Bc = normalize(Bc);
420		✗	BTc = normalize(BTc);
421		✗	}
422
423		✗	void get_constraints(
424			Mesh* mesh, Attribute<vec3>& B, Attribute<index_t>& R_ff,
425			Attribute<index_t>& constraint
426			) {
427
428			geo_argused(R_ff);
429
430		✗	for(index_t c: mesh->facet_corners) {
431		✗	constraint[c] = CNSTR_NONE;
432			}
433
434		✗	for(index_t c: mesh->facet_corners) {
435		✗	index_t edge_constraints = get_edge_constraints(mesh,c,B);
436		✗	index_t f = c/3;
437			index_t c2 = mesh->facets.next_corner_around_facet(f,c);
438		✗	constraint[c] |= edge_constraints;
439		✗	constraint[c2] |= edge_constraints;
440			}
441
442		✗	return;
443
444			/*
445			// Propagate the constraints: all the corners incident to a
446			// vertex that is itself incident to a constrained edge are
447			// constrained.
448			// Normally, I think this would not be required, since the
449			// (u,v) compatibility constraint + the mutiplity constraint of
450			// one of the (u,v)'s imply that all the corners incident to the
451			// considered vertex should have integer coordinates,
452			// however if I do not set this constraint I observed that it
453			// does not work as expected (note: interestingly, when I
454			// added the wheel compatibility constraint there was a big
455			// improvement, but it did not solve all issues).
456			std::stack<index_t> S;
457			std::vector<bool> is_visited(mesh->facet_corners.nb(),false);
458
459			for(index_t c=0; c<mesh->facet_corners.nb(); ++c) {
460			if(constraint[c] != CNSTR_NONE) {
461			S.push(c);
462			is_visited[c] = true;
463			}
464			}
465
466			while(!S.empty()) {
467			index_t c = S.top();
468			S.pop();
469			index_t f = c/3;
470			index_t cprev = mesh->facets.prev_corner_around_facet(f,c);
471			index_t fneigh = mesh->facet_corners.adjacent_facet(cprev);
472			if(fneigh == NO_FACET) {
473			continue;
474			}
475			index_t eneigh = mesh->facets.find_adjacent(fneigh,f);
476			index_t cneigh =
477			mesh->facets.corners_begin(fneigh) + eneigh;
478
479			if(is_visited[cneigh]) {
480			continue;
481			}
482
483			bool cu = (constraint[c] & CNSTR_U) != 0;
484			bool cv = (constraint[c] & CNSTR_V) != 0;
485
486			index_t Rij = R_ff[cprev];
487
488			// If rotation is 90 degrees or 270 degrees, then
489			// u and v are swapped.
490			if((Rij & 1) != 0) {
491			std::swap(cu,cv);
492			}
493
494			if(cu) {
495			constraint[cneigh] |= CNSTR_U;
496			}
497
498			if(cv) {
499			constraint[cneigh] |= CNSTR_V;
500			}
501
502			is_visited[cneigh]=true;
503			S.push(cneigh);
504			}*/
505			}
506
507		✗	index_t get_edge_constraints(
508			Mesh* mesh, index_t c, Attribute<vec3>& B
509			) {
510			index_t result = 0;
511
512		✗	index_t f = c/3;
513		✗	vec3 N = normalize(Geom::mesh_facet_normal(*mesh,f));
514
515			index_t f2 = mesh->facet_corners.adjacent_facet(c);
516		✗	if(f2 != NO_FACET) {
517		✗	if(
518		✗	::fabs(Geom::mesh_normal_angle(*mesh,c)) * 180.0 / M_PI
519			< 45.0
520			) {
521			return 0;
522			}
523			}
524
525			index_t v1 = mesh->facet_corners.vertex(c);
526			index_t c2 = mesh->facets.next_corner_around_facet(c/3,c);
527			index_t v2 = mesh->facet_corners.vertex(c2);
528			vec3 E = mesh->vertices.point(v2) - mesh->vertices.point(v1);
529		✗	vec3 Bf = normalize(B[f]);
530			vec3 Bfrot = cross(N,Bf);
531
532		✗	double a1 = (Geom::angle(E,Bf)) * 180.0 / M_PI;
533		✗	a1 = std::min(a1, 180.0-a1);
534		✗	if(a1 < 10.0) {
535			result |= GlobalParam2d::Internal::CNSTR_V;
536			}
537
538		✗	double a2 = (Geom::angle(E,Bfrot)) * 180.0 / M_PI;
539		✗	a2 = std::min(a2, 180.0-a2);
540		✗	if(a2 < 10.0) {
541		✗	result |= GlobalParam2d::Internal::CNSTR_U;
542			}
543
544			// geo_assert!(a1 < 10.0 && a2 < 10.0));
545			// Should not occur...
546		✗	if(a1 < 10.0 && a2 < 10.0) {
547			result = 0;
548			}
549
550			return result;
551			}
552
553
554		✗	void snap_tex_coord(double& coord) {
555		✗	double snapped = GEO::round(coord);
556		✗	if(std::fabs(coord - snapped) < 0.05) {
557		✗	coord = snapped;
558			}
559		✗	}
560
561		✗	index_t inverse_R(index_t R) {
562		✗	geo_assert(R < 4);
563			static index_t inverse[4] = {
564			0,
565			3,
566			2,
567			1
568			};
569		✗	return inverse[R];
570			}
571
572
573		✗	void transfer_B_to_vertices(
574			Mesh* mesh,
575			Attribute<vec3>& B, Attribute<vec3>& Bv,
576			Attribute<index_t>& R_fv
577			) {
578		✗	for(index_t v: mesh->vertices) {
579		✗	Bv[v] = vec3(0.0, 0.0, 0.0);
580			}
581		✗	for(index_t c: mesh->facet_corners) {
582			index_t v = mesh->facet_corners.vertex(c);
583		✗	index_t f = c/3;
584		✗	vec3 Bf = normalize(B[f]);
585		✗	vec3 N = normalize(Geom::mesh_facet_normal(*mesh, f));
586		✗	FOR(k, R_fv[c]) {
587			Bf = cross(N,Bf);
588			}
589			Bv[v] += Bf;
590			}
591		✗	for(index_t v: mesh->vertices) {
592		✗	Bv[v] = normalize(Bv[v]);
593			}
594		✗	}
595			} // namespace Internal
596
597		✗	void frame_field(
598			Mesh* mesh, Attribute<vec3>& B,
599			double hard_angle_threshold
600			) {
601			FrameField FF;
602			// We will directly query the field on the facets, no need
603			// for the KD-tree.
604			FF.set_use_spatial_search(false);
605		✗	FF.create_from_surface_mesh(*mesh,false,hard_angle_threshold);
606			const vector<double>& frames = FF.frames();
607		✗	for(index_t f: mesh->facets) {
608		✗	B[f] = vec3(
609			frames[9*f+0],
610			frames[9*f+1],
611			frames[9*f+2]
612			);
613			}
614		✗	}
615
616			} // namespace GlobalParam2d
617			} // namespace OGF
618