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/geometry.h>
41			#include <geogram/numerics/matrix_util.h>
42			#include <geogram/basic/logger.h>
43
44			namespace {
45			using namespace GEO;
46
47			const index_t cot_edges[6][2] = {
48			{ 0, 1 }, { 0, 2 }, { 0, 3 }, { 1, 2 }, { 2, 3 }, { 3, 1 }
49			} ;
50
51			}
52
53			namespace GEO {
54
55		✗	CoTan3D::CoTan3D(vec3 P[4], double* anisotropy_as_xx_yy_zz_xy_yz_xz ) {
56		✗	tetvol = (1. / 6.) * dot(P[3] - P[0], cross(P[1] - P[0], P[2] - P[0]));
57		✗	Matrix<6, double> M;
58		✗	M.load_zero();
59		✗	double isotrope_objective[6] = { 1, 1, 1, 0, 0, 0 };
60		✗	double *RHS = anisotropy_as_xx_yy_zz_xy_yz_xz;
61		✗	if (RHS == nullptr)
62		✗	RHS = isotrope_objective;
63		✗	FOR(e, 6){
64		✗	vec3 geom = P[cot_edges[e][1]] - P[cot_edges[e][0]];
65		✗	M(0, e) = geom.x*geom.x;
66		✗	M(1, e) = geom.y*geom.y;
67		✗	M(2, e) = geom.z*geom.z;
68		✗	M(3, e) = 2. * geom.x*geom.y;
69		✗	M(4, e) = 2. * geom.y*geom.z;
70		✗	M(5, e) = 2. * geom.z*geom.x;
71			}
72		✗	Matrix<6, double> inv;
73		✗	bool invertible = M.compute_inverse(inv);
74		✗	if (!invertible) GEO::Logger::out("HexDom") << "Solve did not work" << std::endl;
75		✗	mult(inv, RHS, w);
76		✗	}
77
78
79		✗	index_t CoTan3D::org(index_t e) { return cot_edges[e][0]; }
80		✗	index_t CoTan3D::dest(index_t e) { return cot_edges[e][1]; }
81		✗	double CoTan3D::coeff(index_t e) { return w[e]; }
82
83		✗	void CoTan3D::check_for_grad(vec3 P[4], vec3 grad){
84			double v[4];
85		✗	FOR(i, 4) v[i] = dot(grad, P[i]);
86		✗	double sum = 0;
87		✗	FOR(e, 6) sum += (v[dest(e)] - v[org(e)])*(v[dest(e)] - v[org(e)])* coeff(e);
88		✗	std::cerr << "NRJ = " << grad.length2()
89		✗	<< "\t\twith cot = " << sum
90		✗	<< "\t\tratio = " << grad.length2() / sum
91		✗	<< std::endl;
92		✗	}
93
94			/*******************************************************************************/
95
96		✗	TrglGradient::TrglGradient(const vec3& p0, const vec3& p1, const vec3& p2) {
97		✗	initialize(p0, p1, p2);
98		✗	}
99
100		✗	TrglGradient::TrglGradient() {
101		✗	}
102
103		✗	void TrglGradient::initialize(const vec3& p0, const vec3& p1, const vec3& p2) {
104
105			// Computing TX[] and TY[],
106			// i.e. the coefficients such that:
107			// | df/dX = TX[0].f(v0) + TX[1].f(v1) + TX[2].f(v2)
108			// | df/dY = TY[0].f(v0) + TY[1].f(v1) + TY[2].f(v2)
109			//
110			// (in other words, these coefficient give the gradient of a property
111			// interpolated in the triangle T = (a0, a1, a2).
112
113			// The equations of the coefficients can be simplified,
114			// the general equations are given as code in comments
115			// marked by the tag [simplified]
116
117		✗	vertex_[0] = p0;
118		✗	vertex_[1] = p1;
119		✗	vertex_[2] = p2;
120
121			// Step1: find an orthonormal basis (X, Y) for the triangle.
122			//----------------------------------------------------------
123
124		✗	vec3 origin;
125		✗	vec3 X;
126		✗	vec3 Y;
127		✗	vec3 Z;
128
129		✗	basis(origin, X, Y, Z); // Rem: origin = v0.
130
131			// Step2: compute the coordinates of the three vertices of the
132			// triangle in this basis.
133			//------------------------------------------------------------
134
135			// [simplified] double x0 = 0.0 ;
136			// [simplified] double y0 = 0.0 ;
137		✗	vec3 V1 = p1 - p0;
138		✗	double x1 = sqrt(dot(V1, V1));
139
140			// [simplified] double y1 = 0.0 ;
141		✗	vec3 V2 = p2 - p0;
142		✗	double x2 = dot(V2, X);
143		✗	double y2 = dot(V2, Y);
144
145
146			// Step3: compute the six coefficients TXi and TYi allowing to
147			// compute the two components of the gradient of a function in
148			// the basis (I,X,Y).
149			//------------------------------------------------------------
150
151			// [simplified]
152			// double d = (x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0) ;
153
154		✗	double d = x1 * y2;
155
156		✗	if (fabs(d) < 1e-10) {
157		✗	Logger::warn("TrglGrad")
158		✗	<< "Attempted gradient computation on a flat triangle"
159		✗	<< std::endl;
160		✗	d = 1.0;
161		✗	is_flat_ = true;
162			} else {
163		✗	is_flat_ = false;
164			}
165
166			// [simplified] TX_[0] = (y1 - y2) / d ;
167			// [simplified] TX_[1] = (y2 - y0) / d ;
168			// [simplified] TX_[2] = (y0 - y1) / d ;
169
170		✗	double xx1 = (x1 == 0.0) ? 1.0 : x1;
171
172		✗	TX_[0] = -1.0 / xx1;
173		✗	TX_[1] = 1.0 / xx1;
174		✗	TX_[2] = 0.0;
175
176			// [simplified] TY_[0] = (x2 - x1) / d ;
177			// [simplified] TY_[1] = (x0 - x2) / d ;
178			// [simplified] TY_[2] = (x1 - x0) / d ;
179
180		✗	double yy2 = (y2 == 0.0) ? 1.0 : y2;
181
182		✗	TY_[0] = x2 / d - 1.0 / yy2;
183		✗	TY_[1] = -x2 / d;
184		✗	TY_[2] = 1.0 / yy2;
185		✗	}
186
187		✗	void TrglGradient::basis(vec3& origin, vec3& X, vec3& Y, vec3& Z) const {
188		✗	const vec3& p0 = vertex_[0];
189		✗	const vec3& p1 = vertex_[1];
190		✗	const vec3& p2 = vertex_[2];
191
192		✗	X = normalize(p1 - p0);
193		✗	Z = normalize(cross(X, p2 - p0));
194		✗	Y = cross(Z, X);
195		✗	origin = p0;
196		✗	}
197
198		✗	vec3 TrglGradient::gradient_3d(double value0, double value1, double value2) const {
199		✗	double x = TX_[0] * value0 + TX_[1] * value1; // Note: TX_[2] = 0
200		✗	double y = TY_[0] * value0 + TY_[1] * value1 + TY_[2] * value2;
201		✗	vec3 O;
202		✗	vec3 X, Y, Z;
203		✗	basis(O, X, Y, Z);
204		✗	return x*X + y*Y;
205			}
206
207			/**********************************************************************************/
208
209			// _____ _____ ____ _____ _ _ _
210			// | __ \ / ____| /\ |___ \| __ \ | | | | | |
211			// | |__) | | / \ __) | | | | ______ ______ _ __ ___ | |_ ___ ___ _ __ | |_ ___ _ __ ___ __| |
212			// | ___/| | / /\ \ |__ <| | | | |______|______| | '_ \ / _ \| __| / __/ _ \ '_ \| __/ _ \ '__/ _ \/ _` |
213			// | | | |____ / ____ \ ___) | |__| | | | | | (_) | |_ | (_| __/ | | | || __/ | | __/ (_| |
214			// |_| \_____/_/ \_\ |____/|_____/ |_| |_|\___/ \__| \___\___|_| |_|\__\___|_| \___|\__,_|
215			//
216			//
217
218		✗	void UncenteredPCA3D::begin_points(){
219		✗	nb_points_ = 0;
220		✗	sum_weights_ = M_[0] = M_[1] = M_[2] = M_[3] = M_[4] = M_[5] = 0;
221		✗	}
222
223		✗	void UncenteredPCA3D::point(const vec3& p, double weight){
224		✗	M_[0] += weight * p.x*p.x;
225		✗	M_[1] += weight * p.x*p.y; M_[2] += weight * p.y*p.y;
226		✗	M_[3] += weight * p.x*p.z; M_[4] += weight * p.y*p.z; M_[5] += weight * p.z*p.z;
227		✗	nb_points_++;
228		✗	sum_weights_ += weight;
229		✗	}
230
231		✗	void UncenteredPCA3D::end_points(){
232		✗	for (int i = 0; i < 6; i++) M_[i] = M_[i] / sum_weights_;
233		✗	if (M_[0] <= 0) M_[0] = 1.e-30;
234		✗	if (M_[2] <= 0) M_[2] = 1.e-30;
235		✗	if (M_[5] <= 0) M_[5] = 1.e-30;
236			double eigen_vectors[9];
237			// TODO: remove this dependance
238		✗	GEO::MatrixUtil::semi_definite_symmetric_eigen(M_, 3, eigen_vectors, eigen_value);
239
240		✗	for (int i = 0; i < 3; i++) axis[i] = normalize(vec3(eigen_vectors[3 * i + 0], eigen_vectors[3 * i + 1], eigen_vectors[3 * i + 2]));
241		✗	}
242
243
244			/**********************************************************************************/
245			/* Triangle-Triangle and Triangle-Box intersection routines by Tomas Moller, 1997 */
246			/**********************************************************************************/
247
248			// [Bruno] Modernized it a bit and replaced some macros with templates.
249
250		✗	template <class T> inline T FABS(T x) {
251		✗	return T(::fabs(double(x)));
252			}
253
254			/* if USE_EPSILON_TEST is true then we do a check:
255			* if |dv|<EPSILON then dv=0.0;
256			* else no check is done (which is less robust)
257			*/
258			#define USE_EPSILON_TEST
259
260			const double EPSILON = 0.000001;
261
262		✗	template <class T> inline void CROSS(T dest[3], const T v1[3], const T v2[3]) {
263		✗	dest[0] = v1[1] * v2[2] - v1[2] * v2[1];
264		✗	dest[1] = v1[2] * v2[0] - v1[0] * v2[2];
265		✗	dest[2] = v1[0] * v2[1] - v1[1] * v2[0];
266		✗	}
267
268		✗	template <class T> inline T DOT(T v1[3], T v2[3]) {
269		✗	return (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]);
270			}
271
272		✗	template <class T> inline void SUB(T dest[3], const T v1[3], const T v2[3]) {
273		✗	dest[0] = v1[0] - v2[0];
274		✗	dest[1] = v1[1] - v2[1];
275		✗	dest[2] = v1[2] - v2[2];
276		✗	}
277
278			/* sort so that a<=b */
279		✗	template <class T> inline void SORT(T& a, T& b) {
280		✗	if(a>b) {
281		✗	std::swap(a,b);
282			}
283		✗	}
284
285			/* this edge to edge test is based on Franlin Antonio's gem:
286			* "Faster Line Segment Intersection", in Graphics Gems III,
287			* pp. 199-202
288			*/
289			#define EDGE_EDGE_TEST(V0,U0,U1) \
290			Bx = U0[i0] - U1[i0]; \
291			By = U0[i1] - U1[i1]; \
292			Cx = V0[i0] - U0[i0]; \
293			Cy = V0[i1] - U0[i1]; \
294			f = Ay*Bx - Ax*By; \
295			d = By*Cx - Bx*Cy; \
296			if ((f>0 && d >= 0 && d <= f) || (f<0 && d <= 0 && d >= f)) \
297			{ \
298			e = Ax*Cy - Ay*Cx; \
299			if (f>0) \
300			{ \
301			if (e >= 0 && e <= f) return 1; \
302			} \
303			else \
304			{ \
305			if (e <= 0 && e >= f) return 1; \
306			} \
307			}
308
309			#define EDGE_AGAINST_TRI_EDGES(V0,V1,U0,U1,U2) \
310			{ \
311			double Ax, Ay, Bx, By, Cx, Cy, e, d, f; \
312			Ax = V1[i0] - V0[i0]; \
313			Ay = V1[i1] - V0[i1]; \
314			/* test edge U0,U1 against V0,V1 */ \
315			EDGE_EDGE_TEST(V0, U0, U1); \
316			/* test edge U1,U2 against V0,V1 */ \
317			EDGE_EDGE_TEST(V0, U1, U2); \
318			/* test edge U2,U1 against V0,V1 */ \
319			EDGE_EDGE_TEST(V0, U2, U0); \
320			}
321
322			#define POINT_IN_TRI(V0,U0,U1,U2) \
323			{ \
324			double a, b, c, d0, d1, d2; \
325			/* is T1 completly inside T2? */ \
326			/* check if V0 is inside tri(U0,U1,U2) */ \
327			a = U1[i1] - U0[i1]; \
328			b = -(U1[i0] - U0[i0]); \
329			c = -a*U0[i0] - b*U0[i1]; \
330			d0 = a*V0[i0] + b*V0[i1] + c; \
331			\
332			a = U2[i1] - U1[i1]; \
333			b = -(U2[i0] - U1[i0]); \
334			c = -a*U1[i0] - b*U1[i1]; \
335			d1 = a*V0[i0] + b*V0[i1] + c; \
336			\
337			a = U0[i1] - U2[i1]; \
338			b = -(U0[i0] - U2[i0]); \
339			c = -a*U2[i0] - b*U2[i1]; \
340			d2 = a*V0[i0] + b*V0[i1] + c; \
341			if (d0*d1>0.0) \
342			{ \
343			if (d0*d2>0.0) return 1; \
344			} \
345			}
346
347		✗	static int coplanar_tri_tri(
348			double N[3], double V0[3], double V1[3], double V2[3],
349			double U0[3], double U1[3], double U2[3]
350			) {
351			double A[3];
352			short i0, i1;
353			/* first project onto an axis-aligned plane, that maximizes the area */
354			/* of the triangles, compute indices: i0,i1. */
355		✗	A[0] = FABS(N[0]);
356		✗	A[1] = FABS(N[1]);
357		✗	A[2] = FABS(N[2]);
358		✗	if (A[0]>A[1])
359			{
360		✗	if (A[0]>A[2])
361			{
362		✗	i0 = 1; /* A[0] is greatest */
363		✗	i1 = 2;
364			}
365			else
366			{
367		✗	i0 = 0; /* A[2] is greatest */
368		✗	i1 = 1;
369			}
370			}
371			else /* A[0]<=A[1] */
372			{
373		✗	if (A[2]>A[1])
374			{
375		✗	i0 = 0; /* A[2] is greatest */
376		✗	i1 = 1;
377			}
378			else
379			{
380		✗	i0 = 0; /* A[1] is greatest */
381		✗	i1 = 2;
382			}
383			}
384
385			/* test all edges of triangle 1 against the edges of triangle 2 */
386		✗	EDGE_AGAINST_TRI_EDGES(V0, V1, U0, U1, U2);
387		✗	EDGE_AGAINST_TRI_EDGES(V1, V2, U0, U1, U2);
388		✗	EDGE_AGAINST_TRI_EDGES(V2, V0, U0, U1, U2);
389
390			/* finally, test if tri1 is totally contained in tri2 or vice versa */
391		✗	POINT_IN_TRI(V0, U0, U1, U2);
392		✗	POINT_IN_TRI(U0, V0, V1, V2);
393
394		✗	return 0;
395			}
396
397
398
399			#define NEWCOMPUTE_INTERVALS(VV0,VV1,VV2,D0,D1,D2,D0D1,D0D2,A,B,C,X0,X1) \
400			{ \
401			if (D0D1>0.0) \
402			{ \
403			/* here we know that D0D2<=0.0 */ \
404			/* that is D0, D1 are on the same side, D2 on the other or on the plane */ \
405			A = VV2; B = (VV0 - VV2)*D2; C = (VV1 - VV2)*D2; X0 = D2 - D0; X1 = D2 - D1; \
406			} \
407			else if (D0D2>0.0) \
408			{ \
409			/* here we know that d0d1<=0.0 */ \
410			A = VV1; B = (VV0 - VV1)*D1; C = (VV2 - VV1)*D1; X0 = D1 - D0; X1 = D1 - D2; \
411			} \
412			else if (D1*D2>0.0 || D0 != 0.0) \
413			{ \
414			/* here we know that d0d1<=0.0 or that D0!=0.0 */ \
415			A = VV0; B = (VV1 - VV0)*D0; C = (VV2 - VV0)*D0; X0 = D0 - D1; X1 = D0 - D2; \
416			} \
417			else if (D1 != 0.0) \
418			{ \
419			A = VV1; B = (VV0 - VV1)*D1; C = (VV2 - VV1)*D1; X0 = D1 - D0; X1 = D1 - D2; \
420			} \
421			else if (D2 != 0.0) \
422			{ \
423			A = VV2; B = (VV0 - VV2)*D2; C = (VV1 - VV2)*D2; X0 = D2 - D0; X1 = D2 - D1; \
424			} \
425			else \
426			{ \
427			/* triangles are coplanar */ \
428			return coplanar_tri_tri(N1, V0, V1, V2, U0, U1, U2); \
429			} \
430			}
431
432		✗	int NoDivTriTriIsect(
433			double V0[3], double V1[3], double V2[3],
434			double U0[3], double U1[3], double U2[3]
435			) {
436			double E1[3], E2[3];
437			double N1[3], N2[3], d1, d2;
438			double du0, du1, du2, dv0, dv1, dv2;
439			double D[3];
440			double isect1[2], isect2[2];
441			double du0du1, du0du2, dv0dv1, dv0dv2;
442			short index;
443			double vp0, vp1, vp2;
444			double up0, up1, up2;
445			double bb, cc, max;
446
447			/* compute plane equation of triangle(V0,V1,V2) */
448		✗	SUB(E1, V1, V0);
449		✗	SUB(E2, V2, V0);
450		✗	CROSS(N1, E1, E2);
451		✗	d1 = -DOT(N1, V0);
452			/* plane equation 1: N1.X+d1=0 */
453
454			/* put U0,U1,U2 into plane equation 1 to compute signed distances to the plane*/
455		✗	du0 = DOT(N1, U0) + d1;
456		✗	du1 = DOT(N1, U1) + d1;
457		✗	du2 = DOT(N1, U2) + d1;
458
459			/* coplanarity robustness check */
460			#ifdef USE_EPSILON_TEST
461		✗	if (FABS(du0)<EPSILON) du0 = 0.0;
462		✗	if (FABS(du1)<EPSILON) du1 = 0.0;
463		✗	if (FABS(du2)<EPSILON) du2 = 0.0;
464			#endif
465		✗	du0du1 = du0*du1;
466		✗	du0du2 = du0*du2;
467
468		✗	if (du0du1>0.0 && du0du2>0.0) /* same sign on all of them + not equal 0 ? */
469		✗	return 0; /* no intersection occurs */
470
471			/* compute plane of triangle (U0,U1,U2) */
472		✗	SUB(E1, U1, U0);
473		✗	SUB(E2, U2, U0);
474		✗	CROSS(N2, E1, E2);
475		✗	d2 = -DOT(N2, U0);
476			/* plane equation 2: N2.X+d2=0 */
477
478			/* put V0,V1,V2 into plane equation 2 */
479		✗	dv0 = DOT(N2, V0) + d2;
480		✗	dv1 = DOT(N2, V1) + d2;
481		✗	dv2 = DOT(N2, V2) + d2;
482
483			#ifdef USE_EPSILON_TEST
484		✗	if (FABS(dv0)<EPSILON) dv0 = 0.0;
485		✗	if (FABS(dv1)<EPSILON) dv1 = 0.0;
486		✗	if (FABS(dv2)<EPSILON) dv2 = 0.0;
487			#endif
488
489		✗	dv0dv1 = dv0*dv1;
490		✗	dv0dv2 = dv0*dv2;
491
492		✗	if (dv0dv1>0.0 && dv0dv2>0.0) /* same sign on all of them + not equal 0 ? */
493		✗	return 0; /* no intersection occurs */
494
495			/* compute direction of intersection line */
496		✗	CROSS(D, N1, N2);
497
498			/* compute and index to the largest component of D */
499		✗	max = (double)FABS(D[0]);
500		✗	index = 0;
501		✗	bb = (double)FABS(D[1]);
502		✗	cc = (double)FABS(D[2]);
503		✗	if (bb>max) {max = bb; index = 1;}
504		✗	if (cc>max) {max = cc; index = 2;}
505
506			/* this is the simplified projection onto L*/
507		✗	vp0 = V0[index];
508		✗	vp1 = V1[index];
509		✗	vp2 = V2[index];
510
511		✗	up0 = U0[index];
512		✗	up1 = U1[index];
513		✗	up2 = U2[index];
514
515			/* compute interval for triangle 1 */
516			double a, b, c, x0, x1;
517		✗	NEWCOMPUTE_INTERVALS(vp0, vp1, vp2, dv0, dv1, dv2, dv0dv1, dv0dv2, a, b, c, x0, x1);
518
519			/* compute interval for triangle 2 */
520			double d, e, f, y0, y1;
521		✗	NEWCOMPUTE_INTERVALS(up0, up1, up2, du0, du1, du2, du0du1, du0du2, d, e, f, y0, y1);
522
523			double xx, yy, xxyy, tmp;
524		✗	xx = x0*x1;
525		✗	yy = y0*y1;
526		✗	xxyy = xx*yy;
527
528		✗	tmp = a*xxyy;
529		✗	isect1[0] = tmp + b*x1*yy;
530		✗	isect1[1] = tmp + c*x0*yy;
531
532		✗	tmp = d*xxyy;
533		✗	isect2[0] = tmp + e*xx*y1;
534		✗	isect2[1] = tmp + f*xx*y0;
535
536		✗	SORT(isect1[0], isect1[1]);
537		✗	SORT(isect2[0], isect2[1]);
538
539		✗	if (isect1[1]<isect2[0] || isect2[1]<isect1[0]) return 0;
540		✗	return 1;
541			}
542
543			/*******************************************************************************/
544
545			#define X 0
546			#define Y 1
547			#define Z 2
548
549		✗	template <class T> inline void FINDMINMAX(T x0, T x1, T x2, T& min, T& max) {
550		✗	min = max = x0;
551		✗	if(x1<min) min=x1;
552		✗	if(x1>max) max=x1;
553		✗	if(x2<min) min=x2;
554		✗	if(x2>max) max=x2;
555		✗	}
556
557		✗	static int planeBoxOverlap(float normal[3], float d, float maxbox[3]) {
558			int q;
559			float vmin[3], vmax[3];
560		✗	for (q = X; q <= Z; q++)
561			{
562		✗	if (normal[q]>0.0f)
563			{
564		✗	vmin[q] = -maxbox[q];
565		✗	vmax[q] = maxbox[q];
566			}
567			else
568			{
569		✗	vmin[q] = maxbox[q];
570		✗	vmax[q] = -maxbox[q];
571			}
572			}
573		✗	if (DOT(normal, vmin) + d>0.0f) return 0;
574		✗	if (DOT(normal, vmax) + d >= 0.0f) return 1;
575
576		✗	return 0;
577			}
578
579			/*======================== X-tests ========================*/
580			#define AXISTEST_X01(a, b, fa, fb) \
581			p0 = a*v0[Y] - b*v0[Z]; \
582			p2 = a*v2[Y] - b*v2[Z]; \
583			if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
584			rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z]; \
585			if(min>rad || max<-rad) return 0;
586
587			#define AXISTEST_X2(a, b, fa, fb) \
588			p0 = a*v0[Y] - b*v0[Z]; \
589			p1 = a*v1[Y] - b*v1[Z]; \
590			if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
591			rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z]; \
592			if(min>rad || max<-rad) return 0;
593
594			/*======================== Y-tests ========================*/
595			#define AXISTEST_Y02(a, b, fa, fb) \
596			p0 = -a*v0[X] + b*v0[Z]; \
597			p2 = -a*v2[X] + b*v2[Z]; \
598			if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
599			rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z]; \
600			if(min>rad || max<-rad) return 0;
601
602			#define AXISTEST_Y1(a, b, fa, fb) \
603			p0 = -a*v0[X] + b*v0[Z]; \
604			p1 = -a*v1[X] + b*v1[Z]; \
605			if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
606			rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z]; \
607			if(min>rad || max<-rad) return 0;
608
609			/*======================== Z-tests ========================*/
610
611			#define AXISTEST_Z12(a, b, fa, fb) \
612			p1 = a*v1[X] - b*v1[Y]; \
613			p2 = a*v2[X] - b*v2[Y]; \
614			if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;} \
615			rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y]; \
616			if(min>rad || max<-rad) return 0;
617
618			#define AXISTEST_Z0(a, b, fa, fb) \
619			p0 = a*v0[X] - b*v0[Y]; \
620			p1 = a*v1[X] - b*v1[Y]; \
621			if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
622			rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y]; \
623			if(min>rad || max<-rad) return 0;
624
625
626		✗	int triBoxOverlap(
627			float boxcenter[3], float boxhalfsize[3], float triverts[3][3]
628			) {
629
630			/* use separating axis theorem to test overlap between triangle and box */
631			/* need to test for overlap in these directions: */
632			/* 1) the {x,y,z}-directions (actually, since we use the AABB of the triangle */
633			/* we do not even need to test these) */
634			/* 2) normal of the triangle */
635			/* 3) crossproduct(edge from tri, {x,y,z}-directin) */
636			/* this gives 3x3=9 more tests */
637			float v0[3], v1[3], v2[3];
638			float min, max, d, p0, p1, p2, rad, fex, fey, fez;
639			float normal[3], e0[3], e1[3], e2[3];
640
641			/* This is the fastest branch on Sun */
642			/* move everything so that the boxcenter is in (0,0,0) */
643		✗	SUB(v0, triverts[0], boxcenter);
644		✗	SUB(v1, triverts[1], boxcenter);
645		✗	SUB(v2, triverts[2], boxcenter);
646
647			/* compute triangle edges */
648		✗	SUB(e0, v1, v0); /* tri edge 0 */
649		✗	SUB(e1, v2, v1); /* tri edge 1 */
650		✗	SUB(e2, v0, v2); /* tri edge 2 */
651
652			/* Bullet 3: */
653			/* test the 9 tests first (this was faster) */
654		✗	fex = FABS(e0[X]);
655		✗	fey = FABS(e0[Y]);
656		✗	fez = FABS(e0[Z]);
657		✗	AXISTEST_X01(e0[Z], e0[Y], fez, fey);
658		✗	AXISTEST_Y02(e0[Z], e0[X], fez, fex);
659		✗	AXISTEST_Z12(e0[Y], e0[X], fey, fex);
660
661		✗	fex = FABS(e1[X]);
662		✗	fey = FABS(e1[Y]);
663		✗	fez = FABS(e1[Z]);
664		✗	AXISTEST_X01(e1[Z], e1[Y], fez, fey);
665		✗	AXISTEST_Y02(e1[Z], e1[X], fez, fex);
666		✗	AXISTEST_Z0(e1[Y], e1[X], fey, fex);
667
668		✗	fex = FABS(e2[X]);
669		✗	fey = FABS(e2[Y]);
670		✗	fez = FABS(e2[Z]);
671		✗	AXISTEST_X2(e2[Z], e2[Y], fez, fey);
672		✗	AXISTEST_Y1(e2[Z], e2[X], fez, fex);
673		✗	AXISTEST_Z12(e2[Y], e2[X], fey, fex);
674
675			/* Bullet 1: */
676			/* first test overlap in the {x,y,z}-directions */
677			/* find min, max of the triangle each direction, and test for overlap in */
678			/* that direction -- this is equivalent to testing a minimal AABB around */
679			/* the triangle against the AABB */
680
681			/* test in X-direction */
682		✗	FINDMINMAX(v0[X], v1[X], v2[X], min, max);
683		✗	if (min>boxhalfsize[X] || max<-boxhalfsize[X]) return 0;
684
685			/* test in Y-direction */
686		✗	FINDMINMAX(v0[Y], v1[Y], v2[Y], min, max);
687		✗	if (min>boxhalfsize[Y] || max<-boxhalfsize[Y]) return 0;
688
689			/* test in Z-direction */
690		✗	FINDMINMAX(v0[Z], v1[Z], v2[Z], min, max);
691		✗	if (min>boxhalfsize[Z] || max<-boxhalfsize[Z]) return 0;
692
693			/* Bullet 2: */
694			/* test if the box intersects the plane of the triangle */
695			/* compute plane equation of triangle: normal*x+d=0 */
696		✗	CROSS(normal, e0, e1);
697		✗	d = -DOT(normal, v0); /* plane eq: normal.x+d=0 */
698		✗	if (!planeBoxOverlap(normal, d, boxhalfsize)) return 0;
699
700		✗	return 1; /* box and triangle overlaps */
701			}
702
703			/*******************************************************************************/
704
705
706			}
707