loam_math/
hyperbolic.rs

1//! Poincaré-ball points satisfy `|p| < 1`; isometries act on the Lorentz embedding.
2
3use std::borrow::Cow;
4
5use glam::{Mat4, Quat, Vec3, Vec4};
6use serde::{Deserialize, Serialize};
7
8use crate::space::{IsometryGroup, Space, WgslAccuracy, WgslSpace};
9
10pub const POINCARE_R2_MAX: f32 = 1.0 - 1e-7;
11
12pub const H3_MAX_ARC: f32 = 17.5;
13
14/// H³ hits at least `H3_DEPTH_SEPARATION` apart stay ordered within this hyperbolic distance of an eye at most `H3_EYE_CHART_REACH` from the chart origin, and `valid_point` refuses points farther than this from the origin.
15pub const H3_DEPTH_ENVELOPE: f32 = 6.0;
16pub const H3_DEPTH_SEPARATION: f32 = 0.05;
17pub const H3_EYE_CHART_REACH: f32 = 1.0;
18
19/// False past the radius where the chart's f32 arithmetic is clamped.
20pub fn in_poincare_ball(p: Vec3) -> bool {
21    p.is_finite() && p.length_squared() <= POINCARE_R2_MAX
22}
23
24fn clamp_to_ball(p: Vec3) -> Vec3 {
25    let r2 = p.length_squared();
26    if r2 <= POINCARE_R2_MAX {
27        return p;
28    }
29    #[cfg(debug_assertions)]
30    tracing::warn!("HyperbolicH3: point outside Poincaré ball clamped (|p|²={r2:.4})");
31    let mut q = p * (POINCARE_R2_MAX.sqrt() / r2.sqrt());
32    while q.length_squared() > POINCARE_R2_MAX {
33        q *= 1.0 - f32::EPSILON;
34    }
35    q
36}
37
38/// SO⁺(3,1) acting on `(x, y, z, w)`, with `w` time-like.
39#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
40pub struct Iso3H {
41    pub matrix: Mat4,
42}
43
44impl Iso3H {
45    pub const IDENTITY: Self = Self {
46        matrix: Mat4::IDENTITY,
47    };
48
49    /// Fixes the ball origin.
50    pub fn from_rotation(rotation: Quat) -> Self {
51        let r = glam::Mat3::from_quat(rotation);
52        let c0 = r.col(0);
53        let c1 = r.col(1);
54        let c2 = r.col(2);
55        Self {
56            matrix: Mat4::from_cols(
57                Vec4::new(c0.x, c0.y, c0.z, 0.0),
58                Vec4::new(c1.x, c1.y, c1.z, 0.0),
59                Vec4::new(c2.x, c2.y, c2.z, 0.0),
60                Vec4::new(0.0, 0.0, 0.0, 1.0),
61            ),
62        }
63    }
64
65    /// Maps the ball origin to `target`, with rapidity capped at the chart boundary.
66    pub fn from_translation(target: Vec3) -> Self {
67        let r2 = target.length_squared();
68        if r2 < 1e-14 {
69            return Self::IDENTITY;
70        }
71        let r = r2.sqrt();
72        let dir = target / r;
73
74        let rapidity = 2.0 * artanh(r.min(POINCARE_R2_MAX.sqrt()));
75        let ch = rapidity.cosh();
76        let sh = rapidity.sinh();
77        let k = ch - 1.0;
78        let (dx, dy, dz) = (dir.x, dir.y, dir.z);
79        Self {
80            matrix: Mat4::from_cols(
81                Vec4::new(1.0 + k * dx * dx, k * dx * dy, k * dx * dz, sh * dx),
82                Vec4::new(k * dy * dx, 1.0 + k * dy * dy, k * dy * dz, sh * dy),
83                Vec4::new(k * dz * dx, k * dz * dy, 1.0 + k * dz * dz, sh * dz),
84                Vec4::new(sh * dx, sh * dy, sh * dz, ch),
85            ),
86        }
87    }
88}
89
90#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
91pub struct HyperbolicH3;
92
93impl Space for HyperbolicH3 {
94    type Point = Vec3;
95    type Vector = Vec3;
96    type Frame = Iso3H;
97
98    fn frame_at(&self, at: Vec3) -> Iso3H {
99        Iso3H::from_translation(at)
100    }
101
102    fn distance(&self, a: Vec3, b: Vec3) -> f32 {
103        let a = clamp_to_ball(a);
104        let b = clamp_to_ball(b);
105        let n = mobius_add(-a, b).length();
106        2.0 * artanh(n.min(POINCARE_R2_MAX.sqrt()))
107    }
108
109    fn exp(&self, at: Vec3, v: Vec3) -> Vec3 {
110        let n = v.length();
111        if n < 1e-7 {
112            return at;
113        }
114        let at = clamp_to_ball(at);
115        let lambda = 2.0 / (1.0 - at.length_squared());
116        let dir = v / n;
117        let scale = (lambda * n * 0.5).tanh();
118        mobius_add(at, scale * dir)
119    }
120
121    fn log(&self, from: Vec3, to: Vec3) -> Vec3 {
122        let from = clamp_to_ball(from);
123        let to = clamp_to_ball(to);
124        let d = mobius_add(-from, to);
125        let n = d.length();
126        if n < 1e-7 {
127            return Vec3::ZERO;
128        }
129        let lambda = 2.0 / (1.0 - from.length_squared());
130        let mag = (2.0 / lambda) * artanh(n.min(POINCARE_R2_MAX.sqrt()));
131        mag * d / n
132    }
133
134    fn parallel_transport(&self, from: Vec3, to: Vec3, v: Vec3) -> Vec3 {
135        let from = clamp_to_ball(from);
136        let to = clamp_to_ball(to);
137        let conformal = (1.0 - to.length_squared()) / (1.0 - from.length_squared());
138        conformal * gyr_apply(to, -from, v)
139    }
140
141    fn chart_envelope(&self) -> f32 {
142        H3_DEPTH_ENVELOPE
143    }
144
145    fn valid_point(&self, p: Vec3) -> bool {
146        in_poincare_ball(p) && self.distance(Vec3::ZERO, p) <= H3_DEPTH_ENVELOPE
147    }
148}
149
150impl IsometryGroup for HyperbolicH3 {
151    type Iso = Iso3H;
152
153    fn iso_identity(&self) -> Iso3H {
154        Iso3H::IDENTITY
155    }
156
157    fn iso_compose(&self, a: Iso3H, b: Iso3H) -> Iso3H {
158        Iso3H {
159            matrix: a.matrix * b.matrix,
160        }
161    }
162
163    fn iso_inverse(&self, a: Iso3H) -> Iso3H {
164        let mt = a.matrix.transpose().to_cols_array_2d();
165        let mut out = [[0.0f32; 4]; 4];
166        for col in 0..4 {
167            for row in 0..4 {
168                let sign = if (row == 3) ^ (col == 3) { -1.0 } else { 1.0 };
169                out[col][row] = sign * mt[col][row];
170            }
171        }
172        Iso3H {
173            matrix: Mat4::from_cols_array_2d(&out),
174        }
175    }
176
177    fn iso_apply(&self, iso: Iso3H, p: Vec3) -> Vec3 {
178        let h = poincare_to_hyperboloid(p);
179        let h2 = iso.matrix * h;
180        hyperboloid_to_poincare(h2)
181    }
182
183    fn iso_transport(&self, iso: Iso3H, at: Vec3, v: Vec3) -> Vec3 {
184        let at = clamp_to_ball(at);
185        let h = poincare_to_hyperboloid(at);
186        let dh = poincare_to_hyperboloid_tangent(at, v);
187        hyperboloid_to_poincare_tangent(iso.matrix * h, iso.matrix * dh)
188    }
189}
190
191impl WgslSpace for HyperbolicH3 {
192    fn wgsl_impl(&self) -> Cow<'static, str> {
193        Cow::Owned(format!(
194            "const LOAM_MAX_ARC: f32 = {H3_MAX_ARC:?};{WGSL_IMPL}"
195        ))
196    }
197
198    fn wgsl_accuracy(&self) -> WgslAccuracy {
199        WgslAccuracy::Exact
200    }
201}
202
203const WGSL_IMPL: &str = r#"
204// loam-math :: HyperbolicH3 (v0 Space WGSL ABI)
205const LOAM_H3_R2_MAX: f32 = 0.9999999;
206const LOAM_H3_GYR_N2_MIN: f32 = 1e-20;
207
208fn loam_artanh(x: f32) -> f32 {
209    return 0.5 * log((1.0 + x) / (1.0 - x));
210}
211
212fn loam_clamp_to_ball(p: vec3<f32>) -> vec3<f32> {
213    let r2 = dot(p, p);
214    if (r2 <= LOAM_H3_R2_MAX) {
215        return p;
216    }
217    var q = p * (sqrt(LOAM_H3_R2_MAX) / sqrt(r2));
218    // Normalization can round back onto the singular unit sphere.
219    while (dot(q, q) > LOAM_H3_R2_MAX) {
220        q *= 1.0 - 1.1920929e-7;
221    }
222    return q;
223}
224
225fn loam_mobius_add(a: vec3<f32>, b: vec3<f32>) -> vec3<f32> {
226    let ab = dot(a, b);
227    let aa = dot(a, a);
228    let bb = dot(b, b);
229    let num = (1.0 + 2.0 * ab + bb) * a + (1.0 - aa) * b;
230    let den = 1.0 + 2.0 * ab + aa * bb;
231    if (abs(den) < 1e-12) {
232        return vec3<f32>(0.0, 0.0, 0.0);
233    }
234    return num / den;
235}
236
237fn loam_gyr_apply(a: vec3<f32>, b: vec3<f32>, v: vec3<f32>) -> vec3<f32> {
238    // Ungar, From Möbius to Gyrogroups, 2008, §4, Def. 4.
239    let scalar = 1.0 + dot(a, b);
240    let axis = -cross(a, b);
241    let norm2 = scalar * scalar + dot(axis, axis);
242
243    if (norm2 < LOAM_H3_GYR_N2_MIN) {
244        return v;
245    }
246    return v + (2.0 / norm2) * (scalar * cross(axis, v) + cross(axis, cross(axis, v)));
247}
248
249fn loam_origin_distance(p: vec3<f32>) -> f32 {
250    let r = min(length(p), sqrt(LOAM_H3_R2_MAX));
251    return 2.0 * loam_artanh(r);
252}
253
254fn loam_distance(a: vec3<f32>, b: vec3<f32>) -> f32 {
255
256    let aa = loam_clamp_to_ball(a);
257    let bb = loam_clamp_to_ball(b);
258    let d = loam_mobius_add(-aa, bb);
259    let n = min(length(d), sqrt(LOAM_H3_R2_MAX));
260    return 2.0 * loam_artanh(n);
261}
262
263fn loam_exp(at: vec3<f32>, v: vec3<f32>) -> vec3<f32> {
264    let n = length(v);
265    if (n < 1e-7) { return at; }
266    let p_at = loam_clamp_to_ball(at);
267    let aa = dot(p_at, p_at);
268    let lambda = 2.0 / (1.0 - aa);
269    let dir = v / n;
270    let scale = tanh(lambda * n * 0.5);
271    return loam_clamp_to_ball(loam_mobius_add(p_at, scale * dir));
272}
273
274fn loam_log(p_from: vec3<f32>, p_to: vec3<f32>) -> vec3<f32> {
275    let p_from_clamped = loam_clamp_to_ball(p_from);
276    let p_to_clamped = loam_clamp_to_ball(p_to);
277    let d = loam_mobius_add(-p_from_clamped, p_to_clamped);
278    let n = length(d);
279    if (n < 1e-7) { return vec3<f32>(0.0, 0.0, 0.0); }
280    let aa = dot(p_from_clamped, p_from_clamped);
281    let lambda = 2.0 / (1.0 - aa);
282    let mag = (2.0 / lambda) * loam_artanh(min(n, sqrt(LOAM_H3_R2_MAX)));
283    return mag * d / n;
284}
285
286fn loam_h3_transport(p_from: vec3<f32>, p_to: vec3<f32>, v: vec3<f32>) -> vec3<f32> {
287    let p_from_clamped = loam_clamp_to_ball(p_from);
288    let p_to_clamped = loam_clamp_to_ball(p_to);
289    let conformal = (1.0 - dot(p_to_clamped, p_to_clamped)) / (1.0 - dot(p_from_clamped, p_from_clamped));
290    return conformal * loam_gyr_apply(p_to_clamped, -p_from_clamped, v);
291}
292
293struct LoamGeodesicStep { p: vec3<f32>, v: vec3<f32> }
294
295fn loam_geodesic_step(p: vec3<f32>, v: vec3<f32>, s: f32) -> LoamGeodesicStep {
296    let next = loam_exp(p, v * s);
297    return LoamGeodesicStep(next, loam_h3_transport(p, next, v));
298}
299
300// Cannon, Floyd, Kenyon, Parry, Hyperbolic Geometry, 1997, §7.
301fn loam_poincare_to_hyperboloid(p: vec3<f32>) -> vec4<f32> {
302    let r2 = min(dot(p, p), LOAM_H3_R2_MAX);
303    let den = 1.0 - r2;
304    return vec4<f32>(2.0 * p / den, (1.0 + r2) / den);
305}
306
307fn loam_hyperboloid_to_klein(h: vec4<f32>) -> vec3<f32> {
308    return h.xyz / max(h.w, 1.0);
309}
310"#;
311
312fn artanh(x: f32) -> f32 {
313    0.5 * ((1.0 + x) / (1.0 - x)).ln()
314}
315
316// Ungar, From Möbius to Gyrogroups, 2008, §4, Def. 3.
317fn mobius_add(a: Vec3, b: Vec3) -> Vec3 {
318    let ab = a.dot(b);
319    let aa = a.length_squared();
320    let bb = b.length_squared();
321    let num = (1.0 + 2.0 * ab + bb) * a + (1.0 - aa) * b;
322    let den = 1.0 + 2.0 * ab + aa * bb;
323    if den.abs() < 1e-12 {
324        Vec3::ZERO
325    } else {
326        num / den
327    }
328}
329
330// Ungar, From Möbius to Gyrogroups, 2008, §4, Def. 4.
331fn gyr_apply(a: Vec3, b: Vec3, v: Vec3) -> Vec3 {
332    let scalar = 1.0 + a.dot(b);
333    let axis = -a.cross(b);
334    let norm2 = scalar * scalar + axis.length_squared();
335    v + (2.0 / norm2) * (scalar * axis.cross(v) + axis.cross(axis.cross(v)))
336}
337
338// Cannon, Floyd, Kenyon, Parry, Hyperbolic Geometry, 1997, §7: the Lorentz embedding of the ball.
339pub fn poincare_to_hyperboloid(p: Vec3) -> Vec4 {
340    let r2 = p.length_squared().min(POINCARE_R2_MAX);
341    let den = 1.0 - r2;
342    Vec4::new(
343        2.0 * p.x / den,
344        2.0 * p.y / den,
345        2.0 * p.z / den,
346        (1.0 + r2) / den,
347    )
348}
349
350pub fn hyperboloid_to_poincare(h: Vec4) -> Vec3 {
351    let den = (1.0 + h.w).max(1e-7);
352    Vec3::new(h.x / den, h.y / den, h.z / den)
353}
354
355// Cannon, Floyd, Kenyon, Parry, Hyperbolic Geometry, 1997, §7: Klein is the hyperboloid seen from the origin.
356pub fn hyperboloid_to_klein(h: Vec4) -> Vec3 {
357    Vec3::new(h.x, h.y, h.z) / h.w.max(1.0)
358}
359
360pub fn klein_to_hyperboloid(k: Vec3) -> Vec4 {
361    let w = 1.0 / (1.0 - k.length_squared().min(POINCARE_R2_MAX)).sqrt();
362    Vec4::new(k.x * w, k.y * w, k.z * w, w)
363}
364
365pub fn poincare_to_klein(p: Vec3) -> Vec3 {
366    hyperboloid_to_klein(poincare_to_hyperboloid(p))
367}
368
369pub fn klein_to_poincare(k: Vec3) -> Vec3 {
370    hyperboloid_to_poincare(klein_to_hyperboloid(k))
371}
372
373/// Klein point at hyperbolic distance `t` from the origin along unit `dir`; geodesics from the origin are straight.
374pub fn klein_ray(dir: Vec3, t: f32) -> Vec3 {
375    dir * t.tanh()
376}
377
378fn poincare_to_hyperboloid_tangent(p: Vec3, v: Vec3) -> Vec4 {
379    let r2 = p.length_squared().min(POINCARE_R2_MAX);
380    let den = 1.0 - r2;
381    let radial = 4.0 * p.dot(v) / (den * den);
382    let space = (2.0 / den) * v + radial * p;
383    Vec4::new(space.x, space.y, space.z, radial)
384}
385
386fn hyperboloid_to_poincare_tangent(h: Vec4, dh: Vec4) -> Vec3 {
387    let den = (1.0 + h.w).max(1e-7);
388    let space = Vec3::new(h.x, h.y, h.z);
389    let d_space = Vec3::new(dh.x, dh.y, dh.z);
390    d_space / den - space * (dh.w / (den * den))
391}
392
393#[cfg(test)]
394mod tests {
395    use super::*;
396
397    #[test]
398    fn the_arc_cap_stays_inside_the_saturating_origin_distance() {
399        let farthest = HyperbolicH3.distance(Vec3::ZERO, Vec3::X * POINCARE_R2_MAX.sqrt());
400        assert!(
401            H3_MAX_ARC * 0.92 < farthest,
402            "the escape guard at {} can never fire below {farthest}",
403            H3_MAX_ARC * 0.92
404        );
405        assert!(HyperbolicH3.wgsl_impl().contains(&format!("{H3_MAX_ARC}")));
406    }
407
408    use approx::assert_relative_eq;
409
410    fn h3() -> HyperbolicH3 {
411        HyperbolicH3
412    }
413
414    fn lambda(p: Vec3) -> f32 {
415        2.0 / (1.0 - p.length_squared())
416    }
417
418    #[test]
419    fn iso_translation_moves_origin_to_target() {
420        let s = h3();
421        let target = Vec3::new(0.3, -0.1, 0.2);
422        let iso = Iso3H::from_translation(target);
423        let moved = s.iso_apply(iso, Vec3::ZERO);
424        assert_relative_eq!(moved.x, target.x, epsilon = 1e-5);
425        assert_relative_eq!(moved.y, target.y, epsilon = 1e-5);
426        assert_relative_eq!(moved.z, target.z, epsilon = 1e-5);
427    }
428
429    #[test]
430    fn small_scale_distance_matches_euclidean_via_metric_factor() {
431        let s = h3();
432        let eps = 1e-3;
433        let p = Vec3::new(eps, 0.0, 0.0);
434        assert_relative_eq!(s.distance(Vec3::ZERO, p), 2.0 * eps, epsilon = 1e-6);
435    }
436
437    #[test]
438    fn angle_defect_in_small_triangle_scales_with_area() {
439        let s = h3();
440        let l = 0.05;
441        let v_norm = l * 0.5;
442        let a = Vec3::ZERO;
443        let b = s.exp(a, Vec3::new(v_norm, 0.0, 0.0));
444        let c = s.exp(
445            a,
446            Vec3::new(v_norm * 0.5, v_norm * (3.0_f32).sqrt() * 0.5, 0.0),
447        );
448
449        let angle_at = |p: Vec3, q: Vec3, r: Vec3| -> f32 {
450            let u = s.log(p, q);
451            let w = s.log(p, r);
452            (u.dot(w) / (u.length() * w.length()))
453                .clamp(-1.0, 1.0)
454                .acos()
455        };
456
457        let alpha = angle_at(a, b, c);
458        let beta = angle_at(b, a, c);
459        let gamma = angle_at(c, a, b);
460        let defect = std::f32::consts::PI - (alpha + beta + gamma);
461        let expected_area = (3.0_f32.sqrt() / 4.0) * l * l;
462
463        assert!(
464            defect > 0.0,
465            "hyperbolic triangle should have positive angle defect, got {defect}"
466        );
467        assert_relative_eq!(defect, expected_area, epsilon = 5e-4);
468    }
469
470    fn ball_sweep() -> Vec<Vec3> {
471        let radii = [
472            0.0f32, 0.1, 0.3, 0.5, 0.7, 0.8, 0.9, 0.95, 0.99, 0.999, 0.9999,
473        ];
474        let directions = [
475            Vec3::new(1.0, 0.0, 0.0),
476            Vec3::new(0.0, 1.0, 0.0),
477            Vec3::new(-1.0, 0.0, 0.0),
478            Vec3::new(0.577, 0.577, 0.577),
479            Vec3::new(0.3, -0.9, 0.2),
480        ];
481        let mut out = Vec::new();
482        for r in radii {
483            for d in directions {
484                out.push(d.normalize_or_zero() * r);
485            }
486        }
487        out
488    }
489
490    const SWEEP_TANGENTS: [Vec3; 3] = [
491        Vec3::new(0.06, 0.0, 0.0),
492        Vec3::new(0.0, 0.05, 0.02),
493        Vec3::new(2.0, -1.0, 0.5),
494    ];
495
496    const TWIST_CHORD_TOLERANCE: f32 = 1e-3;
497
498    #[test]
499    fn gyration_matches_ungars_four_addition_definition() {
500        let four_addition = |a: Vec3, b: Vec3, v: Vec3| {
501            let ab = mobius_add(a, b);
502            let bv = mobius_add(b, v);
503            mobius_add(-ab, mobius_add(a, bv))
504        };
505        let cases = [
506            (
507                Vec3::new(0.2, 0.1, -0.05),
508                Vec3::new(-0.1, 0.25, 0.05),
509                Vec3::new(0.03, -0.02, 0.04),
510            ),
511            (
512                Vec3::new(0.3, 0.0, 0.0),
513                Vec3::new(0.0, 0.3, 0.0),
514                Vec3::new(0.0, 0.0, 0.05),
515            ),
516            (
517                Vec3::new(0.05, -0.4, 0.1),
518                Vec3::new(0.15, 0.05, -0.3),
519                Vec3::new(0.06, 0.01, -0.02),
520            ),
521        ];
522        for (a, b, v) in cases {
523            let residual = (gyr_apply(a, b, v) - four_addition(a, b, v)).length();
524
525            assert!(
526                residual <= 1e-5 * v.length(),
527                "gyr[{a:?}, {b:?}] disagrees with its definition by {residual}"
528            );
529        }
530    }
531
532    #[test]
533    fn parallel_transport_is_exact_where_the_clamp_conditions_the_gyration_worst() {
534        let s = h3();
535        let outside = Vec3::new(2.0, 0.0, 0.0);
536        let shell = clamp_to_ball(outside);
537        let scalar = 1.0 + shell.dot(-shell);
538        let bound = (1.0 - POINCARE_R2_MAX) * (1.0 - POINCARE_R2_MAX);
539        assert!(
540            (scalar * scalar - bound).abs() <= 1e-3 * bound,
541            "the clamped shell puts the gyration norm² at {}, not at the \
542             documented bound {bound}",
543            scalar * scalar
544        );
545        for v in SWEEP_TANGENTS {
546            assert_eq!(
547                s.parallel_transport(outside, outside, v),
548                v,
549                "transport at the gyration's conditioning floor moved {v:?}"
550            );
551        }
552    }
553
554    #[test]
555    fn parallel_transport_preserves_the_metric_norm_across_the_whole_ball() {
556        let s = h3();
557        let points = ball_sweep();
558        for &a in &points {
559            for &b in &points {
560                for v in SWEEP_TANGENTS {
561                    let before = lambda(a) * v.length();
562                    let after = lambda(b) * s.parallel_transport(a, b, v).length();
563                    assert!(
564                        (after - before).abs() <= 1e-5 * before,
565                        "transport {a:?} -> {b:?} took the norm of {v:?} \
566                         from {before} to {after}"
567                    );
568                }
569            }
570        }
571    }
572
573    #[test]
574    fn parallel_transport_carries_a_geodesic_tangent_along_its_own_geodesic() {
575        let s = h3();
576        let points = ball_sweep();
577        for &a in &points {
578            for &b in &points {
579                if s.distance(a, b) < 1e-3 {
580                    continue;
581                }
582                let forward = (-s.log(b, a)).normalize();
583                let transported = s.parallel_transport(a, b, s.log(a, b)).normalize();
584
585                assert!(
586                    (transported - forward).length() <= 1e-3,
587                    "transported direction {transported:?} misses the forward \
588                     tangent {forward:?} at {a:?} -> {b:?}"
589                );
590            }
591        }
592    }
593
594    #[test]
595    fn parallel_transport_fixes_the_normal_of_the_geodesic_plane() {
596        let s = h3();
597        let points = ball_sweep();
598        for &a in &points {
599            for &b in &points {
600                let spread = a.normalize_or_zero().cross(b.normalize_or_zero());
601                if spread.length() < 1e-3 {
602                    continue;
603                }
604                let normal = a.cross(b).normalize();
605                let transported = s.parallel_transport(a, b, normal).normalize();
606                assert!(
607                    (transported - normal).length() <= TWIST_CHORD_TOLERANCE,
608                    "transport {a:?} -> {b:?} tilted the plane normal \
609                     {normal:?} to {transported:?}"
610                );
611            }
612        }
613    }
614
615    #[test]
616    fn parallel_transport_along_a_radial_geodesic_fixes_the_orthogonal_plane() {
617        let s = h3();
618        let points = ball_sweep();
619        for &a in &points {
620            for &b in &points {
621                let (ua, ub) = (a.normalize_or_zero(), b.normalize_or_zero());
622                if ua.cross(ub).length() >= 1e-3 || s.distance(a, b) < 1e-3 {
623                    continue;
624                }
625
626                let line = if ua.length_squared() < 0.5 { ub } else { ua };
627                let (w0, w1) = line.any_orthonormal_pair();
628                for w in [w0, w1] {
629                    let transported = s.parallel_transport(a, b, w).normalize();
630                    assert!(
631                        (transported - w).length() <= TWIST_CHORD_TOLERANCE,
632                        "radial transport {a:?} -> {b:?} spun {w:?} to \
633                         {transported:?}"
634                    );
635                }
636            }
637        }
638    }
639
640    #[test]
641    fn parallel_transport_preserves_tangent_orientation() {
642        let s = h3();
643        let points = ball_sweep();
644        for &a in &points {
645            for &b in &points {
646                let frame = glam::Mat3::from_cols(
647                    s.parallel_transport(a, b, Vec3::X),
648                    s.parallel_transport(a, b, Vec3::Y),
649                    s.parallel_transport(a, b, Vec3::Z),
650                );
651                let expected = ((1.0 - b.length_squared()) / (1.0 - a.length_squared())).powi(3);
652                assert!(
653                    (frame.determinant() - expected).abs() <= 1e-4 * expected,
654                    "transport {a:?} -> {b:?} has frame determinant {} against \
655                     the conformal ratio cubed {expected}",
656                    frame.determinant()
657                );
658            }
659        }
660    }
661
662    #[test]
663    fn iso_transport_norm_error_stays_within_the_conformal_factor() {
664        let s = h3();
665        let isos = [
666            Iso3H::from_translation(Vec3::new(0.15, 0.0, 0.0)),
667            Iso3H::from_rotation(Quat::from_rotation_z(0.4)),
668            Iso3H::from_translation(Vec3::new(-0.05, 0.2, 0.1)),
669            Iso3H::from_translation(Vec3::new(0.7, -0.2, 0.1)),
670        ];
671        for at in ball_sweep() {
672            let budget = 16.0 * lambda(at) * f32::EPSILON;
673            for iso in isos {
674                let moved = s.iso_apply(iso, at);
675                for v in SWEEP_TANGENTS {
676                    let before = lambda(at) * v.length();
677                    let after = lambda(moved) * s.iso_transport(iso, at, v).length();
678                    assert!(
679                        (after - before).abs() <= budget * before,
680                        "iso_transport at {at:?} took the norm of {v:?} \
681                         from {before} to {after}, past the {budget} budget"
682                    );
683                }
684            }
685        }
686    }
687
688    #[test]
689    fn a_point_past_the_declared_envelope_is_refused() {
690        let s = h3();
691        let at = |d: f32| Vec3::X * (0.5 * d).tanh();
692        assert!(s.valid_point(at(H3_DEPTH_ENVELOPE - 1.0)));
693        assert!(!s.valid_point(at(H3_DEPTH_ENVELOPE + 1.0)));
694        assert!(!s.valid_point(Vec3::X));
695        assert!(!s.valid_point(Vec3::splat(f32::NAN)));
696    }
697
698    #[test]
699    fn klein_chart_through_the_embedding_matches_its_closed_form() {
700        for p in ball_sweep() {
701            let expected = 2.0 * p / (1.0 + p.length_squared());
702            let k = poincare_to_klein(p);
703            assert!(
704                (k - expected).length() <= 1e-6,
705                "Klein of {p:?} is {k:?}, closed form {expected:?}"
706            );
707        }
708    }
709
710    #[test]
711    fn klein_ray_parameter_is_hyperbolic_arc_length() {
712        let s = h3();
713        let dir = Vec3::new(0.6, -0.48, 0.64);
714        for t in [0.1_f32, 1.0, 3.0, 6.0] {
715            let p = klein_to_poincare(klein_ray(dir, t));
716            assert_relative_eq!(s.distance(Vec3::ZERO, p), t, epsilon = 2e-4 * t);
717        }
718    }
719
720    #[test]
721    fn geodesics_are_straight_in_the_klein_chart() {
722        let s = h3();
723        let at = Vec3::new(0.3, -0.2, 0.5);
724        let v = Vec3::new(0.02, 0.05, -0.03);
725        let k0 = poincare_to_klein(at);
726        let chord = (poincare_to_klein(s.exp(at, 20.0 * v)) - k0).normalize();
727        for i in 1..20 {
728            let k = poincare_to_klein(s.exp(at, i as f32 * v));
729            let off = (k - k0).cross(chord).length();
730            assert!(
731                off <= 1e-5,
732                "geodesic sample {i} leaves the Klein chord by {off}"
733            );
734        }
735    }
736
737    #[test]
738    fn poincare_hyperboloid_round_trip() {
739        let p = Vec3::new(0.2, -0.3, 0.1);
740        let h = poincare_to_hyperboloid(p);
741
742        let lorentz = -h.x * h.x - h.y * h.y - h.z * h.z + h.w * h.w;
743        assert_relative_eq!(lorentz, 1.0, epsilon = 1e-5);
744        let p2 = hyperboloid_to_poincare(h);
745        assert_relative_eq!(p2.x, p.x, epsilon = 1e-6);
746        assert_relative_eq!(p2.y, p.y, epsilon = 1e-6);
747        assert_relative_eq!(p2.z, p.z, epsilon = 1e-6);
748    }
749}
loam_math/hyperbolic.rs

loam_math/
hyperbolic.rs
