loam_render/raymarch/scene/
mod.rs

1//! Emit contract, shared by the 3D and 4D paths: every baked constant goes
2//! through `literal::wgsl_f32`, so parsing the emitted literal recovers the
3//! exact input bits and the emitter contributes no floor to CPU/GPU parity.
4
5mod literal;
6pub mod primitive;
7pub mod primitive4;
8pub mod scene3;
9pub mod scene4;
10
11pub use loam_shape::Shape;
12pub use primitive::Primitive;
13pub use primitive4::Primitive4;
14pub use scene3::{PrimitiveKind, Scene, SceneNode};
15pub use scene4::{
16    Scene4, SceneNode4, PRIM_KIND_HALFSPACE4D, PRIM_KIND_HYPERSPHERE4D, PRIM_KIND_OTHER,
17};
18
19/// Returned by shapes with no closed-form SDF in the emitted dimension.
20pub const SENTINEL_DISTANCE: f32 = 1e9;
21
22#[cfg(test)]
23mod tests {
24    use super::*;
25    use glam::Vec3;
26
27    fn one_shape_per_kind() -> Vec<Shape> {
28        use glam::{Vec2, Vec4};
29        vec![
30            Shape::sphere_at(Vec3::new(0.05, -0.02, 0.03), 0.25),
31            Shape::HalfSpace {
32                normal: Vec3::Y,
33                offset: -0.5,
34            },
35            Shape::HalfSpace4D {
36                normal: Vec4::Y,
37                offset: -0.25,
38            },
39            Shape::Box3 {
40                half_extents: Vec3::new(0.4, 0.3, 0.2),
41            },
42            Shape::Polygon2D {
43                vertices: vec![Vec2::ZERO, Vec2::X, Vec2::Y],
44            },
45            Shape::ConvexPolytope3D {
46                vertices: vec![Vec3::ZERO, Vec3::X, Vec3::Y, Vec3::Z],
47            },
48            Shape::ConvexPolytope4D {
49                vertices: vec![Vec4::ZERO, Vec4::X, Vec4::Y, Vec4::Z, Vec4::W],
50            },
51            Shape::HyperSphere4D {
52                center: Vec4::new(0.1, 0.0, -0.1, 0.2),
53                radius: 0.3,
54            },
55        ]
56    }
57
58    #[test]
59    fn unsupported_primitives_keep_the_distance_sentinel() {
60        use loam_math::{
61            BlendedSpace, EuclideanR3, HyperbolicH3, LinearBlendX, Space, SphericalS3,
62        };
63        fn check<S: Space<Point = Vec3, Vector = Vec3>>(space: &S, halfspace_supported: bool) {
64            let unsupported = [
65                false,
66                !halfspace_supported,
67                true,
68                false,
69                true,
70                true,
71                true,
72                true,
73            ];
74            for (shape, absent) in one_shape_per_kind().into_iter().zip(unsupported) {
75                assert_eq!(
76                    shape.eval(space, Vec3::new(0.11, -0.07, 0.13)) == SENTINEL_DISTANCE,
77                    absent,
78                    "{:?}",
79                    shape.kind()
80                );
81            }
82        }
83        check(&EuclideanR3, true);
84        check(&HyperbolicH3, false);
85        check(&SphericalS3, false);
86        check(
87            &BlendedSpace::new(
88                EuclideanR3,
89                HyperbolicH3,
90                LinearBlendX::new(-0.5, 0.5).unwrap(),
91            ),
92            false,
93        );
94        for (shape, absent) in one_shape_per_kind()
95            .into_iter()
96            .zip([true, true, false, true, true, true, true, false])
97        {
98            assert_eq!(
99                shape.eval_4d(glam::Vec4::new(0.11, -0.07, 0.13, 0.05)) == SENTINEL_DISTANCE,
100                absent,
101                "{:?}",
102                shape.kind()
103            );
104        }
105    }
106
107    #[test]
108    fn halfspace_eval_follows_the_chart_flatness_gate() {
109        use loam_math::{EuclideanR3, HyperbolicH3};
110        let plane = Shape::HalfSpace {
111            normal: Vec3::Y,
112            offset: -0.5,
113        };
114        let p = Vec3::new(0.0, 0.25, 0.0);
115        assert!((plane.eval(&EuclideanR3, p) - 0.75).abs() < 1e-6);
116        assert_eq!(plane.eval(&HyperbolicH3, p), SENTINEL_DISTANCE);
117        let torus = loam_math::FlatTorus3::cube(2.0);
118        assert_eq!(plane.eval(&torus, p), SENTINEL_DISTANCE);
119        assert_eq!(plane.eval(&torus, p + Vec3::Y * 2.0), SENTINEL_DISTANCE);
120    }
121
122    fn deterministic_pair_samples(seed: u32, count: usize, extent: f32) -> Vec<(Vec3, Vec3)> {
123        let mut state = seed;
124        let mut next_f32 = || {
125            state ^= state << 13;
126            state ^= state >> 17;
127            state ^= state << 5;
128            (state as f32 / u32::MAX as f32) * 2.0 - 1.0
129        };
130        (0..count)
131            .map(|_| {
132                let a = Vec3::new(
133                    next_f32() * extent,
134                    next_f32() * extent,
135                    next_f32() * extent,
136                );
137                let b = Vec3::new(
138                    next_f32() * extent,
139                    next_f32() * extent,
140                    next_f32() * extent,
141                );
142                (a, b)
143            })
144            .collect()
145    }
146
147    fn deterministic_samples(seed: u32, count: usize, extent: f32) -> Vec<Vec3> {
148        let mut state = seed;
149        let mut next_f32 = || {
150            state ^= state << 13;
151            state ^= state >> 17;
152            state ^= state << 5;
153            (state as f32 / u32::MAX as f32) * 2.0 - 1.0
154        };
155        (0..count)
156            .map(|_| {
157                Vec3::new(
158                    next_f32() * extent,
159                    next_f32() * extent,
160                    next_f32() * extent,
161                )
162            })
163            .collect()
164    }
165
166    fn assert_lipschitz_1_under_space_metric<S, F>(label: &str, space: &S, sdf: F, extent: f32)
167    where
168        S: loam_math::Space<Point = Vec3, Vector = Vec3>,
169        F: Fn(Vec3) -> f32,
170    {
171        for (a, b) in deterministic_pair_samples(0xABCD_1234, 256, extent) {
172            let separation = space.distance(a, b);
173            if separation < 1e-6 {
174                continue;
175            }
176            let delta = (sdf(a) - sdf(b)).abs();
177            assert!(
178                delta <= separation * (1.0 + 1e-5),
179                "{label}: |sdf({a:?}) - sdf({b:?})| = {delta} exceeds d(a, b) = {separation}",
180            );
181        }
182    }
183
184    #[test]
185    fn scene_eval_is_lipschitz_1_under_the_space_metric() {
186        use loam_math::{EuclideanR3, HyperbolicH3, SphericalS3};
187        let scene = Scene::new(
188            SceneNode::sphere(Vec3::new(0.1, 0.0, 0.0), 0.2)
189                .smooth_union(SceneNode::cube(0.15), 0.06)
190                .union(SceneNode::sphere(Vec3::new(-0.2, 0.1, 0.0), 0.12))
191                .subtract(SceneNode::sphere(Vec3::new(0.0, 0.2, 0.0), 0.08))
192                .intersect(SceneNode::plane(Vec3::Y, -0.6)),
193        );
194
195        assert_lipschitz_1_under_space_metric(
196            "E3",
197            &EuclideanR3,
198            |p| scene.eval(&EuclideanR3, p),
199            1.0,
200        );
201        let scene = Scene::new(
202            SceneNode::sphere(Vec3::new(0.1, 0.0, 0.0), 0.2)
203                .smooth_union(SceneNode::sphere(Vec3::new(-0.2, 0.1, 0.0), 0.12), 0.06)
204                .subtract(SceneNode::sphere(Vec3::new(0.0, 0.2, 0.0), 0.08)),
205        );
206        assert_lipschitz_1_under_space_metric(
207            "H3",
208            &HyperbolicH3,
209            |p| scene.eval(&HyperbolicH3, p),
210            0.3,
211        );
212        assert_lipschitz_1_under_space_metric(
213            "S3",
214            &SphericalS3,
215            |p| scene.eval(&SphericalS3, p),
216            0.3,
217        );
218    }
219
220    #[test]
221    fn sphere_eval_is_zero_on_the_geodesic_surface_in_every_space() {
222        use loam_math::{EuclideanR3, HyperbolicH3, Space, SphericalS3};
223
224        fn check<S: Space<Point = Vec3, Vector = Vec3>>(space: &S, label: &str) {
225            let center = Vec3::new(0.05, -0.03, 0.02);
226            let radius = 0.2_f32;
227            let shape = Shape::sphere_at(center, radius);
228            assert!(
229                (shape.eval(space, center) + radius).abs() < 1e-6,
230                "{label}: centre must read -radius",
231            );
232            for direction in [Vec3::X, Vec3::Y, Vec3::Z, -Vec3::X, Vec3::ONE.normalize()] {
233                let probe = direction * 0.1;
234                let probe_arc = space.distance(center, space.exp(center, probe));
235                let at_arc = |arc: f32| space.exp(center, probe * (arc / probe_arc));
236
237                assert!(
238                    shape.eval(space, at_arc(radius)).abs() < 1e-5,
239                    "{label}: the geodesic sphere of radius r must be the zero set",
240                );
241                assert!(
242                    shape.eval(space, at_arc(radius * 2.0)) > 0.0,
243                    "{label}: sign outside",
244                );
245                assert!(
246                    shape.eval(space, at_arc(radius * 0.5)) < 0.0,
247                    "{label}: sign inside",
248                );
249            }
250        }
251
252        check(&EuclideanR3, "E3");
253        check(&HyperbolicH3, "H3");
254        check(&SphericalS3, "S3");
255    }
256
257    #[test]
258    fn box_and_halfspace_eval_zero_sets_and_signs_in_e3() {
259        use loam_math::EuclideanR3;
260        let half_extents = Vec3::splat(0.4);
261        let box3 = Shape::Box3 { half_extents };
262        assert!((box3.eval(&EuclideanR3, Vec3::ZERO) + 0.4).abs() < 1e-6);
263        assert!(box3.eval(&EuclideanR3, Vec3::new(0.4, 0.0, 0.0)).abs() < 1e-6);
264        assert!((box3.eval(&EuclideanR3, Vec3::new(1.4, 0.0, 0.0)) - 1.0).abs() < 1e-5);
265        let corner_offset = Vec3::splat(1.0);
266        assert!(
267            (box3.eval(&EuclideanR3, half_extents + corner_offset) - corner_offset.length()).abs()
268                < 1e-5
269        );
270
271        let plane = Shape::HalfSpace {
272            normal: Vec3::Y,
273            offset: -0.5,
274        };
275        assert!(plane.eval(&EuclideanR3, Vec3::new(0.0, -0.5, 0.0)).abs() < 1e-6);
276        assert!((plane.eval(&EuclideanR3, Vec3::Y) - 1.5).abs() < 1e-5);
277        assert!(plane.eval(&EuclideanR3, Vec3::new(0.0, -1.0, 0.0)) < 0.0);
278    }
279
280    #[test]
281    fn union_commutes_and_difference_does_not() {
282        use loam_math::EuclideanR3;
283        let left = SceneNode::sphere(Vec3::new(-0.1, 0.0, 0.0), 0.3);
284        let right = SceneNode::sphere(Vec3::new(0.1, 0.0, 0.0), 0.3);
285        let union_lr = Scene::new(left.clone().union(right.clone()));
286        let union_rl = Scene::new(right.clone().union(left.clone()));
287        let diff_lr = Scene::new(left.clone().subtract(right.clone()));
288        let diff_rl = Scene::new(right.subtract(left));
289
290        let mut asymmetric = 0usize;
291        for p in deterministic_samples(0x0BAD_F00D, 256, 0.6) {
292            assert_eq!(
293                union_lr.eval(&EuclideanR3, p),
294                union_rl.eval(&EuclideanR3, p)
295            );
296            if diff_lr.eval(&EuclideanR3, p) != diff_rl.eval(&EuclideanR3, p) {
297                asymmetric += 1;
298            }
299        }
300        assert!(
301            asymmetric > 0,
302            "A minus B and B minus A must differ somewhere in the sampled volume",
303        );
304    }
305
306    #[test]
307    fn smooth_union_underestimates_min_and_converges_to_it() {
308        use loam_math::EuclideanR3;
309        let left = SceneNode::sphere(Vec3::new(-0.2, 0.0, 0.0), 0.25);
310        let right = SceneNode::sphere(Vec3::new(0.2, 0.0, 0.0), 0.25);
311        let hard = Scene::new(left.clone().union(right.clone()));
312        let samples = deterministic_samples(0x51DF_00D5, 512, 0.7);
313
314        for k in [0.25_f32, 0.05, 1e-3] {
315            let soft = Scene::new(left.clone().smooth_union(right.clone(), k));
316            for &p in &samples {
317                let smooth = soft.eval(&EuclideanR3, p);
318                let sharp = hard.eval(&EuclideanR3, p);
319                assert!(
320                    smooth <= sharp + 1e-6,
321                    "k={k}: smooth_union {smooth} must not exceed min {sharp} at {p:?}",
322                );
323                assert!(
324                    sharp - smooth <= k * 0.25 + 1e-6,
325                    "k={k}: gap {} exceeds the k/4 worst case at {p:?}",
326                    sharp - smooth,
327                );
328            }
329        }
330    }
331
332    #[test]
333    fn smooth_union_matches_min_outside_the_blend_band() {
334        use loam_math::EuclideanR3;
335        let left = SceneNode::sphere(Vec3::new(-0.5, 0.0, 0.0), 0.2);
336        let right = SceneNode::sphere(Vec3::new(0.5, 0.0, 0.0), 0.2);
337        let soft = Scene::new(left.clone().smooth_union(right.clone(), 0.02));
338        let hard = Scene::new(left.union(right));
339
340        let p = Vec3::new(-0.5, 0.0, 0.0);
341        assert_eq!(soft.eval(&EuclideanR3, p), hard.eval(&EuclideanR3, p));
342    }
343
344    #[test]
345    fn hyperslice_eval_tracks_distance_and_kind_through_combinators() {
346        use glam::Vec4;
347        use scene4::{PRIM_KIND_HALFSPACE4D, PRIM_KIND_HYPERSPHERE4D, PRIM_KIND_OTHER};
348
349        let ball = SceneNode4::hypersphere(Vec4::ZERO, 0.5);
350        let floor = SceneNode4::halfspace(Vec4::Y, -0.4);
351
352        let union = Scene4::new(ball.clone().union(floor.clone()));
353        let (dist, kind) = union.eval_at(Vec3::new(0.0, 0.6, 0.0), 0.0, true);
354        assert!((dist - 0.1).abs() < 1e-6);
355        assert_eq!(kind, PRIM_KIND_HYPERSPHERE4D);
356        let (_, kind) = union.eval_at(Vec3::new(2.0, -0.3, 0.0), 0.0, true);
357        assert_eq!(kind, PRIM_KIND_HALFSPACE4D);
358
359        let intersection = Scene4::new(ball.clone().intersect(floor.clone()));
360        let (dist, kind) = intersection.eval_at(Vec3::new(0.0, 0.6, 0.0), 0.0, true);
361        assert!((dist - 1.0).abs() < 1e-6);
362        assert_eq!(kind, PRIM_KIND_HALFSPACE4D);
363
364        let difference = Scene4::new(ball.subtract(floor));
365        let (_, kind) = difference.eval_at(Vec3::ZERO, 0.0, true);
366        assert_eq!(kind, PRIM_KIND_OTHER);
367    }
368
369    #[test]
370    fn hyperslice_radius_shrinks_with_the_slice_coordinate() {
371        use glam::Vec4;
372        let radius = 0.5_f32;
373        let scene = Scene4::new(SceneNode4::hypersphere(Vec4::ZERO, radius));
374        for w in [0.0_f32, 0.2, 0.4] {
375            let sliced_radius = (radius * radius - w * w).sqrt();
376            assert!(
377                scene
378                    .eval(Vec3::new(sliced_radius, 0.0, 0.0), w, true)
379                    .abs()
380                    < 1e-6,
381                "w={w}: sliced surface point must read zero",
382            );
383        }
384        assert!(scene.eval(Vec3::ZERO, 0.6, true) > 0.0);
385    }
386
387    #[test]
388    fn hyperslice_gate_off_returns_the_sentinel_for_halfspaces_only() {
389        use glam::Vec4;
390        let scene = Scene4::new(
391            SceneNode4::hypersphere(Vec4::ZERO, 0.5).union(SceneNode4::halfspace(Vec4::Y, -0.4)),
392        );
393        let just_above_floor = Vec3::new(3.0, -0.39, 0.0);
394        assert!(scene.eval(just_above_floor, 0.0, true) < 0.02);
395        let ball_only = Scene4::new(SceneNode4::hypersphere(Vec4::ZERO, 0.5));
396        assert_eq!(
397            scene.eval(just_above_floor, 0.0, false),
398            ball_only.eval(just_above_floor, 0.0, true),
399        );
400    }
401
402    #[test]
403    fn scene4_eval_is_lipschitz_1_in_flat_r4() {
404        use glam::Vec4;
405        let scene = Scene4::new(
406            SceneNode4::hypersphere(Vec4::new(0.1, 0.0, -0.1, 0.05), 0.4)
407                .union(SceneNode4::halfspace(Vec4::Y, -0.5))
408                .subtract(SceneNode4::hypersphere(Vec4::new(0.3, 0.0, 0.0, 0.0), 0.15)),
409        );
410        let mut state: u32 = 0x5555_3333;
411        let mut next_f32 = || {
412            state ^= state << 13;
413            state ^= state >> 17;
414            state ^= state << 5;
415            (state as f32 / u32::MAX as f32) * 2.0 - 1.0
416        };
417        for _ in 0..256 {
418            let a = Vec4::new(next_f32(), next_f32(), next_f32(), next_f32()) * 2.0;
419            let b = Vec4::new(next_f32(), next_f32(), next_f32(), next_f32()) * 2.0;
420            let separation = (a - b).length();
421            if separation < 1e-6 {
422                continue;
423            }
424            let delta =
425                (scene.eval(a.truncate(), a.w, true) - scene.eval(b.truncate(), b.w, true)).abs();
426            assert!(
427                delta <= separation * (1.0 + 1e-5),
428                "|sdf({a:?}) - sdf({b:?})| = {delta} exceeds |a - b| = {separation}",
429            );
430        }
431    }
432
433    const BEYOND_ABSTRACT_INT: f32 = 1.0e19;
434
435    fn assert_naga_accepts(source: &str) {
436        let module = naga::front::wgsl::parse_str(source)
437            .unwrap_or_else(|e| panic!("WGSL parse failed: {e}\n--- source ---\n{source}"));
438        naga::valid::Validator::new(
439            naga::valid::ValidationFlags::all(),
440            naga::valid::Capabilities::empty(),
441        )
442        .validate(&module)
443        .unwrap_or_else(|e| panic!("WGSL validation failed: {e:?}\n--- source ---\n{source}"));
444    }
445
446    #[test]
447    fn scene3_beyond_abstract_int_range_emits_wgsl_naga_accepts() {
448        use loam_math::{EuclideanR3, WgslSpace};
449        let magnitude = BEYOND_ABSTRACT_INT;
450        let scene = Scene::new(
451            SceneNode::sphere(Vec3::new(magnitude, -magnitude, 0.0), magnitude)
452                .union(SceneNode::plane(Vec3::Y, -magnitude))
453                .smooth_union(SceneNode::box_(Vec3::splat(magnitude)), magnitude),
454        );
455        let probe = format!(
456            "{prelude}\n{scene}\n\
457             @compute @workgroup_size(1) fn main() {{\n\
458             \t_ = loam_scene_sdf(vec3<f32>(0.0));\n\
459             }}\n",
460            prelude = EuclideanR3.wgsl_impl(),
461            scene = scene.to_wgsl(&EuclideanR3),
462        );
463        assert_naga_accepts(&probe);
464    }
465
466    #[test]
467    fn scene4_beyond_abstract_int_range_emits_wgsl_naga_accepts() {
468        use glam::Vec4;
469        let magnitude = BEYOND_ABSTRACT_INT;
470        let scene = Scene4::new(
471            SceneNode4::hypersphere(Vec4::splat(magnitude), magnitude)
472                .union(SceneNode4::halfspace(Vec4::Y, -magnitude)),
473        );
474        let native = format!(
475            "{scene}\n\
476             @compute @workgroup_size(1) fn main() {{\n\
477             \t_ = loam_scene_sdf_4d(vec4<f32>(0.0));\n\
478             }}\n",
479            scene = scene.to_wgsl_4d(),
480        );
481        assert_naga_accepts(&native);
482
483        let hyperslice = format!(
484            "{scene}\n\
485             @compute @workgroup_size(1) fn main() {{\n\
486             \t_ = loam_scene_sdf(vec3<f32>(0.0));\n\
487             \t_ = loam_scene_max_t(vec3<f32>(0.0), vec3<f32>(0.0, -1.0, 0.0));\n\
488             }}\n",
489            scene = scene.to_hyperslice_wgsl("0.0"),
490        );
491        assert_naga_accepts(&hyperslice);
492    }
493
494    #[test]
495    #[should_panic(expected = "non-finite")]
496    fn scene3_rejects_a_non_finite_constant() {
497        use loam_math::EuclideanR3;
498        let scene = Scene::new(SceneNode::sphere(Vec3::ZERO, f32::INFINITY));
499        let _ = scene.to_wgsl(&EuclideanR3);
500    }
501
502    #[test]
503    #[should_panic(expected = "non-finite")]
504    fn scene4_rejects_a_non_finite_constant() {
505        use glam::Vec4;
506        let scene = Scene4::new(SceneNode4::halfspace(Vec4::Y, f32::NAN));
507        let _ = scene.to_wgsl_4d();
508    }
509}
loam_render/raymarch/scene/mod.rs

loam_render/raymarch/scene/
mod.rs
