Last couple of analytics. Finished Milestone 1

crates/resolved/tests/analytic.rs

@@ -0,0 +1,327 @@
+//! Rendered results against closed-form answers, through the public renderer, on the
+//! GPU and on the CPU (lavapipe).
+
+mod common;
+
+use resolved::bodies::{Appearance, MOON_RADIUS};
+use resolved::gpu::AdapterChoice;
+use resolved::model::Part;
+use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3};
+use resolved::radiometry::{Radiometry, AU};
+use resolved::{CameraConfig, Material, Mesh, Pose, Renderer, RendererOptions, Scene, Target};
+
+const ADAPTERS: [AdapterChoice; 2] = [AdapterChoice::Gpu, AdapterChoice::Cpu];
+
+fn renderer(adapter: &AdapterChoice) -> Renderer {
+    let options = RendererOptions {
+        adapter: adapter.clone(),
+        read_radiance: true,
+        ..Default::default()
+    };
+    Renderer::new(CameraConfig::example(), options).unwrap()
+}
+
+/// A camera at `position` whose boresight is `boresight`, with image up as close to
+/// inertial +Z as possible.
+fn camera(position: Vector3<f64>, boresight: Vector3<f64>) -> Pose {
+    Pose::looking_at(position, position + boresight, Vector3::z()).unwrap()
+}
+
+/// The Moon on +Y at its mean distance, Earth at the origin, the Sun wherever given.
+const MOON: Vector3<f64> = Vector3::new(0.0, 3.844e8, 0.0);
+
+fn scene(sun: Vector3<f64>, targets: Vec<Target>) -> Scene {
+    Scene {
+        earth_attitude: UnitQuaternion::identity(),
+        sun,
+        moon: MOON,
+        targets,
+    }
+}
+
+/// The radiance-weighted solid angle a Lambertian sphere of unit albedo and unit
+/// irradiance presents to a camera: the integral over its lit, visible surface of
+/// (cos i / pi) cos e / r^2 dA, for incidence angle i, emission angle e and distance r.
+/// Midpoint quadrature in f64 on the CPU, independent of the renderer.
+fn lambert_sphere_integral(
+    center: Vector3<f64>,
+    radius: f64,
+    sun_dir: Vector3<f64>,
+    camera: Vector3<f64>,
+) -> f64 {
+    let (n_theta, n_phi) = (800, 1600);
+    let (d_theta, d_phi) = (
+        std::f64::consts::PI / n_theta as f64,
+        std::f64::consts::TAU / n_phi as f64,
+    );
+    let mut total = 0.0;
+    for i in 0..n_theta {
+        let theta = (i as f64 + 0.5) * d_theta;
+        let area = radius * radius * theta.sin() * d_theta * d_phi;
+        for j in 0..n_phi {
+            let phi = (j as f64 + 0.5) * d_phi;
+            let n = Vector3::new(
+                theta.sin() * phi.cos(),
+                theta.sin() * phi.sin(),
+                theta.cos(),
+            );
+            let cos_i = n.dot(&sun_dir);
+            let to_camera = camera - (center + radius * n);
+            let r = to_camera.norm();
+            let cos_e = n.dot(&to_camera) / r;
+            if cos_i > 0.0 && cos_e > 0.0 {
+                total += cos_i / std::f64::consts::PI * cos_e * area / (r * r);
+            }
+        }
+    }
+    total
+}
+
+#[test]
+fn far_field_limit_of_the_sphere_integral() {
+    // Far from the sphere, the integral approaches the closed-form Lambert phase curve
+    // (2/3) (R/D)^2 Phi(alpha), Phi(alpha) = (sin alpha + (pi - alpha) cos alpha) / pi.
+    // Checks the quadrature the rendering test relies on.
+    let (radius, distance) = (1.0, 1.0e4);
+    for phase_deg in [0.0f64, 60.0, 120.0] {
+        let phase = phase_deg.to_radians();
+        let sun_dir = Vector3::new(phase.cos(), -phase.sin(), 0.0);
+        let integral =
+            lambert_sphere_integral(Vector3::zeros(), radius, sun_dir, Vector3::x() * distance);
+        let phi =
+            (phase.sin() + (std::f64::consts::PI - phase) * phase.cos()) / std::f64::consts::PI;
+        let closed = 2.0 / 3.0 * (radius / distance).powi(2) * phi;
+        assert!(
+            // The quadrature resolves the thin crescent at 120 degrees to about 2e-4.
+            (integral / closed - 1.0).abs() < 5e-4,
+            "phase {phase_deg}: {integral:e} vs {closed:e}"
+        );
+    }
+}
+
+#[test]
+fn moon_disk_brightness_matches_the_sphere_integral() {
+    // Summed over the image, each pixel's channel radiance times the solid angle it
+    // covers, (p/f)^2 cos^3(theta) for a pixel theta off axis, must equal the same
+    // integral taken over the sphere's surface. At 40,000 km the Moon (R/D = 0.043) is
+    // close enough that the far-field phase curve is 3 % low: the near side of the Moon
+    // is nearer the camera than its center.
+    let _devices = common::device_lock();
+    let config = CameraConfig::example();
+    let k = config.intrinsics();
+    let r = Radiometry::new(&config);
+    let albedo = r.effective_albedo(&Appearance::default().moon_albedo);
+    // Off to the Moon's +X side, so Earth is 84 degrees off axis.
+    let distance = 4.0e7;
+    let position = MOON + Vector3::new(distance, 0.0, 0.0);
+    let solid_angle = |u: usize, v: usize| {
+        let cos = k.ray(&Vector2::new(u as f64 + 0.5, v as f64 + 0.5)).z;
+        cos.powi(3) / (k.fx * k.fy)
+    };
+    // The disk is about 395 px in radius around the image center; sum over a box
+    // around it rather than the whole image.
+    let (u0, u1, v0, v1) = (960 - 420, 960 + 420, 540 - 420, 540 + 420);
+    let phases = [0.0f64, 60.0, 120.0].map(|deg| {
+        // The Sun in the X-Y plane, `deg` from the camera direction seen from the Moon.
+        let phase = deg.to_radians();
+        let sun_dir = Vector3::new(phase.cos(), -phase.sin(), 0.0);
+        (
+            deg,
+            sun_dir,
+            lambert_sphere_integral(MOON, MOON_RADIUS, sun_dir, position),
+        )
+    });
+    for adapter in &ADAPTERS {
+        let mut renderer = renderer(adapter);
+        for (phase_deg, sun_dir, geometry) in phases {
+            let sun = MOON + sun_dir * AU;
+            let frame = renderer
+                .submit(&scene(sun, vec![]), &camera(position, -Vector3::x()))
+                .unwrap()
+                .wait()
+                .unwrap();
+            let radiance = frame.radiance.unwrap();
+            let mut rendered = [0.0f64; 3];
+            for v in v0..v1 {
+                for u in u0..u1 {
+                    let omega = solid_angle(u, v);
+                    for (c, sum) in rendered.iter_mut().enumerate() {
+                        *sum += radiance[v * 1920 + u][c] as f64 * omega;
+                    }
+                }
+            }
+            assert_eq!(
+                radiance[v0 * 1920 + u0],
+                [0.0; 4],
+                "the box contains the disk"
+            );
+            for (c, a) in albedo.iter().enumerate() {
+                // Irradiance at the Moon is the solar flux at 1 AU; the Sun is 1 AU away.
+                let expected = a * r.solar_flux[c] * geometry;
+                let error = rendered[c] / expected - 1.0;
+                // Measured: -0.0005 % at 0 deg, -0.005 % at 60, -0.031 % at 120, on both
+                // devices; the crescent's edge pixels carry the error at 120.
+                assert!(
+                    error.abs() < 1e-3,
+                    "{} phase {phase_deg} channel {c}: {:e} vs {expected:e} ({:+.3} %)",
+                    renderer.adapter_info().name,
+                    rendered[c],
+                    100.0 * error
+                );
+            }
+        }
+    }
+}
+
+/// Whether the ray through the center of pixel (u, v) of the example camera at
+/// `camera` hits a sphere of `radius` at `center`, computed in f64 on the CPU.
+fn ray_hits(camera: &Pose, u: f64, v: f64, center: Vector3<f64>, radius: f64) -> bool {
+    let k = CameraConfig::example().intrinsics();
+    let d = camera.attitude * k.ray(&Vector2::new(u, v));
+    let o = center - camera.position;
+    let b = o.dot(&d);
+    (o - b * d).norm() < radius && b > 0.0
+}
+
+#[test]
+fn moon_limb_lands_to_a_hundredth_of_a_pixel() {
+    // Turn the camera so the Moon's limb crosses the center row at column 960.5 + delta,
+    // then ask whether pixel (960, 540), whose center is at 960.5, sees the Moon.
+    let _devices = common::device_lock();
+    let k = CameraConfig::example().intrinsics();
+    let distance = 4.0e7;
+    let position = MOON + Vector3::new(distance, 0.0, 0.0);
+    let limb = (MOON_RADIUS / distance).asin();
+    let sun = MOON + Vector3::x() * AU;
+    for adapter in &ADAPTERS {
+        let mut renderer = renderer(adapter);
+        for delta in [-0.01, 0.01] {
+            // Looking along -X with +Z up, image right is +Y. Turning the boresight from
+            // the Moon's center toward +Y by the limb angle less the angle of 0.5 + delta
+            // pixels puts the right limb at column 960.5 + delta.
+            let beta = limb - ((0.5 + delta) / k.fx).atan();
+            let boresight = Vector3::new(-beta.cos(), beta.sin(), 0.0);
+            let pose = camera(position, boresight);
+            let expected = ray_hits(&pose, 960.5, 540.5, MOON, MOON_RADIUS);
+            assert_eq!(expected, delta > 0.0, "the limb is where the setup intends");
+            let frame = renderer
+                .submit(&scene(sun, vec![]), &pose)
+                .unwrap()
+                .wait()
+                .unwrap();
+            let lit = frame.radiance.unwrap()[540 * 1920 + 960][1] > 0.0;
+            assert_eq!(
+                lit,
+                expected,
+                "{} delta {delta}",
+                renderer.adapter_info().name
+            );
+        }
+    }
+}
+
+#[test]
+fn mesh_edge_lands_to_a_hundredth_of_a_pixel() {
+    // A 4 m square panel 10 m ahead whose left edge projects to column 960.5 + delta.
+    // Pixel (960, 540) is covered exactly when the edge is left of its center.
+    let _devices = common::device_lock();
+    let k = CameraConfig::example().intrinsics();
+    let position = Vector3::new(7.0e6, 0.0, 0.0);
+    // Looking along +Y, sideways from Earth, which is 90 degrees off axis.
+    let pose = camera(position, Vector3::y());
+    let right = pose.attitude * Vector3::x();
+    let sun = position - Vector3::y() * AU;
+    for adapter in &ADAPTERS {
+        let mut renderer = renderer(adapter);
+        for delta in [-0.01, 0.01] {
+            let edge = (0.5 + delta) / k.fx * 10.0;
+            let corners = [
+                [edge, -2.0],
+                [edge + 4.0, -2.0],
+                [edge + 4.0, 2.0],
+                [edge, 2.0],
+            ];
+            // Panel coordinates: x along image right, y along image down, at 10 m.
+            let down = pose.attitude * Vector3::y();
+            let at = |[x, y]: [f64; 2]| {
+                (right * x + down * y + Vector3::y() * 10.0)
+                    .cast::<f32>()
+                    .into()
+            };
+            let panel = Mesh {
+                positions: corners.map(at).to_vec(),
+                normals: vec![[0.0, -1.0, 0.0]; 4],
+                part_of: vec![0; 4],
+                indices: vec![0, 1, 2, 0, 2, 3],
+                parts: vec![Part {
+                    name: "panel".into(),
+                    material: 0,
+                }],
+                materials: vec![Material::gray("gray", 0.5)],
+            };
+            let model = renderer.load_model(&panel).unwrap();
+            let target = Target {
+                model,
+                pose: Pose::new(position, UnitQuaternion::identity()),
+            };
+            let frame = renderer
+                .submit(&scene(sun, vec![target]), &pose)
+                .unwrap()
+                .wait()
+                .unwrap();
+            let covered = frame.radiance.unwrap()[540 * 1920 + 960][1] > 0.0;
+            assert_eq!(
+                covered,
+                delta < 0.0,
+                "{} delta {delta}",
+                renderer.adapter_info().name
+            );
+        }
+    }
+}
+
+#[test]
+fn target_in_umbra_reads_the_bias() {
+    // A panel deep in Earth's shadow, turned to face both the Sun (+X) and the camera
+    // (-Y), so without the eclipse it would be lit at 45 degrees. In the umbra every
+    // pixel it covers reads exactly the 64 DN bias.
+    let _devices = common::device_lock();
+    let position = Vector3::new(-7.0e6, 0.0, 0.0);
+    let sun = Vector3::x() * AU;
+    let pose = camera(position - Vector3::y() * 30.0, Vector3::y());
+    for adapter in &ADAPTERS {
+        let mut renderer = renderer(adapter);
+        // In the plane spanned by (1, 1, 0) and Z, facing (1, -1, 0).
+        let h = 2.0 * std::f32::consts::FRAC_1_SQRT_2;
+        let n = std::f32::consts::FRAC_1_SQRT_2;
+        let panel = Mesh {
+            positions: vec![[-h, -h, -2.0], [h, h, -2.0], [h, h, 2.0], [-h, -h, 2.0]],
+            normals: vec![[n, -n, 0.0]; 4],
+            part_of: vec![0; 4],
+            indices: vec![0, 1, 2, 0, 2, 3],
+            parts: vec![Part {
+                name: "panel".into(),
+                material: 0,
+            }],
+            materials: vec![Material::gray("gray", 0.5)],
+        };
+        let model = renderer.load_model(&panel).unwrap();
+        let target = Target {
+            model,
+            pose: Pose::new(position, UnitQuaternion::identity()),
+        };
+        let frame = renderer
+            .submit(&scene(sun, vec![target]), &pose)
+            .unwrap()
+            .wait()
+            .unwrap();
+        for c in 0..3 {
+            assert_eq!(
+                frame.dn_at(960, 540, c),
+                64,
+                "{}",
+                renderer.adapter_info().name
+            );
+        }
+    }
+}