diff --git a/crates/resolved-cli/src/main.rs b/crates/resolved-cli/src/main.rs index f48e9b7..876b625 100644 --- a/crates/resolved-cli/src/main.rs +++ b/crates/resolved-cli/src/main.rs @@ -11,8 +11,8 @@ use std::time::Instant; use clap::{Args, Parser, Subcommand}; use resolved::gpu::AdapterChoice; use resolved::{ - CameraTier, EarthDetail, EarthTexture, Fidelity, Frame, LightingTier, Look, MoonTexture, - Preview, Renderer, RendererOptions, StarCatalog, StarPhotometry, + CameraTier, CloudMap, EarthDetail, EarthTexture, Fidelity, Frame, LightingTier, Look, + MoonTexture, Preview, Renderer, RendererOptions, StarCatalog, StarPhotometry, }; use crate::scene_file::SceneFile; @@ -105,8 +105,9 @@ struct RenderArgs { #[arg(long, value_name = "FRAMES")] timings: Option, /// Directory of data the tiers need (earth/bmng-2004-MM.jpg from - /// scripts/fetch-earth.sh, moon/wac-hapke-643.png from scripts/fetch-moon.py, - /// stars/stars.bin from scripts/fetch-stars.py). + /// scripts/fetch-earth.sh, clouds/blue-marble.png from scripts/fetch-clouds.sh, + /// moon/wac-hapke-643.png from scripts/fetch-moon.py, stars/stars.bin from + /// scripts/fetch-stars.py). /// Without it: RESOLVED_DATA, then ./data, then the repository's data/. #[arg(long, value_name = "DIR")] data: Option, @@ -244,6 +245,17 @@ fn render(args: &RenderArgs) -> Result<(), String> { } else { None }; + let clouds = if args.lighting >= LightingTier::L1 && file.appearance.clouds != Some(false) { + let image = data_file( + args.data.as_deref(), + "clouds/blue-marble.png", + "fetch-clouds.sh", + "lighting L1", + )?; + Some(Arc::new(CloudMap::load(&image).map_err(|e| e.to_string())?)) + } else { + None + }; let stars = if args.lighting >= LightingTier::L1 && args.camera >= CameraTier::C1 { let catalog = data_file( args.data.as_deref(), @@ -276,6 +288,7 @@ fn render(args: &RenderArgs) -> Result<(), String> { earth_texture, earth_detail: earth_detail.clone(), moon_texture, + clouds, appearance: file.appearance(), read_radiance: false, profile: timings.is_some(), diff --git a/crates/resolved-cli/src/scene_file.rs b/crates/resolved-cli/src/scene_file.rs index 994399f..fedf51b 100644 --- a/crates/resolved-cli/src/scene_file.rs +++ b/crates/resolved-cli/src/scene_file.rs @@ -150,6 +150,9 @@ pub struct AppearanceFile { /// Which month (1-12) Earth's surface shows at lighting L1, for scenes without an /// epoch; with one, the epoch's month is the default. pub earth_month: Option, + /// Whether Earth has clouds at lighting L1 (data/clouds/blue-marble.png, from + /// scripts/fetch-clouds.sh); true if left out. + pub clouds: Option, } /// One spacecraft. diff --git a/crates/resolved-cli/tests/cli.rs b/crates/resolved-cli/tests/cli.rs index 9fab5f1..88317e8 100644 --- a/crates/resolved-cli/tests/cli.rs +++ b/crates/resolved-cli/tests/cli.rs @@ -44,6 +44,11 @@ fn gray_earth() -> PathBuf { .save(dir.join(format!("earth/bmng-2004-{month:02}.jpg"))) .unwrap(); } + // Clear sky: no clouds. + std::fs::create_dir_all(dir.join("clouds")).unwrap(); + image::GrayImage::from_pixel(64, 32, image::Luma([0])) + .save(dir.join("clouds/blue-marble.png")) + .unwrap(); // A uniform Moon: I/F 0.06 (16-bit values of I/F x 100,000). std::fs::create_dir_all(dir.join("moon")).unwrap(); image::ImageBuffer::, _>::from_pixel(64, 32, image::Luma([6000])) @@ -254,6 +259,51 @@ fn lighting_l1_names_a_missing_moon_map() { ); } +#[test] +fn lighting_l1_names_a_missing_cloud_map_unless_clouds_are_off() { + // Earth's images and the Moon's map but no clouds: the error names the file and the + // script, and a scene that turns clouds off renders without them. + let dir = out_dir().join("no-clouds"); + for name in ["earth/bmng-2004-07.jpg", "moon/wac-hapke-643.png"] { + let to = dir.join(name); + std::fs::create_dir_all(to.parent().unwrap()).unwrap(); + std::fs::copy(gray_earth().join(name), to).unwrap(); + } + let render = |scene: &str| { + resolved(&[ + "render", + "--scene", + scene, + "--lighting", + "L1", + "--data", + dir.to_str().unwrap(), + "--preview", + out_dir().join("no-clouds.png").to_str().unwrap(), + ]) + }; + let result = render("examples/scenes/boxwing.toml"); + assert!(!result.status.success()); + let stderr = String::from_utf8_lossy(&result.stderr); + assert!( + stderr.contains("clouds/blue-marble.png") && stderr.contains("fetch-clouds.sh"), + "{stderr}" + ); + let root = Path::new(env!("CARGO_MANIFEST_DIR")).join("../.."); + let scene = std::fs::read_to_string(root.join("examples/scenes/boxwing.toml")).unwrap(); + let scene = scene.replace("earth_month = 7", "earth_month = 7\nclouds = false"); + let models = root.join("examples/models/"); + let scene = scene.replace("../models/", models.to_str().unwrap()); + let path = out_dir().join("boxwing-no-clouds.toml"); + std::fs::write(&path, scene).unwrap(); + let result = render(path.to_str().unwrap()); + assert!( + result.status.success(), + "{}", + String::from_utf8_lossy(&result.stderr) + ); +} + #[test] fn errors_are_reported() { let missing = resolved(&["render", "--scene", "no/such.toml", "--preview", "x.png"]); diff --git a/crates/resolved/src/clouds.rs b/crates/resolved/src/clouds.rs new file mode 100644 index 0000000..6988ce3 --- /dev/null +++ b/crates/resolved/src/clouds.rs @@ -0,0 +1,179 @@ +//! Clouds over Earth, as a thin layer at a fixed height (`shaders/clouds.wgsl`). +//! +//! The map is equirectangular like Earth's ([`crate::EarthTexture`]): longitude -180 to +//! +180 degrees left to right, latitude +90 at the top. Each pixel holds the cloud's +//! reflectance there, 0 for clear sky, as in the cloud layer of NASA's Blue Marble +//! (`scripts/fetch-clouds.sh`). The layer is gray and reflects like a Lambertian +//! surface, and what it does not reflect it passes on (clouds absorb almost no visible +//! light): a ground seen through it, or lit through it, takes 1 less its reflectance. + +use std::path::Path; + +use half::f16; + +use crate::earth::parallel_rows; +use crate::Error; + +/// The layer's height above the ground, in meters. Real cloud tops range from about +/// 1 to 15 km; one height stands for them all. +pub const CLOUD_HEIGHT: f64 = 5e3; + +/// The constant in the cloud's reflectance for optical thickness tau, +/// R = k tau / (1 + k tau): sqrt(3) / 2 (1 - g) for scattering with asymmetry +/// g = 0.85, as in Lacis and Hansen (1974), *J. Atmos. Sci.* 31, 118. +pub const REFLECTANCE_SCALE: f64 = 0.13; + +/// The reflectance of a cloud of optical thickness `tau`. +pub fn reflectance(tau: f64) -> f64 { + let x = REFLECTANCE_SCALE * tau; + x / (1.0 + x) +} + +/// The optical thickness of a cloud of reflectance `r` (below 1). +pub fn optical_thickness(r: f64) -> f64 { + r / (REFLECTANCE_SCALE * (1.0 - r)) +} + +/// The constants `shaders/clouds.wgsl` expects, then the file itself. +pub(crate) fn wgsl() -> String { + format!( + "const CLOUD_HEIGHT: f32 = {:?};\n", + (CLOUD_HEIGHT / 1e3) as f32 + ) + include_str!("shaders/clouds.wgsl") +} + +/// An equirectangular map of the clouds' reflectance. +#[derive(Clone)] +pub struct CloudMap { + width: u32, + height: u32, + /// Reflectance times 255, rows from the top. + values: Vec, +} + +impl std::fmt::Debug for CloudMap { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("CloudMap") + .field("width", &self.width) + .field("height", &self.height) + .finish_non_exhaustive() + } +} + +impl CloudMap { + /// Loads an 8-bit grayscale image (PNG or JPEG) of reflectance times 255, as + /// `scripts/fetch-clouds.sh` writes. + pub fn load(path: &Path) -> Result { + let image = + image::open(path).map_err(|e| Error::Texture(format!("{}: {e}", path.display())))?; + let gray = image.into_luma8(); + Self::from_values(gray.width(), gray.height(), gray.into_raw()) + } + + /// A map from reflectance times 255, rows from the top. + pub fn from_values(width: u32, height: u32, values: Vec) -> Result { + if width == 0 || height == 0 || values.len() != width as usize * height as usize { + return Err(Error::Texture(format!( + "{width} x {height} needs {} values, got {}", + width as usize * height as usize, + values.len() + ))); + } + Ok(Self { + width, + height, + values, + }) + } + + /// Width in pixels. + pub fn width(&self) -> u32 { + self.width + } + + /// Height in pixels. + pub fn height(&self) -> u32 { + self.height + } + + /// The mean reflectance over the sphere, weighting each row by the area it covers. + pub fn mean(&self) -> f64 { + let (mut sum, mut weight) = (0.0, 0.0); + for (y, row) in self.values.chunks(self.width as usize).enumerate() { + let lat = std::f64::consts::PI * (0.5 - (y as f64 + 0.5) / self.height as f64); + let w = lat.cos(); + sum += w * row.iter().map(|&v| f64::from(v)).sum::() / 255.0; + weight += w * row.len() as f64; + } + sum / weight + } + + /// Every mip level from full size down to 1 x 1, each averaging 2 x 2 blocks of the + /// one above, as reflectance in f16. Averaging reflectance, not optical thickness, + /// keeps a far view's brightness. + pub(crate) fn mip_chain(&self) -> Vec<(u32, u32, Vec)> { + let first: Vec = self + .values + .iter() + .map(|&v| f16::from_f32(f32::from(v) / 255.0)) + .collect(); + let mut levels = vec![(self.width, self.height, first)]; + loop { + let (w, h, above) = levels.last().expect("at least one level"); + let (w, h) = (*w, *h); + if w == 1 && h == 1 { + break; + } + let (nw, nh) = ((w / 2).max(1), (h / 2).max(1)); + let mut below = vec![f16::ZERO; nw as usize * nh as usize]; + parallel_rows(&mut below, nw as usize, |y, row| { + let y = y as u32; + for (x, out) in row.iter_mut().enumerate() { + let x = x as u32; + let mut sum = 0.0; + for (dx, dy) in [(0, 0), (1, 0), (0, 1), (1, 1)] { + let sx = (2 * x + dx).min(w - 1); + let sy = (2 * y + dy).min(h - 1); + sum += above[(sy * w + sx) as usize].to_f32(); + } + *out = f16::from_f32(sum / 4.0); + } + }); + levels.push((nw, nh, below)); + } + levels + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn reflectance_and_thickness_invert() { + for tau in [0.5, 5.0, 20.0, 80.0] { + assert_relative_eq!( + optical_thickness(reflectance(tau)), + tau, + max_relative = 1e-12 + ); + } + // Thin cloud reflects about 0.13 tau; thick cloud most of the light. + assert_relative_eq!(reflectance(0.1), 0.01283, max_relative = 1e-3); + assert!(reflectance(50.0) > 0.85); + } + + #[test] + fn mips_average_reflectance() { + let map = CloudMap::from_values(4, 2, vec![0, 255, 0, 255, 255, 255, 0, 0]).unwrap(); + let levels = map.mip_chain(); + assert_eq!(levels.len(), 3); + let (w, h, top) = &levels[1]; + assert_eq!((*w, *h), (2, 1)); + assert_eq!(top[0].to_f32(), 0.75); + assert_eq!(top[1].to_f32(), 0.25); + assert_eq!(levels[2].2[0].to_f32(), 0.5); + assert!(CloudMap::from_values(4, 2, vec![0; 7]).is_err()); + } +} diff --git a/crates/resolved/src/gpu/bodies.rs b/crates/resolved/src/gpu/bodies.rs index 46233b1..0e5b73d 100644 --- a/crates/resolved/src/gpu/bodies.rs +++ b/crates/resolved/src/gpu/bodies.rs @@ -18,6 +18,7 @@ use crate::atmosphere::{ }; use crate::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS, MOON_RADIUS}; use crate::camera::Intrinsics; +use crate::clouds::CloudMap; use crate::earth::{DetailFormat, DetailLayout, EarthDetail, EarthTexture}; use crate::model::rgb_spectrum; use crate::moon::{MoonTexture, NORMALIZED_INCIDENCE, PHASE_SCALE}; @@ -362,6 +363,12 @@ pub struct BodiesPass { /// Earth's shadow around the Moon ([`EarthShadow`]), for the atmosphere's variant: /// its range, then the light's log per entry, in f32. earth_shadow: wgpu::Buffer, + /// The clouds' reflectance (a 1 x 1 clear placeholder until a map is set), and + /// whether a map is set, which compiles the cloud layer into the pipelines. + clouds: wgpu::TextureView, + clouded: bool, + /// The cloud map's height in texels. + cloud_height: u32, bind_group: wgpu::BindGroup, } @@ -388,7 +395,7 @@ impl BodiesPass { /// buffer. pub fn new(ctx: &GpuContext, samples: u32, area: bool) -> Self { let device = &ctx.device; - let pipeline = compile(ctx, samples, area, false); + let pipeline = compile(ctx, samples, area, false, false); let uniforms = device.create_buffer(&wgpu::BufferDescriptor { label: Some("bodies uniforms"), size: std::mem::size_of::() as u64, @@ -447,7 +454,9 @@ impl BodiesPass { &atmosphere, (&detail, &detail_sampler), None, + None, ); + let clouds = upload_clouds(ctx, &CloudMap::from_values(1, 1, vec![0]).expect("1 x 1")); Self { pipeline, bending: None, @@ -465,6 +474,9 @@ impl BodiesPass { detail_sampler, grazing, earth_shadow, + clouds, + clouded: false, + cloud_height: 1, bind_group, } } @@ -476,6 +488,25 @@ impl BodiesPass { self.rebind(ctx); } + /// Uploads the clouds' map with its mip chain and draws the cloud layer over Earth + /// in later frames, where Earth is seen through its atmosphere. + pub fn set_cloud_map(&mut self, ctx: &GpuContext, map: &CloudMap) { + self.clouds = upload_clouds(ctx, map); + self.clouded = true; + self.cloud_height = map.height(); + self.pipeline = compile(ctx, self.samples, self.area, false, true); + if let Some((pipeline, _)) = &mut self.bending { + *pipeline = compile(ctx, self.samples, self.area, true, true); + } + self.rebind(ctx); + } + + /// The clouds' map (a clear 1 x 1 placeholder until one is set), and its height in + /// texels if one is set, for the Earthshine pass. + pub fn clouds(&self) -> (&wgpu::TextureView, Option) { + (&self.clouds, self.clouded.then_some(self.cloud_height)) + } + /// Uploads the Moon's map with its mip chain and binds it for later frames. pub fn set_moon_texture(&mut self, ctx: &GpuContext, texture: &MoonTexture) { self.moon = upload_moon(ctx, texture); @@ -597,7 +628,7 @@ impl BodiesPass { .collect(); ctx.queue .write_buffer(&self.grazing, 0, bytemuck::cast_slice(&values)); - let pipeline = compile(ctx, self.samples, self.area, true); + let pipeline = compile(ctx, self.samples, self.area, true, self.clouded); let group = self.bind_to(ctx, &pipeline, true); self.bending = Some((pipeline, group)); self.rebind(ctx); @@ -652,6 +683,7 @@ impl BodiesPass { &self.atmosphere, (&self.detail, &self.detail_sampler), bending.then_some((&self.grazing, &self.earth_shadow)), + self.clouded.then_some(&self.clouds), ) } @@ -777,6 +809,20 @@ fn upload_earth(ctx: &GpuContext, texture: &EarthTexture) -> wgpu::TextureView { ) } +/// Uploads the clouds' map with its mip chain as one channel of f16 reflectance. +fn upload_clouds(ctx: &GpuContext, map: &CloudMap) -> wgpu::TextureView { + let levels = map.mip_chain(); + upload( + ctx, + "clouds", + wgpu::TextureFormat::R16Float, + 2, + levels + .iter() + .map(|(w, h, p)| (*w, *h, bytemuck::cast_slice(p))), + ) +} + /// Uploads the Moon's map with its mip chain as one channel of f16 reflectance. fn upload_moon(ctx: &GpuContext, texture: &MoonTexture) -> wgpu::TextureView { let levels = texture.mip_chain(); @@ -840,9 +886,15 @@ fn upload<'a>( } /// The bodies pipeline for mesh targets of `samples` samples per pixel, casting grids of -/// rays over pixels on bodies' edges if `area`, and with the Moon and the Sun seen -/// through Earth's atmosphere if `bending`. -fn compile(ctx: &GpuContext, samples: u32, area: bool, bending: bool) -> wgpu::RenderPipeline { +/// rays over pixels on bodies' edges if `area`, with the Moon and the Sun seen through +/// Earth's atmosphere if `bending`, and with the cloud layer if `clouds`. +fn compile( + ctx: &GpuContext, + samples: u32, + area: bool, + bending: bool, + clouds: bool, +) -> wgpu::RenderPipeline { let device = &ctx.device; let (radiance, depth) = if samples > 1 { ( @@ -873,8 +925,13 @@ fn compile(ctx: &GpuContext, samples: u32, area: bool, bending: bool) -> wgpu::R include_str!("../shaders/bending.wgsl") } else { STRAIGHT + } + + &if clouds { + crate::clouds::wgsl() + include_str!("../shaders/cloudy.wgsl") + } else { + CLEAR.to_owned() }) - .into(), + .into(), ), }); device.create_render_pipeline(&wgpu::RenderPipelineDescriptor { @@ -908,6 +965,24 @@ fn compile(ctx: &GpuContext, samples: u32, area: bool, bending: bool) -> wgpu::R }) } +/// `through_clouds` where there are no clouds: the ground lit and seen through the air +/// alone. +const CLEAR: &str = " +fn through_clouds( + camera: vec3, + ground: vec3, + sun: vec3, + albedo: vec4, + lit: f32, + at: MapLookup, +) -> vec4 { + let air = atmosphere_to_ground(camera, normalize(ground - camera), ground, sun); + let light = atmosphere_ground_light(ground, sun); + let surface = u.earth_irradiance * albedo * (light[0] * lit + light[1]); + return surface * air.transmittance + PI * u.earth_irradiance * air.scattered; +} +"; + /// `through_air` and `earth_shadow` where no body is seen or lit through the /// atmosphere: every ray straight from the camera and whole, and Earth's shadow its hard /// edge's. @@ -934,6 +1009,7 @@ fn bind( atmosphere: &AtmosphereTextures, detail: (&[wgpu::TextureView; DETAIL_GROUPS], &wgpu::Sampler), air: Option<(&wgpu::Buffer, &wgpu::Buffer)>, + clouds: Option<&wgpu::TextureView>, ) -> wgpu::BindGroup { let entries = [ wgpu::BindGroupEntry { @@ -986,6 +1062,10 @@ fn bind( }, ] })) + .chain(clouds.map(|view| wgpu::BindGroupEntry { + binding: 17, + resource: wgpu::BindingResource::TextureView(view), + })) .collect(); ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { label: Some("bodies"), diff --git a/crates/resolved/src/gpu/earthshine.rs b/crates/resolved/src/gpu/earthshine.rs index c33f2f9..b8821c4 100644 --- a/crates/resolved/src/gpu/earthshine.rs +++ b/crates/resolved/src/gpu/earthshine.rs @@ -110,6 +110,8 @@ struct Globals { grid: [u32; 4], moon_color: [f32; 4], moon_params: [f32; 4], + // x: 1 if the cloud layer is drawn. y: meters per texel of its map in latitude. + clouds: [f32; 4], } /// The Earthshine pipelines and the per-target buffers and cube maps they fill. @@ -187,6 +189,14 @@ impl EarthshinePass { 8, wgpu::ShaderStages::COMPUTE, )) + .chain([entry( + 13, + wgpu::BindingType::Texture { + sample_type: wgpu::TextureSampleType::Float { filterable: true }, + view_dimension: wgpu::TextureViewDimension::D2, + multisampled: false, + }, + )]) .collect(); let layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { label: Some("earthshine"), @@ -202,6 +212,7 @@ impl EarthshinePass { source: wgpu::ShaderSource::Wgsl( (include_str!("../shaders/common.wgsl").to_owned() + &atmosphere::wgsl(8) + + &crate::clouds::wgsl() + include_str!("../shaders/earthshine.wgsl")) .into(), ), @@ -320,6 +331,7 @@ impl EarthshinePass { ); let (texture, sampler, height) = bodies.earth_texture(); let moon = bodies.moon_texture(); + let (clouds, cloud_height) = bodies.clouds(); let globals = Globals { earth_albedo: radiometry.earth_albedo, earth_rgb: radiometry.earth_rgb, @@ -338,6 +350,14 @@ impl EarthshinePass { 0.0, 0.0, ], + clouds: cloud_height.map_or([0.0; 4], |height| { + [ + 1.0, + (std::f64::consts::PI * EARTH_EQUATORIAL_RADIUS / height as f64) as f32, + 0.0, + 0.0, + ] + }), }; ctx.queue .write_buffer(&self.globals, 0, bytemuck::bytes_of(&globals)); @@ -384,6 +404,10 @@ impl EarthshinePass { let entries: Vec<_> = entries .into_iter() .chain(bodies.atmosphere().entries(8)) + .chain([wgpu::BindGroupEntry { + binding: 13, + resource: wgpu::BindingResource::TextureView(clouds), + }]) .collect(); let group = ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { label: Some("earthshine"), diff --git a/crates/resolved/src/lib.rs b/crates/resolved/src/lib.rs index aa7f1e7..1c7ecbd 100644 --- a/crates/resolved/src/lib.rs +++ b/crates/resolved/src/lib.rs @@ -20,6 +20,7 @@ pub use wgpu; pub mod atmosphere; pub mod bodies; pub mod camera; +pub mod clouds; pub mod earth; pub mod eclipse; #[cfg(feature = "astrokit")] @@ -44,6 +45,7 @@ pub mod sun; pub use camera::{ Blooming, BrownConrady, CameraConfig, Channel, Detector, Intrinsics, Psf, PsfImage, }; +pub use clouds::CloudMap; pub use earth::{EarthDetail, EarthTexture}; pub use error::Error; pub use fidelity::{CameraTier, Fidelity, LightingTier}; diff --git a/crates/resolved/src/render.rs b/crates/resolved/src/render.rs index 704566f..18e0218 100644 --- a/crates/resolved/src/render.rs +++ b/crates/resolved/src/render.rs @@ -16,6 +16,7 @@ use crate::camera::{Blooming, CameraConfig, Intrinsics, Psf}; use std::sync::Arc; use crate::atmosphere::{AtmosphereTables, EarthShadow}; +use crate::clouds::CloudMap; use crate::earth::{EarthDetail, EarthTexture}; use crate::eclipse::{ moon_near_shadow, sun_visible_fraction, sun_visible_light, sun_visible_through_air, @@ -60,6 +61,9 @@ pub struct RendererOptions { pub earth_detail: Option>, /// The Moon's surface, needed from lighting tier L1 up; ignored below it. pub moon_texture: Option>, + /// Clouds over Earth, drawn from lighting tier L1 where Earth is seen through its + /// atmosphere; optional. Blue Marble's 16384 x 8192 map takes 360 MB of GPU memory. + pub clouds: Option>, /// How Earth and the Moon reflect light. pub appearance: Appearance, /// Also read back the channel radiance image each frame, as the pixels receive it @@ -461,6 +465,9 @@ impl Renderer { if let Some(texture) = &options.moon_texture { renderer.set_moon_texture(texture); } + if let Some(map) = &options.clouds { + renderer.set_cloud_map(map); + } Ok(renderer) } @@ -528,6 +535,20 @@ impl Renderer { self.body_radiometry.moon_textured = true; } + /// Uploads the clouds' map with its mip chain and draws the cloud layer over Earth + /// in later frames, as [`RendererOptions::clouds`] does at setup. + /// + /// The layer is drawn only where Earth is seen through its atmosphere, from lighting + /// tier L1 with [`crate::bodies::Appearance::earth_atmosphere`]; otherwise this does + /// nothing. + pub fn set_cloud_map(&mut self, map: &CloudMap) { + if !self.body_radiometry.atmosphere { + log::warn!("clouds are drawn only with Earth's atmosphere; ignoring them"); + return; + } + self.bodies.set_cloud_map(&self.ctx, map); + } + /// Changes the exposure time, in seconds, from the next frame submitted on: what an /// auto-exposure loop does between frames. Only the conversion from radiance to /// electrons depends on it, so this costs nothing beyond a few multiplications. diff --git a/crates/resolved/src/shaders/atmosphere.wgsl b/crates/resolved/src/shaders/atmosphere.wgsl index bf33c21..fc10d26 100644 --- a/crates/resolved/src/shaders/atmosphere.wgsl +++ b/crates/resolved/src/shaders/atmosphere.wgsl @@ -186,44 +186,71 @@ fn atmosphere_sky(camera: vec3, view: vec3, sun: vec3) -> Atmosph return light; } -// Along a ray from `camera` in unit direction `view` to the point `ground` on the -// ground, with the Sun along unit `sun`: the scattering along the whole ray from where -// it enters the atmosphere, less that beyond the ground point seen through the -// transmittance to it (Bruneton's `GetSkyRadianceToPoint`). -fn atmosphere_to_ground( - camera: vec3, - view: vec3, - ground: vec3, - sun: vec3, -) -> AtmosphereLight { +// Where a ray from `camera` in unit direction `view` toward the ground enters the +// atmosphere, with the Sun along unit `sun`: the point, its parameters and the +// scattering from there along the whole ray, for `atmosphere_to_point`. +struct AtmosphereEntry { + x: vec3, + // r, mu, mu_s, nu at `x`. + rmu: vec4, + rayleigh: vec4, + mie: vec4, +} + +fn atmosphere_entering(camera: vec3, view: vec3, sun: vec3) -> AtmosphereEntry { let entry = max(atmosphere_entry(camera, view), 0.0); let x = camera + entry * view; let r = clamp(length(x), ATMOSPHERE_BOTTOM, ATMOSPHERE_TOP); let mu = clamp(dot(x, view) / r, -1.0, 1.0); let mu_s = clamp(dot(x, sun) / r, -1.0, 1.0); let nu = clamp(dot(view, sun), -1.0, 1.0); - let d = length(ground - x); + let uvwz = atmosphere_scattering_uvwz(r, mu, mu_s, nu, true); + return AtmosphereEntry( + x, + vec4(r, mu, mu_s, nu), + atmosphere_lookup(atmosphere_rayleigh_table, uvwz), + atmosphere_lookup(atmosphere_mie_table, uvwz), + ); +} + +// Along the ray of `entry` to the point `point` on it: the scattering along the whole +// ray from where it enters the atmosphere, less that beyond the point seen through the +// transmittance to it (Bruneton's `GetSkyRadianceToPoint`). +fn atmosphere_to_point(entry: AtmosphereEntry, point: vec3) -> AtmosphereLight { + let r = entry.rmu.x; + let mu = entry.rmu.y; + let mu_s = entry.rmu.z; + let nu = entry.rmu.w; + let d = length(point - entry.x); var light: AtmosphereLight; light.transmittance = atmosphere_transmittance(r, mu, d, true); - let uvwz = atmosphere_scattering_uvwz(r, mu, mu_s, nu, true); - let rayleigh = atmosphere_lookup(atmosphere_rayleigh_table, uvwz); - let mie = atmosphere_lookup(atmosphere_mie_table, uvwz); let r_p = clamp(sqrt(d * d + 2.0 * r * mu * d + r * r), ATMOSPHERE_BOTTOM, ATMOSPHERE_TOP); let mu_p = clamp((r * mu + d) / r_p, -1.0, 1.0); let mu_s_p = clamp((r * mu_s + d * nu) / r_p, -1.0, 1.0); let uvwz_p = atmosphere_scattering_uvwz(r_p, mu_p, mu_s_p, nu, true); let rayleigh_p = atmosphere_lookup(atmosphere_rayleigh_table, uvwz_p); let mie_p = atmosphere_lookup(atmosphere_mie_table, uvwz_p); - let rayleigh_between = max(rayleigh - light.transmittance * rayleigh_p, vec4(0.0)); + let rayleigh_between = max(entry.rayleigh - light.transmittance * rayleigh_p, vec4(0.0)); // Bruneton's guard against artifacts in the Mie term with the Sun below the // horizon, where it is all but zero anyway. - let mie_between = max(mie - light.transmittance * mie_p, vec4(0.0)) + let mie_between = max(entry.mie - light.transmittance * mie_p, vec4(0.0)) * smoothstep(0.0, 0.01, mu_s); light.scattered = rayleigh_between * atmosphere_rayleigh_phase(nu) + mie_between * atmosphere_mie_phase(nu); return light; } +// Along a ray from `camera` in unit direction `view` to the point `ground` on the +// ground, with the Sun along unit `sun`. +fn atmosphere_to_ground( + camera: vec3, + view: vec3, + ground: vec3, + sun: vec3, +) -> AtmosphereLight { + return atmosphere_to_point(atmosphere_entering(camera, view, sun), ground); +} + // Sunlight on the ground at `point` with the Sun along unit `sun`, as a fraction of // the Sun's irradiance: x per channel the transmittance toward the Sun (times the // fraction of its disk above the horizon), y the sky's irradiance on a horizontal diff --git a/crates/resolved/src/shaders/bodies.wgsl b/crates/resolved/src/shaders/bodies.wgsl index 52e1d4b..3f17bd3 100644 --- a/crates/resolved/src/shaders/bodies.wgsl +++ b/crates/resolved/src/shaders/bodies.wgsl @@ -10,8 +10,8 @@ // // The code that includes this file first declares SAMPLES, AREA, the types of the mesh // targets, MeshRadiance and MeshDepth (multisampled when SAMPLES > 1), and -// SAMPLE_POSITIONS, where in a pixel each sample lies, and adds `through_air` and -// `earth_shadow` (bodies.rs). +// SAMPLE_POSITIONS, where in a pixel each sample lies, and adds `through_air`, +// `earth_shadow` and `through_clouds` (bodies.rs). struct Uniforms { // fx, fy, cx, cy in pixels. @@ -238,11 +238,7 @@ fn earth_radiance(p: vec3, n_fixed: vec3, at: MapLookup) -> vec4 } let camera = camera_stretched(); let ground = earth_stretched_km(transpose(u.earth_axes) * (p - u.earth_center.xyz), u.earth_shape.x); - let sun = u.earth_sun_stretched.xyz; - let air = atmosphere_to_ground(camera, normalize(ground - camera), ground, sun); - let light = atmosphere_ground_light(ground, sun); - let surface = u.earth_irradiance * albedo * (light[0] * lit + light[1]); - return (surface * air.transmittance + PI * u.earth_irradiance * air.scattered) * eclipse; + return through_clouds(camera, ground, u.earth_sun_stretched.xyz, albedo, lit, at) * eclipse; } // The light Earth's atmosphere scatters toward the camera along `d` where the ray does diff --git a/crates/resolved/src/shaders/clouds.wgsl b/crates/resolved/src/shaders/clouds.wgsl new file mode 100644 index 0000000..7fdfcfe --- /dev/null +++ b/crates/resolved/src/shaders/clouds.wgsl @@ -0,0 +1,55 @@ +// The cloud layer (clouds.rs), shared by bodies.wgsl and earthshine.wgsl, which read the +// map: a sphere CLOUD_HEIGHT km above the ground in the atmosphere's stretched frame, +// where Earth is a sphere, with lengths in km as atmosphere.wgsl takes them. The code +// that includes this file declares CLOUD_HEIGHT and includes atmosphere.wgsl first. + +const CLOUD_TOP: f32 = ATMOSPHERE_BOTTOM + CLOUD_HEIGHT; + +// Where the line from `point`, below the layer, along unit `dir` passes up through it. +fn cloud_crossing(point: vec3, dir: vec3) -> vec3 { + let b = dot(point, dir); + let c = dot(point, point) - CLOUD_TOP * CLOUD_TOP; + return point + (sqrt(max(b * b - c, 0.0)) - b) * dir; +} + +// Where in the map (equirectangular, longitude wrapped to 0..1) the layer's `point` +// lies, Earth's equatorial over polar radius being `flattening`: by the ellipsoid's +// normal there, as the ground's map is read. +fn cloud_uv(point: vec3, flattening: f32) -> vec2 { + let uv = equirectangular(normalize(point * vec3(1.0, 1.0, flattening))); + return vec2(fract(uv.x), uv.y); +} + +// What the camera at `camera` sees of the ground point `ground` under the layer, with +// the Sun along unit `sun`, per channel, in units of the Sun's irradiance over pi. +// The ground, of albedo `albedo` and with the Sun at cosine `cos_sun` to its normal, +// is lit through the cloud of reflectance `shade` where the sunlight crossed the +// layer, sky light included. The line of sight crosses the layer at `top`, through a +// cloud of reflectance `cover`: the cloud reflects the sunlight that reaches it there, +// and passes on 1 - `cover` of the light from below, the ground and the air under the +// layer. The air above the layer adds its own light and dims what is below it. +fn under_clouds( + camera: vec3, + ground: vec3, + sun: vec3, + albedo: vec4, + cos_sun: f32, + top: vec3, + cover: f32, + shade: f32, +) -> vec4 { + let entry = atmosphere_entering(camera, normalize(ground - camera), sun); + let air = atmosphere_to_point(entry, ground); + let light = atmosphere_ground_light(ground, sun); + let seen = albedo * (light[0] * cos_sun + light[1]) * ((1.0 - shade) * air.transmittance) + + PI * air.scattered; + if cover <= 0.0 { + return seen; + } + let air_top = atmosphere_to_point(entry, top); + let light_top = atmosphere_ground_light(top, sun); + let cos_top = max(dot(top, sun) / length(top), 0.0); + let above = PI * air_top.scattered; + let cloud = cover * (light_top[0] * cos_top + light_top[1]) * air_top.transmittance; + return above + cloud + (1.0 - cover) * max(seen - above, vec4(0.0)); +} diff --git a/crates/resolved/src/shaders/cloudy.wgsl b/crates/resolved/src/shaders/cloudy.wgsl new file mode 100644 index 0000000..e980e73 --- /dev/null +++ b/crates/resolved/src/shaders/cloudy.wgsl @@ -0,0 +1,30 @@ +// Earth's ground under the cloud layer: `through_clouds` for bodies.wgsl, in the +// variants of the bodies pass used once a cloud map is set (bodies.rs), with +// clouds.wgsl. + +// The clouds' reflectance (clouds.rs, `CloudMap`), equirectangular like Earth's map. +@group(0) @binding(17) var cloud_map: texture_2d; + +// The clouds' reflectance at the layer's `point`, read with the ground's footprint +// `at` (the two maps share their coordinates). +fn read_cloud(point: vec3, at: MapLookup) -> f32 { + let uv = cloud_uv(point, u.earth_shape.x); + return textureSampleGrad(cloud_map, earth_sampler, uv, at.dx, at.dy).r; +} + +// The channel radiance of ground point `ground` (stretched frame, km) of albedo +// `albedo`, with the Sun along `sun` at cosine `lit` to its normal, seen from `camera` +// through the air and the clouds, and lit through them. +fn through_clouds( + camera: vec3, + ground: vec3, + sun: vec3, + albedo: vec4, + lit: f32, + at: MapLookup, +) -> vec4 { + let top = cloud_crossing(ground, normalize(camera - ground)); + let cover = read_cloud(top, at); + let shade = read_cloud(cloud_crossing(ground, sun), at); + return u.earth_irradiance * under_clouds(camera, ground, sun, albedo, lit, top, cover, shade); +} diff --git a/crates/resolved/src/shaders/earthshine.wgsl b/crates/resolved/src/shaders/earthshine.wgsl index ae69c77..d585fe1 100644 --- a/crates/resolved/src/shaders/earthshine.wgsl +++ b/crates/resolved/src/shaders/earthshine.wgsl @@ -61,6 +61,9 @@ struct Globals { moon_color: vec4, // x: 1 if the Moon has a map. y: meters per texel of the map in latitude. moon_params: vec4, + // x: 1 if the cloud layer is drawn (with the atmosphere). y: meters per texel of + // its map in latitude. + clouds: vec4, } struct Sample { @@ -78,6 +81,8 @@ struct Sample { @group(0) @binding(5) var cube: texture_storage_2d_array; @group(0) @binding(6) var sun_disk: SunDisk; @group(0) @binding(7) var moon_texture: texture_2d; +// The clouds' reflectance (clouds.rs, `CloudMap`); a clear 1 x 1 placeholder without. +@group(0) @binding(13) var cloud_map: texture_2d; // Two unit vectors perpendicular to unit `a` and to each other. fn basis(a: vec3) -> mat2x3 { @@ -129,13 +134,29 @@ fn moon_albedo(t: TargetLight, d: vec3, tm: f32, n: vec3) -> vec4 // Earth's radiance, in units of sunlight / pi, where the ray from the target along `d` // meets the ground at distance `te`, with albedo `albedo` and the Sun at cosine // `cos_sun` to the ground's normal: lit by the sunlight that reaches the ground and by -// the sky, and seen through the air, which adds the sunlight it scatters (as the bodies -// pass draws it). -fn through_air(t: TargetLight, d: vec3, te: f32, albedo: vec4, cos_sun: f32) -> vec4 { +// the sky, and seen through the air, which adds the sunlight it scatters, and through +// the clouds if they are drawn (as the bodies pass draws it). The cloud map is read at +// the mip whose texels match `footprint` meters. +fn through_air( + t: TargetLight, + d: vec3, + te: f32, + albedo: vec4, + cos_sun: f32, + footprint: f32, +) -> vec4 { let to_fixed = transpose(t.earth_axes); let camera = earth_stretched_km(to_fixed * (-t.earth.center.xyz), g.params.z); let ground = earth_stretched_km(to_fixed * (te * d - t.earth.center.xyz), g.params.z); let sun = normalize(to_fixed * t.earth.sun.xyz * vec3(1.0, 1.0, g.params.z)); + if g.clouds.x > 0.5 { + let level = max(log2(footprint / g.clouds.y), 0.0); + let top = cloud_crossing(ground, normalize(camera - ground)); + let shade = cloud_crossing(ground, sun); + let cover = textureSampleLevel(cloud_map, earth_sampler, cloud_uv(top, g.params.z), level).r; + let shaded = textureSampleLevel(cloud_map, earth_sampler, cloud_uv(shade, g.params.z), level).r; + return under_clouds(camera, ground, sun, albedo, cos_sun, top, cover, shaded); + } let air = atmosphere_to_ground(camera, normalize(ground - camera), ground, sun); let light = atmosphere_ground_light(ground, sun); return albedo * (light[0] * cos_sun + light[1]) * air.transmittance + PI * air.scattered; @@ -185,7 +206,9 @@ fn sources(@builtin(global_invocation_id) id: vec3) { let lit = max(dot(normal, source.sun.xyz), 0.0) * visible; var light = albedo * lit; if !is_moon && g.params.w > 0.5 && te > 0.0 && !(tm > 0.0 && tm < te) { - light = through_air(t, d, te, albedo, max(dot(normal, source.sun.xyz), 0.0)) * visible; + let footprint = te * sqrt(source.cone.y / max(dot(normal, -d), 0.05)); + let cos_sun = max(dot(normal, source.sun.xyz), 0.0); + light = through_air(t, d, te, albedo, cos_sun, footprint) * visible; } var out: Sample; out.light = light * (source.sun.w * source.cone.y); diff --git a/crates/resolved/tests/gpu/atmosphere.rs b/crates/resolved/tests/gpu/atmosphere.rs index 8a41289..a7798df 100644 --- a/crates/resolved/tests/gpu/atmosphere.rs +++ b/crates/resolved/tests/gpu/atmosphere.rs @@ -22,7 +22,7 @@ const ADAPTERS: [AdapterChoice; 2] = [AdapterChoice::Gpu, AdapterChoice::Cpu]; /// Earth's surface: one sRGB gray level, 128, everywhere. const GRAY: u8 = 128; -fn renderer(adapter: &AdapterChoice) -> Renderer { +pub(crate) fn renderer(adapter: &AdapterChoice) -> Renderer { let pixels = [GRAY, GRAY, GRAY, 255].repeat(2); let options = RendererOptions { adapter: adapter.clone(), @@ -41,7 +41,7 @@ fn renderer(adapter: &AdapterChoice) -> Renderer { /// Earth unturned, so Earth-fixed axes are inertial, the Sun along `sun` at 1 AU, and the /// camera at `camera` looking at `at`. Everything here is in the equatorial plane, where /// Earth's flattening does not enter. -fn view(sun: Vector3, camera: Vector3, at: Vector3) -> (Scene, Pose) { +pub(crate) fn view(sun: Vector3, camera: Vector3, at: Vector3) -> (Scene, Pose) { let scene = Scene { earth_attitude: UnitQuaternion::identity(), sun: sun.normalize() * AU, @@ -57,7 +57,7 @@ fn view(sun: Vector3, camera: Vector3, at: Vector3) -> (Scene, Po /// symmetric about it: to second order, the radiance along the boresight. (One pixel's /// is half a pixel off, which at the limb moves the tangent point by 140 m and the glow /// by 1.7 %.) -fn center(frame: &resolved::Frame) -> [f32; 4] { +pub(crate) fn center(frame: &resolved::Frame) -> [f32; 4] { let radiance = frame.radiance.as_ref().unwrap(); let mut sum = [0.0f32; 4]; for (x, y) in [(959, 539), (960, 539), (959, 540), (960, 540)] { @@ -69,7 +69,7 @@ fn center(frame: &resolved::Frame) -> [f32; 4] { } /// Per channel: the Sun's irradiance at 1 AU over pi, and the gray surface's albedo. -fn constants() -> (Radiometry, [f64; 3], [f64; 3]) { +pub(crate) fn constants() -> (Radiometry, [f64; 3], [f64; 3]) { let r = Radiometry::new(&CameraConfig::example()); let v = f64::from(GRAY) / 255.0; let linear = ((v + 0.055) / 1.055).powf(2.4); @@ -79,7 +79,7 @@ fn constants() -> (Radiometry, [f64; 3], [f64; 3]) { } /// Where the ray from `c` along unit `v` meets a sphere of radius `radius` first. -fn hit(c: Vector3, v: Vector3, radius: f64) -> f64 { +pub(crate) fn hit(c: Vector3, v: Vector3, radius: f64) -> f64 { let b = c.dot(&v); let q = (c - b * v).norm(); -b - ((radius - q) * (radius + q)).sqrt() @@ -93,13 +93,28 @@ fn ground_reference( view: Vector3, sun: Vector3, ) -> [f64; 3] { - let (_, irradiance, albedo) = constants(); + let albedo = constants().2; + surface_reference(tables, camera, view, sun, BOTTOM_RADIUS, albedo) +} + +/// The same for a Lambertian surface of `albedo` per channel on the sphere of `radius` +/// (meters), lit by the sunlight and sky light that reach that height: per channel, the +/// channel radiance the camera sees. +pub(crate) fn surface_reference( + tables: &AtmosphereTables, + camera: Vector3, + view: Vector3, + sun: Vector3, + radius: f64, + albedo: [f64; 3], +) -> [f64; 3] { + let irradiance = constants().1; let x = camera + hit(camera, view, TOP_RADIUS) * view; - let ground = camera + hit(camera, view, BOTTOM_RADIUS) * view; + let ground = camera + hit(camera, view, radius) * view; let r = TOP_RADIUS; let (mu, mu_s, nu) = (x.dot(&view) / r, x.dot(&sun) / r, view.dot(&sun)); let d = (ground - x).norm(); - let r_d = BOTTOM_RADIUS; + let r_d = radius; let mu_d = (r * mu + d) / r_d; let mu_s_d = (r * mu_s + d * nu) / r_d; let (near, far) = ( @@ -109,16 +124,16 @@ fn ground_reference( let s = tables.scattering(r, mu, mu_s, nu, true); let s_ground = tables.scattering(r_d, mu_d, mu_s_d, nu, true); let mu_g = ground.dot(&sun) / ground.norm(); - let to_sun = tables.transmittance(BOTTOM_RADIUS, mu_g); - let sky = tables.sky_irradiance(BOTTOM_RADIUS, mu_g); + let to_sun = tables.transmittance(radius, mu_g); + let sky = tables.sky_irradiance(radius, mu_g); std::array::from_fn(|k| { let through = (near[k] / far[k]).min(1.0); let between = (s[k] - through * s_ground[k]).max(0.0); - irradiance[k] * (albedo[k] * (to_sun[k] * mu_g + sky[k]) * through + PI * between) + irradiance[k] * (albedo[k] * (to_sun[k] * mu_g.max(0.0) + sky[k]) * through + PI * between) }) } -fn assert_close(got: [f32; 4], want: [f64; 3], tolerance: f64, context: &str) { +pub(crate) fn assert_close(got: [f32; 4], want: [f64; 3], tolerance: f64, context: &str) { for k in 0..3 { let error = f64::from(got[k]) / want[k] - 1.0; assert!( diff --git a/crates/resolved/tests/gpu/clouds.rs b/crates/resolved/tests/gpu/clouds.rs new file mode 100644 index 0000000..d252cd5 --- /dev/null +++ b/crates/resolved/tests/gpu/clouds.rs @@ -0,0 +1,226 @@ +//! The cloud layer at lighting L1 (`clouds.wgsl`): against the atmosphere's tables in +//! f64, a clear map against none, where a cloud's shadow falls, and the light clouds +//! send onto a target, on the GPU and on the CPU (lavapipe). + +use crate::atmosphere::{assert_close, center, constants, renderer, surface_reference, view}; +use crate::common; +use crate::earthshine::{rendered, Case}; + +use half::f16; +use resolved::atmosphere::{AtmosphereTables, BOTTOM_RADIUS}; +use resolved::clouds::CLOUD_HEIGHT; +use resolved::gpu::AdapterChoice; +use resolved::nalgebra::{UnitQuaternion, Vector3}; +use resolved::radiometry::AU; +use resolved::{CloudMap, Pose, Renderer, Scene}; + +const ADAPTERS: [AdapterChoice; 2] = [AdapterChoice::Gpu, AdapterChoice::Cpu]; + +/// `renderer` (gray Earth, uniform Moon) with clouds from `map`, if given. +fn cloudy(adapter: &AdapterChoice, map: Option) -> Renderer { + let mut renderer = renderer(adapter); + if let Some(map) = map { + renderer.set_cloud_map(&map); + } + renderer +} + +/// Clouds of one value everywhere, reflectance `value` / 255. +fn uniform(value: u8) -> CloudMap { + CloudMap::from_values(2, 1, vec![value; 2]).unwrap() +} + +/// The reflectance the GPU reads for map value `value`: in f16. +fn stored(value: u8) -> f64 { + f64::from(f16::from_f32(f32::from(value) / 255.0)) +} + +#[test] +fn a_uniform_layer_matches_the_tables() { + // As the_ground_through_the_air_matches_the_tables: the ground at 0 N 0 E seen 60 + // degrees from its zenith from 20,000 km, the Sun 30 degrees from its zenith on the + // other side. The cloud of reflectance R where the line of sight crosses the layer + // reflects as a Lambertian surface of albedo R at its height, and passes on 1 - R of + // what is below it: the ground, lit through 1 - R of cloud, and the air between. + let _devices = common::device_lock(); + let tables = AtmosphereTables::new(&constants().0); + let ground = Vector3::new(BOTTOM_RADIUS, 0.0, 0.0); + let toward = Vector3::new(60f64.to_radians().cos(), 60f64.to_radians().sin(), 0.0); + let camera = ground + 2e7 * toward; + let sun = Vector3::new(30f64.to_radians().cos(), -(30f64.to_radians().sin()), 0.0); + let albedo = constants().2; + let top = BOTTOM_RADIUS + CLOUD_HEIGHT; + let reference = + |radius, albedo| surface_reference(&tables, camera, -toward, sun, radius, albedo); + for adapter in &ADAPTERS { + for value in [0, 64, 192, 255] { + let r = stored(value); + let cloud = reference(top, [r; 3]); + let air_above = reference(top, [0.0; 3]); + let below = reference(BOTTOM_RADIUS, albedo.map(|a| a * (1.0 - r))); + let want: [f64; 3] = + std::array::from_fn(|k| cloud[k] + (1.0 - r) * (below[k] - air_above[k])); + let mut renderer = cloudy(adapter, Some(uniform(value))); + let (scene, pose) = view(sun, camera, ground); + let frame = renderer.submit(&scene, &pose).unwrap().wait().unwrap(); + let name = format!("{} at {value}", renderer.adapter_info().name); + // Measured: within 3.3e-4 on both devices, with the tables in f16. + assert_close(center(&frame), want, 1e-3, &name); + } + } +} + +#[test] +fn a_clear_map_changes_nothing() { + // The ground seen down through clear sky, with and without the layer. Lavapipe + // agrees to f32 rounding (measured 2.4e-7). The RTX 3060 filters texture reads in + // steps of 1/256 between texels (earth_detail.rs), and where the two pipelines' + // coordinates for the atmosphere's tables differ in the last bit, a read can land + // on the next step: measured 6.7e-3 at most, in 0.3 % of pixels, and 1.8e-6 on + // average. + let _devices = common::device_lock(); + let ground = Vector3::new(BOTTOM_RADIUS, 0.0, 0.0); + let camera = ground * 3.0; + let sun = Vector3::new(1.0, 1.0, 0.3).normalize(); + for adapter in &ADAPTERS { + let mut plain = cloudy(adapter, None); + let mut clear = cloudy(adapter, Some(uniform(0))); + let name = plain.adapter_info().name.clone(); + let (scene, pose) = view(sun, camera, ground); + let a = plain + .submit(&scene, &pose) + .unwrap() + .wait() + .unwrap() + .radiance + .unwrap(); + let b = clear + .submit(&scene, &pose) + .unwrap() + .wait() + .unwrap() + .radiance + .unwrap(); + let errors: Vec = a + .iter() + .zip(&b) + .flat_map(|(p, q)| (0..3).map(move |k| (p[k], q[k]))) + .filter(|(p, _)| *p > 0.0) + .map(|(p, q)| f64::from((q / p - 1.0).abs())) + .collect(); + let worst = errors.iter().fold(0.0, |m: f64, &e| m.max(e)); + let mean = errors.iter().sum::() / errors.len() as f64; + assert!(errors.len() > 1_000_000, "{name}: {} values", errors.len()); + assert!( + worst < 1e-2 && mean < 1e-5, + "{name}: {worst} at most, {mean} on average" + ); + } +} + +#[test] +fn a_cloud_casts_its_shadow_down_sun() { + // One opaque cloud, 2 x 2 texels of an 8192 x 4096 map (9.8 km) centered on 0 N 0 E, + // with the Sun 10 degrees above the horizon to the east, seen from 1000 km straight + // above. Its shadow falls where the line toward the Sun from the ground crosses the + // cloud: 28.7 km to the west. The image's centroid of the cloud (brighter than the + // clear ground) and of its shadow (darker) must land where those points project. + let (w, h) = (8192u32, 4096u32); + let mut values = vec![0u8; (w * h) as usize]; + for y in [h / 2 - 1, h / 2] { + for x in [w / 2 - 1, w / 2] { + values[(y * w + x) as usize] = 255; + } + } + let map = CloudMap::from_values(w, h, values).unwrap(); + let elevation = 10f64.to_radians(); + let sun = Vector3::new(elevation.sin(), elevation.cos(), 0.0); + let top = Vector3::new(BOTTOM_RADIUS + CLOUD_HEIGHT, 0.0, 0.0); + // The ground point whose line toward the Sun crosses the layer at `top`. + let along = top.dot(&sun); + let t = along - (along * along - (top.norm_squared() - BOTTOM_RADIUS * BOTTOM_RADIUS)).sqrt(); + let shadow = top - t * sun; + assert!((shadow.y + 28.7e3).abs() < 0.1e3, "{shadow:?}"); + let camera = Vector3::new(BOTTOM_RADIUS + 1e6, 0.0, 0.0); + let pose = Pose::looking_at(camera, Vector3::zeros(), Vector3::z()).unwrap(); + let scene = Scene { + earth_attitude: UnitQuaternion::identity(), + sun: AU * sun, + moon: Vector3::new(0.0, 0.0, -3.844e8), + moon_attitude: UnitQuaternion::identity(), + targets: vec![], + velocities: Default::default(), + }; + let k = common::example_intrinsics(); + let project = |p: Vector3| { + k.project(&(pose.attitude.inverse() * (p - camera))) + .unwrap() + }; + let (cloud_px, shadow_px) = (project(top), project(shadow)); + let _devices = common::device_lock(); + for adapter in &ADAPTERS { + let mut renderer = cloudy(adapter, Some(map.clone())); + let name = renderer.adapter_info().name.clone(); + let radiance = renderer + .submit(&scene, &pose) + .unwrap() + .wait() + .unwrap() + .radiance + .unwrap(); + common::save("clouds-shadow", renderer.adapter_info(), &radiance); + // Red, which the air scatters least, against the clear ground at the image's + // top-left corner. + let clear = f64::from(radiance[0][0]); + let centroid = |inside: &dyn Fn(f64) -> bool| { + let (mut sum, mut n) = (Vector3::zeros(), 0.0); + for (i, p) in radiance.iter().enumerate() { + if inside(f64::from(p[0]) / clear) { + sum += Vector3::new((i % 1920) as f64 + 0.5, (i / 1920) as f64 + 0.5, 0.0); + n += 1.0; + } + } + assert!(n > 100.0, "{name}: {n} pixels"); + sum / n + }; + let bright = centroid(&|x| x > 2.0); + let dark = centroid(&|x| x < 0.5); + for (got, want, what) in [(bright, cloud_px, "cloud"), (dark, shadow_px, "shadow")] { + let miss = ((got.x - want.x).powi(2) + (got.y - want.y).powi(2)).sqrt(); + // Measured: 0.41 px at most for the cloud, 0.12 for the shadow, 261 px away. + assert!(miss < 1.0, "{name} {what}: at {got:?}, want {want:?}"); + } + } +} + +#[test] +fn clouds_brighten_earthshine() { + // A plate 700 km over 0 N 0 E facing straight down, the Sun overhead, under a + // uniform layer of reflectance 0.5 and with none. The light Earth sends onto the + // plate grows about as Earth seen straight down from there does: measured 2.61, + // 2.44 and 2.14 times in red, green and blue, against 2.61, 2.50 and 2.29 straight + // down. Less in blue, where more of the light from the slant views around the + // plate's nadir is the air's, which the clouds do not change. + let _devices = common::device_lock(); + let target = Vector3::new(BOTTOM_RADIUS + 7e5, 0.0, 0.0); + let case = Case::new("clouds", target, -Vector3::x()); + let down = Pose::looking_at(target, Vector3::zeros(), Vector3::z()).unwrap(); + for adapter in &ADAPTERS { + let mut ratios = vec![]; + for value in [0, 128] { + let mut renderer = cloudy(adapter, Some(uniform(value))); + let light = rendered(&mut renderer, &[case])[0]; + let (scene, _) = view(Vector3::x(), target, Vector3::zeros()); + let frame = renderer.submit(&scene, &down).unwrap().wait().unwrap(); + ratios.push((light, center(&frame))); + } + for k in 0..3 { + let shine = ratios[1].0[k] / ratios[0].0[k]; + let nadir = f64::from(ratios[1].1[k] / ratios[0].1[k]); + assert!( + shine > 2.0 && (shine / nadir - 1.0).abs() < 0.08, + "{adapter:?} {k}: {shine}, {nadir}" + ); + } + } +} diff --git a/crates/resolved/tests/gpu/earthshine.rs b/crates/resolved/tests/gpu/earthshine.rs index 96196f9..7ab028b 100644 --- a/crates/resolved/tests/gpu/earthshine.rs +++ b/crates/resolved/tests/gpu/earthshine.rs @@ -123,7 +123,7 @@ fn srgb_to_linear(level: u8) -> f64 { /// One plate, its normal and the bodies around it. #[derive(Clone, Copy, Debug)] -struct Case { +pub(crate) struct Case { name: &'static str, target: Vector3, normal: Vector3, @@ -137,7 +137,7 @@ struct Case { } impl Case { - fn new(name: &'static str, target: Vector3, normal: Vector3) -> Self { + pub(crate) fn new(name: &'static str, target: Vector3, normal: Vector3) -> Self { Self { name, target, @@ -353,7 +353,7 @@ fn quadrature(case: &Case, surfaces: &Surfaces, rings: usize, sectors: usize) -> /// The light from the bodies on each case's plate as rendered, per channel, in units of /// the Sun's irradiance there: the center pixel's radiance over a Sun-facing plate's, /// less the direct sunlight. -fn rendered(renderer: &mut Renderer, cases: &[Case]) -> Vec<[f64; 3]> { +pub(crate) fn rendered(renderer: &mut Renderer, cases: &[Case]) -> Vec<[f64; 3]> { let model = renderer.load_model(&plate()).unwrap(); let radiometry = Radiometry::new(&CameraConfig::example()); cases diff --git a/crates/resolved/tests/gpu/main.rs b/crates/resolved/tests/gpu/main.rs index 56db604..3ff6770 100644 --- a/crates/resolved/tests/gpu/main.rs +++ b/crates/resolved/tests/gpu/main.rs @@ -9,6 +9,7 @@ mod apparent; mod atmosphere; mod blooming; mod bodies; +mod clouds; mod convolve; mod distortion; mod earth_detail; diff --git a/docs/optics.md b/docs/optics.md index d9c11ad..19d232f 100644 --- a/docs/optics.md +++ b/docs/optics.md @@ -999,8 +999,8 @@ This gets the overall brightness and the phase (how much of the disk is lit) right. It misses some things that are visible in real images: - Earth's oceans, land and clouds differ in brightness by a factor of - several. From lighting L1, Earth's surface comes from an image (next - section); clouds are not modeled yet. + several. From lighting L1, Earth's surface comes from an image and its + clouds from another (next sections). - The full Moon looks almost uniformly bright out to its edge. A Lambertian sphere instead darkens toward the edge, because the light there arrives and leaves at grazing angles. Lunar soil scatters light back toward its @@ -1234,12 +1234,125 @@ Limitations: - **Sunlight reaching the target**, and the Sun and Moon seen past the limb, pass through the atmosphere as "Sunset through the atmosphere" describes; targets, stars and Earth's own surface seen past the limb do not. -- **Clouds** are not drawn: Blue Marble is a cloud-free mosaic. +- **Clouds** are a separate layer, the next section. - **The Moon's shadow** in a solar eclipse dims the air along the line of sight as it dims the ground point behind it. - **Light from Earth onto targets** (Earthshine) includes the atmosphere over the ground, but not the glow beyond Earth's edge. +### Clouds + +Code: `clouds.rs`, `shaders/clouds.wgsl`, `shaders/cloudy.wgsl`, +`shaders/earthshine.wgsl`. From lighting L1, with the atmosphere. + +Blue Marble's surface is a cloud-free mosaic, but at any moment about two +thirds of Earth is under cloud, and clouds are the brightest thing on it. +Without them Earth looks too dark and too clear, and Earthshine on a target is +too low. + +The renderer draws clouds as a thin layer 5 km above the ground. Real cloud +tops run from about 1 to 15 km; one height stands for them all. The layer's map +is the cloud layer of NASA's original Blue Marble (2002): Terra MODIS clouds +of 2001 July 29, with thermal infrared over the poles, 1 km per pixel, resized +to 16384 x 8192 like the Earth maps (`scripts/fetch-clouds.sh`). The same +clouds show whatever the date. + +Each pixel's value is read as the cloud's reflectance there, from 0 for clear +sky to 1. Clouds absorb almost no visible light, so what a cloud does not +reflect it passes on: 1 minus its reflectance. A thick cloud is bright and +hides the ground. A thin one is dim and lets the ground show through. How +reflectance goes with a cloud's optical thickness tau (how many times light +would be scattered crossing it straight) is the two-stream result of Lacis and +Hansen (1974): + + R = 0.13 tau / (1 + 0.13 tau) + +The 0.13 is sqrt(3) / 2 times (1 - g), with g = 0.85 the cloud droplets' strong +forward scattering. So tau = 5 reflects 39 %, tau = 20 reflects 72 %, and tau = +50 reflects 87 %. Over Earth's area the map's mean reflectance is 0.21, and 62 % +of the area has more than 0.05. The map stores reflectance, so its mip levels +average reflectance, which keeps a far view of Earth as bright as a near one. + +For each pixel that sees the ground, the shader: + +1. finds where the line of sight crosses the layer, and reads the cloud there. + In an oblique view that is not the cloud straight above the ground point, + so clouds show parallax against the ground; +2. finds where the sunlight reaching the ground point crossed the layer, and + dims the ground's sunlight and sky light by that cloud. That is the cloud's + shadow, offset from the cloud away from the Sun: 5 km times the cotangent of + the Sun's elevation, 29 km with the Sun 10 degrees up; +3. adds the cloud's own light: sunlight and sky light at 5 km, which pass + through less air than at the ground, times the cloud's reflectance; +4. adds what is below the cloud, the ground and the air between it and the + layer, times 1 minus the cloud's reflectance. Then it applies the air above + the layer: its glow added, and its transmittance on what is under it. + +So high clouds are less hazed than the ground, the way they look from orbit. +The ray's entry into the atmosphere is looked up once and shared by the cloud +and the ground below it. + +The Earthshine pass samples the same layer, so a target over cloud receives +more light. A plate 700 km up facing straight down, the Sun overhead, under +clouds of reflectance 0.5 everywhere, receives 2.6, 2.4 and 2.1 times its +cloud-free Earthshine in red, green and blue. + +How bright Earth is as a whole tests the reading of the map's values as +reflectance. With clouds, rendered at phase angles from 0 to 180 degrees and +integrated over them, Earth reflects 0.27 to 0.29 of the sunlight in green +(the Bond albedo), from four sides in July. Without clouds it reflects 0.13 to +0.16. Measured, Earth's Bond albedo is about 0.30 (CERES, over the whole solar +spectrum). The geometric albedo, the full disk's brightness seen from the +Sun's direction, is 0.19 to 0.21 in green here, against 0.43 for the real +Earth, because clouds and the ocean send back toward the Sun more light than a +Lambertian surface does. + +Tests (`tests/gpu/clouds.rs`), on the GPU and lavapipe: + +- A uniform layer over gray ground, for four reflectances from clear to opaque, + against the atmosphere's tables in f64 at the cloud's height and the + ground's (within 3.3e-4). +- A clear map against no map. Lavapipe agrees to 2.4e-7. The RTX 3060 differs + by up to 6.7e-3 in 0.3 % of pixels, 1.8e-6 on average: its texture filtering + steps by 1/256 between texels, and a coordinate that differs in its last bit + between the two pipelines can land on the next step. +- One opaque cloud with the Sun 10 degrees up, seen from 1000 km straight + above: the cloud and its shadow 28.7 km to the west land within 0.41 and + 0.12 px of where they project, 261 px apart. +- The Earthshine above, against the ratio for Earth seen straight down (within + 8 %; less in blue, where more of the light from the slant views is the + air's). + +Unit tests in `clouds.rs` check the reflectance formula and the mip averaging. +The command-line tests check that a missing map names `scripts/fetch-clouds.sh`, +and that `[appearance] clouds = false` renders without one. + +Cost: the map takes 360 MB of GPU memory. It is compiled into the bodies pass +only once set, so frames without clouds cost what they did. With Earth filling +1920 x 1080 at camera C1, the bodies pass on the RTX 3060 goes from 2.4 to 3.1 ms +over mostly clear Egypt, and from 2.5 to 3.7 ms over the cloudy Pacific, +7.2 ms a frame end to end. On lavapipe, the Pacific view goes from 49 to 64 ms. +Earthshine: unchanged (0.045 ms). + +Limitations: + +- **One height**, and no thickness. Clouds stand neither above the limb nor + above one another, and cast no shadows on each other. Volumetric clouds are + planned for lighting L2. +- **Gray, Lambertian, and not absorbing.** Real clouds reflect more back toward + the Sun and forward along the light, and show a glory and a bright rim at the + terminator that a Lambertian layer does not. +- **No light between the cloud and the ground**: the ground's light reflected + back down by the cloud's base is left out. It matters over snow and ice. +- **The sky light under a cloud** is dimmed by the cloud toward the Sun, not by + the clouds over the whole sky. +- **One day's clouds**, whatever the scene's date. Maps by date are planned. +- **The map** is 2.45 km per pixel. Small cumulus are blurred into a gray + veil. The composite also has straight seams where its pieces meet, at + high latitudes, which show as clouds that end on a straight line (east of + Greenland in the solar eclipse example). +- **Without the atmosphere** (`Appearance::earth_atmosphere` off), no clouds. + ### The Moon's surface Code: `moon.rs`, `shaders/bodies.wgsl`; data from `scripts/fetch-moon.py`. From @@ -1931,7 +2044,7 @@ Limitations: and its error is largest for strongly colored surfaces. - **Lambertian bodies.** Earth and the Moon reflect as Lambertian surfaces here, as they do in the camera's view, with the Moon's map scaled by its - phase curve as seen from the target. Earth has no clouds yet. + phase curve as seen from the target. Clouds are Lambertian too. ## Where the Sun and Moon appear @@ -2103,7 +2216,8 @@ is a run of the script with `--gaia 16`, at twice the size. past the reach of camera C2's grid - Non-Lambertian materials (BRDFs), glints, and the full Moon's edge as bright as its center -- Clouds on Earth +- Clouds as volumes (lighting L2): thickness, the shadows they cast on each + other, and clouds standing above the limb - Earth's surface finer than Blue Marble's 460 m: mosaics of tens of meters (Landsat, Sentinel-2) need tiles streamed in as a view needs them - The atmosphere seen from inside it (a camera below 100 km) diff --git a/readme.md b/readme.md index 62f40fa..cbde4ed 100644 --- a/readme.md +++ b/readme.md @@ -129,7 +129,10 @@ At lighting L1, Earth's surface comes from NASA's Blue Marble, which from `[appearance] earth_month` in a scene without one. `scripts/fetch-earth.sh --sharp 07` also makes July's sharp map (460 m per texel instead of 2.45 km; 330 MB to download and about 3 minutes to make 4.9 GB), which the CLI and the movie use -when it is there. The Moon's surface +when it is there. Clouds over Earth come from the Blue Marble's cloud map (Terra +MODIS, 2001 July 29), which `scripts/fetch-clouds.sh` downloads (420 MB) and +resizes into `data/clouds/blue-marble.png`; `[appearance] clouds = false` leaves +them out. The Moon's surface comes from the LRO Camera's normalized-reflectance mosaic, which `scripts/fetch-moon.py` downloads (1.2 GB) and assembles (with numpy and Pillow) into `data/moon/wac-hapke-643.png`. The Moon is turned by the scene's epoch, or diff --git a/scripts/fetch-clouds.sh b/scripts/fetch-clouds.sh new file mode 100755 index 0000000..2a30b16 --- /dev/null +++ b/scripts/fetch-clouds.sh @@ -0,0 +1,52 @@ +#!/bin/sh +# Downloads NASA's Blue Marble cloud map and resizes it for the renderer. +# +# scripts/fetch-clouds.sh +# +# Source: the cloud layer of NASA's Blue Marble (2002), from +# https://science.nasa.gov/earth/earth-observatory/the-blue-marble-true-color-global-imagery-at-1km-resolution/: +# Terra MODIS clouds of 2001 July 29 (day 210), with thermal infrared over the poles, +# as two 21600 x 21600 PNGs, the western and eastern hemispheres (1 km per pixel). +# Produced by Reto Stöckli, NASA Goddard Space Flight Center; NASA asks that users +# credit "NASA Goddard Space Flight Center". +# +# The two halves are joined and resized to 16384 x 8192 (2.45 km per pixel, the size +# of the Earth maps), one gray channel, into data/clouds/blue-marble.png. The +# renderer reads a pixel's value as the cloud's reflectance there (docs/optics.md, +# "Clouds"); resizing averages those values as they are. Uses ImageMagick (`magick`). +# +# The originals go to data/clouds/source/ and are checked by sha256. Both +# directories are gitignored. + +set -eu + +base='https://eoimages.gsfc.nasa.gov/images/imagerecords/57000/57747' +root="$(cd "$(dirname "$0")/.." && pwd)/data/clouds" +mkdir -p "$root/source" + +# file name, sha256 +sums=$(cat <<'EOF' +cloud.W.2001210.21600x21600.png d7a6b03bb54b8ba56f72af6d1d49f191a25c12e911c43739c900e893d10eb814 +cloud.E.2001210.21600x21600.png 11eb1aaf8e6056b9aebc0de4894dc5ad8a21cbe8738c6ddc92e88701eb7e1b99 +EOF +) + +echo "$sums" | while read -r name sum; do + original="$root/source/$name" + if [ ! -f "$original" ]; then + echo "fetching $name" + curl -fsSL "$base/$name" -o "$original.part" + mv "$original.part" "$original" + fi + echo "$sum $original" | sha256sum -c --quiet - +done + +map="$root/blue-marble.png" +if [ ! -f "$map" ]; then + echo "resizing" + magick "$root/source/cloud.W.2001210.21600x21600.png" \ + "$root/source/cloud.E.2001210.21600x21600.png" +append \ + -colorspace Gray -filter Lanczos -resize 16384x8192! \ + +repage -depth 8 "$map.part.png" + mv "$map.part.png" "$map" +fi