From 7276d0965fb32593735f0ebebf8e432d235dc2c9 Mon Sep 17 00:00:00 2001 From: Nobyuki Umetani Date: Wed, 1 Jan 2025 10:16:04 +0900 Subject: [PATCH] tri3 ray intersection --- del-geo-core/src/tri3.rs | 186 +++++++++++++++++++++++++++++++++ del-geo-core/src/vec3.rs | 25 +++-- del-geo-cpp-headers/src/vec3.h | 10 +- 3 files changed, 207 insertions(+), 14 deletions(-) diff --git a/del-geo-core/src/tri3.rs b/del-geo-core/src/tri3.rs index 0d6cae5..bafd7d2 100644 --- a/del-geo-core/src/tri3.rs +++ b/del-geo-core/src/tri3.rs @@ -168,6 +168,192 @@ where Some(t) } +pub struct RayTriangleIntersectionData { + 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( + org: &[T; 3], + dir: &[T; 3], + p0: &[T; 3], + p1: &[T; 3], + p2: &[T; 3], +) -> Option<(T, T, T, RayTriangleIntersectionData)> +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_( + d_t: Real, + d_u: Real, + d_v: Real, + data: &RayTriangleIntersectionData, +) -> ([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::::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_::(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(p0: &[T; 3], p1: &[T; 3], p2: &[T; 3], q: &[T; 3]) -> [T; 3] where T: num_traits::Float + MulAssign, diff --git a/del-geo-core/src/vec3.rs b/del-geo-core/src/vec3.rs index 127d2ea..6421efb 100644 --- a/del-geo-core/src/vec3.rs +++ b/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 where @@ -25,7 +24,7 @@ where impl Vec3 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(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(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(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(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(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(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) diff --git a/del-geo-cpp-headers/src/vec3.h b/del-geo-cpp-headers/src/vec3.h index 08df770..34fd9d4 100644 --- a/del-geo-cpp-headers/src/vec3.h +++ b/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 +{ + return { + a * x[0] + y[0], + a * x[1] + y[1], + a * x[2] + y[2], + }; +} } \ No newline at end of file