From 33f986044c8f67db9854301936fed3ee84e6891a Mon Sep 17 00:00:00 2001 From: Connor Johnstone Date: Wed, 14 Oct 2026 20:34:12 -0400 Subject: [PATCH] Shadows, self and nearby craft --- crates/resolved-cli/src/main.rs | 4 +- crates/resolved/src/fidelity.rs | 6 +- crates/resolved/src/gpu/mesh.rs | 144 ++++++++- crates/resolved/src/gpu/mod.rs | 1 + crates/resolved/src/gpu/profile.rs | 44 ++- crates/resolved/src/gpu/shadow.rs | 267 ++++++++++++++++ crates/resolved/src/lib.rs | 1 + crates/resolved/src/render.rs | 80 ++++- crates/resolved/src/shaders/mesh.wgsl | 78 ++++- crates/resolved/src/shaders/shadow.wgsl | 16 + crates/resolved/src/shadow.rs | 188 ++++++++++++ crates/resolved/tests/mesh.rs | 12 +- crates/resolved/tests/shadows.rs | 384 ++++++++++++++++++++++++ docs/optics.md | 85 +++++- 14 files changed, 1273 insertions(+), 37 deletions(-) create mode 100644 crates/resolved/src/gpu/shadow.rs create mode 100644 crates/resolved/src/shaders/shadow.wgsl create mode 100644 crates/resolved/src/shadow.rs create mode 100644 crates/resolved/tests/shadows.rs diff --git a/crates/resolved-cli/src/main.rs b/crates/resolved-cli/src/main.rs index ff56f96..64d6e92 100644 --- a/crates/resolved-cli/src/main.rs +++ b/crates/resolved-cli/src/main.rs @@ -156,6 +156,7 @@ fn render(args: &RenderArgs) -> Result<(), String> { appearance: file.appearance(), read_radiance: false, profile: timings.is_some(), + shadow_map_size: None, }; let start = Instant::now(); let mut renderer = Renderer::new(camera, options).map_err(|e| e.to_string())?; @@ -217,7 +218,8 @@ fn render(args: &RenderArgs) -> Result<(), String> { v[v.len() / 2] }; println!( - "gpu bodies {:.3} ms, mesh {:.3} ms, sensor {:.3} ms, first to last {:.3} ms (medians)", + "gpu shadow {:.3} ms, bodies {:.3} ms, mesh {:.3} ms, sensor {:.3} ms, first to last {:.3} ms (medians)", + median(|t| t.shadow), median(|t| t.bodies), median(|t| t.mesh), median(|t| t.sensor), diff --git a/crates/resolved/src/fidelity.rs b/crates/resolved/src/fidelity.rs index 1d16300..b75e03d 100644 --- a/crates/resolved/src/fidelity.rs +++ b/crates/resolved/src/fidelity.rs @@ -19,8 +19,10 @@ pub enum LightingTier { #[default] L0, /// L0, plus Earth's surface from an image (`RendererOptions::earth_texture`, NASA's - /// Blue Marble) with a stand-in for the atmosphere's blue. The rest of M2 (shadows, - /// Earthshine, the Moon's albedo map) joins this tier as it lands. + /// Blue Marble) with a stand-in for the atmosphere's blue, and shadows that targets + /// cast on themselves and on each other, with the penumbra the Sun's disk gives + /// them. The rest of M2 (Earthshine, the Moon's albedo map) joins this tier as it + /// lands. L1, } diff --git a/crates/resolved/src/gpu/mesh.rs b/crates/resolved/src/gpu/mesh.rs index 0ba7742..bb07be7 100644 --- a/crates/resolved/src/gpu/mesh.rs +++ b/crates/resolved/src/gpu/mesh.rs @@ -4,7 +4,7 @@ //! per-target uniform block ([`MeshDraw`]) that places it and says how much sunlight //! reaches it. -use nalgebra::{Matrix3, Vector3}; +use nalgebra::{Matrix3, Matrix4, Vector3}; use super::bodies::MAX_CHANNELS; use super::profile::PassTimestamps; @@ -18,12 +18,25 @@ use crate::Error; #[repr(C)] #[derive(Clone, Copy, Debug, bytemuck::Pod, bytemuck::Zeroable)] -struct Vertex { +pub(crate) struct Vertex { position: [f32; 3], normal: [f32; 3], part: u32, } +impl Vertex { + /// The vertex buffer layout: position, normal, part. + pub(crate) fn layout() -> wgpu::VertexBufferLayout<'static> { + const ATTRIBUTES: [wgpu::VertexAttribute; 3] = + wgpu::vertex_attr_array![0 => Float32x3, 1 => Float32x3, 2 => Uint32]; + wgpu::VertexBufferLayout { + array_stride: std::mem::size_of::() as u64, + step_mode: wgpu::VertexStepMode::Vertex, + attributes: &ATTRIBUTES, + } + } +} + /// A model on the GPU, ready to draw. #[derive(Debug)] pub struct GpuMesh { @@ -31,6 +44,40 @@ pub struct GpuMesh { indices: wgpu::Buffer, index_count: u32, parts: wgpu::BindGroup, + /// Bounding sphere center in the body frame, in meters. + pub center: Vector3, + /// Bounding sphere radius, in meters. + pub radius: f64, +} + +impl GpuMesh { + /// Binds the vertex and index buffers and draws every triangle. + pub(crate) fn draw(&self, pass: &mut wgpu::RenderPass) { + pass.set_vertex_buffer(0, self.vertices.slice(..)); + pass.set_index_buffer(self.indices.slice(..), wgpu::IndexFormat::Uint32); + pass.draw_indexed(0..self.index_count, 0, 0..1); + } +} + +/// Where a target's fragments look up their shadow: its own shadow map. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct ShadowDraw { + /// Body frame to the map's light coordinates ([`crate::shadow::ShadowPlan`]). + pub light_from_body: Matrix4, + /// The map's array layer. + pub layer: u32, + /// Meters per shadow-map texel. + pub texel: f64, + /// The map's depth range, in meters. + pub depth_range: f64, + /// Texels across the map. + pub size: u32, +} + +/// The tangent of the Sun's angular diameter at 1 AU: how fast a shadow's edge blurs +/// with distance from what casts it. 0.93 % of the distance. +fn sun_spread(sun_distance: f64) -> f64 { + (2.0 * (crate::radiometry::SUN_RADIUS / sun_distance).asin()).tan() } /// One target's uniform block; the layout matches `Draw` in `mesh.wgsl`. @@ -43,6 +90,11 @@ pub struct MeshDraw { translation: [f32; 4], sun: [f32; 4], params: [f32; 4], + light_from_body: [[f32; 4]; 4], + // Layer, enabled (0 or 1), meters per texel, depth range in meters. + shadow: [f32; 4], + // Tangent of the Sun's angular diameter, normal offset in meters, map size in texels. + shadow_more: [f32; 4], } impl MeshDraw { @@ -55,6 +107,7 @@ impl MeshDraw { near: f64, sun: Vector3, sun_visible: f64, + shadow: Option, ) -> Self { let to_camera = camera.attitude.inverse(); let rotation: Matrix3 = (to_camera * target.attitude) @@ -80,6 +133,23 @@ impl MeshDraw { 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], + light_from_body: shadow.map_or([[0.0; 4]; 4], |s| { + let m = s.light_from_body.cast::(); + std::array::from_fn(|c| [m[(0, c)], m[(1, c)], m[(2, c)], m[(3, c)]]) + }), + shadow: shadow.map_or([0.0; 4], |s| { + [s.layer as f32, 1.0, s.texel as f32, s.depth_range as f32] + }), + // Offsetting the lookup a texel and a half along the normal keeps a surface + // from shadowing itself where it is sampled at a slant. + shadow_more: shadow.map_or([0.0; 4], |s| { + [ + sun_spread(to_sun.norm()) as f32, + (1.5 * s.texel) as f32, + s.size as f32, + 0.0, + ] + }), } } } @@ -90,6 +160,9 @@ pub struct MeshPass { pipeline: wgpu::RenderPipeline, draw_layout: wgpu::BindGroupLayout, parts_layout: wgpu::BindGroupLayout, + shadow_layout: wgpu::BindGroupLayout, + /// Bound when there are no shadow maps: a 1 x 1 x 1 array the shader does not read. + no_shadows: wgpu::BindGroup, draws: wgpu::Buffer, draws_group: wgpu::BindGroup, stride: u64, @@ -126,11 +199,40 @@ impl MeshPass { count: None, }], }); + let shadow_layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { + label: Some("mesh shadows"), + entries: &[ + wgpu::BindGroupLayoutEntry { + binding: 0, + visibility: wgpu::ShaderStages::FRAGMENT, + ty: wgpu::BindingType::Texture { + sample_type: wgpu::TextureSampleType::Depth, + view_dimension: wgpu::TextureViewDimension::D2Array, + multisampled: false, + }, + count: None, + }, + wgpu::BindGroupLayoutEntry { + binding: 1, + visibility: wgpu::ShaderStages::FRAGMENT, + ty: wgpu::BindingType::Sampler(wgpu::SamplerBindingType::Comparison), + count: None, + }, + ], + }); let layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor { label: Some("mesh"), - bind_group_layouts: &[Some(&draw_layout), Some(&parts_layout)], + bind_group_layouts: &[ + Some(&draw_layout), + Some(&parts_layout), + Some(&shadow_layout), + ], immediate_size: 0, }); + let no_shadows = { + let mut placeholder = super::shadow::ShadowPass::new(ctx, 1); + placeholder.prepare(ctx, 1, &shadow_layout) + }; let module = device.create_shader_module(wgpu::ShaderModuleDescriptor { label: Some("mesh"), source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/mesh.wgsl").into()), @@ -142,11 +244,7 @@ impl MeshPass { 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], - })], + buffers: &[Some(Vertex::layout())], }, // No culling: thin panels are seen from both sides. primitive: wgpu::PrimitiveState::default(), @@ -174,6 +272,8 @@ impl MeshPass { pipeline, draw_layout, parts_layout, + shadow_layout, + no_shadows, draws, draws_group, stride, @@ -234,6 +334,18 @@ impl MeshPass { }) .collect(); let parts: Vec<[f32; 4]> = mesh.parts.iter().map(|p| lambert[p.material]).collect(); + let (mut lo, mut hi) = (Vector3::repeat(f64::MAX), Vector3::repeat(f64::MIN)); + for p in &mesh.positions { + let p = Vector3::from(*p).cast::(); + lo = lo.inf(&p); + hi = hi.sup(&p); + } + let center = (lo + hi) / 2.0; + let radius = mesh + .positions + .iter() + .map(|p| (Vector3::from(*p).cast::() - center).norm()) + .fold(0.0, f64::max); let vertices: Vec = (0..mesh.positions.len()) .map(|i| Vertex { position: mesh.positions[i], @@ -266,6 +378,8 @@ impl MeshPass { wgpu::BufferUsages::INDEX, ), index_count: mesh.indices.len() as u32, + center, + radius, parts: ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { label: Some("mesh parts"), layout: &self.parts_layout, @@ -277,14 +391,21 @@ impl MeshPass { }) } + /// The layout of the shadow-map bind group the shader samples (group 2). + pub fn shadow_layout(&self) -> &wgpu::BindGroupLayout { + &self.shadow_layout + } + /// Uploads each target's uniforms and records the pass into `encoder`, drawing over - /// what earlier passes left in `targets`. + /// what earlier passes left in `targets`. `shadows` is the shadow maps' bind group + /// ([`super::shadow::ShadowPass::prepare`]), or `None` when no draw uses one. pub fn encode( &mut self, ctx: &GpuContext, encoder: &mut wgpu::CommandEncoder, targets: &FrameTargets, draws: &[(&GpuMesh, MeshDraw)], + shadows: Option<&wgpu::BindGroup>, timestamps: Option, ) { if draws.len() > self.capacity { @@ -324,12 +445,11 @@ impl MeshPass { ..Default::default() }); pass.set_pipeline(&self.pipeline); + pass.set_bind_group(2, shadows.unwrap_or(&self.no_shadows), &[]); 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); + mesh.draw(&mut pass); } } } diff --git a/crates/resolved/src/gpu/mod.rs b/crates/resolved/src/gpu/mod.rs index cc4a9e0..b3c477d 100644 --- a/crates/resolved/src/gpu/mod.rs +++ b/crates/resolved/src/gpu/mod.rs @@ -5,6 +5,7 @@ pub mod context; pub mod mesh; pub mod profile; pub mod sensor; +pub mod shadow; pub mod targets; pub use context::{adapters, AdapterChoice, GpuContext}; diff --git a/crates/resolved/src/gpu/profile.rs b/crates/resolved/src/gpu/profile.rs index d0ef516..33e24cf 100644 --- a/crates/resolved/src/gpu/profile.rs +++ b/crates/resolved/src/gpu/profile.rs @@ -8,31 +8,52 @@ use super::GpuContext; /// The passes timed, in the order they run. -pub const PASSES: [&str; 3] = ["bodies", "mesh", "sensor"]; +pub const PASSES: [&str; 4] = ["shadow", "bodies", "mesh", "sensor"]; /// Where one pass writes its begin and end timestamps. #[derive(Clone, Copy, Debug)] pub struct PassTimestamps<'a> { set: &'a wgpu::QuerySet, begin: u32, + write_begin: bool, + write_end: bool, } impl<'a> PassTimestamps<'a> { + /// The same pass's timestamps, keeping only the begin, the end, both or neither: + /// for one logical pass made of several GPU passes. + pub fn split(self, begin: bool, end: bool) -> Self { + Self { + write_begin: self.write_begin && begin, + write_end: self.write_end && end, + ..self + } + } + + fn indices(&self) -> (Option, Option) { + ( + self.write_begin.then_some(self.begin), + self.write_end.then_some(self.begin + 1), + ) + } + /// For a render pass descriptor. pub fn render(self) -> wgpu::RenderPassTimestampWrites<'a> { + let (begin, end) = self.indices(); wgpu::RenderPassTimestampWrites { query_set: self.set, - beginning_of_pass_write_index: Some(self.begin), - end_of_pass_write_index: Some(self.begin + 1), + beginning_of_pass_write_index: begin, + end_of_pass_write_index: end, } } /// For a compute pass descriptor. pub fn compute(self) -> wgpu::ComputePassTimestampWrites<'a> { + let (begin, end) = self.indices(); wgpu::ComputePassTimestampWrites { query_set: self.set, - beginning_of_pass_write_index: Some(self.begin), - end_of_pass_write_index: Some(self.begin + 1), + beginning_of_pass_write_index: begin, + end_of_pass_write_index: end, } } } @@ -40,6 +61,8 @@ impl<'a> PassTimestamps<'a> { /// GPU time spent in each pass of one frame, in milliseconds. #[derive(Clone, Copy, Debug, PartialEq)] pub struct PassTimings { + /// Shadow maps (zero work below lighting L1). + pub shadow: f64, /// Earth, Moon and Sun. pub bodies: f64, /// Target meshes. @@ -105,6 +128,8 @@ impl Profiler { PassTimestamps { set: &self.set, begin: 2 * index as u32, + write_begin: true, + write_end: true, } } @@ -118,10 +143,11 @@ impl Profiler { pub fn timings(&self, ticks: &[u64]) -> PassTimings { let ms = |a: u64, b: u64| b.wrapping_sub(a) as f64 * self.period * 1e-6; PassTimings { - bodies: ms(ticks[0], ticks[1]), - mesh: ms(ticks[2], ticks[3]), - sensor: ms(ticks[4], ticks[5]), - total: ms(ticks[0], ticks[5]), + shadow: ms(ticks[0], ticks[1]), + bodies: ms(ticks[2], ticks[3]), + mesh: ms(ticks[4], ticks[5]), + sensor: ms(ticks[6], ticks[7]), + total: ms(ticks[0], ticks[7]), } } } diff --git a/crates/resolved/src/gpu/shadow.rs b/crates/resolved/src/gpu/shadow.rs new file mode 100644 index 0000000..d1ec2d5 --- /dev/null +++ b/crates/resolved/src/gpu/shadow.rs @@ -0,0 +1,267 @@ +//! The pass that renders shadow maps: each target's depth as seen from the Sun. +//! +//! One layer of a `Depth32Float` array per target, cleared to 1 (farthest from the Sun) +//! and drawn with every caster the planner chose for it ([`crate::shadow::plan`]). The +//! mesh pass then looks each fragment up in its target's layer. + +use nalgebra::Matrix4; + +use super::mesh::{GpuMesh, Vertex}; +use super::profile::PassTimestamps; +use super::GpuContext; + +/// The shadow maps' depth format. +pub const SHADOW_FORMAT: wgpu::TextureFormat = wgpu::TextureFormat::Depth32Float; + +/// One caster drawn into one layer. +#[derive(Clone, Copy, Debug)] +pub struct CasterDraw<'a> { + /// Which target's map (array layer). + pub layer: u32, + /// The caster's model. + pub mesh: &'a GpuMesh, + /// Body frame to the layer's light coordinates. + pub light_from_body: Matrix4, +} + +/// The shadow pipeline, the map array and the uniforms for each caster draw. +#[derive(Debug)] +pub struct ShadowPass { + pipeline: wgpu::RenderPipeline, + caster_layout: wgpu::BindGroupLayout, + casters: wgpu::Buffer, + casters_group: wgpu::BindGroup, + caster_capacity: usize, + stride: u64, + size: u32, + maps: Option<(wgpu::Texture, u32)>, + sampler: wgpu::Sampler, +} + +impl ShadowPass { + /// Compiles the shadow pipeline for maps of `size` x `size` texels. + pub fn new(ctx: &GpuContext, size: u32) -> Self { + let device = &ctx.device; + let caster_layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { + label: Some("shadow caster"), + entries: &[wgpu::BindGroupLayoutEntry { + binding: 0, + visibility: wgpu::ShaderStages::VERTEX, + ty: wgpu::BindingType::Buffer { + ty: wgpu::BufferBindingType::Uniform, + has_dynamic_offset: true, + min_binding_size: wgpu::BufferSize::new(64), + }, + count: None, + }], + }); + let layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor { + label: Some("shadow"), + bind_group_layouts: &[Some(&caster_layout)], + immediate_size: 0, + }); + let module = device.create_shader_module(wgpu::ShaderModuleDescriptor { + label: Some("shadow"), + source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/shadow.wgsl").into()), + }); + let pipeline = device.create_render_pipeline(&wgpu::RenderPipelineDescriptor { + label: Some("shadow"), + layout: Some(&layout), + vertex: wgpu::VertexState { + module: &module, + entry_point: Some("vs_main"), + compilation_options: Default::default(), + buffers: &[Some(Vertex::layout())], + }, + // No culling: thin panels cast shadows whichever way they face. + primitive: wgpu::PrimitiveState::default(), + depth_stencil: Some(wgpu::DepthStencilState { + format: SHADOW_FORMAT, + depth_write_enabled: Some(true), + depth_compare: Some(wgpu::CompareFunction::Less), + stencil: Default::default(), + // Pushes stored depth away from the Sun on slanted surfaces, so a + // surface does not shadow itself where its own depth is sampled. + bias: wgpu::DepthBiasState { + constant: 2, + slope_scale: 2.0, + clamp: 0.0, + }, + }), + multisample: Default::default(), + fragment: None, + multiview_mask: None, + cache: None, + }); + let alignment = device.limits().min_uniform_buffer_offset_alignment as u64; + let stride = 64u64.next_multiple_of(alignment); + let (casters, casters_group) = Self::caster_slots(ctx, &caster_layout, stride, 1); + let sampler = device.create_sampler(&wgpu::SamplerDescriptor { + label: Some("shadow compare"), + mag_filter: wgpu::FilterMode::Linear, + min_filter: wgpu::FilterMode::Linear, + compare: Some(wgpu::CompareFunction::LessEqual), + ..Default::default() + }); + Self { + pipeline, + caster_layout, + casters, + casters_group, + caster_capacity: 1, + stride, + size, + maps: None, + sampler, + } + } + + fn caster_slots( + ctx: &GpuContext, + layout: &wgpu::BindGroupLayout, + stride: u64, + capacity: usize, + ) -> (wgpu::Buffer, wgpu::BindGroup) { + let buffer = ctx.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("shadow casters"), + 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("shadow casters"), + layout, + entries: &[wgpu::BindGroupEntry { + binding: 0, + resource: wgpu::BindingResource::Buffer(wgpu::BufferBinding { + buffer: &buffer, + offset: 0, + size: wgpu::BufferSize::new(64), + }), + }], + }); + (buffer, group) + } + + /// Texels across each map. + pub fn size(&self) -> u32 { + self.size + } + + /// Makes room for `layers` maps, reallocating the array if it is too small, and + /// returns the bind group the mesh pass samples them through (`layout` is its + /// shadow group layout). Allocation happens only when the target count grows. + pub fn prepare( + &mut self, + ctx: &GpuContext, + layers: u32, + layout: &wgpu::BindGroupLayout, + ) -> wgpu::BindGroup { + let have = self.maps.as_ref().map_or(0, |(_, n)| *n); + if have < layers { + let count = layers.next_power_of_two(); + let texture = ctx.device.create_texture(&wgpu::TextureDescriptor { + label: Some("shadow maps"), + size: wgpu::Extent3d { + width: self.size, + height: self.size, + depth_or_array_layers: count, + }, + mip_level_count: 1, + sample_count: 1, + dimension: wgpu::TextureDimension::D2, + format: SHADOW_FORMAT, + usage: wgpu::TextureUsages::RENDER_ATTACHMENT + | wgpu::TextureUsages::TEXTURE_BINDING, + view_formats: &[], + }); + self.maps = Some((texture, count)); + } + let (texture, _) = self.maps.as_ref().expect("allocated above"); + let view = texture.create_view(&wgpu::TextureViewDescriptor { + dimension: Some(wgpu::TextureViewDimension::D2Array), + ..Default::default() + }); + ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("shadow maps"), + layout, + entries: &[ + wgpu::BindGroupEntry { + binding: 0, + resource: wgpu::BindingResource::TextureView(&view), + }, + wgpu::BindGroupEntry { + binding: 1, + resource: wgpu::BindingResource::Sampler(&self.sampler), + }, + ], + }) + } + + /// Records one render pass per map layer in `0..layers`, drawing `draws` into their + /// layers. `timestamps` brackets all of them; with no layers, an empty pass keeps the + /// timestamps written. + pub fn encode( + &mut self, + ctx: &GpuContext, + encoder: &mut wgpu::CommandEncoder, + layers: u32, + draws: &[CasterDraw], + timestamps: Option, + ) { + if draws.len() > self.caster_capacity { + self.caster_capacity = draws.len().next_power_of_two(); + (self.casters, self.casters_group) = + Self::caster_slots(ctx, &self.caster_layout, self.stride, self.caster_capacity); + } + for (i, draw) in draws.iter().enumerate() { + let m = draw.light_from_body; + let columns: [[f32; 4]; 4] = + std::array::from_fn(|c| [m[(0, c)], m[(1, c)], m[(2, c)], m[(3, c)]]); + ctx.queue.write_buffer( + &self.casters, + i as u64 * self.stride, + bytemuck::cast_slice(&columns), + ); + } + let texture = match &self.maps { + Some((texture, _)) if layers > 0 => texture, + _ => { + // Nothing to draw: an empty pass still writes the timestamps. + drop(encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("shadow (none)"), + timestamp_writes: timestamps.map(PassTimestamps::compute), + })); + return; + } + }; + for layer in 0..layers { + let view = texture.create_view(&wgpu::TextureViewDescriptor { + dimension: Some(wgpu::TextureViewDimension::D2), + base_array_layer: layer, + array_layer_count: Some(1), + ..Default::default() + }); + let stamps = timestamps.map(|t| t.split(layer == 0, layer + 1 == layers)); + let mut pass = encoder.begin_render_pass(&wgpu::RenderPassDescriptor { + label: Some("shadow"), + color_attachments: &[], + depth_stencil_attachment: Some(wgpu::RenderPassDepthStencilAttachment { + view: &view, + depth_ops: Some(wgpu::Operations { + load: wgpu::LoadOp::Clear(1.0), + store: wgpu::StoreOp::Store, + }), + stencil_ops: None, + }), + timestamp_writes: stamps.map(PassTimestamps::render), + ..Default::default() + }); + pass.set_pipeline(&self.pipeline); + for (i, draw) in draws.iter().enumerate().filter(|(_, d)| d.layer == layer) { + pass.set_bind_group(0, &self.casters_group, &[(i as u64 * self.stride) as u32]); + draw.mesh.draw(&mut pass); + } + } + } +} diff --git a/crates/resolved/src/lib.rs b/crates/resolved/src/lib.rs index dcd405b..23c2a7d 100644 --- a/crates/resolved/src/lib.rs +++ b/crates/resolved/src/lib.rs @@ -30,6 +30,7 @@ pub mod model; pub mod radiometry; pub mod render; pub mod scene; +pub mod shadow; mod solar_data; pub mod spectrum; diff --git a/crates/resolved/src/render.rs b/crates/resolved/src/render.rs index 8a0a3ee..2eae6a2 100644 --- a/crates/resolved/src/render.rs +++ b/crates/resolved/src/render.rs @@ -18,14 +18,16 @@ use crate::earth::EarthTexture; use crate::eclipse::sun_visible_fraction; use crate::fidelity::{Fidelity, LightingTier}; use crate::gpu::bodies::{BodiesPass, BodiesUniforms, BodyRadiometry}; -use crate::gpu::mesh::{GpuMesh, MeshDraw, MeshPass}; +use crate::gpu::mesh::{GpuMesh, MeshDraw, MeshPass, ShadowDraw}; use crate::gpu::profile::{PassTimings, Profiler}; use crate::gpu::sensor::{SensorPass, SensorUniforms}; +use crate::gpu::shadow::{CasterDraw, ShadowPass}; use crate::gpu::targets::FrameTargets; use crate::gpu::{AdapterChoice, GpuContext}; use crate::model::{Mesh, ModelId}; use crate::radiometry::Radiometry; use crate::scene::{Pose, Scene}; +use crate::shadow::{self, Bounded}; use crate::Error; /// Choices made once, when the renderer is created. @@ -45,6 +47,10 @@ pub struct RendererOptions { /// Time each pass on the GPU and report it in [`Frame::timings`]. Needs an adapter /// with timestamp queries; without one, timings are `None`. pub profile: bool, + /// Texels across each target's shadow map at lighting L1 and up; `None` for 4096. + /// Each map is 4 bytes per texel (64 MB at 4096) and covers its target's bounding + /// sphere, so a 16 m model gets 4 mm texels. + pub shadow_map_size: Option, } /// One rendered image. @@ -88,6 +94,7 @@ pub struct Renderer { targets: FrameTargets, bodies: BodiesPass, meshes: MeshPass, + shadows: ShadowPass, sensor: SensorPass, sensor_uniforms: SensorUniforms, models: Vec, @@ -140,6 +147,7 @@ impl Renderer { let mut renderer = Self { bodies: BodiesPass::new(&ctx), meshes: MeshPass::new(&ctx), + shadows: ShadowPass::new(&ctx, options.shadow_map_size.unwrap_or(4096)), sensor: SensorPass::new(&ctx), sensor_uniforms: SensorUniforms::new(&camera, &radiometry), ctx, @@ -226,12 +234,58 @@ impl Renderer { ]; let bodies = BodiesUniforms::new(scene, &camera, &intrinsics, near, &self.body_radiometry); - let draws: Vec<(&GpuMesh, MeshDraw)> = scene + // Shadow maps, from lighting L1: one per target, with the targets that can + // shadow it drawn in. + let bounded: Vec = scene .targets .iter() .map(|t| { + let mesh = &self.models[t.model.0]; + Bounded { + pose: t.pose, + center: mesh.center, + radius: mesh.radius.max(1e-3), + } + }) + .collect(); + let plans = if self.fidelity.lighting >= LightingTier::L1 { + shadow::plan(&bounded, scene.sun) + } else { + Vec::new() + }; + let size = self.shadows.size(); + let casters: Vec = plans + .iter() + .enumerate() + .flat_map(|(layer, plan)| plan.casters.iter().map(move |&j| (layer, plan, j))) + .map(|(layer, plan, j)| CasterDraw { + layer: layer as u32, + mesh: &self.models[scene.targets[j].model.0], + light_from_body: plan.light_from_body(&bounded[j]).cast(), + }) + .collect(); + let draws: Vec<(&GpuMesh, MeshDraw)> = scene + .targets + .iter() + .enumerate() + .map(|(i, t)| { let visible = sun_visible_fraction(t.pose.position, scene.sun, &occulters); - let draw = MeshDraw::new(&t.pose, &camera, &intrinsics, near, scene.sun, visible); + let shadow = plans.get(i).map(|plan| ShadowDraw { + light_from_body: plan.light_from_body(&bounded[i]), + layer: i as u32, + texel: 2.0 * plan.half_width / size as f64, + depth_range: plan.depth_range(), + size, + }); + let draw = MeshDraw::new( + &t.pose, + &camera, + &intrinsics, + near, + scene.sun, + visible, + shadow, + ); (&self.models[t.model.0], draw) }) .collect(); @@ -239,16 +293,28 @@ impl Renderer { let ctx = &self.ctx; let mut encoder = ctx.device.create_command_encoder(&Default::default()); let stamps = |i| self.profiler.as_ref().map(|p| p.pass(i)); + let shadow_group = (!plans.is_empty()).then(|| { + self.shadows + .prepare(ctx, plans.len() as u32, self.meshes.shadow_layout()) + }); + self.shadows + .encode(ctx, &mut encoder, plans.len() as u32, &casters, stamps(0)); self.bodies - .encode(ctx, &mut encoder, &self.targets, &bodies, stamps(0)); - self.meshes - .encode(ctx, &mut encoder, &self.targets, &draws, stamps(1)); + .encode(ctx, &mut encoder, &self.targets, &bodies, stamps(1)); + self.meshes.encode( + ctx, + &mut encoder, + &self.targets, + &draws, + shadow_group.as_ref(), + stamps(2), + ); self.sensor.encode( ctx, &mut encoder, &self.targets, &self.sensor_uniforms, - stamps(2), + stamps(3), ); if let Some(profiler) = &self.profiler { profiler.resolve(&mut encoder); diff --git a/crates/resolved/src/shaders/mesh.wgsl b/crates/resolved/src/shaders/mesh.wgsl index 47f2e2e..aa0ea83 100644 --- a/crates/resolved/src/shaders/mesh.wgsl +++ b/crates/resolved/src/shaders/mesh.wgsl @@ -14,11 +14,19 @@ struct Draw { sun: vec4, // x: Sun visibility (eclipse) / (Sun distance in AU)^2. params: vec4, + // Body frame to this target's shadow map (shadow.rs). + light_from_body: mat4x4, + // Layer, enabled (0 or 1), meters per texel, depth range in meters. + shadow: vec4, + // Tangent of the Sun's angular diameter, normal offset in meters, map size in texels. + shadow_more: 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>; +@group(2) @binding(0) var shadow_maps: texture_depth_2d_array; +@group(2) @binding(1) var shadow_compare: sampler_comparison; struct VertexIn { @location(0) position: vec3, @@ -31,19 +39,80 @@ struct VertexOut { @location(0) camera_position: vec3, @location(1) normal: vec3, @location(2) @interpolate(flat) part: u32, + // Shadow-map coordinates: x and y in [-1, 1], depth in [0, 1]. + @location(3) light: vec3, } @vertex fn vs_main(v: VertexIn) -> VertexOut { let p = draw.camera_from_body * v.position + draw.translation.xyz; + let n = draw.camera_from_body * v.normal; 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.normal = n; out.part = v.part; + // Look the shadow up slightly off the surface, on its sunward side, so a surface + // does not shadow itself where the map samples it at a slant. + let toward_sun = select(-1.0, 1.0, dot(n, draw.sun.xyz) >= 0.0); + let offset = v.normal * toward_sun * draw.shadow_more.y; + out.light = (draw.light_from_body * vec4(v.position + offset, 1.0)).xyz; return out; } +// The fraction of the Sun this point sees past the target's own geometry and its +// neighbors, from the shadow map, with a penumbra as wide as the Sun's disk makes it. +// +// Percentage-closer soft shadows: first the mean depth of whatever blocks the Sun near +// this point, then a comparison filter whose width is the penumbra that blocker casts, +// (receiver - blocker distance) x the Sun's angular diameter. Both use a regular grid +// of samples rather than a random pattern, so the result has no noise. +fn sunlit(light: vec3) -> f32 { + if draw.shadow.y < 0.5 { + return 1.0; + } + let layer = i32(draw.shadow.x); + let texel = draw.shadow.z; + let range = draw.shadow.w; + let spread = draw.shadow_more.x; + let size = draw.shadow_more.z; + let uv = vec2(light.x * 0.5 + 0.5, 0.5 - light.y * 0.5); + // Two texels of depth margin against comparing a surface with itself. + let depth = light.z - 2.0 * texel / range; + + // Blockers can only be between this point and the map's sunward edge; anything + // there casting onto this point lies within that distance times the spread. + let search = clamp(light.z * range * spread / texel, 1.0, 32.0); + var blockers = 0.0; + var blocker_depth = 0.0; + for (var i = 0; i < 5; i++) { + for (var j = 0; j < 5; j++) { + let offset = (vec2(f32(i), f32(j)) * 0.5 - 1.0) * search; + let texel_xy = clamp(vec2(floor(uv * size + offset)), vec2(0), vec2(i32(size) - 1)); + let d = textureLoad(shadow_maps, texel_xy, layer, 0); + if d < depth { + blockers += 1.0; + blocker_depth += d; + } + } + } + if blockers == 0.0 { + return 1.0; + } + let penumbra = (depth - blocker_depth / blockers) * range * spread; + // Half the penumbra either side of the edge, and never under half a texel, which + // also smooths the map's own stair-steps. + let radius = max(0.5 * penumbra / texel, 0.5) / size; + var lit = 0.0; + for (var i = 0; i < 6; i++) { + for (var j = 0; j < 6; j++) { + let offset = ((vec2(f32(i), f32(j)) + 0.5) / 3.0 - 1.0) * radius; + lit += textureSampleCompareLevel(shadow_maps, shadow_compare, uv + offset, layer, depth); + } + } + return lit / 36.0; +} + @fragment fn fs_main(in: VertexOut) -> @location(0) vec4 { var n = normalize(in.normal); @@ -52,5 +121,10 @@ fn fs_main(in: VertexOut) -> @location(0) vec4 { 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; + let lit = max(dot(n, draw.sun.xyz), 0.0); + var shade = 0.0; + if lit > 0.0 { + shade = sunlit(in.light); + } + return parts[in.part] * lit * shade * draw.params.x; } diff --git a/crates/resolved/src/shaders/shadow.wgsl b/crates/resolved/src/shaders/shadow.wgsl new file mode 100644 index 0000000..ba4f9ac --- /dev/null +++ b/crates/resolved/src/shaders/shadow.wgsl @@ -0,0 +1,16 @@ +// Depth from the Sun: draws one caster into one target's shadow map. +// +// The matrix takes body-frame vertices straight to the map's light coordinates +// (shadow.rs): x and y across the map in [-1, 1], depth along the Sun's rays in [0, 1], +// 0 nearest the Sun. No fragment stage: only depth is written. + +struct Caster { + light_from_body: mat4x4, +} + +@group(0) @binding(0) var caster: Caster; + +@vertex +fn vs_main(@location(0) position: vec3) -> @builtin(position) vec4 { + return caster.light_from_body * vec4(position, 1.0); +} diff --git a/crates/resolved/src/shadow.rs b/crates/resolved/src/shadow.rs new file mode 100644 index 0000000..665cff9 --- /dev/null +++ b/crates/resolved/src/shadow.rs @@ -0,0 +1,188 @@ +//! Planning shadow maps: which targets can shadow which, and how each map is framed. +//! +//! Each target gets an orthographic shadow map looking along the direction from the +//! target to the Sun (the Sun is far enough that its rays are parallel across a +//! spacecraft). The map covers the target's bounding sphere, so every point of the +//! target lands in it, and holds the depth, along the Sun's direction, of whatever is +//! nearest the Sun: the target itself or another target sunward of it. +//! +//! Light coordinates of a point P for the map of a target with bounding center c and +//! radius r, along axes u, v (across the Sun's rays) and w (toward the Sun): +//! +//! ```text +//! x = (P - c) . u / r in [-1, 1] across the map +//! y = (P - c) . v / r +//! depth = (near - (P - c) . w) / (near - far) 0 nearest the Sun, 1 farthest +//! ``` +//! +//! with `near` the sunward extent of the farthest caster and `far = -r`. + +use nalgebra::{Matrix4, Vector3}; + +use crate::scene::Pose; + +/// A target as the planner sees it: where it is and the sphere around its model. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Bounded { + /// The target's pose. + pub pose: Pose, + /// Bounding sphere center in the model's body frame, in meters. + pub center: Vector3, + /// Bounding sphere radius, in meters. + pub radius: f64, +} + +impl Bounded { + /// The bounding sphere's center in the inertial frame. + pub fn world_center(&self) -> Vector3 { + self.pose.position + self.pose.attitude * self.center + } +} + +/// How one target's shadow map is framed and which targets are drawn into it. +#[derive(Clone, Debug, PartialEq)] +pub struct ShadowPlan { + /// The receiver's bounding center, inertial. + pub center: Vector3, + /// Half the map's width, in meters: the receiver's bounding radius. + pub half_width: f64, + /// Unit vectors across the map (u, v) and toward the Sun (w), inertial. + pub axes: [Vector3; 3], + /// Distance from `center` toward the Sun of the map's near plane, in meters. + pub near: f64, + /// The same for the far plane (negative: away from the Sun). + pub far: f64, + /// Indices of the targets drawn into this map, the receiver included. + pub casters: Vec, +} + +impl ShadowPlan { + /// The map's depth range in meters. + pub fn depth_range(&self) -> f64 { + self.near - self.far + } + + /// The matrix taking a point in `target`'s body frame to this map's light + /// coordinates (x, y in [-1, 1], depth in [0, 1], w = 1). + pub fn light_from_body(&self, target: &Bounded) -> Matrix4 { + let [u, v, w] = self.axes; + let range = self.depth_range(); + // Rows act on (P - c); P = R p + t for body point p. + let rows = [u / self.half_width, v / self.half_width, -w / range]; + let rotation = target.pose.attitude.to_rotation_matrix(); + let offset = target.pose.position - self.center; + let mut m = Matrix4::identity(); + for (i, row) in rows.iter().enumerate() { + let r = rotation.matrix().transpose() * row; + m[(i, 0)] = r.x; + m[(i, 1)] = r.y; + m[(i, 2)] = r.z; + m[(i, 3)] = row.dot(&offset); + } + // depth = (near - along) / range: the -along / range is in the row; add near / range. + m[(2, 3)] += self.near / range; + m + } +} + +/// One shadow map per target, lit from the Sun at `sun` (inertial, meters). +/// +/// A target shadows another when its bounding sphere reaches the cylinder swept by the +/// other's sphere toward the Sun. +pub fn plan(targets: &[Bounded], sun: Vector3) -> Vec { + targets + .iter() + .enumerate() + .map(|(i, receiver)| { + let center = receiver.world_center(); + let w = (sun - center).normalize(); + let u = w + .cross(&Vector3::z()) + .try_normalize(1e-6) + .unwrap_or_else(|| w.cross(&Vector3::x()).normalize()); + let v = w.cross(&u); + let r = receiver.radius; + let mut near = r; + let mut casters = vec![i]; + for (j, other) in targets.iter().enumerate().filter(|(j, _)| *j != i) { + let d = other.world_center() - center; + let along = d.dot(&w); + let across = (d - along * w).norm(); + let reach = r + other.radius; + if across < reach && along > -reach { + casters.push(j); + near = near.max(along + other.radius); + } + } + ShadowPlan { + center, + half_width: r, + axes: [u, v, w], + near, + far: -r, + casters, + } + }) + .collect() +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + use nalgebra::{UnitQuaternion, Vector4}; + + fn at(x: f64, y: f64, z: f64, radius: f64) -> Bounded { + Bounded { + pose: Pose::new(Vector3::new(x, y, z), UnitQuaternion::identity()), + center: Vector3::zeros(), + radius, + } + } + + const SUN: Vector3 = Vector3::new(1.496e11, 0.0, 0.0); + + #[test] + fn sunward_neighbor_casts_and_side_neighbor_does_not() { + // Target 0 at the origin; 1 is 20 m sunward (+X); 2 is 20 m off to the side. + let targets = [ + at(0.0, 0.0, 0.0, 5.0), + at(20.0, 2.0, 0.0, 3.0), + at(0.0, 20.0, 0.0, 3.0), + ]; + let plans = plan(&targets, SUN); + assert_eq!(plans[0].casters, [0, 1]); + assert_relative_eq!(plans[0].near, 23.0, epsilon = 1e-6); + assert_relative_eq!(plans[0].far, -5.0); + // Nothing is sunward of target 1. + assert_eq!(plans[1].casters, [1]); + // Target 2 is beside 0, not sunward of it: it does not shadow 0, nor 0 it. + assert_eq!(plans[2].casters, [2]); + } + + #[test] + fn light_coordinates_span_the_bounding_sphere() { + let target = Bounded { + pose: Pose::new( + Vector3::new(7.0e6, 1.0, 2.0), + UnitQuaternion::from_euler_angles(0.3, -0.5, 1.1), + ), + center: Vector3::new(0.5, -1.0, 2.0), + radius: 4.0, + }; + let plan = &plan(&[target], SUN)[0]; + let m = plan.light_from_body(&target); + let light = |body: Vector3| m * Vector4::new(body.x, body.y, body.z, 1.0); + // The bounding center maps to the middle of the map, at the depth of `center`. + let c = light(target.center); + assert_relative_eq!(c.xy(), Vector3::zeros().xy(), epsilon = 1e-9); + assert_relative_eq!(c.z, plan.near / plan.depth_range(), epsilon = 1e-9); + // A body point 4 m toward the Sun is at depth 0 (the near plane); 4 m away, 1. + let toward = target.pose.attitude.inverse() * plan.axes[2]; + assert_relative_eq!(light(target.center + 4.0 * toward).z, 0.0, epsilon = 1e-9); + assert_relative_eq!(light(target.center - 4.0 * toward).z, 1.0, epsilon = 1e-9); + // 4 m across along u is the map's edge. + let across = target.pose.attitude.inverse() * plan.axes[0]; + assert_relative_eq!(light(target.center + 4.0 * across).x, 1.0, epsilon = 1e-9); + } +} diff --git a/crates/resolved/tests/mesh.rs b/crates/resolved/tests/mesh.rs index afd92e9..b29a8ba 100644 --- a/crates/resolved/tests/mesh.rs +++ b/crates/resolved/tests/mesh.rs @@ -112,13 +112,21 @@ fn render(ctx: &GpuContext, name: &str, scene: &Scene, camera: &Pose, mesh: &Mes ); ( &gpu_mesh, - MeshDraw::new(&t.pose, camera, &k, config.near_clip, scene.sun, visible), + MeshDraw::new( + &t.pose, + camera, + &k, + config.near_clip, + scene.sun, + visible, + None, + ), ) }) .collect(); let mut encoder = ctx.device.create_command_encoder(&Default::default()); bodies.encode(ctx, &mut encoder, &targets, &uniforms, None); - meshes.encode(ctx, &mut encoder, &targets, &draws, None); + meshes.encode(ctx, &mut encoder, &targets, &draws, None, None); ctx.queue.submit([encoder.finish()]); let radiance = targets.read_radiance(ctx); common::save(name, &ctx.info, &radiance); diff --git a/crates/resolved/tests/shadows.rs b/crates/resolved/tests/shadows.rs new file mode 100644 index 0000000..881d354 --- /dev/null +++ b/crates/resolved/tests/shadows.rs @@ -0,0 +1,384 @@ +//! Shadows on targets at lighting L1, on the GPU and on the CPU (lavapipe): where they +//! fall, how wide their penumbra is, and that surfaces do not shadow themselves. + +mod common; + +use std::io::Cursor; +use std::sync::Arc; + +use resolved::gpu::AdapterChoice; +use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3}; +use resolved::radiometry::{Radiometry, AU}; +use resolved::{ + CameraConfig, EarthTexture, Fidelity, LightingTier, LoadOptions, Material, Mesh, Pose, + Renderer, RendererOptions, Scene, Target, +}; + +const ALBEDO: f64 = 0.5; +const TARGET: Vector3 = Vector3::new(7.0e6, 0.0, 0.0); +/// Straight up from the target: body +X faces the Sun with the identity attitude. +const SUN: Vector3 = Vector3::new(AU, 0.0, 0.0); +/// Half the width of the plate the shadows fall on, in meters. +const PLATE: f32 = 4.0; +/// The Sun's angular diameter seen from 1 AU, as a tangent: 2 x 696,000 km / 1 AU. +const SUN_SPREAD: f64 = 0.0093; + +/// Four-cornered faces, each split into two triangles, as one gray mesh. +fn quads(name: &str, faces: &[[[f32; 3]; 4]]) -> Mesh { + let mut s = format!("solid {name}\n"); + for f in faces { + 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 [x, y, z] in tri { + s += &format!("vertex {x} {y} {z}\n"); + } + s += "endloop\nendfacet\n"; + } + } + s += &format!("endsolid {name}\n"); + let options = LoadOptions { + default_material: Material::gray("gray", ALBEDO), + ..Default::default() + }; + Mesh::from_stl(&mut Cursor::new(s), name, &options).unwrap() +} + +/// A square in the body YZ plane at height `x`, facing +X. +fn square(x: f32, half: f32) -> [[f32; 3]; 4] { + [ + [x, -half, -half], + [x, half, -half], + [x, half, half], + [x, -half, half], + ] +} + +/// The plate alone. +fn plate() -> Mesh { + quads("plate", &[square(0.0, PLATE)]) +} + +/// The plate with a 2 m square occluder `height` meters above its center. +fn plate_under(height: f32) -> Mesh { + quads("plate", &[square(0.0, PLATE), square(height, 1.0)]) +} + +/// A UV sphere of radius 3 m, flat-faceted: 48 segments around, 24 from pole to pole. +fn sphere() -> Mesh { + let (around, down) = (48, 24); + let point = |i: usize, j: usize| { + let lon = i as f32 / around as f32 * std::f32::consts::TAU; + let lat = (j as f32 / down as f32 - 0.5) * std::f32::consts::PI; + [lat.cos() * lon.cos(), lat.cos() * lon.sin(), lat.sin()].map(|c| 3.0 * c) + }; + let faces: Vec<_> = (0..around) + .flat_map(|i| (0..down).map(move |j| (i, j))) + .map(|(i, j)| { + [ + point(i, j), + point(i + 1, j), + point(i + 1, j + 1), + point(i, j + 1), + ] + }) + .collect(); + quads("sphere", &faces) +} + +fn renderer(adapter: &AdapterChoice, lighting: LightingTier) -> Renderer { + let gray = EarthTexture::from_rgba8(1, 1, vec![128, 128, 128, 255]).unwrap(); + let options = RendererOptions { + adapter: adapter.clone(), + fidelity: Fidelity { + lighting, + ..Default::default() + }, + earth_texture: Some(Arc::new(gray)), + read_radiance: true, + ..Default::default() + }; + Renderer::new(CameraConfig::example(), options).unwrap() +} + +/// 30 m from the target along `direction` (target body frame), looking at it. +fn observer(direction: Vector3, up: Vector3) -> Pose { + Pose::looking_at(TARGET + 30.0 * direction.normalize(), TARGET, up).unwrap() +} + +/// Looking down at the plate 50 degrees from overhead, so an occluder above the plate +/// appears well clear of the shadow it casts. +fn slanted() -> Pose { + let a = 50f64.to_radians(); + observer(Vector3::new(a.cos(), 0.0, a.sin()), Vector3::z()) +} + +/// 10 m from the occluder's shadow edge (body z = 1 m), 40 degrees from overhead, +/// where a pixel spans 1.4 mm of plate across the edge, a seventh of a 1 m occluder's +/// penumbra. +fn close_to_edge() -> Pose { + let a = 40f64.to_radians(); + let edge = TARGET + Vector3::z(); + let eye = edge + 10.0 * Vector3::new(a.cos(), 0.0, a.sin()); + Pose::looking_at(eye, edge, Vector3::z()).unwrap() +} + +fn scene(targets: Vec) -> Scene { + Scene { + earth_attitude: UnitQuaternion::identity(), + sun: SUN, + moon: Vector3::new(0.0, -3.844e8, 0.0), + targets, + } +} + +fn at(model: resolved::ModelId, offset: Vector3) -> Target { + Target { + model, + pose: Pose::new(TARGET + offset, UnitQuaternion::identity()), + } +} + +/// Green-channel radiance of a Lambertian surface at `cos_i` to the Sun, at 1 AU. +fn lit(cos_i: f64) -> f64 { + Radiometry::new(&CameraConfig::example()).lambert_radiance( + 1, + ALBEDO, + cos_i, + (SUN - TARGET).norm(), + ) +} + +/// The plate's (y, z), body frame, and the green radiance there, for each pixel of +/// image column 960, the column through the plate's center line. The plate sits at the +/// identity attitude on `TARGET`. +fn column(observer: &Pose, radiance: &[[f32; 4]]) -> Vec<(f64, f64)> { + let intrinsics = CameraConfig::example().intrinsics(); + let mut samples: Vec<(f64, f64)> = (0..1080) + .filter_map(|row| { + let ray = observer.attitude * intrinsics.ray(&Vector2::new(960.5, row as f64 + 0.5)); + // Where the ray meets the plate's plane, body x = 0. + let t = (TARGET.x - observer.position.x) / ray.x; + let p = observer.position + t * ray - TARGET; + (p.y.abs() < PLATE as f64 && p.z.abs() < PLATE as f64) + .then(|| (p.z, radiance[row * 1920 + 960][1] as f64)) + }) + .collect(); + samples.sort_by(|a, b| a.0.total_cmp(&b.0)); + samples +} + +/// Where, going up through `samples` from `from`, the radiance first reaches `level`, +/// interpolated between pixels. +fn crossing(samples: &[(f64, f64)], from: f64, level: f64) -> f64 { + let i = samples + .iter() + .position(|&(z, r)| z > from && r >= level) + .expect("never reaches the level"); + let ((z0, r0), (z1, r1)) = (samples[i - 1], samples[i]); + z0 + (level - r0) / (r1 - r0) * (z1 - z0) +} + +#[test] +fn occluder_shadows_the_plate_below_it() { + let _devices = common::device_lock(); + let full = lit(1.0); + for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] { + let mut renderer = renderer(&adapter, LightingTier::L1); + let name = renderer.adapter_info().name.clone(); + let model = renderer.load_model(&plate_under(2.0)).unwrap(); + let frame = renderer + .submit(&scene(vec![at(model, Vector3::zeros())]), &slanted()) + .unwrap() + .wait() + .unwrap(); + let radiance = frame.radiance.unwrap(); + common::save("shadow-occluder", renderer.adapter_info(), &radiance); + let samples = column(&slanted(), &radiance); + for &(z, r) in &samples { + // Under the occluder, 2 cm clear of the penumbra: no direct sunlight at all. + if z.abs() < 0.97 { + assert_eq!(r, 0.0, "{name}: z = {z}"); + } + // Beyond it, and short of where the occluder hides the plate from the + // camera: plain Lambert. + if z > 1.03 { + assert!( + (r - full).abs() < 1e-4 * full, + "{name}: z = {z}: {r} vs {full}" + ); + } + } + // The edge is the occluder's edge, z = 1 m, to within a texel of the 4096 map + // across the plate's 11.3 m bounding sphere (2.8 mm). + let edge = crossing(&samples, 0.0, 0.5 * full); + assert!((edge - 1.0).abs() < 2.8e-3, "{name}: edge at {edge}"); + } +} + +#[test] +fn penumbra_widens_with_occluder_height() { + let _devices = common::device_lock(); + let full = lit(1.0); + for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] { + let mut renderer = renderer(&adapter, LightingTier::L1); + let name = renderer.adapter_info().name.clone(); + for height in [1.0, 2.0, 3.0] { + let model = renderer.load_model(&plate_under(height)).unwrap(); + let frame = renderer + .submit(&scene(vec![at(model, Vector3::zeros())]), &close_to_edge()) + .unwrap() + .wait() + .unwrap(); + let samples = column(&close_to_edge(), &frame.radiance.unwrap()); + // The Sun's disk, seen past a straight edge, uncovers linearly in the + // filter this models, so 10 % to 90 % is 0.8 of the full width. + let width = + (crossing(&samples, 0.0, 0.9 * full) - crossing(&samples, 0.0, 0.1 * full)) / 0.8; + let want = height as f64 * SUN_SPREAD; + assert!( + (width - want).abs() < 0.05 * want, + "{name}: occluder {height} m up: penumbra {width} m, want {want} m" + ); + } + } +} + +/// Pixels lit at L0 but darker at L1, among those at least 10 % of `full`: a surface +/// that nothing shadows must not darken (shadow acne). +fn darkened(l0: &[[f32; 4]], l1: &[[f32; 4]], full: f64) -> usize { + l0.iter() + .zip(l1) + .filter(|(a, b)| a[1] as f64 > 0.1 * full && b[1] < 0.999 * a[1]) + .count() +} + +#[test] +fn unshadowed_surfaces_do_not_shadow_themselves() { + let _devices = common::device_lock(); + for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] { + let mut renderers = + [LightingTier::L0, LightingTier::L1].map(|lighting| renderer(&adapter, lighting)); + let name = renderers[1].adapter_info().name.clone(); + // The plate tipped 80 and 87 degrees from the Sun, viewed face on, and a sphere + // of flat facets, whose terminator crosses facets at every angle to the Sun, + // seen against space so only the target can differ. + let a = |deg: f64| deg.to_radians(); + let tipped = |deg: f64| { + let attitude = UnitQuaternion::from_axis_angle(&Vector3::y_axis(), -a(deg)); + let normal = attitude * Vector3::x(); + ( + plate(), + attitude, + observer(normal, Vector3::y()), + format!("plate at {deg}"), + ) + }; + let cases = [ + tipped(80.0), + tipped(87.0), + ( + sphere(), + UnitQuaternion::identity(), + observer(Vector3::new(-0.3, 1.0, 0.3), Vector3::z()), + "sphere".to_string(), + ), + ]; + for (mesh, attitude, observer, case) in cases { + let frames = renderers.each_mut().map(|renderer| { + let model = renderer.load_model(&mesh).unwrap(); + let target = Target { + model, + pose: Pose::new(TARGET, attitude), + }; + renderer + .submit(&scene(vec![target]), &observer) + .unwrap() + .wait() + .unwrap() + .radiance + .unwrap() + }); + let file = format!("shadow-acne-{}", case.replace(' ', "-")); + common::save(&file, renderers[1].adapter_info(), &frames[1]); + let count = darkened(&frames[0], &frames[1], lit(1.0)); + assert_eq!(count, 0, "{name}: {case}: {count} pixels darkened"); + } + } +} + +#[test] +fn sunward_target_shadows_and_side_target_does_not() { + let _devices = common::device_lock(); + let full = lit(1.0); + for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] { + let mut renderer = renderer(&adapter, LightingTier::L1); + let name = renderer.adapter_info().name.clone(); + let plate = renderer.load_model(&plate()).unwrap(); + let panel = renderer + .load_model(&quads("panel", &[square(0.0, 1.0)])) + .unwrap(); + // A 2 m panel tipped 45 degrees, so its +Z half is the nearer the Sun; that + // half's shadow covers plate z from 0 to 0.71 m. + let tipped = |offset: Vector3| Target { + model: panel, + pose: Pose::new( + TARGET + offset, + UnitQuaternion::from_axis_angle(&Vector3::y_axis(), 45f64.to_radians()), + ), + }; + let mut at_z = |targets, want: f64| { + let frame = renderer + .submit(&scene(targets), &slanted()) + .unwrap() + .wait() + .unwrap(); + let samples = column(&slanted(), &frame.radiance.unwrap()); + let (_, r) = *samples + .iter() + .min_by(|a, b| (a.0 - want).abs().total_cmp(&(b.0 - want).abs())) + .unwrap(); + r + }; + // The panel 6 m sunward of the plate's center, listed before or after the + // plate, shadows it on both sides of center. + let sunward = Vector3::new(6.0, 0.0, 0.0); + for z in [-0.4, 0.4] { + let before = at_z(vec![tipped(sunward), at(plate, Vector3::zeros())], z); + let after = at_z(vec![at(plate, Vector3::zeros()), tipped(sunward)], z); + assert_eq!((before, after), (0.0, 0.0), "{name}: z = {z}"); + } + // The same panel 12 m to the side does not. + let side = at_z( + vec![ + at(plate, Vector3::zeros()), + tipped(Vector3::new(0.0, -12.0, 0.0)), + ], + 0.0, + ); + assert!( + (side - full).abs() < 1e-4 * full, + "{name}: {side} vs {full}" + ); + } +} + +#[test] +fn lighting_l0_draws_no_shadows() { + let _devices = common::device_lock(); + let full = lit(1.0); + for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] { + let mut renderer = renderer(&adapter, LightingTier::L0); + let name = renderer.adapter_info().name.clone(); + let model = renderer.load_model(&plate_under(2.0)).unwrap(); + let frame = renderer + .submit(&scene(vec![at(model, Vector3::zeros())]), &slanted()) + .unwrap() + .wait() + .unwrap(); + for (z, r) in column(&slanted(), &frame.radiance.unwrap()) { + if z > -1.3 { + assert!((r - full).abs() < 1e-4 * full, "{name}: z = {z}: {r}"); + } + } + } +} diff --git a/docs/optics.md b/docs/optics.md index 2cb72ca..39d69ed 100644 --- a/docs/optics.md +++ b/docs/optics.md @@ -369,7 +369,8 @@ 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. +only approximate those, for example with the shadow maps lighting L1 uses +(below). ### One sample per pixel @@ -580,6 +581,87 @@ 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. +### Shadows on the target + +Code: `shadow.rs`, `gpu/shadow.rs`, `shaders/mesh.wgsl`. From lighting L1. + +A solar array shades the bus behind it, an antenna dish shades its feed, +and a second spacecraft close by can shade the first. Rasterization draws +each triangle on its own, so by itself it cannot know that something else +stands between a point and the Sun. The renderer finds out with a shadow +map. + +- **Shadow map.** Before the camera's view is drawn, the target is drawn + once more as the Sun sees it, looking straight down the Sun's rays. That + drawing keeps only depth: for each texel (a pixel of the map), how far + along the Sun's direction the first surface is. The Sun is far enough + away that its rays are parallel across a spacecraft, so the view has no + perspective, and a map 4096 texels across covers the target's bounding + sphere. For Landsat 8, whose sphere is 16.6 m across, one texel is 4 mm. +- **Lookup.** When the camera's view shades a point, it finds that point + in the map. If the map holds a surface nearer the Sun than the point, the + point is in shadow. + +Each target gets its own map. Another target is drawn into that map when it +is sunward of this one and close enough across the Sun's rays to shade it. +A target off to the side is left out, since it cannot cast a shadow here. + +The Sun is a disk half a degree across, so a shadow edge is not sharp. A +point near the edge sees part of the disk past the blocker, the same +penumbra as in an eclipse. The width of the penumbra is the distance from +blocker to surface times the tangent of the Sun's angular diameter, 0.93 % +at 1 AU. A solar array 2 m above the bus casts an edge 19 mm wide. One 10 m +away casts an edge 9 cm wide. The renderer reproduces this with +percentage-closer soft shadows: + +1. Search the map around the point for blockers, and average their depth. + That gives the blocker-to-surface distance, and from it the penumbra + width. +2. Compare the point's depth against a 6 x 6 grid of map texels spread over + that width, and take the fraction that are not blocked as the share of + the Sun the point sees. + +Both steps use fixed grids rather than random samples, so the image has no +noise and the same scene always gives the same pixels. + +A surface can wrongly shadow itself, because the map only stores its depth +once per texel. Where the surface is tilted to the Sun, depth changes +across a texel, and half of each texel looks like it is behind itself. The +result is a pattern of dark stripes called shadow acne. The renderer +prevents it three ways. The stored depth is pushed away from the Sun in +proportion to the surface's tilt. The lookup is made from a point 1.5 +texels off the surface, toward the Sun. And a blocker has to be more than 2 +texels nearer the Sun to count. + +Tests (`tests/shadows.rs`) check, on the GPU and lavapipe: + +- A plate under a 2 m square: the shadow reads zero direct light, and its + edge falls where the square's edge projects, to within one texel. +- The measured penumbra width is within 2.2 % of 0.93 % of the occluder's + height, at heights of 1, 2 and 3 m. +- No pixel darkens compared with L0 on a plate tilted 80 and 87 degrees + from the Sun, or on a faceted sphere. +- A target sunward of another shades it; one to the side does not. + +On the RTX 3060, shadows cost 0.14 ms of GPU time per frame for Landsat 8 +(85,000 triangles) and 0.34 ms for ICESat-2 (752,000), plus about 0.15 ms +in the camera view's lookups. Maps 2048 texels across save about 0.15 ms. + +Limitations: + +- **Penumbra shape.** The grid averages a square of the map, so brightness + rises in a straight line across the penumbra. A real disk-shaped Sun + rises along an S-shaped curve: the full width is the same, but the middle + 80 % of the change is 0.69 of the width rather than 0.8. +- **Resolution.** A detail smaller than a texel casts no reliable shadow. A + gap narrower than about 2 texels between a blocker and the surface it + shades (a few mm on Landsat 8) lets light through, because of the acne + margins above. +- **One sample per pixel.** Shadows are sampled at the pixel center, like + everything else, until M2's pixel-area sampling. +- **Bodies.** Earth and the Moon do not use shadow maps. Their shadows on + the target are the eclipse factor below, applied to the whole target. + ### Eclipses Code: `eclipse.rs`. @@ -641,7 +723,6 @@ velocity as an input, which the scene does not carry yet. | 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 | | Pixel-area sampling (anti-aliasing): graded edges, thin parts | M2 | | Lens distortion | M3 |