Skip to main content

loam_math/
euclidean_r4.rs

1use std::borrow::Cow;
2
3use glam::Vec4;
4use serde::{Deserialize, Serialize};
5
6use crate::bivector::{Rotor, Rotor4};
7use crate::space::{IsometryGroup, Space, WgslAccuracy, WgslSpace, FLAT_CHART_MAX_ARC};
8
9/// Rigid motion of R⁴; callers must keep the rotor unit length.
10#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
11pub struct Iso4Flat {
12    pub rotation: Rotor4,
13    pub translation: Vec4,
14}
15
16impl Iso4Flat {
17    pub const IDENTITY: Self = Self {
18        rotation: Rotor4::IDENTITY,
19        translation: Vec4::ZERO,
20    };
21
22    pub fn from_rotation(rotation: Rotor4) -> Self {
23        Self {
24            rotation,
25            translation: Vec4::ZERO,
26        }
27    }
28
29    pub fn from_translation(translation: Vec4) -> Self {
30        Self {
31            rotation: Rotor4::IDENTITY,
32            translation,
33        }
34    }
35}
36
37impl Default for Iso4Flat {
38    fn default() -> Self {
39        Self::IDENTITY
40    }
41}
42
43#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
44pub struct EuclideanR4;
45
46impl Space for EuclideanR4 {
47    type Point = Vec4;
48    type Vector = Vec4;
49    type Frame = Rotor4;
50
51    fn frame_at(&self, _at: Vec4) -> Rotor4 {
52        Rotor4::IDENTITY
53    }
54
55    fn distance(&self, a: Vec4, b: Vec4) -> f32 {
56        (a - b).length()
57    }
58
59    fn exp(&self, at: Vec4, v: Vec4) -> Vec4 {
60        at + v
61    }
62
63    fn log(&self, from: Vec4, to: Vec4) -> Vec4 {
64        to - from
65    }
66
67    fn parallel_transport(&self, _from: Vec4, _to: Vec4, v: Vec4) -> Vec4 {
68        v
69    }
70
71    fn is_chart_flat(&self) -> bool {
72        true
73    }
74
75    fn chart_envelope(&self) -> f32 {
76        f32::INFINITY
77    }
78
79    fn valid_point(&self, p: Vec4) -> bool {
80        p.is_finite()
81    }
82}
83
84impl IsometryGroup for EuclideanR4 {
85    type Iso = Iso4Flat;
86
87    fn iso_identity(&self) -> Iso4Flat {
88        Iso4Flat::IDENTITY
89    }
90
91    fn iso_compose(&self, a: Iso4Flat, b: Iso4Flat) -> Iso4Flat {
92        Iso4Flat {
93            rotation: b.rotation * a.rotation,
94            translation: a.rotation.apply(b.translation) + a.translation,
95        }
96    }
97
98    fn iso_inverse(&self, a: Iso4Flat) -> Iso4Flat {
99        let inv_rot = a.rotation.inverse();
100        Iso4Flat {
101            rotation: inv_rot,
102            translation: inv_rot.apply(-a.translation),
103        }
104    }
105
106    fn iso_apply(&self, iso: Iso4Flat, p: Vec4) -> Vec4 {
107        iso.rotation.apply(p) + iso.translation
108    }
109
110    fn iso_transport(&self, iso: Iso4Flat, _at: Vec4, v: Vec4) -> Vec4 {
111        iso.rotation.apply(v)
112    }
113}
114
115impl WgslSpace for EuclideanR4 {
116    fn wgsl_impl(&self) -> Cow<'static, str> {
117        Cow::Owned(format!(
118            "const LOAM_MAX_ARC: f32 = {FLAT_CHART_MAX_ARC:?};{WGSL_IMPL}"
119        ))
120    }
121
122    fn wgsl_accuracy(&self) -> WgslAccuracy {
123        WgslAccuracy::Exact
124    }
125}
126
127const WGSL_IMPL: &str = r#"
128// loam-math :: EuclideanR4 (v0 Space WGSL ABI)
129fn loam_distance(a: vec4<f32>, b: vec4<f32>) -> f32 { return length(a - b); }
130fn loam_origin_distance(p: vec4<f32>) -> f32 { return length(p); }
131fn loam_exp(at: vec4<f32>, v: vec4<f32>) -> vec4<f32> { return at + v; }
132fn loam_log(p_from: vec4<f32>, p_to: vec4<f32>) -> vec4<f32> { return p_to - p_from; }
133struct LoamGeodesicStep { p: vec4<f32>, v: vec4<f32> }
134fn loam_geodesic_step(p: vec4<f32>, v: vec4<f32>, s: f32) -> LoamGeodesicStep {
135    return LoamGeodesicStep(loam_exp(p, v * s), v);
136}
137"#;
138
139#[cfg(test)]
140mod tests {
141    use super::*;
142    use approx::assert_relative_eq;
143
144    fn r4() -> EuclideanR4 {
145        EuclideanR4
146    }
147
148    #[test]
149    fn iso_transport_ignores_translation() {
150        let s = r4();
151        let iso = Iso4Flat::from_translation(Vec4::new(100.0, -50.0, 7.0, 12.0));
152        let v = Vec4::new(1.0, 0.0, 0.0, 0.0);
153        let transported = s.iso_transport(iso, Vec4::ZERO, v);
154        assert_relative_eq!(transported.x, v.x);
155        assert_relative_eq!(transported.y, v.y);
156        assert_relative_eq!(transported.z, v.z);
157        assert_relative_eq!(transported.w, v.w);
158    }
159}