spherical_cow/util/mod.rs
1//! Useful helper functions such as a fast ray casting method and volume finder for use with arbitrary shaped triangular meshes.
2
3use nalgebra::{Matrix, Point3, Vector3};
4
5/// Find the baycentric coordinates `(u,v)` and distance `t` given three triangle veriticies `vert0`, `vert1`, `vert2` and the
6/// unit vector `dir` (`D`) in the direction of a ray `R(t) = O + tD` such that `R(t)` is equivalent to a point `T(u,v)` on
7/// triangle `T`. Here, `T(u,v) = (1-u-v)V_0+uV_1+uV_2`. If distance `t` is less than the distance of our sphere from the
8/// origin (`sphere_dist`), then add one to the count.
9pub fn ray_intersection_count(
10 triangles: &[(Point3<f32>, Point3<f32>, Point3<f32>)],
11 dir: Vector3<f32>,
12 o_dist: f32,
13) -> i32 {
14 triangles
15 .iter()
16 .map(|triangle| {
17 let &(vert0, vert1, vert2) = triangle;
18 let edge1 = vert1 - vert0;
19 let edge2 = vert2 - vert0;
20 let pvec = Matrix::cross(&dir, &edge2);
21 let det = Matrix::dot(&edge1, &pvec);
22 if det > -1e-6 && det < 1e-6 {
23 return 0;
24 }
25 let inv_det = 1. / det;
26 let tvec = Point3::origin() - vert0;
27 let u = Matrix::dot(&tvec, &pvec) * inv_det;
28 if !(0. ..=1.).contains(&u) {
29 return 0;
30 }
31 let qvec = Matrix::cross(&tvec, &edge1);
32 let v = Matrix::dot(&dir, &qvec) * inv_det;
33 if v < 0. || u + v > 1. {
34 return 0;
35 }
36 // Our ray R(t) = O + tD can now be calculated
37 let t = Matrix::dot(&edge2, &qvec) * inv_det;
38 if t.is_sign_positive() && t < o_dist {
39 1
40 } else {
41 0
42 }
43 })
44 .sum()
45}
46
47/// Identify the volume of a trimesh which contains the cartesian origin (0, 0, 0).
48/// Consider the tetrahedron created by any triangle and the origin OABC and sum thier volumes.
49pub fn trimesh_volume(triangles: &[(Point3<f32>, Point3<f32>, Point3<f32>)]) -> f32 {
50 let sixth = 1. / 6.;
51 triangles
52 .iter()
53 .map(|triangle| {
54 let &(a, b, c) = triangle;
55 sixth
56 * (-c.x * b.y * a.z + b.x * c.y * a.z + c.x * a.y * b.z
57 - a.x * c.y * b.z
58 - b.x * a.y * c.z
59 + a.x * b.y * c.z)
60 })
61 .sum()
62}