diff --git a/Cargo.lock b/Cargo.lock index f41acda..7847574 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -963,6 +963,7 @@ dependencies = [ "approx", "bytemuck", "gltf", + "image", "log", "nalgebra", "pollster", diff --git a/crates/resolved/Cargo.toml b/crates/resolved/Cargo.toml index 44f4125..0e1e34e 100644 --- a/crates/resolved/Cargo.toml +++ b/crates/resolved/Cargo.toml @@ -23,3 +23,4 @@ stl_io = "0.11" [dev-dependencies] approx = "0.5" +image = { version = "0.25", default-features = false, features = ["png"] } diff --git a/crates/resolved/src/eclipse.rs b/crates/resolved/src/eclipse.rs new file mode 100644 index 0000000..cac6e66 --- /dev/null +++ b/crates/resolved/src/eclipse.rs @@ -0,0 +1,153 @@ +//! How much of the Sun a point can see past Earth and the Moon. +//! +//! The Sun is a disk 0.53 degrees across, so a body passing in front of it dims the +//! light gradually: a penumbra where part of the disk is covered, an umbra where all of +//! it is, and an annular case where a small occulter sits inside the disk. All three +//! come from the overlap of two circles on the sky, the same model as astrokit's +//! `eclipse::illumination` (Montenbruck and Gill, *Satellite Orbits*, 3.4.2). + +use nalgebra::Vector3; + +use crate::radiometry::SUN_RADIUS; + +/// The fraction of the solar disk visible from `point`, between 0 and 1, with each +/// occulter given as (center, radius). Positions are in meters in any common frame. +/// +/// Occulters are assumed not to overlap each other on the sky, which holds for Earth +/// and the Moon except during the rare few minutes of a solar eclipse seen from orbit. +pub fn sun_visible_fraction( + point: Vector3, + sun: Vector3, + occulters: &[(Vector3, f64)], +) -> f64 { + let covered: f64 = occulters + .iter() + .map(|(center, radius)| covered_fraction(point, sun, *center, *radius)) + .sum(); + (1.0 - covered).clamp(0.0, 1.0) +} + +/// The fraction of the solar disk one occulter covers, seen from `point`. +fn covered_fraction( + point: Vector3, + sun: Vector3, + center: Vector3, + radius: f64, +) -> f64 { + let to_sun = sun - point; + let to_body = center - point; + let a = (SUN_RADIUS / to_sun.norm()).min(1.0).asin(); + let b = (radius / to_body.norm()).min(1.0).asin(); + // atan2 keeps precision near zero separation, where acos of a dot product loses it. + let c = to_sun.cross(&to_body).norm().atan2(to_sun.dot(&to_body)); + if c >= a + b { + 0.0 + } else if c <= b - a { + 1.0 + } else if c <= a - b { + (b / a).powi(2) + } else { + // A lens-shaped overlap of the two disks. + let x = (c * c + a * a - b * b) / (2.0 * c); + let y = (a * a - x * x).max(0.0).sqrt(); + let area = a * a * (x / a).clamp(-1.0, 1.0).acos() + + b * b * ((c - x) / b).clamp(-1.0, 1.0).acos() + - c * y; + area / (std::f64::consts::PI * a * a) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::bodies::EARTH_EQUATORIAL_RADIUS; + use crate::radiometry::AU; + use approx::assert_relative_eq; + + const RE: f64 = EARTH_EQUATORIAL_RADIUS; + + fn earth() -> [(Vector3, f64); 1] { + [(Vector3::zeros(), RE)] + } + + /// Independent check: the fraction of points spread evenly over the solar disk + /// (perpendicular to the line of sight) whose line to `point` misses Earth. + fn sampled(point: Vector3, sun: Vector3) -> f64 { + let axis = (sun - point).normalize(); + let u = axis.cross(&Vector3::z()).normalize(); + let v = axis.cross(&u); + let n = 150; + let (mut inside, mut visible) = (0, 0); + for i in 0..n { + for j in 0..n { + let (s, t) = ( + 2.0 * (i as f64 + 0.5) / n as f64 - 1.0, + 2.0 * (j as f64 + 0.5) / n as f64 - 1.0, + ); + if s * s + t * t > 1.0 { + continue; + } + inside += 1; + let target = sun + SUN_RADIUS * (s * u + t * v); + let d = (target - point).normalize(); + // Does the ray from `point` toward `target` hit Earth? + let b = point.dot(&d); + let q = point - b * d; + let blocked = b < 0.0 && q.norm() < RE; + visible += usize::from(!blocked); + } + } + visible as f64 / inside as f64 + } + + #[test] + fn sunlit_and_umbra() { + let sun = Vector3::new(AU, 0.0, 0.0); + assert_eq!( + sun_visible_fraction(Vector3::new(7.0e6, 0.0, 0.0), sun, &earth()), + 1.0 + ); + assert_eq!( + sun_visible_fraction(Vector3::new(-7.0e6, 0.0, 0.0), sun, &earth()), + 0.0 + ); + } + + #[test] + fn penumbra_matches_sampling_the_disk() { + // 7000 km behind Earth's center, moved sideways through the shadow edge. The + // penumbra there runs from about 32 km inside Earth's radius to 33 km outside. + let sun = Vector3::new(AU, 0.0, 0.0); + let mut seen = vec![]; + for dy in [-20e3, -10e3, -5e3, 0.0, 5e3, 10e3, 20e3] { + let point = Vector3::new(-7.0e6, RE + dy, 0.0); + let f = sun_visible_fraction(point, sun, &earth()); + assert_relative_eq!(f, sampled(point, sun), epsilon = 0.005); + seen.push(f); + } + assert!( + seen.windows(2).all(|w| w[0] <= w[1]), + "brightens outward: {seen:?}" + ); + assert!( + seen[0] < 0.5 && seen[6] > 0.5, + "spans the penumbra: {seen:?}" + ); + } + + #[test] + fn annular_takes_a_bite() { + // A body small enough to sit inside the solar disk: covered fraction is the + // ratio of apparent areas. + let sun = Vector3::new(AU, 0.0, 0.0); + let point = Vector3::zeros(); + let body = (Vector3::new(1.0e9, 0.0, 0.0), 1.0e6); + let a = (SUN_RADIUS / AU).asin(); + let b = (1.0e6f64 / 1.0e9).asin(); + assert_relative_eq!( + sun_visible_fraction(point, sun, &[body]), + 1.0 - (b / a).powi(2), + max_relative = 1e-12 + ); + } +} diff --git a/crates/resolved/src/gpu/mesh.rs b/crates/resolved/src/gpu/mesh.rs new file mode 100644 index 0000000..054e36c --- /dev/null +++ b/crates/resolved/src/gpu/mesh.rs @@ -0,0 +1,332 @@ +//! The pass that draws target meshes over the bodies. +//! +//! A model is uploaded once ([`MeshPass::upload`]) and drawn every frame with a small +//! per-target uniform block ([`MeshDraw`]) that places it and says how much sunlight +//! reaches it. + +use nalgebra::{Matrix3, Vector3}; + +use super::bodies::MAX_CHANNELS; +use super::targets::{FrameTargets, DEPTH_COMPARE, DEPTH_FORMAT, RADIANCE_FORMAT}; +use super::GpuContext; +use crate::camera::Intrinsics; +use crate::model::Mesh; +use crate::radiometry::{Radiometry, AU}; +use crate::scene::Pose; +use crate::Error; + +#[repr(C)] +#[derive(Clone, Copy, Debug, bytemuck::Pod, bytemuck::Zeroable)] +struct Vertex { + position: [f32; 3], + normal: [f32; 3], + part: u32, +} + +/// A model on the GPU, ready to draw. +#[derive(Debug)] +pub struct GpuMesh { + vertices: wgpu::Buffer, + indices: wgpu::Buffer, + index_count: u32, + parts: wgpu::BindGroup, +} + +/// One target's uniform block; the layout matches `Draw` in `mesh.wgsl`. +#[repr(C)] +#[derive(Clone, Copy, Debug, PartialEq, bytemuck::Pod, bytemuck::Zeroable)] +pub struct MeshDraw { + clip_from_camera: [[f32; 4]; 4], + // mat3x3 columns, each padded to 16 bytes. + camera_from_body: [[f32; 4]; 3], + translation: [f32; 4], + sun: [f32; 4], + params: [f32; 4], +} + +impl MeshDraw { + /// The uniforms to draw a model at `target` seen from `camera`, lit by the Sun at + /// `sun` (inertial, meters) with fraction `sun_visible` of its disk unobstructed. + pub fn new( + target: &Pose, + camera: &Pose, + intrinsics: &Intrinsics, + near: f64, + sun: Vector3, + sun_visible: f64, + ) -> Self { + let to_camera = camera.attitude.inverse(); + let rotation: Matrix3 = (to_camera * target.attitude) + .to_rotation_matrix() + .matrix() + .cast(); + let translation = (to_camera * (target.position - camera.position)).cast::(); + let to_sun = sun - target.position; + let sun_dir = (to_camera * to_sun.normalize()).cast::(); + let au = to_sun.norm() / AU; + let clip = intrinsics.clip_from_camera(near).cast::(); + let col3 = |i: usize| { + let c = rotation.column(i); + [c.x, c.y, c.z, 0.0] + }; + let col4 = |i: usize| { + let c = clip.column(i); + [c.x, c.y, c.z, c.w] + }; + Self { + clip_from_camera: [col4(0), col4(1), col4(2), col4(3)], + camera_from_body: [col3(0), col3(1), col3(2)], + translation: [translation.x, translation.y, translation.z, 0.0], + sun: [sun_dir.x, sun_dir.y, sun_dir.z, 0.0], + params: [(sun_visible / (au * au)) as f32, 0.0, 0.0, 0.0], + } + } +} + +/// The compiled mesh pipeline and a uniform buffer with one slot per target. +#[derive(Debug)] +pub struct MeshPass { + pipeline: wgpu::RenderPipeline, + draw_layout: wgpu::BindGroupLayout, + parts_layout: wgpu::BindGroupLayout, + draws: wgpu::Buffer, + draws_group: wgpu::BindGroup, + stride: u64, + capacity: usize, +} + +impl MeshPass { + /// Compiles the shader and allocates room for one target. + pub fn new(ctx: &GpuContext) -> Self { + let device = &ctx.device; + let draw_layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { + label: Some("mesh draw"), + entries: &[wgpu::BindGroupLayoutEntry { + binding: 0, + visibility: wgpu::ShaderStages::VERTEX_FRAGMENT, + ty: wgpu::BindingType::Buffer { + ty: wgpu::BufferBindingType::Uniform, + has_dynamic_offset: true, + min_binding_size: wgpu::BufferSize::new(std::mem::size_of::() as u64), + }, + count: None, + }], + }); + let parts_layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { + label: Some("mesh parts"), + entries: &[wgpu::BindGroupLayoutEntry { + binding: 0, + visibility: wgpu::ShaderStages::FRAGMENT, + ty: wgpu::BindingType::Buffer { + ty: wgpu::BufferBindingType::Storage { read_only: true }, + has_dynamic_offset: false, + min_binding_size: None, + }, + count: None, + }], + }); + let layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor { + label: Some("mesh"), + bind_group_layouts: &[Some(&draw_layout), Some(&parts_layout)], + immediate_size: 0, + }); + let module = device.create_shader_module(wgpu::ShaderModuleDescriptor { + label: Some("mesh"), + source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/mesh.wgsl").into()), + }); + let pipeline = device.create_render_pipeline(&wgpu::RenderPipelineDescriptor { + label: Some("mesh"), + layout: Some(&layout), + vertex: wgpu::VertexState { + module: &module, + entry_point: Some("vs_main"), + compilation_options: Default::default(), + buffers: &[Some(wgpu::VertexBufferLayout { + array_stride: std::mem::size_of::() as u64, + step_mode: wgpu::VertexStepMode::Vertex, + attributes: &wgpu::vertex_attr_array![0 => Float32x3, 1 => Float32x3, 2 => Uint32], + })], + }, + // No culling: thin panels are seen from both sides. + 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 alignment = device.limits().min_uniform_buffer_offset_alignment as u64; + let stride = (std::mem::size_of::() as u64).next_multiple_of(alignment); + let (draws, draws_group) = Self::draw_slots(ctx, &draw_layout, stride, 1); + Self { + pipeline, + draw_layout, + parts_layout, + draws, + draws_group, + stride, + capacity: 1, + } + } + + fn draw_slots( + ctx: &GpuContext, + layout: &wgpu::BindGroupLayout, + stride: u64, + capacity: usize, + ) -> (wgpu::Buffer, wgpu::BindGroup) { + let buffer = ctx.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("mesh draws"), + size: stride * capacity as u64, + usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }); + let group = ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("mesh draws"), + layout, + entries: &[wgpu::BindGroupEntry { + binding: 0, + resource: wgpu::BindingResource::Buffer(wgpu::BufferBinding { + buffer: &buffer, + offset: 0, + size: wgpu::BufferSize::new(std::mem::size_of::() as u64), + }), + }], + }); + (buffer, group) + } + + /// Uploads `mesh`, folding each part's material through `radiometry`'s channels. + pub fn upload( + &self, + ctx: &GpuContext, + mesh: &Mesh, + radiometry: &Radiometry, + ) -> Result { + use wgpu::util::DeviceExt; + 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: Vec<[f32; 4]> = mesh + .materials + .iter() + .map(|m| { + let mut out = [0.0; 4]; + for (k, a) in radiometry.effective_albedo(&m.albedo).iter().enumerate() { + out[k] = radiometry.lambert_radiance(k, *a, 1.0, AU) as f32; + } + out + }) + .collect(); + let parts: Vec<[f32; 4]> = mesh.parts.iter().map(|p| lambert[p.material]).collect(); + let vertices: Vec = (0..mesh.positions.len()) + .map(|i| Vertex { + position: mesh.positions[i], + normal: mesh.normals[i], + part: mesh.part_of[i], + }) + .collect(); + let buffer = |label, contents: &[u8], usage| { + ctx.device + .create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some(label), + contents, + usage, + }) + }; + let parts = buffer( + "mesh parts", + bytemuck::cast_slice(&parts), + wgpu::BufferUsages::STORAGE, + ); + Ok(GpuMesh { + vertices: buffer( + "mesh vertices", + bytemuck::cast_slice(&vertices), + wgpu::BufferUsages::VERTEX, + ), + indices: buffer( + "mesh indices", + bytemuck::cast_slice(&mesh.indices), + wgpu::BufferUsages::INDEX, + ), + index_count: mesh.indices.len() as u32, + parts: ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("mesh parts"), + layout: &self.parts_layout, + entries: &[wgpu::BindGroupEntry { + binding: 0, + resource: parts.as_entire_binding(), + }], + }), + }) + } + + /// Uploads each target's uniforms and records the pass into `encoder`, drawing over + /// what earlier passes left in `targets`. + pub fn encode( + &mut self, + ctx: &GpuContext, + encoder: &mut wgpu::CommandEncoder, + targets: &FrameTargets, + draws: &[(&GpuMesh, MeshDraw)], + ) { + if draws.len() > self.capacity { + self.capacity = draws.len().next_power_of_two(); + (self.draws, self.draws_group) = + Self::draw_slots(ctx, &self.draw_layout, self.stride, self.capacity); + } + for (i, (_, draw)) in draws.iter().enumerate() { + ctx.queue.write_buffer( + &self.draws, + i as u64 * self.stride, + bytemuck::bytes_of(draw), + ); + } + 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("mesh"), + color_attachments: &[Some(wgpu::RenderPassColorAttachment { + view: &radiance, + depth_slice: None, + resolve_target: None, + ops: wgpu::Operations { + load: wgpu::LoadOp::Load, + store: wgpu::StoreOp::Store, + }, + })], + depth_stencil_attachment: Some(wgpu::RenderPassDepthStencilAttachment { + view: &depth, + depth_ops: Some(wgpu::Operations { + load: wgpu::LoadOp::Load, + store: wgpu::StoreOp::Store, + }), + stencil_ops: None, + }), + ..Default::default() + }); + pass.set_pipeline(&self.pipeline); + for (i, (mesh, _)) in draws.iter().enumerate() { + pass.set_bind_group(0, &self.draws_group, &[(i as u64 * self.stride) as u32]); + pass.set_bind_group(1, &mesh.parts, &[]); + pass.set_vertex_buffer(0, mesh.vertices.slice(..)); + pass.set_index_buffer(mesh.indices.slice(..), wgpu::IndexFormat::Uint32); + pass.draw_indexed(0..mesh.index_count, 0, 0..1); + } + } +} diff --git a/crates/resolved/src/gpu/mod.rs b/crates/resolved/src/gpu/mod.rs index 7ce9e77..1bd3bc6 100644 --- a/crates/resolved/src/gpu/mod.rs +++ b/crates/resolved/src/gpu/mod.rs @@ -2,6 +2,7 @@ pub mod bodies; pub mod context; +pub mod mesh; pub mod targets; pub use context::{adapters, AdapterChoice, GpuContext}; diff --git a/crates/resolved/src/lib.rs b/crates/resolved/src/lib.rs index 4f525a6..290d6b5 100644 --- a/crates/resolved/src/lib.rs +++ b/crates/resolved/src/lib.rs @@ -19,6 +19,7 @@ pub use wgpu; pub mod bodies; pub mod camera; +pub mod eclipse; pub mod error; pub mod gpu; pub mod model; diff --git a/crates/resolved/src/shaders/mesh.wgsl b/crates/resolved/src/shaders/mesh.wgsl new file mode 100644 index 0000000..47f2e2e --- /dev/null +++ b/crates/resolved/src/shaders/mesh.wgsl @@ -0,0 +1,56 @@ +// A target spacecraft mesh, rasterized over the bodies. +// +// Vertices arrive in the model's body frame and are placed in the camera frame by a +// rigid transform computed in f64 on the CPU. Depth uses the same reversed-Z, near / z +// projection the bodies pass wrote, so the depth test sorts mesh and bodies correctly. + +struct Draw { + clip_from_camera: mat4x4, + // Columns: the model's axes in the camera frame. + camera_from_body: mat3x3, + // xyz: the model origin in the camera frame. + translation: vec4, + // xyz: unit vector from the target toward the Sun, camera frame. + sun: vec4, + // x: Sun visibility (eclipse) / (Sun distance in AU)^2. + params: vec4, +} + +@group(0) @binding(0) var draw: Draw; +// Per part: channel radiance of a surface facing the Sun at 1 AU (albedo S_k / pi). +@group(1) @binding(0) var parts: array>; + +struct VertexIn { + @location(0) position: vec3, + @location(1) normal: vec3, + @location(2) part: u32, +} + +struct VertexOut { + @builtin(position) clip: vec4, + @location(0) camera_position: vec3, + @location(1) normal: vec3, + @location(2) @interpolate(flat) part: u32, +} + +@vertex +fn vs_main(v: VertexIn) -> VertexOut { + let p = draw.camera_from_body * v.position + draw.translation.xyz; + var out: VertexOut; + out.clip = draw.clip_from_camera * vec4(p, 1.0); + out.camera_position = p; + out.normal = draw.camera_from_body * v.normal; + out.part = v.part; + return out; +} + +@fragment +fn fs_main(in: VertexOut) -> @location(0) vec4 { + var n = normalize(in.normal); + // Lit from whichever side faces the camera, so a zero-thickness panel seen from + // behind shows its back. On a closed mesh, back faces are hidden by depth anyway. + if dot(n, in.camera_position) > 0.0 { + n = -n; + } + return parts[in.part] * max(dot(n, draw.sun.xyz), 0.0) * draw.params.x; +} diff --git a/crates/resolved/tests/bodies.rs b/crates/resolved/tests/bodies.rs index c0fc5a0..1cf5beb 100644 --- a/crates/resolved/tests/bodies.rs +++ b/crates/resolved/tests/bodies.rs @@ -1,6 +1,8 @@ //! The bodies pass against closed-form geometry and radiometry, on the GPU and on the //! CPU (lavapipe). +mod common; + use approx::assert_relative_eq; use resolved::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS, MOON_RADIUS}; use resolved::gpu::bodies::{BodiesPass, BodiesUniforms, BodyRadiometry}; @@ -36,11 +38,17 @@ fn body_radiometry() -> BodyRadiometry { BodyRadiometry::new(&radiometry, &Appearance::default()).unwrap() } -fn render(ctx: &GpuContext, scene: &Scene, camera: &Pose) -> Frame { - render_with(ctx, scene, camera, &CameraConfig::example()) +fn render(ctx: &GpuContext, name: &str, scene: &Scene, camera: &Pose) -> Frame { + render_with(ctx, name, scene, camera, &CameraConfig::example()) } -fn render_with(ctx: &GpuContext, scene: &Scene, camera: &Pose, config: &CameraConfig) -> Frame { +fn render_with( + ctx: &GpuContext, + name: &str, + 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( @@ -53,8 +61,10 @@ fn render_with(ctx: &GpuContext, scene: &Scene, camera: &Pose, config: &CameraCo let mut encoder = ctx.device.create_command_encoder(&Default::default()); pass.encode(ctx, &mut encoder, &targets, &uniforms); ctx.queue.submit([encoder.finish()]); + let radiance = targets.read_radiance(ctx); + common::save(name, ctx, &radiance); Frame { - radiance: targets.read_radiance(ctx), + radiance, depth: targets.read_depth(ctx), width: config.width as usize, } @@ -91,7 +101,7 @@ fn sunlit_earth_center_radiance_and_depth() { 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 frame = render(&ctx, "earth-sunlit", &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); @@ -113,7 +123,7 @@ fn earth_limb_lands_where_projection_says() { .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 frame = render(&ctx, "earth-limb", &scene, &earth_camera()); let hit = (0..frame.width) .filter(|&u| frame.at(u, 540).1 > 0.0) .count(); @@ -131,7 +141,7 @@ fn night_side_is_dark_and_hides_the_sun() { // 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()); + let frame = render(&ctx, "earth-night", &scene, &earth_camera()); assert!( frame.radiance.iter().all(|r| *r == [0.0; 4]), "{}", @@ -153,7 +163,7 @@ fn solar_disk_size_and_radiance() { 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); + let frame = render(&ctx, "sun-disk", &scene, &camera); assert_eq!(frame.at(960, 540).0, sun_radiance); let lit = frame .radiance @@ -177,7 +187,7 @@ fn moon_in_front_of_earth_wins() { .map(|l| l / (au * au) as f32) }; for ctx in contexts() { - let frame = render(&ctx, &scene, &earth_camera()); + let frame = render(&ctx, "moon-in-front-of-earth", &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); @@ -207,7 +217,7 @@ fn limb_follows_the_ellipsoid() { 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 frame = render_with(&ctx, "limb-north", &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)); @@ -219,7 +229,7 @@ fn limb_follows_the_ellipsoid() { // 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 frame = render_with(&ctx, "limb-east", &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)); diff --git a/crates/resolved/tests/common/mod.rs b/crates/resolved/tests/common/mod.rs new file mode 100644 index 0000000..bcc8bd4 --- /dev/null +++ b/crates/resolved/tests/common/mod.rs @@ -0,0 +1,54 @@ +//! Shared by the GPU integration tests: saving what a test rendered, for looking at. +//! +//! Set `RESOLVED_TEST_IMAGES=1` and each test writes PNGs to +//! `target/tmp/test-images/`, one per test case and device. + +use std::path::PathBuf; + +use resolved::gpu::GpuContext; +use resolved::radiometry::{Radiometry, AU}; +use resolved::CameraConfig; + +/// Writes `radiance` (example camera, row-major from the top) as an sRGB PNG named +/// `-.png`, if `RESOLVED_TEST_IMAGES` is set. +/// +/// Channels r, g, b map to red, green, blue, scaled so a white Lambertian surface +/// facing the Sun at 1 AU is full white. Brighter than that (the Sun) clips. This is a +/// display of radiance, not a simulated camera image. +pub fn save(name: &str, ctx: &GpuContext, radiance: &[[f32; 4]]) { + if std::env::var_os("RESOLVED_TEST_IMAGES").is_none() { + return; + } + let config = CameraConfig::example(); + let r = Radiometry::new(&config); + let white: Vec = (0..3) + .map(|k| r.lambert_radiance(k, 1.0, 1.0, AU) as f32) + .collect(); + let srgb = |x: f32| { + let x = x.clamp(0.0, 1.0); + let v = if x <= 0.003_130_8 { + 12.92 * x + } else { + 1.055 * x.powf(1.0 / 2.4) - 0.055 + }; + (v * 255.0).round() as u8 + }; + let pixels: Vec = radiance + .iter() + .flat_map(|p| (0..3).map(move |k| (k, p[k]))) + .map(|(k, v)| srgb(v / white[k])) + .collect(); + let dir = PathBuf::from(env!("CARGO_TARGET_TMPDIR")).join("test-images"); + std::fs::create_dir_all(&dir).unwrap(); + let device = if ctx.is_cpu() { "cpu" } else { "gpu" }; + let path = dir.join(format!("{name}-{device}.png")); + image::save_buffer( + &path, + &pixels, + config.width, + config.height, + image::ExtendedColorType::Rgb8, + ) + .unwrap(); + eprintln!("wrote {}", path.display()); +} diff --git a/crates/resolved/tests/mesh.rs b/crates/resolved/tests/mesh.rs new file mode 100644 index 0000000..f879992 --- /dev/null +++ b/crates/resolved/tests/mesh.rs @@ -0,0 +1,287 @@ +//! The mesh pass over the bodies pass: shading, silhouette, occlusion in both +//! directions and eclipse, on the GPU and on the CPU (lavapipe). + +mod common; + +use std::io::Cursor; + +use approx::assert_relative_eq; +use resolved::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS}; +use resolved::eclipse::sun_visible_fraction; +use resolved::gpu::bodies::{BodiesPass, BodiesUniforms, BodyRadiometry}; +use resolved::gpu::mesh::{MeshDraw, MeshPass}; +use resolved::gpu::targets::FrameTargets; +use resolved::gpu::{AdapterChoice, GpuContext}; +use resolved::model::{LoadOptions, Material, Mesh, Part}; +use resolved::nalgebra::{UnitQuaternion, Vector3}; +use resolved::radiometry::{Radiometry, AU}; +use resolved::{CameraConfig, Pose, Scene, Target}; + +const ALBEDO: f64 = 0.5; + +struct Frame { + radiance: Vec<[f32; 4]>, + depth: Vec, +} + +impl Frame { + fn at(&self, u: usize, v: usize) -> ([f32; 4], f32) { + (self.radiance[v * 1920 + u], self.depth[v * 1920 + u]) + } +} + +fn contexts() -> [GpuContext; 2] { + [AdapterChoice::Gpu, AdapterChoice::Cpu].map(|c| GpuContext::new(&c).unwrap()) +} + +/// A cube with edges of 2 `scale` meters, centered on its origin, gray albedo 0.5. +fn cube(scale: f64) -> Mesh { + let mut s = String::from("solid cube\n"); + let v = |i: usize| [i & 1, (i >> 1) & 1, (i >> 2) & 1].map(|b| 2.0 * b as f32 - 1.0); + for f in [ + [0, 2, 3, 1], + [4, 5, 7, 6], + [0, 1, 5, 4], + [2, 6, 7, 3], + [0, 4, 6, 2], + [1, 3, 7, 5], + ] { + for tri in [[f[0], f[1], f[2]], [f[0], f[2], f[3]]] { + s += "facet normal 0 0 0\nouter loop\n"; + for i in tri { + let [x, y, z] = v(i); + s += &format!("vertex {x} {y} {z}\n"); + } + s += "endloop\nendfacet\n"; + } + } + s += "endsolid cube\n"; + let options = LoadOptions { + scale, + default_material: Material::gray("gray", ALBEDO), + }; + Mesh::from_stl(&mut Cursor::new(s), "cube", &options).unwrap() +} + +/// Channel radiance of a gray-0.5 surface facing the Sun at `sun_distance`. +fn facing_sun(sun_distance: f64) -> [f32; 4] { + let r = Radiometry::new(&CameraConfig::example()); + let mut out = [0.0; 4]; + for (k, o) in out.iter_mut().take(3).enumerate() { + *o = r.lambert_radiance(k, ALBEDO, 1.0, sun_distance) as f32; + } + out +} + +fn scene(sun: Vector3, target: Option) -> Scene { + Scene { + earth_attitude: UnitQuaternion::identity(), + sun, + moon: Vector3::new(0.0, -3.844e8, 0.0), + targets: target.into_iter().map(|pose| Target { pose }).collect(), + } +} + +/// Renders the bodies, then `mesh` at each target in `scene`, with the eclipse factor +/// from Earth. +fn render(ctx: &GpuContext, name: &str, scene: &Scene, camera: &Pose, mesh: &Mesh) -> Frame { + let config = CameraConfig::example(); + let k = config.intrinsics(); + let radiometry = Radiometry::new(&config); + let targets = FrameTargets::new(ctx, config.width, config.height); + let bodies = BodiesPass::new(ctx); + let mut meshes = MeshPass::new(ctx); + let gpu_mesh = meshes.upload(ctx, mesh, &radiometry).unwrap(); + let body_radiometry = BodyRadiometry::new(&radiometry, &Appearance::default()).unwrap(); + let uniforms = BodiesUniforms::new(scene, camera, &k, config.near_clip, &body_radiometry); + let draws: Vec<_> = scene + .targets + .iter() + .map(|t| { + let visible = sun_visible_fraction( + t.pose.position, + scene.sun, + &[(Vector3::zeros(), EARTH_EQUATORIAL_RADIUS)], + ); + ( + &gpu_mesh, + MeshDraw::new(&t.pose, camera, &k, config.near_clip, scene.sun, visible), + ) + }) + .collect(); + let mut encoder = ctx.device.create_command_encoder(&Default::default()); + bodies.encode(ctx, &mut encoder, &targets, &uniforms); + meshes.encode(ctx, &mut encoder, &targets, &draws); + ctx.queue.submit([encoder.finish()]); + let radiance = targets.read_radiance(ctx); + common::save(name, ctx, &radiance); + Frame { + radiance, + depth: targets.read_depth(ctx), + } +} + +fn identity_at(position: Vector3) -> Pose { + Pose::new(position, UnitQuaternion::identity()) +} + +#[test] +fn sunlit_face_radiance_depth_and_silhouette() { + // A 2 m cube 7000 km from Earth's center on the sunward side, seen from 30 m + // farther out, looking down at Earth. Its +X face faces both the camera and the Sun. + let sun = Vector3::new(AU, 0.0, 0.0); + let target = Vector3::new(7.0e6, 0.0, 0.0); + let camera = + Pose::looking_at(target + Vector3::new(30.0, 0.0, 0.0), target, Vector3::z()).unwrap(); + let scene = scene(sun, Some(identity_at(target))); + let config = CameraConfig::example(); + let k = config.intrinsics(); + let expected = facing_sun((sun - target).norm()); + // Pixel centers within fx * 1 / 29 of the principal point see the face. + let half_width = k.fx / 29.0; + let expected_width = (0..1920) + .filter(|u| ((*u as f64 + 0.5) - k.cx).abs() < half_width) + .count(); + for ctx in contexts() { + let frame = render(&ctx, "cube-sunlit", &scene, &camera, &cube(1.0)); + let (radiance, depth) = frame.at(960, 540); + for c in 0..3 { + assert_relative_eq!(radiance[c], expected[c], max_relative = 1e-4); + } + assert_relative_eq!(depth as f64, config.near_clip / 29.0, max_relative = 1e-5); + // Earth fills the rest of the view, much farther away. + let width = (0..1920).filter(|&u| frame.at(u, 540).1 > 1e-5).count(); + assert!( + width.abs_diff(expected_width) <= 2, + "{}: {width} vs {expected_width}", + ctx.info.name + ); + let (earth, earth_depth) = frame.at(100, 540); + assert!(earth[1] > 0.0 && earth_depth > 0.0 && earth_depth < 1e-7); + } +} + +#[test] +fn earth_hides_a_target_behind_it_but_not_in_front() { + // A 100 km cube (about 5 px across from 200,000 km), lit from +Y so all of it that + // could show is sunlit. + let sun = Vector3::new(0.0, AU, 0.0); + let camera = Pose::looking_at( + Vector3::new(2.0e8, 0.0, 0.0), + Vector3::zeros(), + Vector3::z(), + ) + .unwrap(); + let big = cube(5.0e4); + for ctx in contexts() { + let empty = render(&ctx, "cube-scene-empty", &scene(sun, None), &camera, &big); + let behind = render( + &ctx, + "cube-behind-earth", + &scene(sun, Some(identity_at(Vector3::new(-7.0e6, 0.0, 0.0)))), + &camera, + &big, + ); + assert!( + behind.radiance == empty.radiance && behind.depth == empty.depth, + "{}", + ctx.info.name + ); + let front = render( + &ctx, + "cube-in-front-of-earth", + &scene(sun, Some(identity_at(Vector3::new(7.0e6, 0.0, 0.0)))), + &camera, + &big, + ); + let changed = (0..front.depth.len()) + .filter(|&i| front.depth[i] != empty.depth[i]) + .count(); + assert!( + (10..100).contains(&changed), + "{}: {changed} px changed", + ctx.info.name + ); + } +} + +#[test] +fn penumbra_dims_and_umbra_darkens() { + // Targets 7000 km behind Earth, one on the shadow axis and one at Earth's radius + // off it (inside the penumbra), each seen from 30 m on its sunward side. + let sun = Vector3::new(AU, 0.0, 0.0); + for ctx in contexts() { + for y in [0.0, EARTH_EQUATORIAL_RADIUS] { + let target = Vector3::new(-7.0e6, y, 0.0); + let camera = + Pose::looking_at(target + Vector3::new(30.0, 0.0, 0.0), target, Vector3::z()) + .unwrap(); + let frame = render( + &ctx, + &format!( + "cube-shadow-{}", + if y == 0.0 { "umbra" } else { "penumbra" } + ), + &scene(sun, Some(identity_at(target))), + &camera, + &cube(1.0), + ); + let visible = + sun_visible_fraction(target, sun, &[(Vector3::zeros(), EARTH_EQUATORIAL_RADIUS)]); + let full = facing_sun((sun - target).norm()); + let (radiance, depth) = frame.at(960, 540); + assert!(depth > 0.0); + for c in 0..3 { + assert_relative_eq!( + radiance[c], + full[c] * visible as f32, + max_relative = 1e-4, + epsilon = 1e-30 + ); + } + if y == 0.0 { + assert_eq!(visible, 0.0); + } else { + assert!( + visible > 0.3 && visible < 0.7, + "penumbra fraction {visible}" + ); + } + } + } +} + +#[test] +fn thin_panel_is_lit_on_the_side_facing_the_camera() { + // One triangle whose stored normal points away from the camera; the Sun is on the + // camera's side. Two-sided shading lights the side the camera sees. + let panel = Mesh { + positions: vec![[0.0, -5.0, -5.0], [0.0, 5.0, -5.0], [0.0, 0.0, 5.0]], + normals: vec![[-1.0, 0.0, 0.0]; 3], + part_of: vec![0; 3], + indices: vec![0, 1, 2], + parts: vec![Part { + name: "panel".into(), + material: 0, + }], + materials: vec![Material::gray("gray", ALBEDO)], + }; + let sun = Vector3::new(AU, 0.0, 0.0); + let target = Vector3::new(7.0e6, 0.0, 0.0); + let camera = + Pose::looking_at(target + Vector3::new(30.0, 0.0, 0.0), target, Vector3::z()).unwrap(); + let expected = facing_sun((sun - target).norm()); + for ctx in contexts() { + let frame = render( + &ctx, + "thin-panel", + &scene(sun, Some(identity_at(target))), + &camera, + &panel, + ); + let (radiance, _) = frame.at(960, 540); + for c in 0..3 { + assert_relative_eq!(radiance[c], expected[c], max_relative = 1e-4); + } + } +} diff --git a/docs/optics.md b/docs/optics.md index a9a4c18..8ef2dee 100644 --- a/docs/optics.md +++ b/docs/optics.md @@ -313,13 +313,45 @@ 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). +## How shapes become pixels + +Two techniques turn 3D shapes into an image. The renderer uses one for the +bodies and the other for the target. + +- **Mesh.** A surface approximated by many small flat triangles that share + corners (vertices). A cube needs 12 triangles. A curved surface needs + many more, and it is still faceted up close. Each vertex stores a position + and usually a normal, the direction the surface faces there. STL and glTF + files are meshes, and so is the target spacecraft. +- **Rasterization.** The GPU's built-in way of drawing a mesh. Each + triangle's corners are projected into the image with the pinhole + equations above. The GPU then finds every pixel whose center falls inside + the projected triangle, and runs a small program (a fragment shader) for + each of those pixels to compute its value. It is fast because it only + touches the pixels each triangle covers, and GPUs have dedicated hardware + for it. It can only draw triangles. +- **Ray casting.** The reverse direction: for each pixel, follow the line + of sight from the camera through the pixel center and solve for where it + first meets a surface. For a sphere or an ellipsoid that is a quadratic + equation with an exact answer, so there are no triangles and no facets. + Every pixel does the work, whether or not it ends up seeing the body. + +The bodies use ray casting, because their shapes have exact equations and +their edges are what limb-based navigation measures. The target uses +rasterization, because a spacecraft model has no equation, only its mesh. +Both write into the same depth image, so they sort against each other +correctly. Ray casting a mesh is also possible, by testing rays against +its triangles. That is what M5's ray tracer does, because following rays +gives exact shadows and light bouncing between surfaces. Rasterization can +only approximate those, for example with the shadow maps planned for M2. + ## 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 +The three bodies are not drawn as triangle meshes; they are ray cast. 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. @@ -398,6 +430,49 @@ measured reflectance spectra (M4 adds those, and BRDFs). STL has no units. CAD tools usually export millimeters, so loading a CAD model usually needs `scale: 0.001`. glTF is always meters. +## Drawing the target + +Code: `gpu/mesh.rs`, `shaders/mesh.wgsl`. + +The target's mesh is rasterized after the bodies, into the same radiance +and depth images. Each pixel the target covers is kept only if the target +is nearer than what is already there. That one rule handles a target in +front of Earth, a target behind Earth's limb, and parts of the target +hiding other parts. + +Each point on the target is shaded like the bodies: albedo times cos(i) +times the sunlight that reaches it. A panel modeled with no thickness, such +as a solar array, has one surface but two sides. Whichever side faces the +camera is the side that is lit. A closed shape never shows its inside, so +this only matters for thin parts. + +### Eclipses + +Code: `eclipse.rs`. + +The Sun is a disk, not a point. So as a spacecraft passes into Earth's +shadow, the light fades over a few seconds instead of switching off: + +- **Umbra.** Earth covers the whole solar disk. No direct sunlight. +- **Penumbra.** Earth covers part of the disk. 7000 km from Earth's + center, on the night side, the penumbra is about 65 km wide; a satellite + crosses it in under 10 seconds. +- **Annular.** A body small enough to sit inside the solar disk blocks a + ring-shaped bite. From orbit this happens with the Moon. + +The renderer finds the visible fraction of the solar disk from the overlap +of two circles on the sky. The circles are the Sun's disk and the +occulter's, each with its apparent angular radius. It multiplies the +target's sunlight by that fraction. The model is the standard one for +satellites (Montenbruck and Gill), the same as astrokit's. A test checks it +against tracing lines from the target to thousands of points spread over +the solar disk. + +The fraction is computed once per target, at its center. The whole target +dims together, which is accurate unless the target is a sizable fraction of +the penumbra's width. Light reaching the target from Earth itself +(Earthshine), which matters most when the target is in shadow, comes in M2. + ## Not modeled yet | Effect | Where | diff --git a/readme.md b/readme.md index d3d067d..4dca84f 100644 --- a/readme.md +++ b/readme.md @@ -61,6 +61,14 @@ them needs a Vulkan driver for the GPU and lavapipe installed. Set `WGPU_BACKEND`, `WGPU_ADAPTER_NAME` and the other standard wgpu environment variables to steer adapter discovery. +To see what the GPU tests render, set `RESOLVED_TEST_IMAGES=1`. Each test then +writes a PNG per device to `target/tmp/test-images/`, scaled so a white +surface facing the Sun is full white: + +```sh +RESOLVED_TEST_IMAGES=1 cargo test --test bodies --test mesh +``` + ## License MIT