diff --git a/crates/resolved/tests/analytic.rs b/crates/resolved/tests/analytic.rs new file mode 100644 index 0000000..f75a55d --- /dev/null +++ b/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, boresight: Vector3) -> 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 = Vector3::new(0.0, 3.844e8, 0.0); + +fn scene(sun: Vector3, targets: Vec) -> 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, + radius: f64, + sun_dir: Vector3, + camera: Vector3, +) -> 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, 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::() + .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 + ); + } + } +} diff --git a/docs/optics.md b/docs/optics.md index c4574c3..1c7f561 100644 --- a/docs/optics.md +++ b/docs/optics.md @@ -433,6 +433,16 @@ the surface normal and the Sun. The night side has i past 90 degrees and is black. The Sun's distance d is taken from each body's center, since the difference across a body is a few parts in 10^5. +The total light a sphere sends a far-away camera has a closed form. A +Lambertian sphere of albedo rho and radius R, at distance D, delivers an +irradiance of (2/3) rho E (R/D)^2 Phi(alpha), where alpha is the phase angle +(Sun-sphere-camera) and Phi(alpha) = (sin alpha + (pi - alpha) cos alpha) / pi. +It falls from 1 at full phase to 0.61 at 60 degrees and 0.11 at 120 degrees. +Up close it is not exact: the near side of the sphere is nearer the camera +than its center, so it counts for more. At 40,000 km from the Moon (R/D = +0.043), the real total is 3 % above the formula. The renderer's image agrees +with the exact surface integral to 0.03 % or better (`tests/analytic.rs`). + This gets the overall brightness and the phase (how much of the disk is lit) right. It misses some things that are visible in real images: diff --git a/hooks/pre-push b/hooks/pre-push index 54fca37..4a5afe8 100755 --- a/hooks/pre-push +++ b/hooks/pre-push @@ -7,8 +7,14 @@ set -e echo "pre-push: cargo test" cargo test --workspace --all-features --quiet -echo "pre-push: cargo test --no-default-features" -cargo test --workspace --no-default-features --quiet +# Without default features: every target must still build cleanly, and the unit +# tests run. The GPU integration tests do not depend on features, so the run above +# covers them. +echo "pre-push: cargo clippy --no-default-features" +cargo clippy --workspace --all-targets --no-default-features --quiet -- -D warnings + +echo "pre-push: cargo test --no-default-features (unit tests)" +cargo test --workspace --no-default-features --lib --bins --quiet echo "pre-push: cargo doc" RUSTDOCFLAGS="-D warnings" cargo doc --workspace --no-deps --all-features --quiet