Skip to main content

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}