Fix tests triangle3d:

- triangle3d::isPointInsideFast now using some epsilon to catch all points on the borders.
- triangle3d::getIntersectionOfPlaneWithLine calculates now with higher precision for more exact results.
- triangle3d::isOnSameSide (used by isPointInside) calculates now with higher precision and uses some epsilon to make it work with larger integers and less floating point troubles. Slightly slower now.


git-svn-id: svn://svn.code.sf.net/p/irrlicht/code/trunk@4183 dfc29bdd-3216-0410-991c-e03cc46cb475

changes.txt

 Changes in 1.8 (??.??.2011)
 
+ - triangle3d::isPointInsideFast now using some epsilon to catch all points on the borders.
+
+ - triangle3d::getIntersectionOfPlaneWithLine calculates now with higher precision for more exact results.
+
+ - triangle3d::isOnSameSide (used by isPointInside) calculates now with higher precision and uses some epsilon to make it work with larger integers and less floating point troubles. Slightly slower now.
+ 
  - new function equalsByUlp to test for spacing between floating point numbers.
 
  - speedup for collada writing.

include/triangle3d.h

@@ -82,22 +82,25 @@ namespace core
 		}
 
 		//! Check if a point is inside the triangle (border-points count also as inside)
-		/** NOTE: When working with T='int' you should prefer isPointInsideFast, as 
-		isPointInside will run into number-overflows already with coordinates in the 3-digit-range.
+		/*
 		\param p Point to test. Assumes that this point is already
 		on the plane of the triangle.
 		\return True if the point is inside the triangle, otherwise false. */
 		bool isPointInside(const vector3d<T>& p) const
 		{
-			return (isOnSameSide(p, pointA, pointB, pointC) &&
- 				isOnSameSide(p, pointB, pointA, pointC) &&
- 				isOnSameSide(p, pointC, pointA, pointB));
+			vector3d<f64> af64((f64)pointA.X, (f64)pointA.Y, (f64)pointA.Z);
+			vector3d<f64> bf64((f64)pointB.X, (f64)pointB.Y, (f64)pointB.Z);
+			vector3d<f64> cf64((f64)pointC.X, (f64)pointC.Y, (f64)pointC.Z);
+			vector3d<f64> pf64((f64)p.X, (f64)p.Y, (f64)p.Z);
+			return (isOnSameSide(pf64, af64, bf64, cf64) &&
+ 				isOnSameSide(pf64, bf64, af64, cf64) &&
+ 				isOnSameSide(pf64, cf64, af64, bf64));
 		}
 
 		//! Check if a point is inside the triangle (border-points count also as inside)
