dualis_core/vector.rs
1//! Vector arithmetic that no single domain owns.
2//!
3//! An orthonormal basis about a direction is needed to sample a hemisphere, to
4//! orient a scatter lobe and to build a local frame for a contact normal.
5//! Reflection about a plane is a light ray off a mirror and a ball off a wall, in
6//! the same three lines. Neither belongs to optics, so neither lives there.
7//!
8//! Refraction does belong to optics โ Snell's law is about refractive indices โ
9//! and is in `dualis-optics` accordingly.
10
11use glam::DVec3;
12
13/// Two unit vectors completing an orthonormal frame with `dir`.
14///
15/// `dir` must be normalised. The helper axis is chosen by which component of
16/// `dir` is small, which keeps the cross product well conditioned โ picking a
17/// fixed helper would give a near-zero cross product whenever `dir` happened to
18/// be parallel to it, and the resulting basis would be numerical noise.
19pub fn basis_for(dir: DVec3) -> (DVec3, DVec3) {
20 let helper = if dir.x.abs() < 0.9 {
21 DVec3::X
22 } else {
23 DVec3::Y
24 };
25 let u = helper.cross(dir).normalize();
26 let v = dir.cross(u);
27 (u, v)
28}
29
30/// Flip `n` so it opposes the incident direction (`d ยท n < 0`).
31///
32/// A surface has two sides and a stored normal points at one of them. Whether a
33/// ray is arriving from the front or the back is a fact about the ray, not the
34/// surface, so the normal is oriented per hit rather than per surface.
35pub fn oriented_against(n: DVec3, d: DVec3) -> DVec3 {
36 if d.dot(n) > 0.0 {
37 -n
38 } else {
39 n
40 }
41}
42
43/// Mirror `d` about the plane with normal `n`. Both should be normalised.
44pub fn reflect(d: DVec3, n: DVec3) -> DVec3 {
45 (d - n * (2.0 * d.dot(n))).normalize()
46}
47
48#[cfg(test)]
49mod tests {
50 use super::*;
51
52 /// The basis is orthonormal and right-handed for every direction, including
53 /// the axis-aligned ones where a badly chosen helper vector would collapse it.
54 #[test]
55 fn the_basis_is_orthonormal_everywhere() {
56 let directions = [
57 DVec3::X,
58 DVec3::Y,
59 DVec3::Z,
60 -DVec3::X,
61 -DVec3::Z,
62 DVec3::new(1.0, 1.0, 1.0).normalize(),
63 DVec3::new(0.9999, 0.01, 0.0).normalize(),
64 DVec3::new(-0.7, 0.2, 0.68).normalize(),
65 ];
66 for d in directions {
67 let (u, v) = basis_for(d);
68 assert!((u.length() - 1.0).abs() < 1e-12, "u not unit for {d}");
69 assert!((v.length() - 1.0).abs() < 1e-12, "v not unit for {d}");
70 assert!(u.dot(d).abs() < 1e-12, "u not perpendicular to {d}");
71 assert!(v.dot(d).abs() < 1e-12, "v not perpendicular to {d}");
72 assert!(u.dot(v).abs() < 1e-12, "u and v not perpendicular for {d}");
73 // Right-handed: u x v = d.
74 assert!((u.cross(v) - d).length() < 1e-12, "left-handed for {d}");
75 }
76 }
77
78 /// Reflection reverses the normal component and keeps the tangential one, so
79 /// the angle in equals the angle out and the reflection of a reflection is the
80 /// original.
81 #[test]
82 fn reflection_preserves_the_angle_and_is_its_own_inverse() {
83 let n = DVec3::new(0.0, 1.0, 0.0);
84 let d = DVec3::new(1.0, -1.0, 0.0).normalize();
85 let r = reflect(d, n);
86 assert!((r - DVec3::new(1.0, 1.0, 0.0).normalize()).length() < 1e-12);
87 // Equal angles either side.
88 assert!((d.dot(n).abs() - r.dot(n).abs()).abs() < 1e-12);
89 // Reflecting the reflection comes back, up to direction.
90 assert!((reflect(r, n) - d).length() < 1e-12);
91 // Normal incidence turns straight around.
92 assert!((reflect(-n, n) - n).length() < 1e-12);
93 }
94
95 #[test]
96 fn a_normal_is_oriented_per_hit_not_per_surface() {
97 let n = DVec3::Z;
98 // Arriving from the front, the normal is left alone.
99 assert_eq!(oriented_against(n, -DVec3::Z), n);
100 // Arriving from behind, it is flipped to face the ray.
101 assert_eq!(oriented_against(n, DVec3::Z), -n);
102 // And either way the result opposes the ray, which is the contract.
103 for d in [DVec3::Z, -DVec3::Z, DVec3::new(0.3, 0.1, 0.9).normalize()] {
104 assert!(oriented_against(n, d).dot(d) <= 0.0);
105 }
106 }
107}