loam_math/
rasterizable.rs

1use glam::{Vec3, Vec4};
2
3use crate::space::Space;
4use crate::{EuclideanR3, EuclideanR4};
5
6const PROJECTION_DENOM_EPSILON: f32 = 1e-4;
7
8/// Floor on `1 - dot(p, pole)` at the projection singularity.
9pub const STEREOGRAPHIC_POLE_EPSILON: f32 = 1e-4;
10
11// do Carmo, Differential Geometry of Curves and Surfaces, 1976, §1.4.
12fn perp_frame(n: Vec4) -> (Vec4, Vec4, Vec4) {
13    let ax = n.x.abs();
14    let ay = n.y.abs();
15    let az = n.z.abs();
16    let aw = n.w.abs();
17    let drop = ax.max(ay).max(az).max(aw);
18    let mut seeds = [Vec4::X, Vec4::Y, Vec4::Z, Vec4::W];
19    let drop_idx = if drop == ax {
20        0
21    } else if drop == ay {
22        1
23    } else if drop == az {
24        2
25    } else {
26        3
27    };
28    seeds[drop_idx] = Vec4::ZERO;
29
30    let mut basis = [Vec4::ZERO; 3];
31    let mut count = 0usize;
32    for s in seeds {
33        if s == Vec4::ZERO {
34            continue;
35        }
36        let mut v = s - s.dot(n) * n;
37        for b in basis.iter().take(count) {
38            v -= v.dot(*b) * *b;
39        }
40        basis[count] = v.normalize();
41        count += 1;
42    }
43    (basis[0], basis[1], basis[2])
44}
45
46// Wikipedia, Stereographic projection.
47pub(crate) fn stereographic_to_r3(p: Vec4, pole: Vec4) -> Vec3 {
48    let dot = p.dot(pole).clamp(-1.0, 1.0);
49    let denom = (1.0 - dot).max(STEREOGRAPHIC_POLE_EPSILON);
50    if pole == Vec4::W {
51        return Vec3::new(p.x, p.y, p.z) / denom;
52    }
53    let perp = p - dot * pole;
54    let scaled = perp / denom;
55    let (e1, e2, e3) = perp_frame(pole);
56    Vec3::new(scaled.dot(e1), scaled.dot(e2), scaled.dot(e3))
57}
58
59/// Unsupported dimension/variant combinations return `None`.
60#[derive(Clone, Copy, Debug, Default, PartialEq)]
61pub enum Projection<const N: usize> {
62    /// Keeps the first three components, padding with zero if needed.
63    #[default]
64    Identity,
65
66    /// Drops the indexed axis; an out-of-range index returns `None`.
67    Orthographic { drop_axis: usize },
68
69    /// Pinhole at `(0, 0, 0, focal_distance)`; a vertex at or beyond the eye in w is refused with `None`.
70    Perspective4D { focal_distance: f32 },
71
72    // Coxeter, Regular Polytopes, 3rd ed., ch. 13.
73    /// Projects onto the cell plane; the normal must be outward and unit, and the eye outside the polytope.
74    Schlegel {
75        cell_normal: Vec4,
76        cell_offset: f32,
77        viewpoint_distance: f32,
78        basis: [Vec4; 3],
79    },
80
81    // Wikipedia, Stereographic projection.
82    /// Projects unit S³ from a unit pole into its perpendicular 3-flat.
83    Stereographic {
84        /// Unit pole; `Vec4::W` gives the usual xyz readout.
85        pole: Vec4,
86    },
87}
88
89impl Projection<4> {
90    pub fn schlegel(cell_normal: Vec4, cell_offset: f32, viewpoint_distance: f32) -> Projection<4> {
91        let (e1, e2, e3) = perp_frame(cell_normal);
92        Self::schlegel_with_basis(cell_normal, cell_offset, viewpoint_distance, [e1, e2, e3])
93    }
94
95    pub fn schlegel_with_basis(
96        cell_normal: Vec4,
97        cell_offset: f32,
98        viewpoint_distance: f32,
99        basis: [Vec4; 3],
100    ) -> Projection<4> {
101        Projection::Schlegel {
102            cell_normal,
103            cell_offset,
104            viewpoint_distance,
105            basis,
106        }
107    }
108}
109
110pub trait RasterizableSpace<const N: usize>: Space {
111    fn point_to_array(p: Self::Point) -> [f32; N];
112
113    fn array_to_point(arr: [f32; N]) -> Self::Point;
114
115    fn project_point(point: Self::Point, projection: &Projection<N>) -> Option<Vec3>;
116
117    /// Emits both endpoints and `samples - 1` interior points when `samples > 1`.
118    fn tessellate_segment(
119        p0: Self::Point,
120        p1: Self::Point,
121        samples: usize,
122        emit: impl FnMut(Self::Point),
123    );
124}
125
126impl RasterizableSpace<3> for EuclideanR3 {
127    fn point_to_array(p: Vec3) -> [f32; 3] {
128        p.to_array()
129    }
130
131    fn array_to_point(arr: [f32; 3]) -> Vec3 {
132        Vec3::from_array(arr)
133    }
134
135    fn project_point(point: Vec3, projection: &Projection<3>) -> Option<Vec3> {
136        match projection {
137            Projection::Identity => Some(point),
138            Projection::Orthographic { drop_axis } => match *drop_axis {
139                0 => Some(Vec3::new(point.y, point.z, 0.0)),
140                1 => Some(Vec3::new(point.x, point.z, 0.0)),
141                2 => Some(Vec3::new(point.x, point.y, 0.0)),
142                _ => None,
143            },
144            Projection::Perspective4D { .. }
145            | Projection::Schlegel { .. }
146            | Projection::Stereographic { .. } => None,
147        }
148    }
149
150    fn tessellate_segment(p0: Vec3, p1: Vec3, samples: usize, mut emit: impl FnMut(Vec3)) {
151        emit(p0);
152        for i in 1..samples {
153            let t = i as f32 / samples as f32;
154            emit(p0.lerp(p1, t));
155        }
156        emit(p1);
157    }
158}
159
160impl RasterizableSpace<4> for EuclideanR4 {
161    fn point_to_array(p: Vec4) -> [f32; 4] {
162        p.to_array()
163    }
164
165    fn array_to_point(arr: [f32; 4]) -> Vec4 {
166        Vec4::from_array(arr)
167    }
168
169    fn project_point(point: Vec4, projection: &Projection<4>) -> Option<Vec3> {
170        match projection {
171            Projection::Identity => Some(Vec3::new(point.x, point.y, point.z)),
172            Projection::Orthographic { drop_axis } => match *drop_axis {
173                0 => Some(Vec3::new(point.y, point.z, point.w)),
174                1 => Some(Vec3::new(point.x, point.z, point.w)),
175                2 => Some(Vec3::new(point.x, point.y, point.w)),
176                3 => Some(Vec3::new(point.x, point.y, point.z)),
177                _ => None,
178            },
179            Projection::Perspective4D { focal_distance } => {
180                let denom = focal_distance - point.w;
181                if denom < PROJECTION_DENOM_EPSILON {
182                    return None;
183                }
184                let scale = focal_distance / denom;
185                Some(Vec3::new(point.x, point.y, point.z) * scale)
186            }
187            Projection::Schlegel {
188                cell_normal,
189                cell_offset,
190                viewpoint_distance,
191                basis,
192            } => {
193                // Coxeter, Regular Polytopes, 3rd ed., ch. 13.
194                let n = *cell_normal;
195                let eye = *viewpoint_distance * n;
196                let n_dot_eye = n.dot(eye);
197
198                let denom = n.dot(point) - n_dot_eye;
199                if denom.abs() < PROJECTION_DENOM_EPSILON {
200                    return None;
201                }
202                let t = (*cell_offset - n_dot_eye) / denom;
203                let result = eye + t * (point - eye);
204                let [e1, e2, e3] = *basis;
205                Some(Vec3::new(result.dot(e1), result.dot(e2), result.dot(e3)))
206            }
207            Projection::Stereographic { pole } => {
208                Some(stereographic_to_r3(point.normalize(), *pole))
209            }
210        }
211    }
212
213    fn tessellate_segment(p0: Vec4, p1: Vec4, samples: usize, mut emit: impl FnMut(Vec4)) {
214        emit(p0);
215        for i in 1..samples {
216            let t = i as f32 / samples as f32;
217            emit(p0.lerp(p1, t));
218        }
219        emit(p1);
220    }
221}
222
223#[cfg(test)]
224mod tests {
225    use super::*;
226    use crate::SphericalS3Embedded;
227    use approx::assert_relative_eq;
228
229    const GOLDEN_TIE_FRAME: Vec3 = Vec3::new(0.5, 0.5, std::f32::consts::FRAC_1_SQRT_2);
230
231    fn project4(p: Vec4, projection: &Projection<4>) -> Vec3 {
232        <EuclideanR4 as RasterizableSpace<4>>::project_point(p, projection)
233            .expect("projection should accept this vertex")
234    }
235
236    #[test]
237    fn stereographic_r4_normalizes_scaled_input() {
238        let proj = Projection::Stereographic { pole: Vec4::W };
239        for p in [
240            Vec4::new(0.3, -0.1, 0.2, -0.5).normalize(),
241            Vec4::new(-0.4, 0.6, 0.1, 0.3).normalize(),
242            Vec4::new(0.0, 0.0, 0.0, -1.0),
243        ] {
244            let unit = project4(p, &proj);
245            for k in [0.25_f32, 1.5, 3.0] {
246                let scaled = project4(k * p, &proj);
247                assert!(
248                    scaled.abs_diff_eq(unit, 1e-5),
249                    "scale {k}: {scaled:?} vs {unit:?}"
250                );
251            }
252        }
253    }
254
255    #[test]
256    fn r3_tessellate_one_sample_appends_endpoints() {
257        let p0 = Vec3::new(0.0, 0.0, 0.0);
258        let p1 = Vec3::new(2.0, 4.0, -6.0);
259        let mut out = Vec::new();
260        <EuclideanR3 as RasterizableSpace<3>>::tessellate_segment(p0, p1, 1, |p| out.push(p));
261        assert_eq!(out.len(), 2);
262        assert_eq!(out[0], p0);
263        assert_eq!(out[1], p1);
264    }
265
266    #[test]
267    fn r3_orthographic_drops_named_axis() {
268        let p = Vec3::new(1.0, 2.0, 3.0);
269        let pj = |drop_axis| {
270            <EuclideanR3 as RasterizableSpace<3>>::project_point(
271                p,
272                &Projection::Orthographic { drop_axis },
273            )
274        };
275        assert_eq!(pj(0), Some(Vec3::new(2.0, 3.0, 0.0)));
276        assert_eq!(pj(1), Some(Vec3::new(1.0, 3.0, 0.0)));
277        assert_eq!(pj(2), Some(Vec3::new(1.0, 2.0, 0.0)));
278        assert_eq!(pj(3), None);
279        assert_eq!(pj(99), None);
280    }
281
282    #[test]
283    fn r4_identity_drops_w() {
284        let p = Vec4::new(1.0, 2.0, 3.0, 4.0);
285        let projected = project4(p, &Projection::Identity);
286        assert_eq!(projected, Vec3::new(1.0, 2.0, 3.0));
287    }
288
289    #[test]
290    fn r4_orthographic_drops_each_axis() {
291        let p = Vec4::new(1.0, 2.0, 3.0, 4.0);
292        let pj = |drop_axis| {
293            <EuclideanR4 as RasterizableSpace<4>>::project_point(
294                p,
295                &Projection::Orthographic { drop_axis },
296            )
297        };
298        assert_eq!(pj(0), Some(Vec3::new(2.0, 3.0, 4.0)));
299        assert_eq!(pj(1), Some(Vec3::new(1.0, 3.0, 4.0)));
300        assert_eq!(pj(2), Some(Vec3::new(1.0, 2.0, 4.0)));
301        assert_eq!(pj(3), Some(Vec3::new(1.0, 2.0, 3.0)));
302        assert_eq!(pj(4), None);
303    }
304
305    #[test]
306    fn r4_perspective4d_w_zero_is_unchanged() {
307        let p = Vec4::new(1.0, 2.0, 3.0, 0.0);
308        let proj = Projection::Perspective4D {
309            focal_distance: 2.0,
310        };
311        let got = project4(p, &proj);
312        assert_eq!(got, Vec3::new(1.0, 2.0, 3.0));
313    }
314
315    #[test]
316    fn r4_perspective4d_cube_within_cube_scaling() {
317        let focal = 2.0;
318        let proj = Projection::Perspective4D {
319            focal_distance: focal,
320        };
321        let near = Vec4::new(0.5, 0.5, 0.5, 0.5);
322        let far = Vec4::new(0.5, 0.5, 0.5, -0.5);
323        let pn = project4(near, &proj);
324        let pf = project4(far, &proj);
325
326        let r_near = (pn.length() / 0.5_f32.mul_add(3.0_f32.sqrt(), 0.0)).abs();
327        let r_far = (pf.length() / 0.5_f32.mul_add(3.0_f32.sqrt(), 0.0)).abs();
328        assert!((r_near - 4.0 / 3.0).abs() < 1e-5, "near scale {r_near}");
329        assert!((r_far - 4.0 / 5.0).abs() < 1e-5, "far scale {r_far}");
330        assert!(pn.length() > pf.length(), "near={pn:?} far={pf:?}");
331    }
332
333    #[test]
334    fn r4_perspective4d_at_viewer_refuses() {
335        let p = Vec4::new(0.1, 0.2, 0.3, 2.0);
336        let proj = Projection::Perspective4D {
337            focal_distance: 2.0,
338        };
339        let got = <EuclideanR4 as RasterizableSpace<4>>::project_point(p, &proj);
340        assert_eq!(got, None, "expected a refusal at the viewer, got {got:?}");
341    }
342
343    #[test]
344    fn r3_perspective4d_refuses() {
345        let p = Vec3::new(1.0, 2.0, 3.0);
346        let proj = Projection::Perspective4D {
347            focal_distance: 2.0,
348        };
349        let got = <EuclideanR3 as RasterizableSpace<3>>::project_point(p, &proj);
350        assert_eq!(got, None);
351    }
352
353    const TESSERACT_VERTS: [Vec4; 16] = [
354        Vec4::new(0.5, 0.5, 0.5, 0.5),
355        Vec4::new(-0.5, 0.5, 0.5, 0.5),
356        Vec4::new(0.5, -0.5, 0.5, 0.5),
357        Vec4::new(-0.5, -0.5, 0.5, 0.5),
358        Vec4::new(0.5, 0.5, -0.5, 0.5),
359        Vec4::new(-0.5, 0.5, -0.5, 0.5),
360        Vec4::new(0.5, -0.5, -0.5, 0.5),
361        Vec4::new(-0.5, -0.5, -0.5, 0.5),
362        Vec4::new(0.5, 0.5, 0.5, -0.5),
363        Vec4::new(-0.5, 0.5, 0.5, -0.5),
364        Vec4::new(0.5, -0.5, 0.5, -0.5),
365        Vec4::new(-0.5, -0.5, 0.5, -0.5),
366        Vec4::new(0.5, 0.5, -0.5, -0.5),
367        Vec4::new(-0.5, 0.5, -0.5, -0.5),
368        Vec4::new(0.5, -0.5, -0.5, -0.5),
369        Vec4::new(-0.5, -0.5, -0.5, -0.5),
370    ];
371
372    #[test]
373    fn schlegel_chosen_cell_renders_undistorted() {
374        let cell_offset = 0.5;
375        let proj = Projection::schlegel(Vec4::W, cell_offset, 1.5 * cell_offset);
376        let cell: Vec<Vec4> = TESSERACT_VERTS.iter().take(8).copied().collect();
377        let projected: Vec<Vec3> = cell.iter().map(|v| project4(*v, &proj)).collect();
378
379        for i in 0..cell.len() {
380            for j in (i + 1)..cell.len() {
381                let orig = (Vec3::new(cell[i].x, cell[i].y, cell[i].z)
382                    - Vec3::new(cell[j].x, cell[j].y, cell[j].z))
383                .length();
384                let got = (projected[i] - projected[j]).length();
385                assert!(
386                    (orig - got).abs() < 1e-5,
387                    "chosen-cell distance v{i}-v{j} should be {orig}, got {got}"
388                );
389            }
390        }
391    }
392
393    #[test]
394    fn schlegel_non_axis_aligned_cell_is_not_flattened() {
395        let verts = [
396            Vec4::X,
397            -Vec4::X,
398            Vec4::Y,
399            -Vec4::Y,
400            Vec4::Z,
401            -Vec4::Z,
402            Vec4::W,
403            -Vec4::W,
404        ];
405
406        let centroid = (Vec4::X + Vec4::Y + Vec4::Z + Vec4::W) / 4.0;
407        let cell_offset = centroid.length();
408        let cell_normal = centroid / cell_offset;
409        let proj = Projection::schlegel(cell_normal, cell_offset, 1.5 * cell_offset);
410
411        let boundary = [Vec4::X, Vec4::Y, Vec4::Z, Vec4::W];
412        let inner = [-Vec4::X, -Vec4::Y, -Vec4::Z, -Vec4::W];
413        let proj_pt = |v: Vec4| project4(v, &proj);
414
415        let boundary_r: Vec<f32> = boundary.iter().map(|&v| proj_pt(v).length()).collect();
416        let inner_r: Vec<f32> = inner.iter().map(|&v| proj_pt(v).length()).collect();
417        let r0 = boundary_r[0];
418        assert!(r0 > 1e-3, "boundary must not collapse to the origin");
419        for r in &boundary_r {
420            assert!(
421                (r - r0).abs() < 1e-5,
422                "boundary tetrahedron must be regular, radii {boundary_r:?}"
423            );
424        }
425        for r in &inner_r {
426            assert!(
427                *r < r0 - 1e-3,
428                "every nested vertex must sit inside the boundary, inner {inner_r:?} vs {r0}"
429            );
430        }
431
432        let proj_boundary: Vec<Vec3> = boundary.iter().map(|&v| proj_pt(v)).collect();
433        let edge_len = 2.0_f32.sqrt();
434        for i in 0..proj_boundary.len() {
435            for j in (i + 1)..proj_boundary.len() {
436                let got = (proj_boundary[i] - proj_boundary[j]).length();
437                assert!(
438                    (got - edge_len).abs() < 1e-5,
439                    "oblique boundary edge {i}-{j} should stay {edge_len}, got {got}"
440                );
441            }
442        }
443
444        let all: Vec<Vec3> = verts.iter().map(|&v| proj_pt(v)).collect();
445        for axis in 0..3 {
446            let comp = |p: Vec3| [p.x, p.y, p.z][axis];
447            let spread = all
448                .iter()
449                .map(|&p| comp(p))
450                .fold(f32::NEG_INFINITY, f32::max)
451                - all.iter().map(|&p| comp(p)).fold(f32::INFINITY, f32::min);
452            assert!(
453                spread > 0.5,
454                "axis {axis} should have real spread, got {spread}"
455            );
456        }
457    }
458
459    #[test]
460    fn schlegel_uses_supplied_basis_for_readout() {
461        let p = Vec4::new(0.5, -0.25, 0.125, 0.5);
462        let xyz = Projection::schlegel_with_basis(Vec4::W, 0.5, 0.75, [Vec4::X, Vec4::Y, Vec4::Z]);
463        let yxz = Projection::schlegel_with_basis(Vec4::W, 0.5, 0.75, [Vec4::Y, Vec4::X, Vec4::Z]);
464
465        let a = project4(p, &xyz);
466        let b = project4(p, &yxz);
467
468        assert_eq!(a, Vec3::new(0.5, -0.25, 0.125));
469        assert_eq!(b, Vec3::new(-0.25, 0.5, 0.125));
470    }
471
472    #[test]
473    fn schlegel_zero_denominator_refuses() {
474        let viewpoint_distance = 0.75;
475        let proj = Projection::schlegel(Vec4::W, 0.5, viewpoint_distance);
476
477        let p = Vec4::new(0.3, -0.2, 0.1, viewpoint_distance);
478        let got = <EuclideanR4 as RasterizableSpace<4>>::project_point(p, &proj);
479        assert_eq!(
480            got, None,
481            "a degenerate denominator should refuse, got {got:?}"
482        );
483    }
484
485    // Wikipedia, Stereographic projection.
486    fn stereo_inverse_w_pole(q: Vec3) -> Vec4 {
487        let s = q.length_squared();
488        let w = (s - 1.0) / (s + 1.0);
489        let xyz = q * (1.0 - w);
490        Vec4::new(xyz.x, xyz.y, xyz.z, w)
491    }
492
493    fn stereo_inverse_general(q: Vec3, pole: Vec4) -> Vec4 {
494        let (e1, e2, e3) = perp_frame(pole);
495        let perp = q.x * e1 + q.y * e2 + q.z * e3;
496
497        let s = perp.length_squared();
498        let dot = (s - 1.0) / (s + 1.0);
499        dot * pole + (1.0 - dot) * perp
500    }
501
502    #[test]
503    fn stereographic_default_pole_is_drop_w_of_scaled() {
504        for p in [
505            Vec4::new(0.5, 0.5, 0.5, 0.5),
506            Vec4::new(-0.5, 0.5, 0.5, -0.5),
507            Vec4::new(-0.6, 0.0, 0.8, 0.0),
508        ] {
509            let got = stereographic_to_r3(p, Vec4::W);
510            let want = Vec3::new(p.x, p.y, p.z) / (1.0 - p.w);
511            assert_eq!(
512                got, want,
513                "fast path must match canonical formula for {p:?}"
514            );
515        }
516    }
517
518    #[test]
519    fn stereographic_inverts_off_pole() {
520        let proj_w = Projection::Stereographic { pole: Vec4::W };
521        for p in [
522            Vec4::new(0.2, 0.1, -0.3, 0.4).normalize(),
523            Vec4::new(-0.5, 0.5, 0.5, -0.5),
524            Vec4::new(0.7, -0.2, 0.1, 0.1).normalize(),
525        ] {
526            let img = project4(p, &proj_w);
527            let back = stereo_inverse_w_pole(img);
528            assert_relative_eq!(back.x, p.x, epsilon = 1e-5);
529            assert_relative_eq!(back.y, p.y, epsilon = 1e-5);
530            assert_relative_eq!(back.z, p.z, epsilon = 1e-5);
531            assert_relative_eq!(back.w, p.w, epsilon = 1e-5);
532        }
533        let pole = Vec4::new(0.1, -0.2, 0.3, 0.9).normalize();
534        let proj_n = Projection::Stereographic { pole };
535        for p in [
536            Vec4::new(0.6, 0.5, -0.2, 0.0).normalize(),
537            Vec4::new(-0.3, 0.4, 0.5, -0.2).normalize(),
538        ] {
539            let img = project4(p, &proj_n);
540            let back = stereo_inverse_general(img, pole);
541            assert_relative_eq!(back.x, p.x, epsilon = 1e-5);
542            assert_relative_eq!(back.y, p.y, epsilon = 1e-5);
543            assert_relative_eq!(back.z, p.z, epsilon = 1e-5);
544            assert_relative_eq!(back.w, p.w, epsilon = 1e-5);
545        }
546    }
547
548    #[test]
549    fn stereographic_pole_denominator_clamped_finite() {
550        for pole in [Vec4::W, Vec4::new(0.5, 0.5, 0.5, 0.5)] {
551            let proj = Projection::Stereographic { pole };
552            let at_pole = project4(pole, &proj);
553            for c in [at_pole.x, at_pole.y, at_pole.z] {
554                assert!(
555                    c.is_finite(),
556                    "pole input must clamp finite, got {at_pole:?}"
557                );
558            }
559
560            let (e1, _, _) = perp_frame(pole);
561            let near = (pole * 0.9999 + e1 * 0.01).normalize();
562            let near_img = project4(near, &proj);
563            for c in [near_img.x, near_img.y, near_img.z] {
564                assert!(
565                    c.is_finite(),
566                    "near-pole input must stay finite, got {near_img:?}"
567                );
568            }
569        }
570    }
571
572    #[test]
573    fn stereographic_antipode_maps_to_origin() {
574        for pole in [Vec4::W, Vec4::new(0.1, -0.2, 0.3, 0.9).normalize()] {
575            let proj = Projection::Stereographic { pole };
576            let got = project4(-pole, &proj);
577            assert_relative_eq!(got.x, 0.0, epsilon = 1e-6);
578            assert_relative_eq!(got.y, 0.0, epsilon = 1e-6);
579            assert_relative_eq!(got.z, 0.0, epsilon = 1e-6);
580        }
581    }
582
583    #[test]
584    fn stereographic_frame_breaks_axis_ties_in_index_order() {
585        let pole = Vec4::new(0.0, 0.0, 1.0, 1.0).normalize();
586        let proj = Projection::Stereographic { pole };
587        let p = Vec4::new(0.5, 0.5, -0.5, 0.5);
588        let first = project4(p, &proj);
589
590        assert_relative_eq!(first.x, GOLDEN_TIE_FRAME.x, epsilon = 1e-6);
591        assert_relative_eq!(first.y, GOLDEN_TIE_FRAME.y, epsilon = 1e-6);
592        assert_relative_eq!(first.z, GOLDEN_TIE_FRAME.z, epsilon = 1e-6);
593    }
594
595    #[test]
596    fn stereographic_frame_orthonormal_for_every_pole() {
597        let mut poles = vec![Vec4::X, Vec4::Y, Vec4::Z, Vec4::W];
598        for axis in [Vec4::X, Vec4::Y, Vec4::Z, Vec4::W] {
599            for k in 1..6 {
600                let t = k as f32 / 6.0;
601                poles.push((Vec4::splat(0.5) * (1.0 - t) + axis * t).normalize());
602            }
603        }
604        for pole in poles {
605            let (e1, e2, e3) = perp_frame(pole);
606            for e in [e1, e2, e3] {
607                assert!(e.is_finite(), "frame must be finite for pole {pole:?}");
608            }
609            assert_relative_eq!(e1.dot(e1), 1.0, epsilon = 1e-5);
610            assert_relative_eq!(e2.dot(e2), 1.0, epsilon = 1e-5);
611            assert_relative_eq!(e3.dot(e3), 1.0, epsilon = 1e-5);
612            assert!(e1.dot(e2).abs() < 1e-5, "e1·e2 for pole {pole:?}");
613            assert!(e1.dot(e3).abs() < 1e-5, "e1·e3 for pole {pole:?}");
614            assert!(e2.dot(e3).abs() < 1e-5, "e2·e3 for pole {pole:?}");
615            assert!(e1.dot(pole).abs() < 1e-5, "e1 ⟂ pole {pole:?}");
616            assert!(e2.dot(pole).abs() < 1e-5, "e2 ⟂ pole {pole:?}");
617            assert!(e3.dot(pole).abs() < 1e-5, "e3 ⟂ pole {pole:?}");
618        }
619    }
620
621    #[test]
622    fn stereographic_is_conformal() {
623        let s = SphericalS3Embedded;
624        let pole = Vec4::W;
625        let proj = Projection::Stereographic { pole };
626        let v = Vec4::new(0.3, -0.1, 0.2, -0.5).normalize();
627        let a = Vec4::new(0.5, 0.4, -0.1, -0.3).normalize();
628        let b = Vec4::new(-0.2, 0.3, 0.6, -0.4).normalize();
629
630        let ta = s.log(v, a);
631        let tb = s.log(v, b);
632        let intrinsic = (ta.dot(tb) / (ta.length() * tb.length()))
633            .clamp(-1.0, 1.0)
634            .acos();
635
636        let step = 1e-3;
637        let pv = project4(v, &proj);
638        let pa = project4(s.exp(v, ta.normalize() * step), &proj);
639        let pb = project4(s.exp(v, tb.normalize() * step), &proj);
640        let da = pa - pv;
641        let db = pb - pv;
642        let projected = (da.dot(db) / (da.length() * db.length()))
643            .clamp(-1.0, 1.0)
644            .acos();
645        assert_relative_eq!(projected, intrinsic, epsilon = 1e-2);
646    }
647
648    #[test]
649    fn r4_tessellate_lerps_all_components() {
650        let p0 = Vec4::new(0.0, 0.0, 0.0, 0.0);
651        let p1 = Vec4::new(4.0, 8.0, 12.0, 16.0);
652        let mut out = Vec::new();
653        <EuclideanR4 as RasterizableSpace<4>>::tessellate_segment(p0, p1, 2, |p| out.push(p));
654        assert_eq!(out.len(), 3);
655        assert_eq!(out[0], p0);
656        assert_eq!(out[1], Vec4::new(2.0, 4.0, 6.0, 8.0));
657        assert_eq!(out[2], p1);
658    }
659}
loam_math/rasterizable.rs

loam_math/
rasterizable.rs
