GCC Code Coverage Report

Line	Branch	Exec	Source
1			/*
2			* Copyright (c) 2000-2025 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/en/bruno-levy-1
32			*
33			* Inria,
34			* Domaine de Voluceau,
35			* 78150 Le Chesnay - Rocquencourt
36			* FRANCE
37			*
38			*/
39
40			#include <geogram/mesh/mesh_minkowski.h>
41			#include <geogram/mesh/mesh.h>
42			#include <geogram/mesh/mesh_geometry.h>
43			#include <geogram/mesh/mesh_convex_hull.h>
44			#include <geogram/mesh/mesh_surface_intersection.h>
45			#include <geogram/mesh/mesh_topology.h>
46			#include <geogram/mesh/mesh_fill_holes.h>
47			#include <geogram/delaunay/delaunay.h>
48			#include <geogram/numerics/predicates.h>
49
50			// Implement an exact version of Hichem Barki's algorithm:
51			// 1) 2D convex hull: implement QuickHull, make a generic version that has
52			// its own orient2d_compare() predicate has a "plug-in"
53
54			namespace {
55			using namespace GEO;
56
57			/**
58			* \brief Convexity of an edge in a 3D mesh
59			* \param[in] M a const reference to the mesh
60			* \param[in] f1 an edge facet
61			* \param[in] le1 a local edge index in \p f1
62			* \param[in] angle_tolerance maximum normal angle deviation for
63			* considering that the two adjacent facets are co-planar
64			* \retval POSITIVE if the edge is convex
65			* \retval ZERO if the two facets adjacent to the edge are coplanar,
66			* up to \p angle_tolerance
67			* \retval NEGATIVE if the edge is concave
68			*/
69		✗	Sign edge_convexity(
70			const Mesh& M, index_t f1, index_t le1,
71			double angle_tolerance = 0.0
72			) {
73			// Tiny tolerance (0.01 degrees) for nearly co-planar facets (like
74			// square facets in CSGBuilder's tesselated spheres).
75			index_t n1 = M.facets.nb_vertices(f1);
76			index_t f2 = M.facets.adjacent(f1,le1);
77			geo_debug_assert(f2 != NO_INDEX);
78			index_t n2 = M.facets.nb_vertices(f2);
79		✗	index_t le2 = M.facets.find_adjacent(f2,f1);
80		✗	vec3 p1 = M.vertices.point(M.facets.vertex(f1,le1));
81		✗	vec3 p2 = M.vertices.point(M.facets.vertex(f1,(le1+1)%n1));
82		✗	vec3 p3 = M.vertices.point(M.facets.vertex(f1,(le1+2)%n1));
83		✗	vec3 p4 = M.vertices.point(M.facets.vertex(f2,(le2+2)%n2));
84		✗	if(angle_tolerance != 0.0) {
85			vec3 N1 = cross(p2-p1, p3-p1);
86			vec3 N2 = cross(p4-p1, p2-p1);
87			double a1 = Geom::angle(N1,N2);
88		✗	if(::fabs(a1) < angle_tolerance) {
89		✗	return ZERO;
90			}
91			}
92			return PCK::orient_3d(p1,p2,p3,p4);
93			}
94
95		✗	bool mesh_is_convex_3d(const Mesh& M, double angle_tolerance = 0.0) {
96			index_t nb_positive = 0;
97			index_t nb_negative = 0;
98		✗	for(index_t f: M.facets) {
99		✗	for(index_t le=0; le<M.facets.nb_vertices(f); ++le) {
100			index_t g = M.facets.adjacent(f,le);
101			geo_debug_assert(g != NO_INDEX);
102		✗	if(f > g) {
103		✗	Sign s = edge_convexity(M,f,le,angle_tolerance);
104		✗	nb_positive += (s > 0);
105		✗	nb_negative += (s < 0);
106			}
107			}
108			}
109		✗	return (nb_positive == 0) || (nb_negative == 0);
110			}
111
112		✗	bool mesh_is_convex_2d(const Mesh& M) {
113			vector<index_t> nxt(M.vertices.nb(), NO_INDEX);
114		✗	for(index_t e: M.edges) {
115			index_t v1 = M.edges.vertex(e,0);
116			index_t v2 = M.edges.vertex(e,1);
117			geo_debug_assert(nxt[v1] == NO_INDEX);
118		✗	nxt[v1] = v2;
119			}
120
121			index_t nb_positive = 0;
122			index_t nb_negative = 0;
123		✗	for(index_t v1: M.vertices) {
124		✗	index_t v2 = nxt[v1];
125		✗	index_t v3 = nxt[v2];
126		✗	Sign s = PCK::orient_2d(
127			M.vertices.point_ptr(v1),
128			M.vertices.point_ptr(v2),
129			M.vertices.point_ptr(v3)
130			);
131		✗	nb_positive += (s > 0);
132		✗	nb_negative += (s < 0);
133			}
134		✗	return (nb_positive == 0) || (nb_negative == 0);
135			}
136
137		✗	void compute_minkowski_sum_convex_convex_3d(
138			Mesh& result, const Mesh& op1, const Mesh& op2
139			) {
140		✗	result.clear();
141		✗	result.vertices.set_dimension(3);
142
143		✗	result.vertices.create_vertices(op1.vertices.nb()*op2.vertices.nb());
144		✗	for(index_t v1: op1.vertices) {
145		✗	for(index_t v2: op2.vertices) {
146		✗	result.vertices.point(v1*op2.vertices.nb()+v2) =
147			op1.vertices.point(v1) + op2.vertices.point(v2) ;
148			}
149			}
150		✗	compute_convex_hull_3d(result);
151		✗	MeshSurfaceIntersection I(result);
152			// TODO: make it work without binding original_facet_id
153			Attribute<index_t> original_facet_id(
154		✗	result.facets.attributes(), "original_facet_id"
155			);
156		✗	for(index_t f: result.facets) {
157		✗	original_facet_id[f] = f;
158			}
159		✗	I.simplify_coplanar_facets();
160		✗	}
161
162			/**
163			* \brief Reference: Hichem Barki Ph.D. thesis, Algo 2 P. 90
164			* A is possibly non-convex
165			* B is convex
166			*/
167			class Minkovski {
168			public:
169		✗	Minkovski(
170			const Mesh& A, const Mesh& B, Mesh& result
171		✗	) : A_(A), B_(B), result_(result) {
172
173		✗	Av_f_.assign(A_.facets.nb(), NO_INDEX);
174		✗	Af_to_Bvcontrib_.resize(A_.facets.nb());
175
176		✗	A_has_borders_ = false;
177			// Store one facet incident to each vertex for vertex_is_elevated()
178		✗	for(index_t f: A_.facets) {
179		✗	for(index_t v: A_.facets.vertices(f)) {
180		✗	Av_f_[v] = f;
181			}
182		✗	for(index_t g: A_.facets.adjacent(f)) {
183		✗	A_has_borders_ = (A_has_borders_ || g == NO_INDEX);
184			}
185			}
186		✗	}
187
188		✗	void compute() {
189		✗	result_.clear();
190		✗	result_.vertices.set_dimension(3);
191
192			// Translated facets
193			vector<index_t> afv;
194		✗	for(index_t af: A_.facets) {
195		✗	afv.resize(0);
196		✗	for(index_t v: A_.facets.vertices(af)) {
197			afv.push_back(v);
198			}
199
200		✗	vec3 ap1 = A_.facets.point(af, 0);
201		✗	vec3 ap2 = A_.facets.point(af, 1);
202		✗	vec3 ap3 = A_.facets.point(af, 2);
203
204			vec3 N = cross(ap2-ap1,ap3-ap1);
205		✗	if(length(N) < 1e-6) {
206		✗	continue;
207			}
208
209			vector<index_t>& af_bvcontrib = Af_to_Bvcontrib_[af];
210		✗	for(index_t bv: B_.vertices) {
211		✗	if(af_bvcontrib.size() == 0) {
212		✗	af_bvcontrib.push_back(bv);
213			} else {
214		✗	vec3 bp = B_.vertices.point(bv);
215		✗	vec3 best_bp_so_far = B_.vertices.point(af_bvcontrib[0]);
216		✗	Sign s = N_dot_compare(
217			ap1, ap2, ap3, bp, best_bp_so_far
218			);
219		✗	if(s >= 0) {
220		✗	if(s > 0) {
221		✗	af_bvcontrib.resize(0);
222			}
223		✗	af_bvcontrib.push_back(bv);
224			}
225			}
226			}
227		✗	minkowski_2d(afv, af_bvcontrib);
228			}
229
230			// Corner facets
231		✗	for(index_t bf: B_.facets) {
232		✗	vec3 bp1 = B_.facets.point(bf, 0);
233		✗	vec3 bp2 = B_.facets.point(bf, 1);
234		✗	vec3 bp3 = B_.facets.point(bf, 2);
235			vec3 bfN = cross(bp2-bp1,bp3-bp1);
236		✗	if(length(bfN) < 1e-6) {
237			continue;
238			}
239		✗	for(index_t av: A_.vertices) {
240		✗	if(A_vertex_is_elevated(av, bp1, bp2, bp3)) {
241		✗	vec3 T = A_.vertices.point(av);
242		✗	index_t deg = B_.facets.nb_vertices(bf);
243		✗	index_t first_v = result_.vertices.create_vertices(deg);
244		✗	index_t new_f = result_.facets.create_polygon(deg);
245		✗	for(index_t lv=0; lv<deg; ++lv) {
246		✗	result_.facets.set_vertex(new_f, lv, first_v+lv);
247		✗	result_.vertices.point(first_v+lv) =
248		✗	B_.facets.point(bf,lv) + T;
249			}
250			}
251			}
252			// TODO: edges of op1 elevated w.r.t. V
253			}
254
255			// Edge facets
256		✗	for_each_edge(A_, [&](const Edge& aE) {
257		✗	if(!A_facets_have_distinct_contributing_vertices(aE.f, aE.g)) {
258			return;
259			}
260		✗	if(edge_convexity(A_, aE.f, aE.le) == POSITIVE) {
261			return;
262			}
263		✗	for_each_edge(B_, [&](const Edge& bE) {
264			// B is almost always a tessellated sphere with
265			// quads made of two triangles. Ignore the diagonal
266			// edges of those quads.
267			//if(Geom::angle(bE.Nf, bE.Ng) < 0.01 * M_PI / 180.0) {
268			// return;
269			//}
270
271		✗	Sign s1 = geo_sgn(dot(aE.dir,bE.Nf));
272			Sign s2 = geo_sgn(dot(aE.dir,bE.Ng));
273
274		✗	if(s1 == s2) {
275			return;
276			}
277
278			vec3 abN = cross(aE.dir,bE.dir);
279			if(
280		✗	double(s1)*dot(cross(aE.Nf,abN),aE.dir) < 0 &&
281		✗	double(s1)*dot(cross(abN,aE.Ng),aE.dir) < 0
282			) {
283		✗	create_quad(
284			aE.p1+bE.p1, aE.p1+bE.p2, aE.p2+bE.p2, aE.p2+bE.p1
285			);
286			}
287			});
288			});
289
290		✗	tessellate_facets(result_, 3);
291		✗	}
292
293
294			protected:
295
296			struct Edge {
297		✗	Edge(
298			const Mesh& mesh_in, index_t f_in, index_t le_in
299		✗	) : mesh(mesh_in), f(f_in), le(le_in) {
300		✗	g = mesh.facets.adjacent(f,le);
301			if(!valid()) {
302			return;
303			}
304		✗	v1 = mesh.facets.vertex(f,le);
305		✗	v2 = mesh.facets.vertex(f,(le+1)%mesh.facets.nb_vertices(f));
306		✗	p1 = mesh.vertices.point(v1);
307		✗	p2 = mesh.vertices.point(v2);
308		✗	dir = p2-p1;
309		✗	Nf = Geom::mesh_facet_normal(mesh,f);
310		✗	Ng = Geom::mesh_facet_normal(mesh,g);
311			}
312			bool valid() {
313		✗	return g != NO_INDEX && f > g;
314			}
315			const Mesh& mesh;
316			index_t f;
317			index_t le;
318			index_t g;
319			index_t v1;
320			index_t v2;
321			vec3 p1;
322			vec3 p2;
323			vec3 dir;
324			vec3 Nf;
325			vec3 Ng;
326			};
327
328		✗	void for_each_edge(
329			const Mesh& M, std::function<void(const Edge& E)> doit
330			) {
331		✗	for(index_t f: M.facets) {
332		✗	for(index_t le=0; le<M.facets.nb_vertices(f); ++le) {
333		✗	Edge E(M,f,le);
334			if(E.valid()) {
335			doit(E);
336			}
337			}
338			}
339		✗	}
340
341		✗	void create_quad(vec3 p1, vec3 p2, vec3 p3, vec3 p4) {
342		✗	index_t first_v =result_.vertices.create_vertices(4);
343		✗	index_t new_f = result_.facets.create_polygon(4);
344		✗	for(index_t lv=0; lv<4; ++lv) {
345		✗	result_.facets.set_vertex(new_f, lv, first_v+lv);
346			}
347		✗	result_.vertices.point(first_v ) = p1;
348		✗	result_.vertices.point(first_v+1) = p2;
349		✗	result_.vertices.point(first_v+2) = p3;
350		✗	result_.vertices.point(first_v+3) = p4;
351		✗	}
352
353		✗	bool A_facets_have_distinct_contributing_vertices(
354			index_t Af1, index_t Af2
355			) const {
356			const vector<index_t>& c1 = Af_to_Bvcontrib_[Af1];
357			const vector<index_t>& c2 = Af_to_Bvcontrib_[Af2];
358		✗	for(index_t v1: c1) {
359		✗	if(std::find(c2.begin(), c2.end(), v1) == c2.end()) {
360			return true;
361			}
362			}
363			return false;
364			}
365
366		✗	bool A_vertex_is_elevated(
367			index_t av, vec3 bp1, vec3 bp2, vec3 bp3
368			) const {
369		✗	vec3 ap1 = A_.vertices.point(av);
370		✗	index_t first_f = Av_f_[av];
371			index_t f = first_f;
372		✗	index_t lv = A_.facets.find_vertex(f,av);
373			do {
374		✗	index_t N = A_.facets.nb_vertices(f);
375		✗	index_t v2 = A_.facets.vertex(f, (lv + 1) % N);
376		✗	vec3 ap2 = A_.vertices.point(v2);
377		✗	Sign s = N_dot_compare(bp1,bp2,bp3,ap1,ap2);
378		✗	if(s < 0) {
379		✗	return false;
380			}
381		✗	f = A_.facets.adjacent(f,lv);
382		✗	lv = A_.facets.find_vertex(f,av);
383		✗	} while(f != first_f);
384			return true;
385			}
386
387		✗	void minkowski_2d(
388			const vector<index_t>& av, const vector<index_t>& bv
389			) {
390			geo_debug_assert(av.size() >= 3);
391			geo_debug_assert(bv.size() > 0);
392		✗	if(bv.size() == 1) {
393			// simple translation
394		✗	vec3 T = B_.vertices.point(bv[0]);
395			index_t N = av.size();
396		✗	index_t first_v = result_.vertices.create_vertices(N);
397		✗	index_t new_f = result_.facets.create_polygon(N);
398		✗	for(index_t lv=0; lv<N; ++lv) {
399		✗	result_.facets.set_vertex(new_f, lv, first_v+lv);
400		✗	result_.vertices.point(first_v+lv) =
401		✗	A_.vertices.point(av[lv]) + T;
402			}
403			} else {
404			// Super ugly: use Delaunay in 2D orthogonal plane to
405			// compute convex hull.
406		✗	vec3 p1 = A_.vertices.point(av[0]);
407		✗	vec3 p2 = A_.vertices.point(av[1]);
408		✗	vec3 p3 = A_.vertices.point(av[2]);
409		✗	vec3 NN = normalize(cross(p2-p1,p3-p1));
410		✗	vec3 U = normalize(Geom::perpendicular(NN));
411			vec3 V = cross(NN,U);
412			vector<vec2> uv;
413			vector<vec3> xyz;
414		✗	uv.reserve(av.size()*bv.size());
415		✗	xyz.reserve(av.size()*bv.size());
416		✗	for(index_t curav: av) {
417		✗	vec3 ap = A_.vertices.point(curav);
418		✗	for(index_t curbv: bv) {
419		✗	vec3 bp = B_.vertices.point(curbv);
420			vec3 abp = ap+bp-p1;
421		✗	uv.emplace_back(dot(abp,U), dot(abp,V));
422		✗	xyz.emplace_back(ap+bp);
423			}
424			}
425			Delaunay_var delaunay = Delaunay::create(
426		✗	coord_index_t(2), "BDEL2d"
427		✗	);
428		✗	delaunay->set_keeps_infinite(true);
429		✗	delaunay->set_vertices(uv.size(), uv[0].data());
430			vector<index_t> nxt(uv.size(), NO_INDEX);
431			index_t first = NO_INDEX;
432		✗	for(index_t t=delaunay->nb_finite_cells();
433		✗	t<delaunay->nb_cells(); ++t) {
434			index_t v1= NO_INDEX, v2=NO_INDEX;
435		✗	for(index_t lv=0; lv<3; ++lv) {
436		✗	if(delaunay->cell_vertex(t,lv) == NO_INDEX) {
437		✗	v1 = delaunay->cell_vertex(t,(lv+1)%3);
438		✗	v2 = delaunay->cell_vertex(t,(lv+2)%3);
439			}
440			}
441			first = v2;
442		✗	nxt[v2] = v1;
443			}
444			vector<vec3> f_xyz;
445			index_t v = first;
446			do {
447			f_xyz.push_back(xyz[v]);
448		✗	v = nxt[v];
449		✗	} while(v != first);
450
451		✗	if(f_xyz.size() < 3) {
452			return;
453			}
454
455			index_t N = f_xyz.size();
456		✗	index_t first_v = result_.vertices.create_vertices(N);
457		✗	index_t new_f = result_.facets.create_polygon(N);
458		✗	for(index_t lv=0; lv<N; ++lv) {
459		✗	result_.facets.set_vertex(new_f, lv, first_v+lv);
460		✗	result_.vertices.point(first_v+lv) = f_xyz[lv];
461			}
462			}
463			}
464
465			/**
466			* \brief Compares the dot product between the normal to a triangle
467			* and two vectors
468			* \param[in] p1 , p2 , p3 the three vertices of the triangle, that
469			* define the normal vector N = cross(p2-p1,p3-p1)
470			* \param[in] q1 , q2 the two points to be compared relative to N
471			* \retval POSITIVE if dot(N,q2) > dot(N,q1)
472			* \retval ZERO if dot(N,q2) = dot(N,q1)
473			* \retval NEGATIVE if dot(N,q2) < dot(N,q1)
474			*/
475		✗	Sign N_dot_compare(
476			const vec3& p1, const vec3& p2, const vec3& p3,
477			const vec3& q1, const vec3& q2
478			) const {
479			// TODO: new specialized predicate
480		✗	return PCK::det_3d(p3-p1, p2-p1, q2-q1);
481			}
482
483			private:
484			const Mesh& A_;
485			const Mesh& B_;
486			bool A_has_borders_;
487			vector<index_t> Av_f_;
488			vector<vector<index_t>> Af_to_Bvcontrib_;
489			Mesh& result_;
490			};
491
492		✗	void compute_minkowski_sum_non_convex_convex_3d(
493			Mesh& result, const Mesh& A, const Mesh& B
494			) {
495			// TODO: why needed ?
496		✗	reorient_connected_components(const_cast<Mesh&>(A));
497		✗	reorient_connected_components(const_cast<Mesh&>(B));
498		✗	Minkovski mink(A,B,result);
499		✗	mink.compute();
500		✗	}
501
502		✗	void compute_minkowski_sum_convex_convex_2d(
503			Mesh& result, const Mesh& op1, const Mesh& op2
504			) {
505		✗	result.clear();
506		✗	result.vertices.set_dimension(2);
507		✗	result.vertices.create_vertices(op1.vertices.nb()*op2.vertices.nb());
508		✗	for(index_t v1: op1.vertices) {
509		✗	for(index_t v2: op2.vertices) {
510		✗	result.vertices.point<2>(v1*op2.vertices.nb()+v2) =
511			op1.vertices.point<2>(v1) + op2.vertices.point<2>(v2) ;
512			}
513			}
514		✗	compute_convex_hull_2d(result);
515		✗	}
516
517
518			}
519
520			namespace GEO {
521
522		✗	void compute_minkowski_sum_3d(
523			Mesh& result, const Mesh& op1, const Mesh& op2
524			) {
525			const double angle_tol = 0.01 * M_PI / 180.0;
526		✗	bool op1_is_convex = mesh_is_convex_3d(op1, angle_tol);
527		✗	bool op2_is_convex = mesh_is_convex_3d(op2, angle_tol);
528
529		✗	if(op1_is_convex && op2_is_convex) {
530		✗	compute_minkowski_sum_convex_convex_3d(result, op1, op2);
531		✗	return;
532			}
533
534		✗	if(!op1_is_convex && op2_is_convex) {
535			// op2 may be not exactly convex -> convexify it
536		✗	Mesh op2_convex;
537		✗	op2_convex.copy(op2);
538		✗	compute_convex_hull_3d(op2_convex);
539		✗	compute_minkowski_sum_non_convex_convex_3d(result, op1, op2_convex);
540			return;
541		✗	}
542
543		✗	if(op1_is_convex && !op2_is_convex) {
544			// op1 may be not exactly convex -> convexify it
545		✗	Mesh op1_convex;
546		✗	op1_convex.copy(op2);
547		✗	compute_convex_hull_3d(op1_convex);
548		✗	compute_minkowski_sum_non_convex_convex_3d(result, op2, op1_convex);
549			return;
550		✗	}
551
552			throw(
553			std::logic_error(
554			"compute_Minkowski_sum_3d: "
555			"non-convex case not implemented yet"
556			)
557		✗	);
558			}
559
560		✗	void compute_minkowski_sum_2d(
561			Mesh& result, const Mesh& op1, const Mesh& op2
562			) {
563		✗	if(mesh_is_convex_2d(op1) && mesh_is_convex_2d(op2)) {
564		✗	compute_minkowski_sum_convex_convex_2d(result, op1, op2);
565		✗	return;
566			}
567
568			throw(
569			std::logic_error(
570			"compute_Minkowski_sum_2d: "
571			"non-convex case not implemented yet"
572			)
573		✗	);
574			}
575
576			}
577