tri3 ray intersection

del-geo-core/src/tri3.rs

@@ -168,6 +168,192 @@ where
     Some(t)
 }
 
+pub struct RayTriangleIntersectionData<Real> {
+    inv_det: Real,
+    n_dot_c: Real,
+    n_dot_dir: Real,
+    r_dot_e2: Real,
+    r_dot_e1: Real,
+    c: [Real; 3],
+    n: [Real; 3],
+    r: [Real; 3],
+    e1: [Real; 3],
+    e2: [Real; 3],
+    pub dir: [Real; 3],
+}
+
+/// ray triangle intersection.
+/// * `dir` - any nonzero vector (not necessary to be a unit vector)
+/// * `t` - ratio of `dir` vector from
+/// * `u` - barycentric coordinate
+/// * `v` - barycentric coordinate
+///
+/// `org + t * dir = (1 - u - v) * p0 + u * p1 + v * p2`
+pub fn ray_triangle_intersection<T>(
+    org: &[T; 3],
+    dir: &[T; 3],
+    p0: &[T; 3],
+    p1: &[T; 3],
+    p2: &[T; 3],
+) -> Option<(T, T, T, RayTriangleIntersectionData<T>)>
+where
+    T: num_traits::Float,
+{
+    use crate::vec3::Vec3;
+    let e1 = p0.sub(p1);
+    let e2 = p2.sub(p0);
+    let n = e1.cross(&e2);
+    let n_dot_dir = n.dot(dir);
+    if n_dot_dir.is_zero() {
+        return None;
+    }
+    let inv_det = T::one() / n_dot_dir;
+    let c = p0.sub(org);
+    let n_dot_c = n.dot(&c);
+    let t_ = inv_det * n_dot_c;
+    if t_ < T::zero() {
+        return None;
+    }
+    let r = dir.cross(&c);
+    //
+    let r_dot_e2 = r.dot(&e2);
+    let u_ = inv_det * r_dot_e2;
+    if u_ < T::zero() {
+        return None;
+    }
+    //
+    let r_dot_e1 = r.dot(&e1);
+    let v_ = inv_det * r_dot_e1;
+    if v_ < T::zero() {
+        return None;
+    }
+
+    if u_ + v_ >= T::one() {
+        return None;
+    }
+    Some((
+        t_,
+        u_,
+        v_,
+        RayTriangleIntersectionData {
+            inv_det,
+            n_dot_c,
+            n_dot_dir,
+            r_dot_e2,
+            r_dot_e1,
+            c,
+            n,
+            r,
+            e1,
+            e2,
+            dir: *dir,
+        },
+    ))
+}
+
+#[test]
+fn test_ray_triangle_intersection() {
+    let p0 = [-13f32, -5., 8.];
+    let p1 = [14., -5., 8.];
+    let p2 = [1., 3., -3.];
+    let origin = [8., 11., 10.];
+    let dir = [-7., -11., -8.];
+    let t1 = intersection_against_line(&p0, &p1, &p2, &origin, &dir).unwrap();
+    let (t0, _u0, _v0, _data) = ray_triangle_intersection(&origin, &dir, &p0, &p1, &p2).unwrap();
+    assert!((t0 - t1).abs() < 1.0e-8);
+}
+
+pub fn dldw_ray_triangle_intersection_<Real>(
+    d_t: Real,
+    d_u: Real,
+    d_v: Real,
+    data: &RayTriangleIntersectionData<Real>,
+) -> ([Real; 3], [Real; 3], [Real; 3])
+where
+    Real: num_traits::Float + Copy,
+{
+    let d_n_dot_c = d_t * data.inv_det;
+    let d_r_dot_e2 = d_u * data.inv_det;
+    let d_r_dot_e1 = d_v * data.inv_det;
+    let d_inv_det = d_t * data.n_dot_c + d_u * data.r_dot_e2 + d_v * data.r_dot_e1;
+    //
+    let mut d_n = data.c.scale(d_n_dot_c);
+    let mut d_c = data.n.scale(d_n_dot_c);
+    let mut d_e2 = data.r.scale(d_r_dot_e2);
+    let mut d_e1 = data.r.scale(d_r_dot_e1);
+    let d_r = data.e2.scale(d_r_dot_e2).add(&data.e1.scale(d_r_dot_e1));
+    //
+    let d_n_dot_dir = -d_inv_det / data.n_dot_dir / data.n_dot_dir;
+    d_n.add_in_place(&data.dir.scale(d_n_dot_dir));
+    d_c.add_in_place(&d_r.cross(&data.dir));
+    d_e2.add_in_place(&d_n.cross(&data.e1));
+    d_e1.add_in_place(&data.e2.cross(&d_n));
+    //
+    let d_p0 = d_e1.sub(&d_e2).add(&d_c);
+    let d_p1 = d_e1.scale(-Real::one());
+    let d_p2 = d_e2;
+    (d_p0, d_p1, d_p2)
+}
+
+#[test]
+fn test_dw_ray_triangle_intersection() {
+    type Real = f64;
+    let p0: [[Real; 3]; 3] = [[-13., -5., 8.], [14., -5., 8.], [1., 3., -3.]];
+
+    // let origin = nalgebra::Vector3::<Real>::new(9., 11., 12.);
+    let origin = [8., 11., 10.];
+    let dir = [1., 0., 2.].sub(&origin);
+
+    let Some((t0, u0, v0, data)) = ray_triangle_intersection(&origin, &dir, &p0[0], &p0[1], &p0[2])
+    else {
+        panic!()
+    };
+
+    let ha = p0[0]
+        .scale(1. - u0 - v0)
+        .add(&p0[1].scale(u0))
+        .add(&p0[2].scale(v0));
+    assert!(crate::tet::volume(&p0[0], &p0[1], &p0[2], &ha).abs() < 1.0e-8);
+    let hb = crate::vec3::axpy(t0, &dir, &origin);
+    assert!(ha.sub(&hb).norm() < 1.0e-6);
+
+    // d_t, d_u, d_u are back-propagated from the loss
+    let d_t = 0.7;
+    let d_u = 1.3;
+    let d_v = 1.1;
+
+    let dp = dldw_ray_triangle_intersection_::<Real>(d_t, d_u, d_v, &data);
+    let dp = [dp.0, dp.1, dp.2];
+
+    let eps = 1.0e-5;
+    for i_node in 0..3 {
+        for i_dim in 0..3 {
+            let p1 = {
+                let mut p1 = p0;
+                p1[i_node][i_dim] += eps;
+                p1
+            };
+            let Some((t1, u1, v1, _data)) =
+                ray_triangle_intersection(&origin, &dir, &p1[0], &p1[1], &p1[2])
+            else {
+                panic!()
+            };
+            let dloss = ((t1 - t0) * d_t + (u1 - u0) * d_u + (v1 - v0) * d_v) / eps;
+            let dloss_analytic = dp[i_node][i_dim];
+            // dbg!(dloss, dloss_analytic);
+            assert!(
+                (dloss - dloss_analytic).abs() < 6.0e-7,
+                "{} {} {}",
+                dloss_analytic,
+                dloss,
+                dloss - dloss_analytic
+            );
+        }
+    }
+}
+
+// ----------------------------
+
 pub fn to_barycentric_coords<T>(p0: &[T; 3], p1: &[T; 3], p2: &[T; 3], q: &[T; 3]) -> [T; 3]
 where
     T: num_traits::Float + MulAssign,