-		/** This method uses a barycentric coordinate system. 
-		It is faster than isPointInside but is more susceptible to floating point rounding 
-		errors. This will especially be noticable when the FPU is in single precision mode 
+		/** This method uses a barycentric coordinate system.
+		It is faster than isPointInside but is more susceptible to floating point rounding
+		errors. This will especially be noticable when the FPU is in single precision mode
 		(which is for example set on default by Direct3D).
 		\param p Point to test. Assumes that this point is already
 		on the plane of the triangle.
@@ -107,21 +110,21 @@ namespace core
 			const vector3d<T> a = pointC - pointA;
 			const vector3d<T> b = pointB - pointA;
 			const vector3d<T> c = p - pointA;
-			
+
 			const f64 dotAA = a.dotProduct( a);
 			const f64 dotAB = a.dotProduct( b);
 			const f64 dotAC = a.dotProduct( c);
 			const f64 dotBB = b.dotProduct( b);
 			const f64 dotBC = b.dotProduct( c);
-			 
+
 			// get coordinates in barycentric coordinate system
-			const f64 invDenom =  1/(dotAA * dotBB - dotAB * dotAB); 
+			const f64 invDenom =  1/(dotAA * dotBB - dotAB * dotAB);
 			const f64 u = (dotBB * dotAC - dotAB * dotBC) * invDenom;
 			const f64 v = (dotAA * dotBC - dotAB * dotAC ) * invDenom;
-		 
+
 			// We count border-points as inside to keep downward compatibility.
-			// That's why we use >= and <= instead of > and < as more commonly seen on the web.
-			return (u >= 0) && (v >= 0) && (u + v <= 1);
+			// Rounding-error also needed for some test-cases.
+			return (u > -ROUNDING_ERROR_f32) && (v >= 0) && (u + v < 1+ROUNDING_ERROR_f32);
 
 		}
 
@@ -166,15 +169,27 @@ namespace core
 		bool getIntersectionOfPlaneWithLine(const vector3d<T>& linePoint,
 			const vector3d<T>& lineVect, vector3d<T>& outIntersection) const
 		{
-			const vector3d<T> normal = getNormal().normalize();
-			T t2;
-
-			if ( core::iszero ( t2 = normal.dotProduct(lineVect) ) )
+			// Work with f64 to get more precise results (makes enough difference to be worth the casts).
+			const vector3d<f64> linePointf64(linePoint.X, linePoint.Y, linePoint.Z);
+			const vector3d<f64> lineVectf64(lineVect.X, lineVect.Y, lineVect.Z);
+			vector3d<f64> outIntersectionf64;
+
+			core::triangle3d<irr::f64> trianglef64(vector3d<f64>((f64)pointA.X, (f64)pointA.Y, (f64)pointA.Z)
+										,vector3d<f64>((f64)pointB.X, (f64)pointB.Y, (f64)pointB.Z)
+										, vector3d<f64>((f64)pointC.X, (f64)pointC.Y, (f64)pointC.Z));
+			const vector3d<irr::f64> normalf64 = trianglef64.getNormal().normalize();
+			f64 t2;
+
+			if ( core::iszero ( t2 = normalf64.dotProduct(lineVectf64) ) )
 				return false;
 
-			T d = pointA.dotProduct(normal);
-			T t = -(normal.dotProduct(linePoint) - d) / t2;
-			outIntersection = linePoint + (lineVect * t);
+			f64 d = trianglef64.pointA.dotProduct(normalf64);
+			f64 t = -(normalf64.dotProduct(linePointf64) - d) / t2;
+			outIntersectionf64 = linePointf64 + (lineVectf64 * t);
+
+			outIntersection.X = (T)outIntersectionf64.X;
+			outIntersection.Y = (T)outIntersectionf64.Y;
+			outIntersection.Z = (T)outIntersectionf64.Z;
 			return true;
 		}
 
@@ -226,13 +241,27 @@ namespace core
 		vector3d<T> pointC;
 
 	private:
-		bool isOnSameSide(const vector3d<T>& p1, const vector3d<T>& p2,
-			const vector3d<T>& a, const vector3d<T>& b) const
+		// Using f64 instead of <T> to avoid integer overflows when T=int (maybe also less floating point troubles).
+		bool isOnSameSide(const vector3d<f64>& p1, const vector3d<f64>& p2,
+			const vector3d<f64>& a, const vector3d<f64>& b) const
 		{
-			vector3d<T> bminusa = b - a;
-			vector3d<T> cp1 = bminusa.crossProduct(p1 - a);
-			vector3d<T> cp2 = bminusa.crossProduct(p2 - a);
-			return (cp1.dotProduct(cp2) >= 0.0f);
+			vector3d<f64> bminusa = b - a;
+			vector3d<f64> cp1 = bminusa.crossProduct(p1 - a);
+			vector3d<f64> cp2 = bminusa.crossProduct(p2 - a);
+			f64 res = cp1.dotProduct(cp2);
+			if ( res < 0 )
+			{
+				// This catches some floating point troubles.
+				// Unfortunately slightly expensive and we don't really know the best epsilon for iszero.
+				vector3d<f64> cp1 = bminusa.normalize().crossProduct((p1 - a).normalize());
+				if ( 	core::iszero(cp1.X, (f64)ROUNDING_ERROR_f32)
+					&& 	core::iszero(cp1.Y, (f64)ROUNDING_ERROR_f32)
+					&& 	core::iszero(cp1.Z, (f64)ROUNDING_ERROR_f32) )
+				{
+					res = 0.f;
+				}
+			}
+			return (res >= 0.0f);
 		}
 	};
 

tests/tests-last-passed-at.txt

-Tests finished. 1 test of 1 passed.
-Compiled as DEBUG
-Test suite pass at GMT Tue Jun 05 19:06:55 2012
-
+Tests finished. 1 test of 1 passed.
+Compiled as DEBUG
+Test suite pass at GMT Mon Jun 11 15:15:12 2012
+

tests/triangle3d.cpp

@@ -7,22 +7,13 @@
 using namespace irr;
 using namespace core;
 
-template<class T>
-static bool isOnSameSide(const vector3d<T>& p1, const vector3d<T>& p2,
-	const vector3d<T>& a, const vector3d<T>& b)
-{
-	vector3d<T> bminusa = b - a;
-	vector3d<T> cp1 = bminusa.crossProduct(p1 - a);
-	vector3d<T> cp2 = bminusa.crossProduct(p2 - a);
-	return (cp1.dotProduct(cp2)+core::ROUNDING_ERROR_f64 >= 0.0f);
-}
-
 template<class T>
 static bool testGetIntersectionWithLine(core::triangle3d<T>& triangle, const core::line3d<T>& ray)
 {
 	bool allExpected=true;
 	const vector3d<T> linevect = ray.getVector().normalize();
 	vector3d<T> intersection;
+	// When we just raise Y parallel to the ray then either all should fail or all succeed (the latter in our case).
 	for (u32 i=0; i<100; ++i)
 	{
 		if (!triangle.getIntersectionOfPlaneWithLine(ray.start, linevect, intersection))
@@ -30,6 +21,7 @@ static bool testGetIntersectionWithLine(core::triangle3d<T>& triangle, const cor
 			allExpected=false;
 			logTestString("triangle3d plane test %d failed\n", i);
 		}
+		//logTestString("intersection: %f %f %f\n", intersection.X, intersection.Y, intersection.Z);
 		if (!triangle.isPointInsideFast(intersection))
 		{
 			allExpected=false;
@@ -39,12 +31,6 @@ static bool testGetIntersectionWithLine(core::triangle3d<T>& triangle, const cor
 		{
 			allExpected=false;
 			logTestString("triangle3d point test %d failed\n", i);
-		        if (!isOnSameSide(intersection, triangle.pointA, triangle.pointB, triangle.pointC))
-				logTestString("triangle3d side1 test %d failed\n", i);
-			if (!isOnSameSide(intersection, triangle.pointB, triangle.pointA, triangle.pointC))
-				logTestString("triangle3d side2 test %d failed\n", i);
-			if (!isOnSameSide(intersection, triangle.pointC, triangle.pointA, triangle.pointB))
-				logTestString("triangle3d side3 test %d failed\n", i);
 		}
 
 		if (!triangle.getIntersectionWithLine(ray.start, linevect, intersection))
@@ -276,41 +262,33 @@ bool testTriangle3d(void)
 {
 	bool allExpected = true;
 
-	/* TODO: disabled for now. I (aka CuteAlien) have by now an example which allows debugging
-	   that problem easier and also found some workaround (which needs an interface change
-	   and a behaviour change and won't get into 1.7 therefore).
 	logTestString("Test getIntersectionWithLine with f32\n");
 	{
 		triangle3df triangle(
-				vector3df(11300.000000f, 129.411758f, 200.000000f),
-				vector3df(11200.000000f, 94.117645f, 300.000000f),
-				vector3df(11300.000000f, 129.411758f, 300.000000f));
+				vector3df(11300.f, 129.411758f, 200.f),
+				vector3df(11200.f, 94.117645f, 300.f),
+				vector3df(11300.f, 129.411758f, 300.f));
 		line3df ray;
-		ray.start = vector3df(11250.000000f, 329.000000f, 250.000000f);
-		ray.end = vector3df(11250.000000, -1000.000000, 250.000000);
+		ray.start = vector3df(11250.f, 329.f, 250.f);
+		ray.end = vector3df(11250.f, -1000.f, 250.f);
 		allExpected &= testGetIntersectionWithLine(triangle, ray);
 	}
 	logTestString("Test getIntersectionWithLine with f64\n");
 	{
 		triangle3d<f64> triangle(
-				vector3d<f64>(11300.000000f, 129.411758f, 200.000000f),
-				vector3d<f64>(11200.000000f, 94.117645f, 300.000000f),
-				vector3d<f64>(11300.000000f, 129.411758f, 300.000000f));
+				vector3d<f64>(11300., 129.411758, 200.),
+				vector3d<f64>(11200., 94.117645, 300.),
+				vector3d<f64>(11300., 129.411758, 300.));
 		line3d<f64> ray;
-		ray.start = vector3d<f64>(11250.000000f, 329.000000f, 250.000000f);
-		ray.end = vector3d<f64>(11250.000000, -1000.000000, 250.000000);
+		ray.start = vector3d<f64>(11250., 329., 250.);
+		ray.end = vector3d<f64>(11250., -1000., 250.);
 		allExpected &= testGetIntersectionWithLine(triangle, ray);
 	}
-	*/
 
-	/* For now we have no solution yet to fix isPointInside for large integers without
-	getting worse floating-point precision at the same time.
-	So instead isPointInsideFast got fixed and should be used for int's.
 	bool testEigen = triangle3di(vector3di(250, 0, 0), vector3di(0, 0, 500), vector3di(500, 0, 500)).isPointInside(vector3di(300,0,300));
 	if ( !testEigen )	// test from Eigen from here: http://irrlicht.sourceforge.net/forum/viewtopic.php?f=7&t=44372&p=254331#p254331
 		logTestString("Test isPointInside fails with integers\n");
 	allExpected &= testEigen;
-	*/
 
 	logTestString("Test isPointInside with f32\n");
 	{