Skip to main content

falcon_library/model3d/
math.rs

1//! Small vector and matrix helpers for the 3D viewer and the thumbnail
2//! renderer. Vectors are `[f32; 3]`; matrices are column-major `[[f32; 4]; 4]`
3//! (`m[column][row]`), as WGSL reads them. The camera space is right-handed
4//! (the camera looks down −Z) and clip depth runs from 0 (near) to 1 (far),
5//! as wgpu expects.
6
7pub type V3 = [f32; 3];
8pub type M4 = [[f32; 4]; 4];
9
10pub fn add(a: V3, b: V3) -> V3 {
11    [a[0] + b[0], a[1] + b[1], a[2] + b[2]]
12}
13
14pub fn sub(a: V3, b: V3) -> V3 {
15    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
16}
17
18pub fn scale(a: V3, s: f32) -> V3 {
19    [a[0] * s, a[1] * s, a[2] * s]
20}
21
22pub fn dot(a: V3, b: V3) -> f32 {
23    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
24}
25
26pub fn cross(a: V3, b: V3) -> V3 {
27    [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]]
28}
29
30pub fn length(a: V3) -> f32 {
31    dot(a, a).sqrt()
32}
33
34/// `a` scaled to length 1; the zero vector stays zero.
35pub fn normalize(a: V3) -> V3 {
36    let l = length(a);
37    if l > 0.0 { scale(a, 1.0 / l) } else { a }
38}
39
40/// `p` turned by `angle` radians about the axis through `origin` along the
41/// unit vector `axis`, counter-clockwise looking down the axis (Rodrigues).
42pub fn rotate_about(p: V3, origin: V3, axis: V3, angle: f32) -> V3 {
43    let v = sub(p, origin);
44    let (s, c) = angle.sin_cos();
45    let r = add(add(scale(v, c), scale(cross(axis, v), s)), scale(axis, dot(axis, v) * (1.0 - c)));
46    add(origin, r)
47}
48
49pub fn mul(a: &M4, b: &M4) -> M4 {
50    let mut m = [[0.0; 4]; 4];
51    for (c, col) in m.iter_mut().enumerate() {
52        for (r, cell) in col.iter_mut().enumerate() {
53            *cell = (0..4).map(|k| a[k][r] * b[c][k]).sum();
54        }
55    }
56    m
57}
58
59/// `m` applied to the point `p`, after the divide by w.
60pub fn transform_point(m: &M4, p: V3) -> V3 {
61    let v = [p[0], p[1], p[2], 1.0];
62    let mut o = [0.0f32; 4];
63    for (r, out) in o.iter_mut().enumerate() {
64        *out = (0..4).map(|k| m[k][r] * v[k]).sum();
65    }
66    let w = if o[3] != 0.0 { o[3] } else { 1.0 };
67    [o[0] / w, o[1] / w, o[2] / w]
68}
69
70/// The view matrix of a camera at `eye` looking at `center` with `up`
71/// pointing up the screen.
72pub fn look_at(eye: V3, center: V3, up: V3) -> M4 {
73    let f = normalize(sub(center, eye));
74    let s = normalize(cross(f, up));
75    let u = cross(s, f);
76    [
77        [s[0], u[0], -f[0], 0.0],
78        [s[1], u[1], -f[1], 0.0],
79        [s[2], u[2], -f[2], 0.0],
80        [-dot(s, eye), -dot(u, eye), dot(f, eye), 1.0],
81    ]
82}
83
84/// A perspective projection: `fovy` is the vertical field of view in
85/// radians, `aspect` width over height.
86pub fn perspective(fovy: f32, aspect: f32, near: f32, far: f32) -> M4 {
87    let f = 1.0 / (fovy / 2.0).tan();
88    [
89        [f / aspect, 0.0, 0.0, 0.0],
90        [0.0, f, 0.0, 0.0],
91        [0.0, 0.0, far / (near - far), -1.0],
92        [0.0, 0.0, near * far / (near - far), 0.0],
93    ]
94}
95
96/// An orthographic projection showing `half_h` units above and below the
97/// centre of the screen.
98pub fn orthographic(half_h: f32, aspect: f32, near: f32, far: f32) -> M4 {
99    let half_w = half_h * aspect;
100    [
101        [1.0 / half_w, 0.0, 0.0, 0.0],
102        [0.0, 1.0 / half_h, 0.0, 0.0],
103        [0.0, 0.0, 1.0 / (near - far), 0.0],
104        [0.0, 0.0, near / (near - far), 1.0],
105    ]
106}
107
108#[cfg(test)]
109mod tests {
110    use super::*;
111
112    fn close(a: V3, b: V3) -> bool {
113        (0..3).all(|i| (a[i] - b[i]).abs() < 1e-4)
114    }
115
116    #[test]
117    fn rotation_turns_counter_clockwise() {
118        let p = rotate_about([1.0, 0.0, 0.0], [0.0; 3], [0.0, 0.0, 1.0], std::f32::consts::FRAC_PI_2);
119        assert!(close(p, [0.0, 1.0, 0.0]));
120        // About an axis through another point.
121        let p = rotate_about([2.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], std::f32::consts::PI);
122        assert!(close(p, [0.0, 0.0, 0.0]));
123    }
124
125    #[test]
126    fn look_at_puts_the_target_in_front() {
127        let v = look_at([0.0, 0.0, 5.0], [0.0; 3], [0.0, 1.0, 0.0]);
128        assert!(close(transform_point(&v, [0.0; 3]), [0.0, 0.0, -5.0]));
129        assert!(close(transform_point(&v, [0.0, 1.0, 0.0]), [0.0, 1.0, -5.0]));
130        assert!(close(transform_point(&v, [1.0, 0.0, 0.0]), [1.0, 0.0, -5.0]));
131    }
132
133    #[test]
134    fn projections_map_near_and_far_to_zero_and_one() {
135        let p = perspective(1.0, 1.5, 0.5, 50.0);
136        assert!((transform_point(&p, [0.0, 0.0, -0.5])[2]).abs() < 1e-5);
137        assert!((transform_point(&p, [0.0, 0.0, -50.0])[2] - 1.0).abs() < 1e-4);
138        // The edge of the field of view lands on the edge of the screen.
139        let y = (0.5f32).tan() * 10.0;
140        assert!((transform_point(&p, [0.0, y, -10.0])[1] - 1.0).abs() < 1e-4);
141        let o = orthographic(2.0, 2.0, 1.0, 9.0);
142        assert!(close(transform_point(&o, [4.0, 2.0, -1.0]), [1.0, 1.0, 0.0]));
143        assert!((transform_point(&o, [0.0, 0.0, -9.0])[2] - 1.0).abs() < 1e-5);
144        let id = mul(&o, &[[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]);
145        assert_eq!(id, o);
146    }
147}