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/hexdom/quads_from_boundary.h>
42			#include <exploragram/hexdom/intersect_tools.h>
43			#include <exploragram/hexdom/polygon.h>
44
45			namespace GEO {
46
47
48			struct QuadsFromBoundry {
49			index_t n;
50			vector<vec2>& pts;
51			double ave;
52			vector<index_t> global_vid[2]; // vertex indices of the sub problem
53			vector<index_t>& quads;
54
55			QuadsFromBoundry(vector<vec2>& p_pts, vector<index_t>& p_quads) : pts(p_pts), quads(p_quads){}
56
57			double angle(int i) {
58			vec2 P[3];
59			FOR(p, 3) P[p] = aupp(i + int(p) - 1, pts);
60			return (180. / M_PI)*atan2(det(P[2] - P[1], P[0] - P[1]), dot(P[2] - P[1], P[0] - P[1]));
61			}
62
63			double angle(index_t i) {
64			return angle(int(i));
65			}
66
67			vec2 point(int i) {
68			return aupp(i , pts);
69			}
70
71			vec2 point(index_t i) {
72			return aupp(i , pts);
73			}
74
75			void init() {
76			n = pts.size();
77			ave = 0;
78			FOR(p, pts.size()) ave += (pts[p] - aupp(p + 1, pts)).length();
79			ave /= double(pts.size());
80			}
81
82
83			bool solve_and_merge_subproblems() {
84			plop("recurs");
85			// solve on two halves
86			vector<vec2> poly[2];
87			FOR(i, 2) FOR(fv, global_vid[i].size()) poly[i].push_back(pts[global_vid[i][fv]]);
88
89			vector<index_t> poly_quad[2];
90
91			FOR(i, 2) {
92			QuadsFromBoundry sub(poly[i], poly_quad[i]);
93			if (!sub.apply()) return false;
94			plop(pts.size());
95			plop(global_vid[i].size());
96			}
97			// add new pts to global
98			FOR(i, 2) for (index_t d = global_vid[i].size(); d < poly[i].size(); d++) {
99			global_vid[i].push_back(pts.size());
100			pts.push_back(poly[i][d]);
101			}
102			FOR(i, 2) FOR(qu, poly_quad[i].size()) quads.push_back(global_vid[i][poly_quad[i][qu]]);
103
104			return true;
105
106			}
107
108
109			bool apply() {
110			init();
111			if (n == 4) { FOR(v, 4) quads.push_back(v); return true; }
112
113
114
115
116			// try to find a quad to close
117			FOR(v, n) {
118			double alpha[2] = { angle(v),angle((v + 1) % n) };
119			double tolerance = 45;
120			if (alpha[0] > 90 - tolerance && alpha[0] < 90 + tolerance
121			&& alpha[1] > 90 - tolerance && alpha[1] < 90 + tolerance) {
122			FOR(i, n) if (i != v && (i != ((v + 1) % n))) global_vid[0].push_back(i);
123			FOR(i, 4) global_vid[1].push_back((v + i + n - 1) % n);
124			return solve_and_merge_subproblems();
125			//goto split_done;
126			}
127			}
128
129
130			// try to find a new vertex to punch
131			FOR(it, 2)
132			FOR(v, n) {
133			double alpha[3] = { angle(v),angle((v + 1) % n),angle((v + 2) % n) };
134			double tolerance = 45;
135			bool can_punch;
136			if (it == 0)
137			can_punch = alpha[0] > 90 - tolerance && alpha[0] < 90 + tolerance
138			&& alpha[1] > 0 - tolerance && alpha[1] < 0 + tolerance
139			&& alpha[2] > 90 - tolerance && alpha[2] < 90 + tolerance
140			;
141			else //(it == 1)
142			can_punch = alpha[0] > 90 - tolerance && alpha[0] < 90 + tolerance;
143			if (can_punch) {// more or less 90 degree ;)
144			vec2 nvP;
145			vec2 vec[2] = { point(v + 1) - point(v),point(v - 1) - point(v)};
146			mat2 mat;
147			FOR(i, 2)FOR(j, 2) mat(i, j) = vec[i][j];
148			mat2 inv = mat.inverse();
149			vec2 b(vec[0].length2(), vec[1].length2());
150			mult(inv, b.data(), nvP.data());
151			nvP = nvP + point(v);
152			nvP = point(v) + vec[0] + vec[1];
153
154			// check that there is no existing vertex
155			bool geometric_issue = false;
156			{
157			vec2 new_quad[4] = { point(v - 1),point(v),point(v + 1),nvP };
158			for (double du = 0; du < .3; du += .2)
159			for (double dv = .8; dv < 1; dv += .2) {
160			vec2 pixel =
161			du*dv*new_quad[0]
162			+ du*(1. - dv)*new_quad[1]
163			+ (1. - du)*dv*new_quad[3]
164			+ (1. - du)*(1. - dv)*new_quad[2];
165			for (double c = 0; c < 1.; c += .2) FOR(vv, n)
166			if (((c*aupp(vv + 1, pts) + (1. - c)*point(vv)) - pixel).length() < .5*ave) geometric_issue = true;
167			}
168			}
169			if (geometric_issue) continue;
170
171
172			FOR(i, n) global_vid[0].push_back((i != v) ? ((i) % n) : pts.size());
173			FOR(i, 3) global_vid[1].push_back((v + i + n - 1) % n);
174			global_vid[1].push_back(pts.size());
175			pts.push_back(nvP);
176			return solve_and_merge_subproblems();
177			//goto split_done;
178			}
179			}
180
181			return true;
182			// split_done:
183			}
184			};
185
186
187			//double angle(vector<vec2>& pts, int i) {
188			// vec2 P[3];
189			// FOR(p, 3) P[p] = aupp(i + p - 1, pts);
190			// return (180. / M_PI)*atan2(det(P[2] - P[1], P[0] - P[1]), dot(P[2] - P[1], P[0] - P[1]));
191			//}
192			//double ave_length(vector<vec2>& pts) {
193			// double res = 0;
194			// FOR(p, pts.size()) res += (pts[p] - aupp(p + 1, pts)).length();
195			// return res/double(pts.size());
196			//}
197
198
199
200			//bool try_quadrangulate_with_punch_vertex(vector<vec2>& P, vector<index_t>& quads) {
201			// int n = P.size();
202			// if (n == 4) { FOR(v, 4) quads.push_back(v); return true; }
203			// vector<index_t> global_vid[2];
204
205
206			// // try to find a quad to close
207			// FOR(v, n) {
208			// double alpha[2] = { angle(P, v),angle(P, (v + 1) % n) };
209			// double tolerance = 45;
210			// if (alpha[0] > 90-tolerance && alpha[0] < 90 + tolerance
211			// && alpha[1] > 90 - tolerance && alpha[1] < 90 + tolerance) {
212			// // check that there is no existing vertex
213			// //bool geometric_issue = false;
214			// //for (double c = 0; c < 1.; c += .2) FOR(vv, n) if (((c*aupp(vv + 1, P) + (1. - c)*P[vv]) - nvP).length() < .5*ave_length(P)) geometric_issue = true;
215			// //if (geometric_issue) continue;
216			// plop("go");
217			// FOR(i, n) if (i != v && (i != ((v + 1) % n))) global_vid[0].push_back(i);
218			// FOR(i, 4) global_vid[1].push_back((v + i + n - 1) % n);
219			// goto split_done;
220			// }
221			// }
222
223
224			// // try to find a new vertex to punch
225			// FOR(it,2)
226			// FOR(v, n) {
227			// double alpha[3] = { angle(P, v),angle(P, (v + 1) % n),angle(P, (v + 2) % n) };
228			// double tolerance = 45;
229			// bool can_punch ;
230			// if (it == 0)
231			// can_punch = alpha[0] > 90 - tolerance && alpha[0] < 90 + tolerance
232			// && alpha[1] > 0 - tolerance && alpha[1] < 0 + tolerance
233			// && alpha[2] > 90 - tolerance && alpha[2] < 90 + tolerance
234			// ;
235			// if (it == 1)
236			// can_punch = alpha[0] > 90 - tolerance && alpha[0] < 90 + tolerance;
237			// if (can_punch) {// more or less 90 degree ;)
238			// vec2 nvP;
239			// vec2 vec[2] = { aupp(v + 1, P) - P[v],aupp(v - 1, P) - P[v] };
240			// mat2 mat;
241			// FOR(i,2)FOR(j,2) mat(i, j) = vec[i][j];
242			// mat2 inv = mat.inverse();
243			// vec2 b (vec[0].length2(), vec[1].length2());
244			// mult(inv, b.data(), nvP.data());
245			// nvP = nvP + P[v];
246			// nvP = P[v] + vec[0] + vec[1];
247
248			// // check that there is no existing vertex
249			// bool geometric_issue = false;
250			// {
251			// double ave = ave_length(P);
252			// vec2 new_quad[4] = { aupp(v - 1, P),aupp(v , P),aupp(v + 1, P),nvP };
253			// for (double du = 0; du < .3; du += .2)
254			// for (double dv = .8; dv < 1; dv += .2) {
255			// vec2 pixel =
256			// du*dv*new_quad[0]
257			// + du*(1. - dv)*new_quad[1]
258			// + (1. - du)*dv*new_quad[3]
259			// + (1. - du)*(1. - dv)*new_quad[2];
260			// for (double c = 0; c < 1.; c += .2) FOR(vv, n) if (((c*aupp(vv + 1, P) + (1. - c)*P[vv]) - pixel).length() < .5*ave) geometric_issue = true;
261			// }
262			// }
263			// if (geometric_issue) continue;
264
265			//
266			// FOR(i, n) global_vid[0].push_back((i != v) ? ((i)%n) : P.size());
267			// FOR(i, 3) global_vid[1].push_back((v+i+n-1)%n);
268			// global_vid[1].push_back(P.size());
269			// P.push_back(nvP);
270			// goto split_done;
271			// }
272			// }
273
274			// return true;
275			// split_done:
276			// // solve on two halves
277			// vector<vec2> poly[2];
278			// FOR(i, 2) FOR(fv, global_vid[i].size()) poly[i].push_back(P[global_vid[i][fv]]);
279
280			// vector<index_t> poly_quad[2];
281			// FOR(i, 2) if (!try_quadrangulate_with_punch_vertex(poly[i],poly_quad[i])) return false;
282			//
283			// // add new pts to global
284			// FOR(i, 2) for (int d = global_vid[i].size(); d < poly[i].size(); d++) {
285			// global_vid[i].push_back(P.size());
286			// P.push_back(poly[i][d]);
287			// }
288			// FOR(i, 2) FOR(qu, poly_quad[i].size()) quads.push_back(global_vid[i][poly_quad[i][qu]]);
289
290			// return true;
291			//}
292
293
294
295
296
297
298		✗	bool try_quadrangulate(vector<vec2>& pts, vector<index_t>& quads) {
299
300		✗	return Poly2d(pts).try_quadrangulate(quads);
301
302			//QuadsFromBoundry sub(pts, quads);
303			//return sub.apply();
304
305			}
306			}
307