diff --git a/crates/resolved/Cargo.toml b/crates/resolved/Cargo.toml index ae46dca..bbdef5b 100644 --- a/crates/resolved/Cargo.toml +++ b/crates/resolved/Cargo.toml @@ -16,7 +16,7 @@ nalgebra = "0.35" thiserror = "2.0" wgpu = "30.0.1" pollster = "1.0" -bytemuck = "1.0" +bytemuck = { version = "1.0", features = ["derive"] } log = "0.4" [dev-dependencies] diff --git a/crates/resolved/src/bodies.rs b/crates/resolved/src/bodies.rs new file mode 100644 index 0000000..0e934b7 --- /dev/null +++ b/crates/resolved/src/bodies.rs @@ -0,0 +1,35 @@ +//! Earth, Moon and Sun: their shapes and how they reflect light. + +use crate::spectrum::Spectrum; + +/// WGS84 equatorial radius, in meters. +pub const EARTH_EQUATORIAL_RADIUS: f64 = 6_378_137.0; +/// WGS84 polar radius, in meters. +pub const EARTH_POLAR_RADIUS: f64 = 6_356_752.314_245; +/// Mean lunar radius, in meters (IAU WGCCRE). +pub const MOON_RADIUS: f64 = 1_737_400.0; + +/// How Earth and the Moon reflect light. +/// +/// Both are treated as Lambertian with one reflectance spectrum over the whole body. +/// That gives the right overall brightness and phase-dependent shading, but not +/// oceans, clouds, maria or the Moon's backscatter peak; those come with albedo maps +/// and better reflectance models in later milestones. +#[derive(Clone, Debug, PartialEq)] +pub struct Appearance { + /// Earth's reflectance. + pub earth_albedo: Spectrum, + /// The Moon's reflectance. + pub moon_albedo: Spectrum, +} + +impl Default for Appearance { + /// Flat spectra: 0.3 for Earth (its Bond albedo) and 0.12 for the Moon (its + /// average visible normal albedo). + fn default() -> Self { + Self { + earth_albedo: Spectrum::constant(0.3), + moon_albedo: Spectrum::constant(0.12), + } + } +} diff --git a/crates/resolved/src/error.rs b/crates/resolved/src/error.rs index a84ce6a..8465f96 100644 --- a/crates/resolved/src/error.rs +++ b/crates/resolved/src/error.rs @@ -17,6 +17,10 @@ pub enum Error { available: Vec, }, + /// A request this version cannot carry out yet. + #[error("unsupported: {0}")] + Unsupported(String), + /// The adapter refused to open a device. #[error("could not open a device: {0}")] RequestDevice(#[from] wgpu::RequestDeviceError), diff --git a/crates/resolved/src/gpu/bodies.rs b/crates/resolved/src/gpu/bodies.rs new file mode 100644 index 0000000..fa0cbfa --- /dev/null +++ b/crates/resolved/src/gpu/bodies.rs @@ -0,0 +1,226 @@ +//! The pass that draws Earth, the Moon and the Sun. +//! +//! It is the first pass of a frame: it clears the targets, then writes radiance and +//! depth wherever a body is hit. Pixels that see empty space keep zero radiance and +//! depth 0, which reversed Z treats as infinitely far. + +use nalgebra::Vector3; + +use super::targets::{FrameTargets, DEPTH_COMPARE, DEPTH_FORMAT, RADIANCE_FORMAT}; +use super::GpuContext; +use crate::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS, MOON_RADIUS}; +use crate::camera::Intrinsics; +use crate::radiometry::{Radiometry, AU, SUN_RADIUS}; +use crate::scene::{Pose, Scene}; +use crate::Error; + +/// Channels one radiance target holds. +pub const MAX_CHANNELS: usize = 4; + +/// Per-channel constants for the bodies that do not change from frame to frame. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct BodyRadiometry { + /// Channel radiance of Earth's surface facing the Sun at 1 AU: albedo S_k / pi. + pub earth_lambert_1au: [f32; 4], + /// The same for the Moon. + pub moon_lambert_1au: [f32; 4], + /// Channel radiance of the solar disk. + pub sun_radiance: [f32; 4], +} + +impl BodyRadiometry { + /// Folds `appearance` through `radiometry`'s channels. Integrates over wavelength, so + /// it belongs in setup, not in the frame loop. + pub fn new(radiometry: &Radiometry, appearance: &Appearance) -> Result { + let n = radiometry.solar_flux.len(); + if n > MAX_CHANNELS { + return Err(Error::Unsupported(format!( + "{n} channels; at most {MAX_CHANNELS} are supported" + ))); + } + let lambert = |albedo: Vec| { + let mut out = [0.0; 4]; + for (k, a) in albedo.iter().enumerate() { + out[k] = radiometry.lambert_radiance(k, *a, 1.0, AU) as f32; + } + out + }; + let mut sun_radiance = [0.0; 4]; + for (k, s) in sun_radiance.iter_mut().take(n).enumerate() { + *s = radiometry.sun_radiance(k) as f32; + } + Ok(Self { + earth_lambert_1au: lambert(radiometry.effective_albedo(&appearance.earth_albedo)), + moon_lambert_1au: lambert(radiometry.effective_albedo(&appearance.moon_albedo)), + sun_radiance, + }) + } +} + +/// The shader's uniform block; the layout matches `Uniforms` in `bodies.wgsl`. +#[repr(C)] +#[derive(Clone, Copy, Debug, PartialEq, bytemuck::Pod, bytemuck::Zeroable)] +pub struct BodiesUniforms { + intrinsics: [f32; 4], + params: [f32; 4], + earth_center: [f32; 4], + // mat3x3 columns, each padded to 16 bytes. + earth_axes: [[f32; 4]; 3], + earth_shape: [f32; 4], + earth_sun: [f32; 4], + earth_lambert: [f32; 4], + moon_center: [f32; 4], + moon_sun: [f32; 4], + moon_lambert: [f32; 4], + sun_center: [f32; 4], + sun_radiance: [f32; 4], +} + +impl BodiesUniforms { + /// The uniforms for `scene` seen from `camera`. Geometry is differenced in f64 + /// before the cast, as in [`Scene::relative_to`]. + pub fn new( + scene: &Scene, + camera: &Pose, + intrinsics: &Intrinsics, + near: f64, + radiometry: &BodyRadiometry, + ) -> Self { + let rel = scene.relative_to(camera); + let to_camera = camera.attitude.inverse(); + let vec4 = |v: Vector3, w: f64| [v.x, v.y, v.z, w as f32]; + let sun_from = |p: Vector3| (to_camera * (scene.sun - p).normalize()).cast::(); + let scaled = |lambert: [f32; 4], distance: f64| { + let au2 = ((distance / AU) * (distance / AU)) as f32; + lambert.map(|l| l / au2) + }; + let axes = rel.earth.rotation.to_rotation_matrix(); + let column = |i: usize| { + let c = axes.matrix().column(i); + [c.x, c.y, c.z, 0.0] + }; + Self { + intrinsics: [intrinsics.fx, intrinsics.fy, intrinsics.cx, intrinsics.cy] + .map(|x| x as f32), + params: [near as f32, 0.0, 0.0, 0.0], + earth_center: vec4(rel.earth.translation, EARTH_EQUATORIAL_RADIUS), + earth_axes: [column(0), column(1), column(2)], + earth_shape: [ + (EARTH_EQUATORIAL_RADIUS / EARTH_POLAR_RADIUS) as f32, + 0.0, + 0.0, + 0.0, + ], + earth_sun: vec4(sun_from(Vector3::zeros()), 0.0), + earth_lambert: scaled(radiometry.earth_lambert_1au, scene.sun.norm()), + moon_center: vec4(rel.moon, MOON_RADIUS), + moon_sun: vec4(sun_from(scene.moon), 0.0), + moon_lambert: scaled(radiometry.moon_lambert_1au, (scene.sun - scene.moon).norm()), + sun_center: vec4(rel.sun, SUN_RADIUS), + sun_radiance: radiometry.sun_radiance, + } + } +} + +/// The compiled bodies pipeline and its uniform buffer. +#[derive(Debug)] +pub struct BodiesPass { + pipeline: wgpu::RenderPipeline, + uniforms: wgpu::Buffer, + bind_group: wgpu::BindGroup, +} + +impl BodiesPass { + /// Compiles the shader and allocates the uniform buffer. + pub fn new(ctx: &GpuContext) -> Self { + let device = &ctx.device; + let module = device.create_shader_module(wgpu::ShaderModuleDescriptor { + label: Some("bodies"), + source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/bodies.wgsl").into()), + }); + let pipeline = device.create_render_pipeline(&wgpu::RenderPipelineDescriptor { + label: Some("bodies"), + layout: None, + vertex: wgpu::VertexState { + module: &module, + entry_point: Some("vs_main"), + compilation_options: Default::default(), + buffers: &[], + }, + primitive: wgpu::PrimitiveState::default(), + depth_stencil: Some(wgpu::DepthStencilState { + format: DEPTH_FORMAT, + depth_write_enabled: Some(true), + depth_compare: Some(DEPTH_COMPARE), + stencil: Default::default(), + bias: Default::default(), + }), + multisample: Default::default(), + fragment: Some(wgpu::FragmentState { + module: &module, + entry_point: Some("fs_main"), + compilation_options: Default::default(), + targets: &[Some(RADIANCE_FORMAT.into())], + }), + multiview_mask: None, + cache: None, + }); + let uniforms = device.create_buffer(&wgpu::BufferDescriptor { + label: Some("bodies uniforms"), + size: std::mem::size_of::() as u64, + usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }); + let bind_group = device.create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("bodies"), + layout: &pipeline.get_bind_group_layout(0), + entries: &[wgpu::BindGroupEntry { + binding: 0, + resource: uniforms.as_entire_binding(), + }], + }); + Self { + pipeline, + uniforms, + bind_group, + } + } + + /// Uploads this frame's uniforms and records the pass into `encoder`. + pub fn encode( + &self, + ctx: &GpuContext, + encoder: &mut wgpu::CommandEncoder, + targets: &FrameTargets, + uniforms: &BodiesUniforms, + ) { + ctx.queue + .write_buffer(&self.uniforms, 0, bytemuck::bytes_of(uniforms)); + let radiance = targets.radiance.create_view(&Default::default()); + let depth = targets.depth.create_view(&Default::default()); + let mut pass = encoder.begin_render_pass(&wgpu::RenderPassDescriptor { + label: Some("bodies"), + color_attachments: &[Some(wgpu::RenderPassColorAttachment { + view: &radiance, + depth_slice: None, + resolve_target: None, + ops: wgpu::Operations { + load: wgpu::LoadOp::Clear(wgpu::Color::TRANSPARENT), + store: wgpu::StoreOp::Store, + }, + })], + depth_stencil_attachment: Some(wgpu::RenderPassDepthStencilAttachment { + view: &depth, + depth_ops: Some(wgpu::Operations { + load: wgpu::LoadOp::Clear(0.0), + store: wgpu::StoreOp::Store, + }), + stencil_ops: None, + }), + ..Default::default() + }); + pass.set_pipeline(&self.pipeline); + pass.set_bind_group(0, &self.bind_group, &[]); + pass.draw(0..3, 0..1); + } +} diff --git a/crates/resolved/src/gpu/mod.rs b/crates/resolved/src/gpu/mod.rs index 5c9e701..7ce9e77 100644 --- a/crates/resolved/src/gpu/mod.rs +++ b/crates/resolved/src/gpu/mod.rs @@ -1,5 +1,7 @@ -//! GPU plumbing: adapter selection, device setup. +//! GPU plumbing: adapter selection, device setup, render targets and passes. +pub mod bodies; pub mod context; +pub mod targets; pub use context::{adapters, AdapterChoice, GpuContext}; diff --git a/crates/resolved/src/gpu/targets.rs b/crates/resolved/src/gpu/targets.rs new file mode 100644 index 0000000..9151138 --- /dev/null +++ b/crates/resolved/src/gpu/targets.rs @@ -0,0 +1,120 @@ +//! The per-frame images passes render into, and reading them back. + +use super::GpuContext; + +/// Channel radiance, one channel per component. Four channels fit in one target; f32 +/// because solar radiance overflows f16 (#21 measures a pre-exposed f16 alternative). +pub const RADIANCE_FORMAT: wgpu::TextureFormat = wgpu::TextureFormat::Rgba32Float; +/// Reversed-Z depth, near / z: cleared to 0 (infinitely far), nearer is greater. +pub const DEPTH_FORMAT: wgpu::TextureFormat = wgpu::TextureFormat::Depth32Float; +/// Depth test for reversed Z. +pub const DEPTH_COMPARE: wgpu::CompareFunction = wgpu::CompareFunction::Greater; + +/// The images one frame renders into, sized for one camera. +#[derive(Debug)] +pub struct FrameTargets { + /// Image width in pixels. + pub width: u32, + /// Image height in pixels. + pub height: u32, + /// Channel radiance, [`RADIANCE_FORMAT`]. + pub radiance: wgpu::Texture, + /// Depth, [`DEPTH_FORMAT`]. + pub depth: wgpu::Texture, +} + +impl FrameTargets { + /// Allocates targets of `width` x `height` pixels. + pub fn new(ctx: &GpuContext, width: u32, height: u32) -> Self { + let texture = |label, format, usage| { + ctx.device.create_texture(&wgpu::TextureDescriptor { + label: Some(label), + size: wgpu::Extent3d { + width, + height, + depth_or_array_layers: 1, + }, + mip_level_count: 1, + sample_count: 1, + dimension: wgpu::TextureDimension::D2, + format, + usage, + view_formats: &[], + }) + }; + let usage = wgpu::TextureUsages::RENDER_ATTACHMENT + | wgpu::TextureUsages::TEXTURE_BINDING + | wgpu::TextureUsages::COPY_SRC; + Self { + width, + height, + radiance: texture("radiance", RADIANCE_FORMAT, usage), + depth: texture("depth", DEPTH_FORMAT, usage), + } + } + + /// Reads the radiance target back, row by row from the top. Slow; for tests and + /// debugging, not the per-frame path. + pub fn read_radiance(&self, ctx: &GpuContext) -> Vec<[f32; 4]> { + read_texture(ctx, &self.radiance, wgpu::TextureAspect::All, 16) + } + + /// Reads the depth target back, row by row from the top. Slow, like + /// [`FrameTargets::read_radiance`]. + pub fn read_depth(&self, ctx: &GpuContext) -> Vec { + read_texture(ctx, &self.depth, wgpu::TextureAspect::DepthOnly, 4) + } +} + +/// Copies a whole texture into a buffer, waits, and unpacks it into texels of +/// `texel_bytes` each. +fn read_texture( + ctx: &GpuContext, + texture: &wgpu::Texture, + aspect: wgpu::TextureAspect, + texel_bytes: u32, +) -> Vec { + let (width, height) = (texture.width(), texture.height()); + let row = width * texel_bytes; + // Buffer rows must start on 256-byte boundaries. + let padded = row.next_multiple_of(wgpu::COPY_BYTES_PER_ROW_ALIGNMENT); + let buffer = ctx.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("readback"), + size: padded as u64 * height as u64, + usage: wgpu::BufferUsages::MAP_READ | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }); + let mut encoder = ctx.device.create_command_encoder(&Default::default()); + encoder.copy_texture_to_buffer( + wgpu::TexelCopyTextureInfo { + texture, + mip_level: 0, + origin: wgpu::Origin3d::ZERO, + aspect, + }, + wgpu::TexelCopyBufferInfo { + buffer: &buffer, + layout: wgpu::TexelCopyBufferLayout { + offset: 0, + bytes_per_row: Some(padded), + rows_per_image: Some(height), + }, + }, + texture.size(), + ); + ctx.queue.submit([encoder.finish()]); + buffer + .slice(..) + .map_async(wgpu::MapMode::Read, |r| r.expect("readback map")); + ctx.device + .poll(wgpu::PollType::wait_indefinitely()) + .expect("readback poll"); + let view = buffer.slice(..).get_mapped_range().expect("readback range"); + view.chunks(padded as usize) + .flat_map(|r| { + bytemuck::cast_slice::(&r[..row as usize]) + .iter() + .copied() + }) + .collect() +} diff --git a/crates/resolved/src/lib.rs b/crates/resolved/src/lib.rs index 086ed74..a857043 100644 --- a/crates/resolved/src/lib.rs +++ b/crates/resolved/src/lib.rs @@ -17,6 +17,7 @@ pub use nalgebra; /// [`gpu::GpuContext`] hands out wgpu types. pub use wgpu; +pub mod bodies; pub mod camera; pub mod error; pub mod gpu; diff --git a/crates/resolved/src/scene.rs b/crates/resolved/src/scene.rs index fc74d80..91603a4 100644 --- a/crates/resolved/src/scene.rs +++ b/crates/resolved/src/scene.rs @@ -30,6 +30,26 @@ impl Pose { Self { position, attitude } } + /// A camera pose at `position` with its boresight toward `target` and the image's + /// up direction (-Y in the camera frame) as close to `up` as possible. + /// + /// Returns `None` if `target` is at `position` or straight along `up`. + pub fn looking_at( + position: Vector3, + target: Vector3, + up: Vector3, + ) -> Option { + let z = (target - position).try_normalize(0.0)?; + let y = -(up - up.dot(&z) * z).try_normalize(1e-12)?; + let x = y.cross(&z); + let rotation = + nalgebra::Rotation3::from_matrix_unchecked(nalgebra::Matrix3::from_columns(&[x, y, z])); + Some(Self::new( + position, + UnitQuaternion::from_rotation_matrix(&rotation), + )) + } + /// The pose of a frame mounted on this body at `mount`, which is expressed in this /// body's frame. /// @@ -150,6 +170,20 @@ mod tests { assert_relative_eq!(camera.attitude, body.attitude, epsilon = 1e-12); } + #[test] + fn looking_at_points_boresight_and_keeps_up_up() { + let position = Vector3::new(7.0e6, 0.0, 0.0); + let camera = Pose::looking_at(position, Vector3::zeros(), Vector3::z()).unwrap(); + let axis = |v: Vector3| camera.attitude * v; + assert_relative_eq!(axis(Vector3::z()), -Vector3::x(), epsilon = 1e-12); + // Image up is camera -Y. + assert_relative_eq!(axis(-Vector3::y()), Vector3::z(), epsilon = 1e-12); + // Image right is boresight x up: (-X) x Z = +Y. + assert_relative_eq!(axis(Vector3::x()), Vector3::y(), epsilon = 1e-12); + assert!(Pose::looking_at(position, position, Vector3::z()).is_none()); + assert!(Pose::looking_at(position, Vector3::zeros(), Vector3::x()).is_none()); + } + #[test] fn boresight_target_lands_on_positive_z() { // Camera at 7000 km on +X, looking back at Earth: camera +Z is inertial -X. diff --git a/crates/resolved/src/shaders/bodies.wgsl b/crates/resolved/src/shaders/bodies.wgsl new file mode 100644 index 0000000..9a7e7f6 --- /dev/null +++ b/crates/resolved/src/shaders/bodies.wgsl @@ -0,0 +1,111 @@ +// Earth, Moon and Sun, intersected analytically per pixel. +// +// Everything is in the camera frame (+Z boresight, +X right, +Y down), in meters, with +// the camera at the origin. Each pixel casts one ray through its center, takes the +// nearest hit, writes that body's channel radiance, and writes depth as near / z so the +// target mesh drawn later occludes and is occluded correctly. + +struct Uniforms { + // fx, fy, cx, cy in pixels. + intrinsics: vec4, + // x: near clip distance. + params: vec4, + // xyz: Earth center. w: equatorial radius. + earth_center: vec4, + // Columns: Earth-fixed axes in the camera frame. + earth_axes: mat3x3, + // x: equatorial / polar radius. + earth_shape: vec4, + // xyz: unit vector from Earth toward the Sun. + earth_sun: vec4, + // Per channel: channel radiance of the surface facing the Sun (albedo S_k / (pi d^2)). + earth_lambert: vec4, + // xyz: Moon center. w: radius. + moon_center: vec4, + moon_sun: vec4, + moon_lambert: vec4, + // xyz: Sun center. w: radius. + sun_center: vec4, + // Per channel: channel radiance of the solar disk. + sun_radiance: vec4, +} + +@group(0) @binding(0) var u: Uniforms; + +struct Out { + @location(0) radiance: vec4, + @builtin(frag_depth) depth: f32, +} + +const MISS: f32 = -1.0; + +@vertex +fn vs_main(@builtin(vertex_index) i: u32) -> @builtin(position) vec4 { + // One triangle that covers the screen. + let p = vec2(f32((i << 1u) & 2u), f32(i & 2u)); + return vec4(p * 2.0 - 1.0, 0.0, 1.0); +} + +// Distance along unit ray `d` from `o` to a sphere of radius `r` at the origin, or MISS. +// +// The perpendicular miss distance comes from the point of closest approach, and +// r^2 - |q|^2 is factored as (r - |q|)(r + |q|). Both avoid subtracting two large +// nearly equal numbers, which in f32 would lose the limb of a body thousands of +// kilometers away. +fn sphere_hit(o: vec3, d: vec3, r: f32) -> f32 { + let b = dot(o, d); + let q = o - b * d; + let lq = length(q); + if lq > r { + return MISS; + } + let t = -b - sqrt((r - lq) * (r + lq)); + return select(MISS, t, t > 0.0); +} + +@fragment +fn fs_main(@builtin(position) pos: vec4) -> Out { + let k = u.intrinsics; + // `pos.xy` is the pixel center, (i + 0.5, j + 0.5). + let d = normalize(vec3((pos.x - k.z) / k.x, (pos.y - k.w) / k.y, 1.0)); + + var best = 3.4e38; + var radiance = vec4(0.0); + + // Earth: in Earth-fixed axes, stretch Z by a / b so the ellipsoid becomes a sphere + // of radius a. Distances along the stretched ray are rescaled by its length. + let to_cam = transpose(u.earth_axes) * (-u.earth_center.xyz); + let dir = transpose(u.earth_axes) * d; + let stretch = vec3(1.0, 1.0, u.earth_shape.x); + let ds = dir * stretch; + let ls = length(ds); + let te = sphere_hit(to_cam * stretch, ds / ls, u.earth_center.w) / ls; + if te > 0.0 && te < best { + best = te; + let p = to_cam + te * dir; + let n_fixed = normalize(p * stretch * stretch); + let n = u.earth_axes * n_fixed; + radiance = u.earth_lambert * max(dot(n, u.earth_sun.xyz), 0.0); + } + + let tm = sphere_hit(-u.moon_center.xyz, d, u.moon_center.w); + if tm > 0.0 && tm < best { + best = tm; + let n = normalize(tm * d - u.moon_center.xyz); + radiance = u.moon_lambert * max(dot(n, u.moon_sun.xyz), 0.0); + } + + let ts = sphere_hit(-u.sun_center.xyz, d, u.sun_center.w); + if ts > 0.0 && ts < best { + best = ts; + radiance = u.sun_radiance; + } + + if best > 3.0e38 { + discard; + } + var out: Out; + out.radiance = radiance; + out.depth = u.params.x / (best * d.z); + return out; +} diff --git a/crates/resolved/tests/bodies.rs b/crates/resolved/tests/bodies.rs new file mode 100644 index 0000000..c0fc5a0 --- /dev/null +++ b/crates/resolved/tests/bodies.rs @@ -0,0 +1,232 @@ +//! The bodies pass against closed-form geometry and radiometry, on the GPU and on the +//! CPU (lavapipe). + +use approx::assert_relative_eq; +use resolved::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS, MOON_RADIUS}; +use resolved::gpu::bodies::{BodiesPass, BodiesUniforms, BodyRadiometry}; +use resolved::gpu::targets::FrameTargets; +use resolved::gpu::{AdapterChoice, GpuContext}; +use resolved::nalgebra::{UnitQuaternion, Vector3}; +use resolved::radiometry::{Radiometry, AU, SUN_RADIUS}; +use resolved::{CameraConfig, Pose, Scene}; + +/// Camera distance from Earth's center: far enough that Earth fits in the example +/// camera's 6.8 degree height. +const RANGE: f64 = 2.0e8; + +struct Frame { + radiance: Vec<[f32; 4]>, + depth: Vec, + width: usize, +} + +impl Frame { + fn at(&self, u: usize, v: usize) -> ([f32; 4], f32) { + let i = v * self.width + u; + (self.radiance[i], self.depth[i]) + } +} + +fn contexts() -> [GpuContext; 2] { + [AdapterChoice::Gpu, AdapterChoice::Cpu].map(|c| GpuContext::new(&c).unwrap()) +} + +fn body_radiometry() -> BodyRadiometry { + let radiometry = Radiometry::new(&CameraConfig::example()); + BodyRadiometry::new(&radiometry, &Appearance::default()).unwrap() +} + +fn render(ctx: &GpuContext, scene: &Scene, camera: &Pose) -> Frame { + render_with(ctx, scene, camera, &CameraConfig::example()) +} + +fn render_with(ctx: &GpuContext, scene: &Scene, camera: &Pose, config: &CameraConfig) -> Frame { + let targets = FrameTargets::new(ctx, config.width, config.height); + let pass = BodiesPass::new(ctx); + let uniforms = BodiesUniforms::new( + scene, + camera, + &config.intrinsics(), + config.near_clip, + &body_radiometry(), + ); + let mut encoder = ctx.device.create_command_encoder(&Default::default()); + pass.encode(ctx, &mut encoder, &targets, &uniforms); + ctx.queue.submit([encoder.finish()]); + Frame { + radiance: targets.read_radiance(ctx), + depth: targets.read_depth(ctx), + width: config.width as usize, + } +} + +fn scene(sun: Vector3, moon: Vector3) -> Scene { + Scene { + earth_attitude: UnitQuaternion::identity(), + sun, + moon, + targets: vec![], + } +} + +/// Camera on +X at `RANGE`, looking at Earth's center, north up. +fn earth_camera() -> Pose { + Pose::looking_at( + Vector3::new(RANGE, 0.0, 0.0), + Vector3::zeros(), + Vector3::z(), + ) + .unwrap() +} + +/// The Moon somewhere out of view. +fn moon_away() -> Vector3 { + Vector3::new(0.0, -3.844e8, 0.0) +} + +#[test] +fn sunlit_earth_center_radiance_and_depth() { + // The Sun is behind the camera, so the point under the boresight faces both. + let scene = scene(Vector3::new(AU, 0.0, 0.0), moon_away()); + let expected = body_radiometry().earth_lambert_1au; + let near = CameraConfig::example().near_clip; + for ctx in contexts() { + let frame = render(&ctx, &scene, &earth_camera()); + let (radiance, depth) = frame.at(960, 540); + for k in 0..3 { + assert_relative_eq!(radiance[k], expected[k], max_relative = 1e-4); + } + assert_eq!(radiance[3], 0.0); + let z = RANGE - EARTH_EQUATORIAL_RADIUS; + assert_relative_eq!(depth as f64, near / z, max_relative = 1e-5); + } +} + +#[test] +fn earth_limb_lands_where_projection_says() { + let scene = scene(Vector3::new(AU, 0.0, 0.0), moon_away()); + let k = CameraConfig::example().intrinsics(); + // Viewed from the equatorial plane, the silhouette's horizontal extent is the + // equatorial circle: half-width fx tan(asin(a / range)) pixels. + let half_width = k.fx * (EARTH_EQUATORIAL_RADIUS / RANGE).asin().tan(); + let expected = (0..1920) + .filter(|u| ((*u as f64 + 0.5) - k.cx).abs() < half_width) + .count(); + for ctx in contexts() { + let frame = render(&ctx, &scene, &earth_camera()); + let hit = (0..frame.width) + .filter(|&u| frame.at(u, 540).1 > 0.0) + .count(); + assert!( + hit.abs_diff(expected) <= 2, + "{}: {hit} px across, expected {expected}", + ctx.info.name + ); + } +} + +#[test] +fn night_side_is_dark_and_hides_the_sun() { + // The Sun directly behind Earth: the visible hemisphere is unlit and the solar disk + // is eclipsed. + let scene = scene(Vector3::new(-AU, 0.0, 0.0), moon_away()); + for ctx in contexts() { + let frame = render(&ctx, &scene, &earth_camera()); + assert!( + frame.radiance.iter().all(|r| *r == [0.0; 4]), + "{}", + ctx.info.name + ); + assert!(frame.at(960, 540).1 > 0.0, "Earth still occupies depth"); + assert_eq!(frame.at(0, 0).1, 0.0, "empty space stays at depth 0"); + } +} + +#[test] +fn solar_disk_size_and_radiance() { + let sun = Vector3::new(0.0, AU, 0.0); + let scene = scene(sun, moon_away()); + let position = Vector3::new(RANGE, 0.0, 0.0); + let camera = Pose::looking_at(position, sun, Vector3::z()).unwrap(); + let k = CameraConfig::example().intrinsics(); + let radius = k.fx * (SUN_RADIUS / (sun - position).norm()).asin().tan(); + let expected = std::f64::consts::PI * radius * radius; + let sun_radiance = body_radiometry().sun_radiance; + for ctx in contexts() { + let frame = render(&ctx, &scene, &camera); + assert_eq!(frame.at(960, 540).0, sun_radiance); + let lit = frame + .radiance + .iter() + .filter(|r| **r == sun_radiance) + .count() as f64; + assert_relative_eq!(lit, expected, max_relative = 0.01); + } +} + +#[test] +fn moon_in_front_of_earth_wins() { + let moon = Vector3::new(RANGE / 2.0, 0.0, 0.0); + let sun = Vector3::new(AU, 0.0, 0.0); + let scene = scene(sun, moon); + let near = CameraConfig::example().near_clip; + let expected = { + let au = (sun - moon).norm() / AU; + body_radiometry() + .moon_lambert_1au + .map(|l| l / (au * au) as f32) + }; + for ctx in contexts() { + let frame = render(&ctx, &scene, &earth_camera()); + let (radiance, depth) = frame.at(960, 540); + for k in 0..3 { + assert_relative_eq!(radiance[k], expected[k], max_relative = 1e-4); + } + let z = RANGE / 2.0 - MOON_RADIUS; + assert_relative_eq!(depth as f64, near / z, max_relative = 1e-5); + // Earth (290 px radius) still shows around the Moon (158 px), and is farther. + let earth = frame.at(960 + 200, 540).1; + assert!(earth > 0.0 && earth < depth); + } +} + +#[test] +fn limb_follows_the_ellipsoid() { + // A 1 m focal length makes a pixel 5.5 urad, so at this range Earth's 21 km of + // flattening is 19 px: a spherical Earth would put the north limb 19 px off. + let mut config = CameraConfig::example(); + config.focal_length = 1.0; + let scene = scene(Vector3::new(AU, 0.0, 0.0), moon_away()); + let position = Vector3::new(RANGE, 0.0, 0.0); + // From a point on the equatorial plane, the line of sight that grazes the meridian + // ellipse x^2/a^2 + z^2/b^2 = 1 touches it at x = a^2 / D, z = b sqrt(1 - a^2/D^2). + let (a, b) = (EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS); + let x = a * a / RANGE; + let north = Vector3::new(x, 0.0, b * (1.0 - (a / RANGE).powi(2)).sqrt()); + let east = Vector3::new(x, a * (1.0 - (a / RANGE).powi(2)).sqrt(), 0.0); + for ctx in contexts() { + // Boresight on the north limb, north up: Earth fills rows below the center. + let camera = Pose::looking_at(position, north, Vector3::z()).unwrap(); + let frame = render_with(&ctx, &scene, &camera, &config); + let column: Vec = (0..1080).map(|v| frame.at(960, v).1 > 0.0).collect(); + let edge = column.iter().position(|hit| *hit).unwrap(); + assert!(column[edge..].iter().all(|hit| *hit)); + assert!( + edge.abs_diff(540) <= 1, + "{}: north limb at row {edge}", + ctx.info.name + ); + + // Boresight on the east limb, north up: Earth is to the left of center. + let camera = Pose::looking_at(position, east, Vector3::z()).unwrap(); + let frame = render_with(&ctx, &scene, &camera, &config); + let row: Vec = (0..1920).map(|u| frame.at(u, 540).1 > 0.0).collect(); + let edge = row.iter().position(|hit| !*hit).unwrap(); + assert!(row[edge..].iter().all(|hit| !*hit)); + assert!( + edge.abs_diff(960) <= 1, + "{}: east limb at column {edge}", + ctx.info.name + ); + } +} diff --git a/docs/optics.md b/docs/optics.md index 1698de1..6c949f2 100644 --- a/docs/optics.md +++ b/docs/optics.md @@ -313,17 +313,75 @@ Picking the exposure time is the same trade a photographer makes. A longer exposure saturates the sunlit target. A shorter one leaves the side of the target lit only by Earthshine lost in noise once noise is modeled (M3). +## Earth, Moon and Sun + +Code: `bodies.rs`, `gpu/bodies.rs`, `shaders/bodies.wgsl`. + +The three bodies are not drawn as triangle meshes. For every pixel, the +shader sends a ray from the camera through the pixel's center and solves +for where it meets each body's exact shape. The nearest hit wins. This +gives an exact edge at any zoom, at the cost of a few arithmetic operations +per pixel. A mesh would show its facets along the limb, the edge of the +disk, which is the feature horizon-based navigation measures. + +- **Earth** is the WGS84 ellipsoid: 6378.137 km equatorial radius, + 6356.752 km polar. The 21 km difference matters at fine pixel scales. A + camera with 5.5 urad pixels at 200,000 km sees the north limb 19 pixels + away from where a spherical Earth would put it. A test checks the limb + against the exact tangent to the ellipsoid to within a pixel. +- **The Moon** is a sphere of 1737.4 km radius. From Earth's distance it + is about 0.52 degrees across, 82 pixels in the example camera. +- **The Sun** is a sphere of 695,700 km radius: 0.53 degrees across, about + 85 pixels in the example camera. Its disk has the solar radiance from the + folding section above, and it is drawn only where no nearer body covers + it. So Earth eclipses it naturally. + +### Shading + +Earth and the Moon are Lambertian, with one albedo each, set by +`Appearance`: 0.3 for Earth and 0.12 for the Moon, flat across wavelength. +Each point's radiance is rho cos(i) S_k / (pi d^2), with i the angle between +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. + +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: + +- Earth's oceans, land and clouds differ in brightness by a factor of + several. An albedo map comes in M2. +- The full Moon looks almost uniformly bright out to its edge. A Lambertian + sphere instead darkens toward the edge, because the light there arrives + and leaves at grazing angles. Lunar soil scatters light back toward its + source much more strongly than Lambert predicts, and a dedicated lunar + reflectance model (M4) fixes that. The Moon's maria (the dark plains of + old basalt lava, the "seas" seen with the naked eye) and its bright + highlands also differ in albedo by about a factor of two (M2 albedo + map). +- Earth's atmosphere softens the limb and adds a thin bright band around it + (M4). Here Earth's edge is a sharp line. + +### Depth from the bodies + +The pass also writes each hit's depth, near / z, as described in the +geometry section. The target spacecraft, drawn after the bodies, uses the +same depth buffer. So a target in front of Earth covers it, and a target +behind Earth's limb is hidden. + ## Not modeled yet | Effect | Where | |---|---| | Light from Earth and the Moon illuminating the target (Earthshine, Moonshine) | M2 | +| Earth's shadow on the Moon (lunar eclipses) | M2 | +| Earthshine on the Moon's night side | M4 | | Shadows the target casts on itself | M2 | | The Sun as a disk that is darker toward its edge (limb darkening) | M2 | | Lens distortion | M3 | | Blur from diffraction and lens imperfections (point spread function) | M3 | | Noise: photon shot noise, read noise, dark current, pixel-to-pixel variation | M3 | | Saturation spilling into neighbors (blooming) | M3 | -| Non-Lambertian materials (BRDFs), glints | M4 | +| Non-Lambertian materials (BRDFs), glints, lunar backscatter | M4 | +| Earth and Moon albedo maps | M2 | | Earth's atmosphere: scattering, the bright limb | M4 | | Motion blur, rolling shutter, stray light | M6 |