GetClosestPointOnTriangle reorders vertices so that A lies on the shortest edge (#383)

This makes the result less sensitive to round off errors. Also added unit test.

.github/workflows/determinism_check.yml

@@ -1,8 +1,8 @@
 name: Determinism Check
 
 env:
-    CONVEX_VS_MESH_HASH: '0xbf4b66a448ae6f68'
-    RAGDOLL_HASH: '0x588d002b30f5b3d'
+    CONVEX_VS_MESH_HASH: '0xd956e5d4fe63b095'
+    RAGDOLL_HASH: '0xb874a5f0c92ba9c3'
 
 on:
   push:

Build/CMakeLists.txt

@@ -181,9 +181,8 @@ elseif ("${CMAKE_SYSTEM_NAME}" STREQUAL "Linux" OR "${CMAKE_SYSTEM_NAME}" STREQU
 		# Do not use -ffast-math since it cannot be turned off in a single compilation unit under clang, see Core.h
 		set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ffp-model=precise")
 
-		# On clang 14 and later we can turn off float contraction through a pragma, older versions need it off always, see Core.h
-		# Note that this does not seem to work on ARM so we turn fp contract off there always
-		if (CMAKE_CXX_COMPILER_VERSION LESS 14 OR CROSS_COMPILE_ARM OR CMAKE_OSX_ARCHITECTURES MATCHES "arm64" OR "${CMAKE_SYSTEM_PROCESSOR}" STREQUAL "aarch64")
+		# On clang 14 and later we can turn off float contraction through a pragma, older versions and deterministic versions need it off always, see Core.h
+		if (CMAKE_CXX_COMPILER_VERSION LESS 14 OR CROSS_PLATFORM_DETERMINISTIC)
 			set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ffp-contract=off")
 		endif()
 	endif()

Jolt/Geometry/ClosestPoint.h

@@ -150,22 +150,24 @@ namespace ClosestPoint
 		// Taken from: Real-Time Collision Detection - Christer Ericson (Section: Closest Point on Triangle to Point)
 		// With p = 0
 
-		// Calculate edges
-		Vec3 ab = inB - inA; 
-		Vec3 ac = inC - inA; 
-		Vec3 bc = inC - inB;
-
 		// The most accurate normal is calculated by using the two shortest edges
 		// See: https://box2d.org/posts/2014/01/troublesome-triangle/
 		// The difference in normals is most pronounced when one edge is much smaller than the others (in which case the other 2 must have roughly the same length).
 		// Therefore we can suffice by just picking the shortest from 2 edges and use that with the 3rd edge to calculate the normal.
-		// We first check which of the edges is shorter
-		UVec4 bc_shorter_than_ac = Vec4::sLess(bc.DotV4(bc), ac.DotV4(ac));
+		// We first check which of the edges is shorter and if bc is shorter than ac then we swap a with c to a is always on the shortest edge
+		UVec4 swap_ac;
+		{
+			Vec3 ac = inC - inA;
+			Vec3 bc = inC - inB;
+			swap_ac = Vec4::sLess(bc.DotV4(bc), ac.DotV4(ac));
+		}
+		Vec3 a = Vec3::sSelect(inA, inC, swap_ac);
+		Vec3 c = Vec3::sSelect(inC, inA, swap_ac);
 
-		// We calculate both normals and then select the one that had the shortest edge for our normal (this avoids branching)
-		Vec3 normal_bc = ab.Cross(bc);
-		Vec3 normal_ac = ab.Cross(ac);
-		Vec3 n = Vec3::sSelect(normal_ac, normal_bc, bc_shorter_than_ac);
+		// Calculate normal
+		Vec3 ab = inB - a;
+		Vec3 ac = c - a;
+		Vec3 n = ab.Cross(ac);
 		float n_len_sq = n.LengthSq();
 
 		// Check degenerate
@@ -203,13 +205,13 @@ namespace ClosestPoint
 		}
 
 		// Check if P in vertex region outside A 
-		Vec3 ap = -inA; 
+		Vec3 ap = -a; 
 		float d1 = ab.Dot(ap); 
 		float d2 = ac.Dot(ap); 
 		if (d1 <= 0.0f && d2 <= 0.0f)
 		{
-			outSet = 0b0001;
-			return inA; // barycentric coordinates (1,0,0)
+			outSet = swap_ac.TestAnyTrue()? 0b0100 : 0b0001;
+			return a; // barycentric coordinates (1,0,0)
 		}
 
 		// Check if P in vertex region outside B 
@@ -223,42 +225,39 @@ namespace ClosestPoint
 		}
 
 		// Check if P in edge region of AB, if so return projection of P onto AB 
