From 6b410963ca6f428b33b618527529ad838372e5aa Mon Sep 17 00:00:00 2001
From: snowyukitty <270071858+snowyukitty@users.noreply.github.com>
Date: Fri, 9 Oct 2026 18:39:37 +0900
Subject: [PATCH] Fix zero-area cancellation in tetrahedron face projection

---
 .../tests/issue_451_convex_hull_raycast.rs    | 61 +++++++++++++++++++
 src/query/point/point_tetrahedron.rs          | 18 ++++--
 2 files changed, 74 insertions(+), 5 deletions(-)
 create mode 100644 crates/parry3d/tests/issue_451_convex_hull_raycast.rs

diff --git a/crates/parry3d/tests/issue_451_convex_hull_raycast.rs b/crates/parry3d/tests/issue_451_convex_hull_raycast.rs
new file mode 100644
index 00000000..b80a1e3d
--- /dev/null
+++ b/crates/parry3d/tests/issue_451_convex_hull_raycast.rs
@@ -0,0 +1,61 @@
+use parry3d::math::Vector;
+use parry3d::query::{PointQueryWithLocation, Ray};
+use parry3d::shape::{SharedShape, Tetrahedron, TetrahedronPointLocation};
+
+#[test]
+fn tetrahedron_distant_face_projection() {
+    let a = Vector::ZERO;
+    let b = Vector::new(1.0, -1.0, 0.0);
+    let c = Vector::new(1.0, 0.0, -1.0);
+    let d = Vector::splat(-1.0);
+    let expected = Vector::new(0.5, -0.25, -0.25);
+    for (face, tetrahedron) in [
+        Tetrahedron::new(a, b, c, d),
+        Tetrahedron::new(a, b, d, c),
+        Tetrahedron::new(a, d, b, c),
+        Tetrahedron::new(d, a, b, c),
+    ]
+    .into_iter()
+    .enumerate()
+    {
+        for distance in [0.25, 4096.0] {
+            let (projection, location) = tetrahedron
+                .project_local_point_and_get_location(expected + Vector::splat(distance), true);
+            assert!(!projection.is_inside);
+            assert!((projection.point - expected).length() < 1.0e-3);
+            let TetrahedronPointLocation::OnFace(index, coordinates) = location else {
+                panic!("unexpected location: {location:?}");
+            };
+            assert_eq!(index as usize, face);
+            for (actual, expected) in coordinates.into_iter().zip([0.5, 0.25, 0.25]) {
+                assert!((actual - expected).abs() < 1.0e-3);
+            }
+        }
+    }
+}
+
+fn hull() -> SharedShape {
+    SharedShape::convex_hull(&[
+        Vector::new(-34.0, -14.0, 29.0),
+        Vector::new(-24.0, -12.0, 29.0),
+        Vector::new(-32.0, -12.0, 29.0),
+        Vector::new(-32.0, -12.0, 2.0),
+        Vector::new(-24.0, -12.0, 2.0),
+    ])
+    .unwrap()
+}
+
+#[test]
+fn convex_hull_nonintersecting_ray() {
+    let ray = Ray::new(
+        Vector::new(-16.0, 8.0, 0.0),
+        Vector::new(-0.5956519, -0.67131174, 0.44106603),
+    );
+    assert!(!hull().intersects_local_ray(&ray, f32::MAX));
+}
+
+#[test]
+fn convex_hull_intersecting_ray() {
+    let ray = Ray::new(Vector::new(-28.0, 0.0, 10.0), -Vector::Y);
+    assert!(hull().intersects_local_ray(&ray, f32::MAX));
+}
diff --git a/src/query/point/point_tetrahedron.rs b/src/query/point/point_tetrahedron.rs
index faa3df65..9c271c30 100644
--- a/src/query/point/point_tetrahedron.rs
+++ b/src/query/point/point_tetrahedron.rs
@@ -268,11 +268,19 @@ impl PointQueryWithLocation for Tetrahedron {
                     // numbers and there are not enough decimal bits to perform
                     // the normalization.
                     let normal = n.try_normalize()?;
-                    let vc = normal.dot(ap.cross(bp));
-                    let va = normal.dot(bp.cross(cp));
-                    let vb = normal.dot(cp.cross(ap));
-
-                    let denom = va + vb + vc;
+                    let mut vc = normal.dot(ap.cross(bp));
+                    let mut va = normal.dot(bp.cross(cp));
+                    let mut vb = normal.dot(cp.cross(ap));
+
+                    let mut denom = va + vb + vc;
+                    if denom == 0.0 {
+                        // Use face edges to avoid cancellation between products of
+                        // nearly parallel vectors when the query point is far away.
+                        vc = normal.dot(ab.cross(ap));
+                        va = normal.dot((c - b).cross(bp));
+                        vb = normal.dot(cp.cross(ac));
+                        denom = va + vb + vc;
+                    }
                     assert!(denom != 0.0);
                     let inv_denom = 1.0 / denom;