loam_render/
hypergimbal.rs

1use glam::{Vec3, Vec4};
2use loam_math::{Bivector, Plane4, Rotor4};
3use loam_shape::LineMesh;
4
5/// A 16-cell cell centre, on no coordinate great circle: every ring is finite and congruent.
6pub const POLE: Vec4 = Vec4::new(0.5, 0.5, 0.5, 0.5);
7
8const MIN_PLANE_INCIDENCE: f32 = 1e-2;
9
10const MIN_POLE_DISTANCE: f32 = 1e-4;
11
12const IMAGE_AXES: [Vec4; 3] = [
13    Vec4::new(0.5, 0.5, -0.5, -0.5),
14    Vec4::new(0.5, -0.5, 0.5, -0.5),
15    Vec4::new(0.5, -0.5, -0.5, 0.5),
16];
17
18fn to_r3(q: Vec4) -> Vec3 {
19    Vec3::new(
20        IMAGE_AXES[0].dot(q),
21        IMAGE_AXES[1].dot(q),
22        IMAGE_AXES[2].dot(q),
23    )
24}
25
26fn plane_axes(plane: Plane4) -> (Vec4, Vec4) {
27    match plane {
28        Plane4::Xy => (Vec4::X, Vec4::Y),
29        Plane4::Xz => (Vec4::X, Vec4::Z),
30        Plane4::Xw => (Vec4::X, Vec4::W),
31        Plane4::Yz => (Vec4::Y, Vec4::Z),
32        Plane4::Yw => (Vec4::Y, Vec4::W),
33        Plane4::Zw => (Vec4::Z, Vec4::W),
34    }
35}
36
37fn wrap_pi(angle: f32) -> f32 {
38    use std::f32::consts::{PI, TAU};
39    let wrapped = angle.rem_euclid(TAU);
40    if wrapped > PI {
41        wrapped - TAU
42    } else {
43        wrapped
44    }
45}
46
47#[derive(Copy, Clone, Debug, PartialEq)]
48pub struct Hypergimbal {
49    pub center: Vec3,
50    pub scale: f32,
51}
52
53impl Hypergimbal {
54    /// `None` within f32 noise of [`POLE`], where the image is unbounded.
55    // Coxeter, Introduction to Geometry (1969), section 6.9.
56    pub fn project(&self, p: Vec4) -> Option<Vec3> {
57        let denom = 1.0 - p.dot(POLE);
58        if denom.abs() < MIN_POLE_DISTANCE {
59            return None;
60        }
61        Some(self.center + to_r3(p) / denom * self.scale)
62    }
63
64    /// All six rings, in [`Plane4::ALL`] order.
65    pub fn rings(&self) -> [Ring; 6] {
66        Plane4::ALL.map(|plane| self.ring(plane))
67    }
68
69    pub fn ring(&self, plane: Plane4) -> Ring {
70        let (a, b) = plane_axes(plane);
71        let alpha = a.dot(POLE);
72        let beta = b.dot(POLE);
73        let pole_overlap = (alpha * alpha + beta * beta).sqrt();
74        let s = (1.0 - pole_overlap * pole_overlap).sqrt();
75        let (a_phase, b_phase) = (
76            (a * alpha + b * beta) / pole_overlap,
77            (b * alpha - a * beta) / pole_overlap,
78        );
79        let u = to_r3((a_phase - POLE * pole_overlap) / s);
80        Ring {
81            plane,
82            center: self.center + u * (self.scale * pole_overlap / s),
83            u,
84            v: to_r3(b_phase),
85            radius: self.scale / s,
86            pole_overlap,
87            phase: beta.atan2(alpha),
88        }
89    }
90
91    /// Nearest ring within `tolerance` of the plane hit; ties break by depth.
92    pub fn pick(&self, ray_origin: Vec3, ray_direction: Vec3, tolerance: f32) -> Option<Ring> {
93        let mut best: Option<(f32, Ring)> = None;
94        for ring in self.rings() {
95            let Some(hit) = ring.ray_hit(ray_origin, ray_direction) else {
96                continue;
97            };
98            if ((hit - ring.center).length() - ring.radius).abs() > tolerance {
99                continue;
100            }
101            let depth = (hit - ray_origin).length();
102            if best.is_none_or(|(nearest, _)| depth < nearest) {
103                best = Some((depth, ring));
104            }
105        }
106        best.map(|(_, ring)| ring)
107    }
108
109    pub fn append_line_mesh(&self, style: &RingStyle, out: &mut LineMesh<3>) {
110        let step = std::f32::consts::TAU / style.segments as f32;
111        for ring in self.rings() {
112            let color = style.colors[ring.plane as usize];
113            let mut prev = ring.point(0.0);
114            for k in 1..=style.segments {
115                let next = ring.point(k as f32 * step);
116                out.segments.push((prev.to_array(), next.to_array()));
117                out.colors.push((color, color));
118                out.widths.push(style.width_px);
119                prev = next;
120            }
121        }
122    }
123}
124
125#[derive(Copy, Clone, Debug, PartialEq)]
126pub struct Ring {
127    pub plane: Plane4,
128    pub center: Vec3,
129    pub u: Vec3,
130    pub v: Vec3,
131    pub radius: f32,
132    pub pole_overlap: f32,
133    pub phase: f32,
134}
135
136impl Ring {
137    pub fn normal(&self) -> Vec3 {
138        self.u.cross(self.v)
139    }
140
141    pub fn point(&self, chi: f32) -> Vec3 {
142        self.center + (self.u * chi.cos() + self.v * chi.sin()) * self.radius
143    }
144
145    pub fn ring_angle(&self, world: Vec3) -> f32 {
146        let offset = world - self.center;
147        offset.dot(self.v).atan2(offset.dot(self.u))
148    }
149
150    // Danby, Fundamentals of Celestial Mechanics (1992), section 6.3.
151    pub fn arc_angle(&self, chi: f32) -> f32 {
152        let s = (1.0 - self.pole_overlap * self.pole_overlap).sqrt();
153        (s * chi.sin()).atan2(chi.cos() + self.pole_overlap) + self.phase
154    }
155
156    pub fn ray_hit(&self, ray_origin: Vec3, ray_direction: Vec3) -> Option<Vec3> {
157        let normal = self.normal();
158        let incidence = ray_direction.dot(normal);
159        if incidence.abs() < MIN_PLANE_INCIDENCE * ray_direction.length() {
160            return None;
161        }
162        let t = (self.center - ray_origin).dot(normal) / incidence;
163        (t > 0.0).then(|| ray_origin + ray_direction * t)
164    }
165
166    pub fn drag_angle(&self, grab: Vec3, cursor: Vec3) -> f32 {
167        wrap_pi(self.arc_angle(self.ring_angle(cursor)) - self.arc_angle(self.ring_angle(grab)))
168    }
169
170    pub fn drag_rotor(&self, grab: Vec3, cursor: Vec3) -> Rotor4 {
171        (self.plane.unit_bivector() * self.drag_angle(grab, cursor)).exp()
172    }
173}
174
175#[derive(Clone, Debug)]
176pub struct RingStyle {
177    pub segments: usize,
178    pub width_px: f32,
179    /// RGBA per plane, in [`Plane4::ALL`] order.
180    pub colors: [[f32; 4]; 6],
181}
182
183impl RingStyle {
184    pub const PAIRED_HUE_COLORS: [[f32; 4]; 6] = [
185        [0.95, 0.30, 0.32, 1.0],
186        [0.35, 0.85, 0.40, 1.0],
187        [0.30, 0.42, 0.75, 1.0],
188        [0.38, 0.62, 0.98, 1.0],
189        [0.24, 0.55, 0.30, 1.0],
190        [0.70, 0.25, 0.28, 1.0],
191    ];
192}
193
194impl Default for RingStyle {
195    fn default() -> Self {
196        Self {
197            segments: 96,
198            width_px: 2.5,
199            colors: Self::PAIRED_HUE_COLORS,
200        }
201    }
202}
203
204#[cfg(test)]
205mod tests {
206    use super::*;
207    use loam_math::Rotor;
208    use std::f32::consts::{PI, TAU};
209
210    const WIDGET: Hypergimbal = Hypergimbal {
211        center: Vec3::new(0.3, -0.7, 1.1),
212        scale: 0.8,
213    };
214
215    fn great_circle_point(plane: Plane4, theta: f32) -> Vec4 {
216        let (a, b) = plane_axes(plane);
217        a * theta.cos() + b * theta.sin()
218    }
219
220    fn assert_close(actual: f32, expected: f32, tol: f32, what: &str) {
221        assert!(
222            (actual - expected).abs() < tol,
223            "{what}: {actual} != {expected}"
224        );
225    }
226
227    #[test]
228    fn ring_is_the_projected_great_circle() {
229        for plane in Plane4::ALL {
230            let ring = WIDGET.ring(plane);
231            for step in 0..24 {
232                let theta = step as f32 / 24.0 * TAU;
233                let projected = WIDGET
234                    .project(great_circle_point(plane, theta))
235                    .expect("no coordinate great circle passes through the pole");
236                assert_close(
237                    (projected - ring.center).length(),
238                    ring.radius,
239                    2e-4,
240                    "off the circle",
241                );
242                let chi = ring.ring_angle(projected);
243                assert!(
244                    (ring.point(chi) - projected).length() < 2e-4,
245                    "{plane:?} θ={theta}: ring.point(χ) != σ(p(θ))"
246                );
247            }
248        }
249    }
250
251    #[test]
252    fn arc_angle_inverts_the_projection_along_each_ring() {
253        for plane in Plane4::ALL {
254            let ring = WIDGET.ring(plane);
255            for step in 0..37 {
256                let theta = -PI + step as f32 / 37.0 * TAU;
257                let projected = WIDGET.project(great_circle_point(plane, theta)).unwrap();
258                let recovered = ring.arc_angle(ring.ring_angle(projected));
259                assert_close(
260                    wrap_pi(recovered - theta),
261                    0.0,
262                    2e-4,
263                    &format!("{plane:?} arc angle at θ={theta}"),
264                );
265            }
266        }
267    }
268
269    #[test]
270    fn drag_yields_the_rotation_its_plane_predicts() {
271        for plane in Plane4::ALL {
272            let ring = WIDGET.ring(plane);
273            for grab_step in 0..8 {
274                let theta = -PI + grab_step as f32 / 8.0 * TAU;
275                for delta in [-1.1_f32, -0.4, 0.05, 0.9, 2.2] {
276                    let start = great_circle_point(plane, theta);
277                    let end = great_circle_point(plane, theta + delta);
278                    let grab = WIDGET.project(start).unwrap();
279                    let cursor = WIDGET.project(end).unwrap();
280
281                    assert_close(
282                        wrap_pi(ring.drag_angle(grab, cursor) - delta),
283                        0.0,
284                        1e-3,
285                        &format!("{plane:?} drag angle from θ={theta} by {delta}"),
286                    );
287
288                    let rotor = ring.drag_rotor(grab, cursor);
289                    assert!(
290                        (rotor.apply(start) - end).length() < 2e-3,
291                        "{plane:?}: rotor does not carry the grabbed point to the cursor"
292                    );
293                }
294            }
295        }
296    }
297
298    #[test]
299    fn drag_across_the_seam_takes_the_short_way() {
300        for plane in Plane4::ALL {
301            let ring = WIDGET.ring(plane);
302            let grab = ring.point(PI - 0.05);
303            let cursor = ring.point(-PI + 0.05);
304            let angle = ring.drag_angle(grab, cursor);
305            assert!(
306                angle.abs() < 0.6,
307                "{plane:?}: seam crossing gave {angle} rad"
308            );
309        }
310    }
311
312    #[test]
313    fn grazing_and_backward_rays_miss_the_ring_plane() {
314        let ring = WIDGET.ring(Plane4::Xy);
315        let eye = ring.center + ring.normal() * 3.0;
316        assert!(ring.ray_hit(eye, -ring.normal()).is_some());
317        assert!(
318            ring.ray_hit(eye, ring.normal()).is_none(),
319            "plane behind the eye must not hit"
320        );
321        assert!(
322            ring.ray_hit(eye, ring.u).is_none(),
323            "parallel ray must miss"
324        );
325        let grazing = (ring.u + ring.normal() * (MIN_PLANE_INCIDENCE * 0.5)).normalize();
326        assert!(
327            ring.ray_hit(eye, grazing).is_none(),
328            "grazing ray must miss"
329        );
330    }
331
332    #[test]
333    fn line_mesh_is_well_formed_and_on_the_rings() {
334        let style = RingStyle {
335            segments: 12,
336            ..RingStyle::default()
337        };
338        let mut mesh = LineMesh::<3>::default();
339        mesh.segments.push(([0.0; 3], [1.0; 3]));
340        mesh.colors.push(([0.0; 4], [0.0; 4]));
341        mesh.widths.push(1.0);
342
343        WIDGET.append_line_mesh(&style, &mut mesh);
344        assert_eq!(mesh.segments.len(), 1 + 6 * style.segments);
345        assert_eq!(mesh.colors.len(), mesh.segments.len());
346        assert_eq!(mesh.widths.len(), mesh.segments.len());
347
348        for (plane, edges) in Plane4::ALL
349            .into_iter()
350            .zip(mesh.segments[1..].chunks_exact(style.segments))
351        {
352            let ring = WIDGET.ring(plane);
353            for (start, end) in edges {
354                for point in [*start, *end] {
355                    let offset = Vec3::from_array(point) - ring.center;
356                    assert!(offset.dot(ring.normal()).abs() < 1e-4);
357                    assert!((offset.length() - ring.radius).abs() < 1e-4);
358                }
359            }
360        }
361        let first = mesh.segments[1].0;
362        let last = mesh.segments[style.segments].1;
363        assert!((Vec3::from_array(first) - Vec3::from_array(last)).length() < 1e-4);
364    }
365
366    #[test]
367    fn the_pole_has_no_image() {
368        assert!(WIDGET.project(POLE).is_none());
369        assert!(WIDGET.project(POLE * 0.99999).is_none());
370        assert!(WIDGET.project(-POLE).is_some());
371    }
372}
loam_render/hypergimbal.rs

loam_render/
hypergimbal.rs