del-geo-core/src/vec3.rs

@@ -1,6 +1,5 @@
 //! methods for 3D vector
-//!
-use std::ops::MulAssign;
+
 /// trait for 3D vector
 pub trait Vec3<Real>
 where
@@ -25,7 +24,7 @@ where
 
 impl<Real> Vec3<Real> for [Real; 3]
 where
-    Real: num_traits::Float + std::ops::MulAssign,
+    Real: num_traits::Float,
 {
     fn normalize(&self) -> Self {
         normalize(self)
@@ -76,7 +75,7 @@ where
 
 pub fn orthogonalize<Real>(u: &[Real; 3], v: &[Real; 3]) -> [Real; 3]
 where
-    Real: num_traits::Float + MulAssign,
+    Real: num_traits::Float,
 {
     let t = u.dot(v) / u.dot(u);
     [v[0] - t * u[0], v[1] - t * u[1], v[2] - t * u[2]]
@@ -210,20 +209,20 @@ where
 /// in-place normalize function
 pub fn normalize_in_place<T>(v: &mut [T; 3]) -> T
 where
-    T: num_traits::Float + std::ops::MulAssign,
+    T: num_traits::Float,
 {
     let l = v.norm();
     let linv = T::one() / l;
-    v[0] *= linv;
-    v[1] *= linv;
-    v[2] *= linv;
+    v[0] = v[0] * linv;
+    v[1] = v[1] * linv;
+    v[2] = v[2] * linv;
     l
 }
 
 /// return normalized 3D vector
 pub fn normalize<T>(v: &[T; 3]) -> [T; 3]
 where
-    T: num_traits::Float + MulAssign,
+    T: num_traits::Float,
 {
     let l = v.norm();
     let linv = T::one() / l;
@@ -277,16 +276,16 @@ where
 
 pub fn scale_in_place<T>(a: &mut [T; 3], s: T)
 where
-    T: MulAssign + Copy,
+    T: num_traits::Float,
 {
     for x in a.iter_mut() {
-        *x *= s;
+        *x = *x * s;
     }
 }
 
 pub fn distance<T>(p0: &[T; 3], p1: &[T; 3]) -> T
 where
-    T: num_traits::Float + MulAssign,
+    T: num_traits::Float,
 {
     let [v0, v1, v2] = p1.sub(p0);
     (v0 * v0 + v1 * v1 + v2 * v2).sqrt()
@@ -332,7 +331,7 @@ where
 
 pub fn mirror_reflection<Real>(v: &[Real; 3], nrm: &[Real; 3]) -> [Real; 3]
 where
-    Real: num_traits::Float + MulAssign,
+    Real: num_traits::Float,
 {
     let a = nrm.dot(v);
     std::array::from_fn(|i| v[i] - nrm[i] * Real::from(2).unwrap() * a)

del-geo-cpp-headers/src/vec3.h

@@ -38,7 +38,15 @@ auto dot(const float* a, const float* b) -> float
     return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
 }
 
-
+__device__
+auto axpy(float a, const float* x, const float* y) -> cuda::std::array<float,3>
+{
+    return {
+        a * x[0] + y[0],
+        a * x[1] + y[1],
+        a * x[2] + y[2],
+    };
+}
 
 
 }