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/intersect_tools.h>
41			#include <exploragram/hexdom/mesh_utils.h>
42			#define FPG_UNCERTAIN_VALUE 0
43			#include <geogram/numerics/predicates/orient3d.h>
44			#include <geogram/mesh/triangle_intersection.h>
45			#include <geogram/mesh/mesh_io.h>
46			namespace {
47			using namespace GEO;
48
49			struct IndexPointCmp {
50			IndexPointCmp(vector<vec3>& p_data, index_t p_dim)
51			: dim(p_dim), data(p_data) {}
52			bool operator()(const index_t A, const index_t B) {
53		✗	return data[A][dim] > data[B][dim];
54			}
55			index_t dim;
56			vector<vec3>& data;
57			};
58			}
59
60			namespace GEO {
61
62
63			const index_t quad_rand_split[2][3] = { { 0, 1, 2 },{ 0, 2, 3 } }; // random triangulation
64			const index_t quad_split[4][3] = { { 0, 1, 2 },{ 0, 2, 3 },{ 1, 2, 3 },{ 1, 3, 0 } }; // convex hull
65			const index_t diamon_split[12][3] = { // surface and inside
66			{ 4, 5, 0 },{ 4, 5, 1 },{ 4, 5, 2 },{ 4, 5, 3 },
67			{ 0, 1, 4 },{ 1, 2, 4 },{ 2, 3, 4 },{ 3, 0, 4 },
68			{ 0, 1, 5 },{ 1, 2, 5 },{ 2, 3, 5 },{ 3, 0, 5 } // orientation doesn't matter
69			};
70
71			const index_t diamon_quad_split[16][3] = { // surface and inside
72			{ 4, 5, 0 },{ 4, 5, 1 },{ 4, 5, 2 },{ 4, 5, 3 },
73			{ 0, 1, 4 },{ 1, 2, 4 },{ 2, 3, 4 },{ 3, 0, 4 },
74			{ 0, 1, 5 },{ 1, 2, 5 },{ 2, 3, 5 },{ 3, 0, 5 }, // orientation doesn't matter
75			{ 0, 1, 2 },{ 0, 2, 3 },{ 1, 2, 3 },{ 1, 3, 0 }
76			};
77
78
79
80
81		✗	bool BBox::intersect(const BBox& b) const {
82		✗	FOR(d, 3) {
83		✗	if (min[d] > b.max[d] || max[d] < b.min[d]) {
84			return false;
85			}
86			}
87			return true;
88			}
89
90		✗	bool BBox::contains(const vec3& v) const {
91		✗	FOR(d, 3) {
92		✗	if (min[d] > v[d] || max[d] < v[d]) {
93			return false;
94			}
95			}
96			return true;
97			}
98
99		✗	bool BBox::is_null() const {
100		✗	FOR(d, 3) {
101		✗	if (max[d] - min[d] < 0) {
102			return true;
103			}
104			}
105			return false;
106			}
107
108		✗	void BBox::add(const BBox& b) {
109		✗	if (b.is_null()) return;
110		✗	add(b.min);
111		✗	add(b.max);
112			}
113
114		✗	void BBox::add(const vec3& P) {
115		✗	FOR(d, 3) {
116		✗	min[d] = std::min(min[d], P[d]);
117		✗	max[d] = std::max(max[d], P[d]);
118			}
119		✗	}
120
121		✗	vec3 BBox::bary() const {
122		✗	return 0.5*(min + max);
123			}
124
125
126			/**********************************************************/
127
128			inline unsigned int mylog2( unsigned int x ) {
129			unsigned int ans = 0 ;
130		✗	while( x>>=1 ) ans++;
131			return ans ;
132			}
133
134
135		✗	void HBoxes::init(vector<BBox>& inboxes) {
136			vector<vec3> G(inboxes.size());
137		✗	tree_pos_to_org.resize(inboxes.size());
138		✗	FOR(p, G.size()) G[p] = inboxes[p].bary();
139		✗	FOR(p, G.size()) tree_pos_to_org[p] = p;
140		✗	sort(G, 0, tree_pos_to_org.size());
141
142		✗	offset = index_t(pow(2.0, 1.0 + mylog2(G.size()))) - 1;
143		✗	tree.resize(offset + G.size());
144		✗	FOR(i, G.size()) tree[offset + i] = inboxes[tree_pos_to_org[i]];
145		✗	for (int i = int(offset) - 1; i >= 0; i--) {
146		✗	for (int son = 2 * i + 1; son < 2 * i + 3; son++)
147		✗	if (son < int(tree.size())) tree[i].add(tree[son]);
148			}
149
150		✗	STAT_nb_visits = 0;
151		✗	STAT_nb_leafs = 0;
152		✗	STAT_nb_requests = 0;
153		✗	}
154
155		✗	void HBoxes::sort(vector<vec3> &G, index_t org, index_t dest) {
156
157			// find the best dim to cut
158			index_t dim = 2;
159			BBox b;
160		✗	for (index_t i = org; i < dest; i++) b.add(G[tree_pos_to_org[i]]);
161		✗	FOR(d, 2) if (b.max[d] - b.min[d] > b.max[dim] - b.min[dim]) dim = d;
162			// sort
163			IndexPointCmp cmp(G, dim);
164		✗	std::sort(tree_pos_to_org.begin() + int(org), tree_pos_to_org.begin() + int(dest), cmp);
165		✗	if (dest - org <= 2) return;
166		✗	index_t m = org + index_t(pow(2.0, int(mylog2(dest - org - 1))));
167		✗	sort(G, org, m);
168		✗	sort(G, m, dest);
169			}
170
171		✗	void HBoxes::intersect(BBox& b, vector<index_t>& primitives, index_t node) {
172		✗	if (node == 0) STAT_nb_requests++;
173		✗	geo_assert(node < tree.size());
174		✗	STAT_nb_visits++;
175		✗	if (!tree[node].intersect(b)) return;
176		✗	if (node >= offset) {
177		✗	STAT_nb_leafs++;
178		✗	primitives.push_back(tree_pos_to_org[node - offset]);
179			} else {
180		✗	for (index_t son = 2 * node + 1; son < 2 * node + 3; son++)
181		✗	if (son < tree.size())
182		✗	intersect(b, primitives, son);
183			}
184			}
185
186			/**********************************************************/
187
188		✗	void DynamicHBoxes::init(vector<BBox>& inboxes) {
189		✗	hbox.init(inboxes);
190			moved.clear();
191			movedbbox.clear();
192		✗	}
193
194		✗	void DynamicHBoxes::intersect(BBox& b, vector<index_t>& primitives) {
195		✗	hbox.intersect(b, primitives);
196		✗	FOR(i, moved.size()) {
197		✗	if (movedbbox[i].intersect(b))
198		✗	primitives.push_back(moved[i]);
199			}
200		✗	}
201
202		✗	void DynamicHBoxes::update_bbox(index_t id, BBox b) {
203		✗	geo_assert(id<hbox.tree_pos_to_org.size());
204		✗	FOR(i, moved.size()) {
205		✗	if (moved[i] == id) {
206		✗	movedbbox[i] = b;
207		✗	return;
208			}
209			}
210		✗	moved.push_back(id);
211		✗	movedbbox.push_back(b);
212			}
213
214
215
216
217		✗	static double tetra_volume(vec3 A, vec3 B, vec3 C, vec3 D) {
218		✗	return dot(cross(B - A, C - A), D - A);
219			}
220
221		✗	double tetra_volume_sign(vec3 A, vec3 B, vec3 C, vec3 D) {
222		✗	double res = tetra_volume(A, B, C, D);
223		✗	if (std::abs(res) > 1e-15) return res;
224		✗	return dot(normalize(cross(normalize(B - A), normalize(C - A))), normalize(D - A));
225			}
226
227		✗	bool same_sign(double a, double b) { return (a > 0) == (b > 0); }
228
229
230
231		✗	vector<BBox> facets_bbox(Mesh* m) {
232			vector<BBox> inboxes(m->facets.nb());
233		✗	FOR(f, m->facets.nb()) {
234			index_t nbv = m->facets.nb_vertices(f);
235		✗	FOR(fv, nbv) inboxes[f].add(X(m)[m->facets.vertex(f, fv)]);
236			}
237		✗	return inboxes;
238			}
239
240
241
242
243		✗	FacetIntersect::FacetIntersect(Mesh* p_m) { m = p_m;
244		✗	inboxes = facets_bbox(m);
245		✗	hb.init(inboxes);
246		✗	}
247
248		✗	static void save_conflict(std::string name, vec3 A0, vec3 B0, vec3 C0, vec3 A1, vec3 B1, vec3 C1) {
249		✗	Mesh conflict;
250			conflict.vertices.create_vertices(6);
251			conflict.facets.create_triangles(2);
252		✗	X(&conflict)[0] = A0; X(&conflict)[1] = B0; X(&conflict)[2] = C0;
253		✗	X(&conflict)[3] = A1; X(&conflict)[4] = B1; X(&conflict)[5] = C1;
254		✗	FOR(f, 2) FOR(lv, 3) conflict.facets.set_vertex(f, lv, 3 * f + lv);
255		✗	mesh_save(conflict, "C:/DATA/debug/" + name + "conflict.geogram");
256		✗	}
257
258		✗	static bool polyintersect_both_triangulation(vector<vec3>& P, vector<vec3>& Q) {
259			bool conflict = false;
260		✗	FOR(trP, 4) {
261		✗	FOR(trQ, 4) {
262		✗	if (trP > 0 && P.size() == 3) continue;
263		✗	if (trQ > 0 && Q.size() == 3) continue;
264			TriangleIsects trash;
265		✗	conflict = conflict || triangles_intersections(
266		✗	P[quad_split[trP][0]], P[quad_split[trP][1]], P[quad_split[trP][2]],
267		✗	Q[quad_split[trQ][0]], Q[quad_split[trQ][1]], Q[quad_split[trQ][2]],
268			trash
269			);
270			FOR(i, trash.size()) {
271			//if (trash[i].first > 2 || trash[i].second > 2) conflict = true;
272			}
273			// static int nb_intersects = 0;
274		✗	if (conflict) {
275			// nb_intersects++;
276			// save_conflict("gna"+ String::to_string(nb_intersects), P[quad_split[trP][0]], P[quad_split[trP][1]], P[quad_split[trP][2]],
277			// Q[quad_split[trQ][0]], Q[quad_split[trQ][1]], Q[quad_split[trQ][2]]);
278		✗	return true;
279			}
280			}
281			}
282			return false;
283			}
284		✗	bool polyintersect(vector<vec3>& P, vector<vec3>& Q) {
285		✗	geo_assert(P.size() == 3 || P.size() == 4);
286		✗	geo_assert(Q.size() == 3 || Q.size() == 4);
287
288			// check for same facet
289		✗	if (P.size() == Q.size()) FOR(off, P.size()) {
290			bool is_same = true;
291		✗	FOR(v, P.size()) {
292		✗	if ((P[v] - Q[(v + off) % P.size()]).length2() != 0) {
293			is_same = false;
294			break;
295			}
296			}
297		✗	if (is_same) return false;
298			}
299		✗	return polyintersect_both_triangulation(P, Q);
300			}
301
302		✗	vector<index_t> FacetIntersect::get_intersections(vector<vec3>& P) {
303			vector<index_t> res;
304			vector<index_t> primitives;
305			BBox request_bbox;
306		✗	FOR(v, P.size()) request_bbox.add(P[v]);
307			request_bbox.dilate(1e-15);
308		✗	hb.intersect(request_bbox, primitives);
309		✗	FOR(i, primitives.size()) {
310		✗	index_t opp_f = primitives[i];
311			vector<vec3> Q;
312		✗	FOR(fv, m->facets.nb_vertices(opp_f))
313		✗	Q.push_back(X(m)[m->facets.vertex(opp_f, fv)]);
314		✗	if (polyintersect(P, Q) || polyintersect(Q,P))
315		✗	res.push_back(opp_f);
316			}
317		✗	return res;
318			}
319		✗	vector<index_t> FacetIntersect::get_intersections(index_t& f) {
320		✗	vector<vec3> verts(m->facets.nb_vertices(f));
321		✗	FOR(v, m->facets.nb_vertices(f)) verts[v] = X(m)[m->facets.vertex(f,v)];
322		✗	return get_intersections(verts);
323			}
324
325
326		✗	vector<index_t> get_intersecting_faces(Mesh* m) {
327			vector<index_t> res;
328		✗	FacetIntersect finter(m);
329		✗	FOR(f, m->facets.nb()) {
330			// if ((f%50) ==0)plop(double(f) / double(m->facets.nb()));
331		✗	vector<index_t> opp = finter.get_intersections(f);
332		✗	FOR(i, opp.size()) {
333		✗	res.push_back(f);
334			res.push_back(opp[i]);
335			}
336			}
337		✗	return res;
338		✗	}
339
340
341
342
343		✗	void check_no_intersecting_faces(Mesh* m, bool allow_duplicated ) {
344		✗	vector<index_t> intersect = get_intersecting_faces(m);
345		✗	if (allow_duplicated) {
346			index_t offpair = 0;
347		✗	while (offpair < intersect.size()) {
348		✗	index_t f0 = intersect[offpair];
349		✗	index_t f1 = intersect[offpair + 1];
350		✗	if (f0 != f1 && (facet_bary(m, f0) - facet_bary(m, f1)).length2() < 1e-15) {
351		✗	FOR(d, 2) std::swap(intersect[offpair + d], intersect[intersect.size() - 2 + d]);
352		✗	FOR(d, 2) intersect.pop_back();
353			}
354		✗	else offpair += 2;
355			}
356			}
357		✗	if (intersect.empty()) return;
358		✗	Attribute<int> intersection(m->facets.attributes(), "intersection");
359		✗	FOR(f, m->facets.nb()) intersection[f] = 0;
360		✗	FOR(i, intersect.size()) intersection[intersect[i]] = 1;
361		✗	intersection[intersect[0]] = 2;
362		✗	intersection[intersect[1]] = 3;
363		✗	FOR(i, 2) plop(m->facets.nb_vertices(intersect[i]));
364		✗	FOR(i, 2) FOR(lv, m->facets.nb_vertices(intersect[i])) plop(m->facets.vertex(intersect[i], lv));
365		✗	FOR(i, 2) FOR(lv, m->facets.nb_vertices(intersect[i])) plop(X(m)[m->facets.vertex(intersect[i], lv)]);
366		✗	mesh_save(*m, "C:/DATA/debug/intersectingsurface.geogram");
367
368		✗	save_conflict("first",
369		✗	X(m)[m->facets.vertex(intersect[0], 0)],
370		✗	X(m)[m->facets.vertex(intersect[0], 1)],
371		✗	X(m)[m->facets.vertex(intersect[0], 2)],
372		✗	X(m)[m->facets.vertex(intersect[1], 0)],
373		✗	X(m)[m->facets.vertex(intersect[1], 1)],
374		✗	X(m)[m->facets.vertex(intersect[1], 2)]
375			);
376		✗	Mesh conflict;
377			conflict.vertices.create_vertices(6);
378			conflict.facets.create_triangles(2);
379		✗	FOR(f, 2) FOR(lv, 3) X(&conflict)[3*f+lv] = X(m)[m->facets.vertex(intersect[f], lv)];
380		✗	FOR(f, 2) FOR(lv, 3) conflict.facets.set_vertex(f, lv, 3 * f + lv);
381		✗	mesh_save(conflict, "C:/DATA/debug/conflict.geogram");
382		✗	geo_assert_not_reached;
383		✗	}
384
385
386			}
387