-		float vc = d1 * d4 - d3 * d2; 
-		if (vc <= 0.0f && d1 >= 0.0f && d3 <= 0.0f) 
+		if (d1 * d4 <= d3 * d2 && d1 >= 0.0f && d3 <= 0.0f) 
 		{ 
 			float v = d1 / (d1 - d3); 
-			outSet = 0b0011;
-			return inA + v * ab; // barycentric coordinates (1-v,v,0) 
+			outSet = swap_ac.TestAnyTrue()? 0b0110 : 0b0011;
+			return a + v * ab; // barycentric coordinates (1-v,v,0) 
 		}
 
 		// Check if P in vertex region outside C 
-		Vec3 cp = -inC; 
+		Vec3 cp = -c; 
 		float d5 = ab.Dot(cp); 
 		float d6 = ac.Dot(cp); 
 		if (d6 >= 0.0f && d5 <= d6) 
 		{
-			outSet = 0b0100;
-			return inC; // barycentric coordinates (0,0,1)
+			outSet = swap_ac.TestAnyTrue()? 0b0001 : 0b0100;
+			return c; // barycentric coordinates (0,0,1)
 		}
 
 		// Check if P in edge region of AC, if so return projection of P onto AC 
-		float vb = d5 * d2 - d1 * d6; 
-		if (vb <= 0.0f && d2 >= 0.0f && d6 <= 0.0f) 
+		if (d5 * d2 <= d1 * d6 && d2 >= 0.0f && d6 <= 0.0f) 
 		{ 
 			float w = d2 / (d2 - d6); 
 			outSet = 0b0101;
-			return inA + w * ac; // barycentric coordinates (1-w,0,w) 
+			return a + w * ac; // barycentric coordinates (1-w,0,w) 
 		}
 
 		// Check if P in edge region of BC, if so return projection of P onto BC 
-		float va = d3 * d6 - d5 * d4;
 		float d4_d3 = d4 - d3;
 		float d5_d6 = d5 - d6;
-		if (va <= 0.0f && d4_d3 >= 0.0f && d5_d6 >= 0.0f) 
+		if (d3 * d6 <= d5 * d4 && d4_d3 >= 0.0f && d5_d6 >= 0.0f) 
 		{ 
 			float w = d4_d3 / (d4_d3 + d5_d6); 
-			outSet = 0b0110;
-			return inB + w * bc; // barycentric coordinates (0,1-w,w) 
+			outSet = swap_ac.TestAnyTrue()? 0b0011 : 0b0110;
+			return inB + w * (c - inB); // barycentric coordinates (0,1-w,w) 
 		}
 
 		// P inside face region.
@@ -268,7 +267,7 @@ namespace ClosestPoint
 		// Note that this way of calculating the closest point is much more accurate than first calculating barycentric coordinates 
 		// and then calculating the closest point based on those coordinates.
 		outSet = 0b0111;
-		return n * (inA + inB + inC).Dot(n) / (3.0f * n_len_sq);
+		return n * (a + inB + c).Dot(n) / (3.0f * n_len_sq);
 	}
 
 	/// Check if the origin is outside the plane of triangle (inA, inB, inC). inD specifies the front side of the plane.

UnitTests/Geometry/ClosestPointTests.cpp

@@ -6,17 +6,88 @@
 
 TEST_SUITE("ClosestPointTests")
 {
-	TEST_CASE("TestNearDegenerateTriangle")
+	// Test closest point from inPoint to triangle (inA, inB, inC)
+	inline static void TestClosestPointToTriangle(Vec3Arg inA, Vec3Arg inB, Vec3Arg inC, Vec3Arg inPoint, Vec3Arg inExpectedClosestPoint, uint32 inExpectedSet)
+	{
+		// Make triangle relative to inPoint so we can get the closest poin to the origin
+		Vec3 a = inA - inPoint;
+		Vec3 b = inB - inPoint;
+		Vec3 c = inC - inPoint;
+
+		// Extract bits for A, B and C
+		uint32 expected_a = inExpectedSet & 1;
+		uint32 expected_b = (inExpectedSet & 2) >> 1;
+		uint32 expected_c = (inExpectedSet & 4) >> 2;
+
+		// Test all permutations of ABC
+		uint32 set = 0;
+		Vec3 closest_point = inPoint + ClosestPoint::GetClosestPointOnTriangle(a, b, c, set);
+		CHECK(set == inExpectedSet);
+		CHECK_APPROX_EQUAL(closest_point, inExpectedClosestPoint, 2.0e-5f);
+
+		closest_point = inPoint + ClosestPoint::GetClosestPointOnTriangle(a, c, b, set);
+		CHECK(set == ((expected_b << 2) | (expected_c << 1) | expected_a));
+		CHECK_APPROX_EQUAL(closest_point, inExpectedClosestPoint, 2.0e-5f);
+
+		closest_point = inPoint + ClosestPoint::GetClosestPointOnTriangle(b, a, c, set);
+		CHECK(set == ((expected_c << 2) | (expected_a << 1) | expected_b));
+		CHECK_APPROX_EQUAL(closest_point, inExpectedClosestPoint, 2.0e-5f);
+
+		closest_point = inPoint + ClosestPoint::GetClosestPointOnTriangle(b, c, a, set);
+		CHECK(set == ((expected_a << 2) | (expected_c << 1) | expected_b));
+		CHECK_APPROX_EQUAL(closest_point, inExpectedClosestPoint, 2.0e-5f);
+
+		closest_point = inPoint + ClosestPoint::GetClosestPointOnTriangle(c, a, b, set);
+		CHECK(set == ((expected_b << 2) | (expected_a << 1) | expected_c));
+		CHECK_APPROX_EQUAL(closest_point, inExpectedClosestPoint, 2.0e-5f);
+
+		closest_point = inPoint + ClosestPoint::GetClosestPointOnTriangle(c, b, a, set);
+		CHECK(set == ((expected_a << 2) | (expected_b << 1) | expected_c));
+		CHECK_APPROX_EQUAL(closest_point, inExpectedClosestPoint, 2.0e-5f);
+	}
+
+	TEST_CASE("TestLongTriangle")
+	{
+		Vec3 a(100, 1, 0);
+		Vec3 b(100, 1, 1);
+		Vec3 c(-100, 1, 0);
+
+		// Test interior
+		TestClosestPointToTriangle(a, b, c, Vec3(0, 0, 0.1f), Vec3(0, 1, 0.1f), 0b0111);
+
+		// Edge AB
+		TestClosestPointToTriangle(a, b, c, Vec3(101, 0, 0.5f), Vec3(100, 1, 0.5f), 0b0011);
+
+		// Edge AC
+		TestClosestPointToTriangle(a, b, c, Vec3(0, 0, -0.1f), Vec3(0, 1, 0), 0b0101);
+
+		// Edge BC
+		Vec3 point_bc(0, 0, 1);
+		Vec3 bc = c - b;
+		Vec3 closest_bc = b + ((point_bc - b).Dot(bc) / bc.LengthSq()) * bc;
+		TestClosestPointToTriangle(a, b, c, point_bc, closest_bc, 0b0110);
+
+		// Vertex A
+		TestClosestPointToTriangle(a, b, c, Vec3(101, 0, -1), a, 0b0001);
+
+		// Vertex B
+		TestClosestPointToTriangle(a, b, c, Vec3(101, 0, 2), b, 0b0010);
+
+		// Vertex C
+		TestClosestPointToTriangle(a, b, c, Vec3(-101, 0, 0), c, 0b0100);
+	}
+	
+	TEST_CASE("TestNearColinearTriangle")
 	{
 		// A very long triangle that is nearly colinear
 		Vec3 a(99.9999847f, 0.946687222f, 99.9999847f);
 		Vec3 b(-100.010002f, 0.977360725f, -100.010002f);
 		Vec3 c(-100.000137f, 0.977310658f, -100.000137f);
 
-		uint32 set;
-		Vec3 p = ClosestPoint::GetClosestPointOnTriangle(a, b, c, set);
+		// Closest point is on edge AC
+		Vec3 ac = c - a;
+		Vec3 expected_closest = a + (-a.Dot(ac) / ac.LengthSq()) * ac;
 
-		CHECK(set == 0x3);
-		CHECK_APPROX_EQUAL(p, Vec3(7.62939453e-05f, 0.962023199f, 7.62939453e-05f));
+		TestClosestPointToTriangle(a, b, c, Vec3::sZero(), expected_closest, 0b0101);
 	}
 }

UnitTests/Physics/CastShapeTests.cpp

@@ -223,10 +223,10 @@ TEST_SUITE("CastShapeTests")
 			CHECK(collector.HadHit());
 			CHECK(collector.mHit.mBodyID2 == bodies.front()->GetID());
 			CHECK_APPROX_EQUAL(collector.mHit.mFraction, 4.0f / 10.0f);
-			CHECK_APPROX_EQUAL(collector.mHit.mPenetrationAxis.Normalized(), Vec3(1, 0, 0), 2.0e-3f);
+			CHECK_APPROX_EQUAL(collector.mHit.mPenetrationAxis.Normalized(), Vec3(1, 0, 0), 4.0e-3f);
 			CHECK_APPROX_EQUAL(collector.mHit.mPenetrationDepth, 0.0f);
-			CHECK_APPROX_EQUAL(collector.mHit.mContactPointOn1, Vec3(0, 0, 0), 1.0e-4f);
-			CHECK_APPROX_EQUAL(collector.mHit.mContactPointOn2, Vec3(0, 0, 0), 1.0e-4f);
+			CHECK_APPROX_EQUAL(collector.mHit.mContactPointOn1, Vec3(0, 0, 0), 2.0e-3f);
+			CHECK_APPROX_EQUAL(collector.mHit.mContactPointOn2, Vec3(0, 0, 0), 2.0e-3f);
 			CHECK(!collector.mHit.mIsBackFaceHit);
 		}