diff --git a/crates/resolved-cli/src/main.rs b/crates/resolved-cli/src/main.rs index d154769..94b3ab9 100644 --- a/crates/resolved-cli/src/main.rs +++ b/crates/resolved-cli/src/main.rs @@ -51,10 +51,10 @@ struct RenderArgs { /// rolls highlights off to white at saturation. --out is unchanged. #[arg(long)] pretty: bool, - /// Lighting fidelity tier: L0, or L1 (adds Earth's surface, shadows on targets, - /// light from Earth and the Moon, eclipses on them, the Sun's limb darkening, and - /// stars, which camera C1 shows; needs Blue Marble data, and with C1 the star - /// catalog). + /// Lighting fidelity tier: L0, or L1 (adds Earth's surface and atmosphere, shadows + /// on targets, light from Earth and the Moon, eclipses on them, the Sun's limb + /// darkening, and stars, which camera C1 shows; needs Blue Marble data, and with C1 + /// the star catalog). #[arg(long, default_value = "L0")] lighting: LightingTier, /// Camera fidelity tier: C0, C1 (adds pixels that collect light over their area, diff --git a/crates/resolved/src/atmosphere.rs b/crates/resolved/src/atmosphere.rs new file mode 100644 index 0000000..83325d6 --- /dev/null +++ b/crates/resolved/src/atmosphere.rs @@ -0,0 +1,1177 @@ +//! Earth's atmosphere: the sunlight it scatters toward a viewer, and how much light it +//! lets through, precomputed as tables that the shaders look up +//! (`shaders/atmosphere.wgsl`). +//! +//! The model is Bruneton's precomputed atmospheric scattering (Bruneton and Neyret +//! 2008), in the form of Bruneton's 2017 implementation +//! (, BSD license): +//! its texture parameterizations, and its Earth atmosphere of Rayleigh scattering by +//! air, Mie scattering and absorption by aerosols, and absorption by ozone. Multiple +//! scattering is approximated as Hillaire (2020) does, from a small table of the light +//! scattered more than once at each altitude and Sun angle, rather than computed order +//! by order: that keeps the setup to a fraction of a second on the CPU. +//! +//! The tables are computed for wavelengths every 20 nm and folded into the camera's +//! channels with sunlight, the optics' transmission and each channel's QE, as the rest +//! of the renderer folds spectra ([`crate::radiometry`]). Scattered light is linear in +//! sunlight at each wavelength, so folding it is exact. Transmittance is not: a channel +//! whose band spans very different transmittances (blue light along a grazing path) is +//! given their sunlight-weighted mean, which `docs/optics.md` discusses. +//! +//! Lengths are in meters here and in kilometers in the shaders, where f32 needs them +//! small. The atmosphere is spherical, of Earth's equatorial radius; the shaders use it +//! in Earth-fixed axes stretched along the polar axis, where Earth's ellipsoid is that +//! sphere ([`crate::bodies`]). + +use std::f64::consts::PI; +use std::sync::{Arc, Mutex}; + +use crate::bodies::EARTH_EQUATORIAL_RADIUS; +use crate::radiometry::{Radiometry, AU, SUN_RADIUS}; + +/// The ground's radius, in meters. +pub const BOTTOM_RADIUS: f64 = EARTH_EQUATORIAL_RADIUS; +/// Height of the top of the atmosphere above the ground, in meters. Above it the air is +/// too thin to matter: 3.7e-6 of the ground's density at Rayleigh's 8 km scale height. +pub const TOP_ALTITUDE: f64 = 100e3; +/// The top of the atmosphere's radius, in meters. +pub const TOP_RADIUS: f64 = BOTTOM_RADIUS + TOP_ALTITUDE; + +/// Rayleigh scattering coefficient of air at the ground at 1 um, per meter; it scales as +/// 1 / lambda^4 (Bruneton's value). +const RAYLEIGH: f64 = 1.24062e-6; +const RAYLEIGH_SCALE_HEIGHT: f64 = 8000.0; +/// Aerosols' optical depth straight up at 550 nm: the global annual mean of the +/// aerosol climatology of Kinne et al. (2013), about 0.12. Real values run from a few +/// hundredths over clean ocean to over 1 in dust and smoke. +const AEROSOL_DEPTH: f64 = 0.12; +/// How aerosols' extinction falls with wavelength, as (wavelength / 550 nm)^-ANGSTROM: +/// Angstrom's own typical value. Fine dust and smoke are higher, sea salt lower. +const ANGSTROM: f64 = 1.3; +const MIE_SCALE_HEIGHT: f64 = 1200.0; +/// Aerosols' extinction at the ground at `nm`, per meter. +fn mie_extinction(nm: f64) -> f64 { + AEROSOL_DEPTH / MIE_SCALE_HEIGHT * (nm / 550.0).powf(-ANGSTROM) +} + +/// The fraction of the aerosols' extinction that is scattering, not absorption. +const MIE_ALBEDO: f64 = 0.9; +/// The asymmetry of the aerosols' Cornette-Shanks phase function. +pub const MIE_G: f64 = 0.8; +/// Ozone's peak number density, per cubic meter: 300 Dobson units spread over a profile +/// whose integral is 15 km. +const OZONE_PEAK: f64 = 300.0 * 2.687e20 / 15000.0; +/// Ozone's absorption cross-section, in m^2, from 360 to 830 nm every 10 nm: the +/// Serdyuchenko et al. (2014) 233 K spectrum averaged over each 10 nm bin, as tabulated +/// in Bruneton's demo. Taken as the 360 nm value below that, and as 0 above 830 nm, +/// where ozone's optical depth is under 6e-4. +#[rustfmt::skip] // Eight to a row, as tabulated. +const OZONE_CROSS_SECTION: [f64; 48] = [ + 1.18e-27, 2.182e-28, 2.818e-28, 6.636e-28, 1.527e-27, 2.763e-27, 5.52e-27, 8.451e-27, + 1.582e-26, 2.316e-26, 3.669e-26, 4.924e-26, 7.752e-26, 9.016e-26, 1.48e-25, 1.602e-25, + 2.139e-25, 2.755e-25, 3.091e-25, 3.5e-25, 4.266e-25, 4.672e-25, 4.398e-25, 4.701e-25, + 5.019e-25, 4.305e-25, 3.74e-25, 3.215e-25, 2.662e-25, 2.238e-25, 1.852e-25, 1.473e-25, + 1.209e-25, 9.423e-26, 7.455e-26, 6.566e-26, 5.105e-26, 4.15e-26, 4.228e-26, 3.237e-26, + 2.451e-26, 2.801e-26, 2.534e-26, 1.624e-26, 1.465e-26, 2.078e-26, 1.383e-26, 7.105e-27, +]; +/// The ground's albedo in the light scattered more than once: the light the ground sends +/// back up into the air. A cloud-free Earth's mean, about Blue Marble's. +const GROUND_ALBEDO: f64 = 0.1; +/// The Sun's angular radius at 1 AU, in radians. +const SUN_ANGULAR_RADIUS: f64 = SUN_RADIUS / AU; +/// The lowest Sun the scattering table covers, as the cosine of its zenith angle at the +/// ground: 120 degrees, deep into twilight. +const MU_S_MIN: f64 = -0.5; +/// Wavelengths the tables are computed at, in nm: every LAMBDA_STEP over this range, less +/// any where no channel responds. The tables change smoothly with wavelength, and a +/// channel is tens of nm wide. +const LAMBDA_RANGE: (u32, u32) = (360, 1100); +const LAMBDA_STEP: usize = 20; + +/// Transmittance table size: (mu, r). +pub const TRANSMITTANCE_SIZE: [usize; 2] = [256, 64]; +/// Scattering table size, (nu, mu_s, mu, r), stored as a 3D texture of +/// (nu * mu_s) x mu x r. +pub const SCATTERING_SIZE: [usize; 4] = [8, 32, 128, 16]; +/// Sky irradiance table size: (mu_s, r). +pub const IRRADIANCE_SIZE: [usize; 2] = [64, 16]; +/// Multiple scattering table size, per side: (mu_s, r). +const MULTIPLE_SIZE: usize = 32; + +/// Steps in the integral along a ray for the transmittance table, weighted for +/// exponential densities. +const TRANSMITTANCE_STEPS: usize = 500; +/// Steps along a ray for the scattering table (Bruneton's), weighted for exponential +/// densities ([`exponential_weights`]). +const SCATTERING_STEPS: usize = 50; +/// Directions and steps per direction for the multiple scattering table (Hillaire's). +const MULTIPLE_DIRECTIONS: usize = 8; +const MULTIPLE_STEPS: usize = 20; + +/// The tables, folded into a camera's channels (up to 4; unused ones are 0). Values are +/// per unit of sunlight: scattered radiance per unit solar irradiance (1/sr), and +/// irradiance as a fraction of the Sun's. +#[derive(Debug)] +pub struct AtmosphereTables { + /// Transmittance from a point to the top of the atmosphere, [`TRANSMITTANCE_SIZE`]. + pub transmittance: Vec<[f32; 4]>, + /// Rayleigh single scattering without its phase function, plus all multiple + /// scattering divided by Rayleigh's phase function, [`SCATTERING_SIZE`]. + pub rayleigh: Vec<[f32; 4]>, + /// Mie single scattering without its phase function, [`SCATTERING_SIZE`]. + pub mie: Vec<[f32; 4]>, + /// Irradiance from the sky on a horizontal surface, as a fraction of the Sun's, + /// [`IRRADIANCE_SIZE`]. + pub irradiance: Vec<[f32; 4]>, +} + +impl AtmosphereTables { + /// The tables for `radiometry`'s channels. Computed once per set of channels in a + /// process and shared: 0.25 s on a 24-core desktop CPU, 1.1 s on 4 of its cores. + pub fn new(radiometry: &Radiometry) -> Arc { + static CACHE: Mutex, Arc)>> = Mutex::new(Vec::new()); + let spectrum = Spectral::new(radiometry); + let key: Vec = spectrum + .weights + .iter() + .flatten() + .chain(&spectrum.lambda) + .map(|x| x.to_bits()) + .collect(); + let mut cache = CACHE.lock().expect("not poisoned"); + if let Some((_, tables)) = cache.iter().find(|(k, _)| *k == key) { + return tables.clone(); + } + let tables = Arc::new(Self::compute(&spectrum)); + cache.push((key, tables.clone())); + tables + } + + fn compute(spectrum: &Spectral) -> Self { + let columns = Columns::new(); + let multiple = Multiple::new(spectrum, &columns); + let (rayleigh, mie) = scattering(spectrum, &columns, &multiple); + let transmittance = columns + .values + .iter() + .map(|c| spectrum.fold(|l| spectrum.transmittance(l, c))) + .collect(); + let mut tables = Self { + transmittance, + rayleigh, + mie, + irradiance: Vec::new(), + }; + tables.irradiance = sky_irradiance(&tables); + tables + } + + /// Per channel, the transmittance from radius `r` to the top of the atmosphere along + /// a ray at cosine `mu` from the zenith, which must not meet the ground. + pub fn transmittance(&self, r: f64, mu: f64) -> [f64; 4] { + let (u, v) = transmittance_uv(r, mu); + bilinear(&self.transmittance, TRANSMITTANCE_SIZE, u, v) + } + + /// Per channel, the light scattered toward a point at radius `r`, looking along a ray + /// at cosine `mu` from its zenith, with the Sun at cosine `mu_s` from the zenith and + /// `nu` from the ray, from the whole ray up to the ground or the top of the + /// atmosphere: radiance per unit solar irradiance, phase functions included. + pub fn scattering(&self, r: f64, mu: f64, mu_s: f64, nu: f64, ground: bool) -> [f64; 4] { + let rayleigh = lookup(&self.rayleigh, r, mu, mu_s, nu, ground); + let mie = lookup(&self.mie, r, mu, mu_s, nu, ground); + let (pr, pm) = (rayleigh_phase(nu), mie_phase(nu)); + std::array::from_fn(|k| rayleigh[k] * pr + mie[k] * pm) + } + + /// Per channel, the sky's irradiance on a horizontal surface at radius `r` with the + /// Sun at cosine `mu_s` from the zenith, as a fraction of the Sun's irradiance. + pub fn sky_irradiance(&self, r: f64, mu_s: f64) -> [f64; 4] { + let (u, v) = irradiance_uv(r, mu_s); + bilinear(&self.irradiance, IRRADIANCE_SIZE, u, v) + } +} + +/// The height of inertial point `p` above the atmosphere's sphere, in meters, measured +/// as the shaders see it: in Earth-fixed axes stretched along the polar axis, where +/// Earth's ellipsoid is a sphere of the equatorial radius. +pub fn altitude(scene: &crate::Scene, p: nalgebra::Vector3) -> f64 { + let fixed = scene.earth_attitude.inverse() * p; + let stretch = BOTTOM_RADIUS / crate::bodies::EARTH_POLAR_RADIUS; + nalgebra::Vector3::new(fixed.x, fixed.y, fixed.z * stretch).norm() - BOTTOM_RADIUS +} + +/// The constants `shaders/atmosphere.wgsl` expects: lengths in kilometers, and the +/// tables' sizes. +pub(crate) fn wgsl_constants() -> String { + let [tw, th] = TRANSMITTANCE_SIZE; + let [nu, mu_s, mu, r] = SCATTERING_SIZE; + let [iw, ih] = IRRADIANCE_SIZE; + format!( + "const ATMOSPHERE_BOTTOM: f32 = {:?};\n\ + const ATMOSPHERE_TOP: f32 = {:?};\n\ + const ATMOSPHERE_MU_S_MIN: f32 = {:?};\n\ + const ATMOSPHERE_MIE_G: f32 = {:?};\n\ + const ATMOSPHERE_SUN_ANGULAR_RADIUS: f32 = {:?};\n\ + const ATMOSPHERE_TRANSMITTANCE_SIZE = vec2({tw}.0, {th}.0);\n\ + const ATMOSPHERE_SCATTERING_SIZE = vec4({nu}.0, {mu_s}.0, {mu}.0, {r}.0);\n\ + const ATMOSPHERE_IRRADIANCE_SIZE = vec2({iw}.0, {ih}.0);\n", + BOTTOM_RADIUS / 1000.0, + TOP_RADIUS / 1000.0, + MU_S_MIN, + MIE_G, + SUN_ANGULAR_RADIUS, + ) +} + +/// The wavelengths the tables are computed at, the atmosphere's coefficients there, and +/// how each folds into the camera's channels. +struct Spectral { + /// Wavelengths, in nm. + lambda: Vec, + /// Rayleigh scattering at the ground, per meter. + rayleigh: Vec, + /// Aerosol extinction at the ground, per meter. + mie: Vec, + /// Ozone absorption at peak density, per meter. + ozone: Vec, + /// Per wavelength, per channel: the weight of that wavelength's value in the + /// channel's, so that a channel's value is the sum of weights times values. + weights: Vec<[f64; 4]>, +} + +impl Spectral { + fn new(radiometry: &Radiometry) -> Self { + let channels = radiometry.solar_flux.len().min(4); + let grid: Vec = (LAMBDA_RANGE.0..=LAMBDA_RANGE.1) + .step_by(LAMBDA_STEP) + .map(f64::from) + .collect(); + // Each wavelength's value stands for a tent from its neighbors to it and down + // again, so the sum is the integral of the linear interpolation; the end tents + // reach out flat to cover the channels' tails. + let tent = |i: usize, l: f64| { + let x = l * 1e9; + let (lo, hi) = (grid[i.saturating_sub(1)], grid[(i + 1).min(grid.len() - 1)]); + let c = grid[i]; + if (i == 0 && x <= c) || (i == grid.len() - 1 && x >= c) { + 1.0 + } else if x >= lo && x <= c && c > lo { + (x - lo) / (c - lo) + } else if x > c && x <= hi && hi > c { + (hi - x) / (hi - c) + } else { + 0.0 + } + }; + let mut lambda = Vec::new(); + let mut weights = Vec::new(); + for (i, &nm) in grid.iter().enumerate() { + let mut w = [0.0; 4]; + for (k, wk) in w.iter_mut().enumerate().take(channels) { + *wk = radiometry.channel_integral(k, |l| tent(i, l)) / radiometry.solar_flux[k]; + } + // Wavelengths where no channel gets a millionth of its light are left out. + if w.iter().any(|&x| x > 1e-6) { + lambda.push(nm); + weights.push(w); + } + } + let rayleigh = lambda + .iter() + .map(|l| RAYLEIGH * (l * 1e-3).powi(-4)) + .collect(); + let mie = lambda.iter().map(|&l| mie_extinction(l)).collect(); + // Ozone's cross-section averaged over each wavelength's tent. + let step = LAMBDA_STEP as f64; + let ozone = lambda + .iter() + .map(|&l| { + let parts = [(-0.5, 0.25), (0.0, 0.5), (0.5, 0.25)]; + OZONE_PEAK + * parts + .iter() + .map(|(at, w)| w * ozone_cross_section(l + at * step)) + .sum::() + }) + .collect(); + Self { + lambda, + rayleigh, + mie, + ozone, + weights, + } + } + + /// Per channel, the sum over wavelengths of `value(wavelength index)` with the + /// channels' weights. + fn fold(&self, value: impl Fn(usize) -> f64) -> [f32; 4] { + let mut out = [0.0f64; 4]; + for (l, w) in self.weights.iter().enumerate() { + let v = value(l); + for k in 0..4 { + out[k] += w[k] * v; + } + } + out.map(|x| x as f32) + } + + /// Optical depth at wavelength index `l` of the column densities `c` (Rayleigh, Mie, + /// ozone, each in meters at ground or peak density). + fn depth(&self, l: usize, c: &[f64; 3]) -> f64 { + self.rayleigh[l] * c[0] + self.mie[l] * c[1] + self.ozone[l] * c[2] + } + + fn transmittance(&self, l: usize, c: &[f64; 3]) -> f64 { + (-self.depth(l, c)).exp() + } +} + +/// Ozone's cross-section at `nm`, linear between the table's bins' centers. +fn ozone_cross_section(nm: f64) -> f64 { + if nm > 830.0 { + return 0.0; + } + let x = ((nm - 365.0) / 10.0).clamp(0.0, 47.0); + let i = (x.floor() as usize).min(46); + let f = x - i as f64; + OZONE_CROSS_SECTION[i] * (1.0 - f) + OZONE_CROSS_SECTION[i + 1] * f +} + +/// Relative densities at `altitude` meters: air, aerosols and ozone. +fn densities(altitude: f64) -> [f64; 3] { + let ozone = if altitude < 25e3 { + (altitude - 10e3) / 15e3 + } else { + (40e3 - altitude) / 15e3 + }; + [ + (-altitude / RAYLEIGH_SCALE_HEIGHT).exp(), + (-altitude / MIE_SCALE_HEIGHT).exp(), + ozone.clamp(0.0, 1.0), + ] +} + +fn safe_sqrt(a: f64) -> f64 { + a.max(0.0).sqrt() +} + +fn clamp_cosine(mu: f64) -> f64 { + mu.clamp(-1.0, 1.0) +} + +fn clamp_radius(r: f64) -> f64 { + r.clamp(BOTTOM_RADIUS, TOP_RADIUS) +} + +/// Distance from radius `r` along a ray at cosine `mu` from the zenith to the top of the +/// atmosphere. +fn distance_to_top(r: f64, mu: f64) -> f64 { + (-r * mu + safe_sqrt(r * r * (mu * mu - 1.0) + TOP_RADIUS * TOP_RADIUS)).max(0.0) +} + +/// The same to the ground, for a ray that meets it. +fn distance_to_bottom(r: f64, mu: f64) -> f64 { + (-r * mu - safe_sqrt(r * r * (mu * mu - 1.0) + BOTTOM_RADIUS * BOTTOM_RADIUS)).max(0.0) +} + +fn intersects_ground(r: f64, mu: f64) -> bool { + mu < 0.0 && r * r * (mu * mu - 1.0) + BOTTOM_RADIUS * BOTTOM_RADIUS >= 0.0 +} + +/// The distance to the top of the atmosphere along a horizontal ray at the ground. +fn horizon_span() -> f64 { + (TOP_RADIUS * TOP_RADIUS - BOTTOM_RADIUS * BOTTOM_RADIUS).sqrt() +} + +/// A value in [0, 1] as a texture coordinate whose texel centers hold the ends. +fn coord_from_unit(x: f64, size: usize) -> f64 { + 0.5 / size as f64 + x * (1.0 - 1.0 / size as f64) +} + +fn unit_from_coord(u: f64, size: usize) -> f64 { + (u - 0.5 / size as f64) / (1.0 - 1.0 / size as f64) +} + +fn transmittance_uv(r: f64, mu: f64) -> (f64, f64) { + let h = horizon_span(); + let rho = safe_sqrt(r * r - BOTTOM_RADIUS * BOTTOM_RADIUS); + let d = distance_to_top(r, mu); + let (d_min, d_max) = (TOP_RADIUS - r, rho + h); + ( + coord_from_unit((d - d_min) / (d_max - d_min), TRANSMITTANCE_SIZE[0]), + coord_from_unit(rho / h, TRANSMITTANCE_SIZE[1]), + ) +} + +fn transmittance_r_mu(u: f64, v: f64) -> (f64, f64) { + let h = horizon_span(); + let x_mu = unit_from_coord(u, TRANSMITTANCE_SIZE[0]); + let rho = h * unit_from_coord(v, TRANSMITTANCE_SIZE[1]); + let r = (rho * rho + BOTTOM_RADIUS * BOTTOM_RADIUS).sqrt(); + let (d_min, d_max) = (TOP_RADIUS - r, rho + h); + let d = d_min + x_mu * (d_max - d_min); + let mu = if d == 0.0 { + 1.0 + } else { + (h * h - rho * rho - d * d) / (2.0 * r * d) + }; + (r, clamp_cosine(mu)) +} + +fn scattering_uvwz(r: f64, mu: f64, mu_s: f64, nu: f64, ground: bool) -> [f64; 4] { + let [_, size_mu_s, size_mu, size_r] = SCATTERING_SIZE; + let h = horizon_span(); + let rho = safe_sqrt(r * r - BOTTOM_RADIUS * BOTTOM_RADIUS); + let u_r = coord_from_unit(rho / h, size_r); + let r_mu = r * mu; + let discriminant = r_mu * r_mu - r * r + BOTTOM_RADIUS * BOTTOM_RADIUS; + let u_mu = if ground { + let d = -r_mu - safe_sqrt(discriminant); + let (d_min, d_max) = (r - BOTTOM_RADIUS, rho); + let x = if d_max == d_min { + 0.0 + } else { + (d - d_min) / (d_max - d_min) + }; + 0.5 - 0.5 * coord_from_unit(x, size_mu / 2) + } else { + let d = -r_mu + safe_sqrt(discriminant + h * h); + let (d_min, d_max) = (TOP_RADIUS - r, rho + h); + 0.5 + 0.5 * coord_from_unit((d - d_min) / (d_max - d_min), size_mu / 2) + }; + let d = distance_to_top(BOTTOM_RADIUS, mu_s); + let (d_min, d_max) = (TOP_RADIUS - BOTTOM_RADIUS, h); + let a = (d - d_min) / (d_max - d_min); + let big_a = (distance_to_top(BOTTOM_RADIUS, MU_S_MIN) - d_min) / (d_max - d_min); + let u_mu_s = coord_from_unit((1.0 - a / big_a).max(0.0) / (1.0 + a), size_mu_s); + [(nu + 1.0) / 2.0, u_mu_s, u_mu, u_r] +} + +/// The parameters of a scattering table texel, from its integer coordinates in the 3D +/// texture: r, mu, mu_s, nu and whether the ray meets the ground. +fn scattering_texel(x: usize, y: usize, z: usize) -> (f64, f64, f64, f64, bool) { + let [size_nu, size_mu_s, size_mu, size_r] = SCATTERING_SIZE; + let frag = [x as f64 + 0.5, y as f64 + 0.5, z as f64 + 0.5]; + let uvwz = [ + (frag[0] / size_mu_s as f64).floor() / (size_nu - 1) as f64, + (frag[0] % size_mu_s as f64) / size_mu_s as f64, + frag[1] / size_mu as f64, + frag[2] / size_r as f64, + ]; + let h = horizon_span(); + let rho = h * unit_from_coord(uvwz[3], size_r); + let r = (rho * rho + BOTTOM_RADIUS * BOTTOM_RADIUS).sqrt(); + let (mu, ground) = if uvwz[2] < 0.5 { + let (d_min, d_max) = (r - BOTTOM_RADIUS, rho); + let d = d_min + (d_max - d_min) * unit_from_coord(1.0 - 2.0 * uvwz[2], size_mu / 2); + let mu = if d == 0.0 { + -1.0 + } else { + clamp_cosine(-(rho * rho + d * d) / (2.0 * r * d)) + }; + (mu, true) + } else { + let (d_min, d_max) = (TOP_RADIUS - r, rho + h); + let d = d_min + (d_max - d_min) * unit_from_coord(2.0 * uvwz[2] - 1.0, size_mu / 2); + let mu = if d == 0.0 { + 1.0 + } else { + clamp_cosine((h * h - rho * rho - d * d) / (2.0 * r * d)) + }; + (mu, false) + }; + let x_mu_s = unit_from_coord(uvwz[1], size_mu_s); + let (d_min, d_max) = (TOP_RADIUS - BOTTOM_RADIUS, h); + let big_a = (distance_to_top(BOTTOM_RADIUS, MU_S_MIN) - d_min) / (d_max - d_min); + let a = (big_a - x_mu_s * big_a) / (1.0 + x_mu_s * big_a); + let d = d_min + a.min(big_a) * (d_max - d_min); + let mu_s = if d == 0.0 { + 1.0 + } else { + clamp_cosine((h * h - d * d) / (2.0 * BOTTOM_RADIUS * d)) + }; + let nu = clamp_cosine(uvwz[0] * 2.0 - 1.0); + // nu's range depends on mu and mu_s. + let spread = ((1.0 - mu * mu) * (1.0 - mu_s * mu_s)).sqrt(); + ( + r, + mu, + mu_s, + nu.clamp(mu * mu_s - spread, mu * mu_s + spread), + ground, + ) +} + +fn irradiance_uv(r: f64, mu_s: f64) -> (f64, f64) { + ( + coord_from_unit(mu_s * 0.5 + 0.5, IRRADIANCE_SIZE[0]), + coord_from_unit( + (r - BOTTOM_RADIUS) / (TOP_RADIUS - BOTTOM_RADIUS), + IRRADIANCE_SIZE[1], + ), + ) +} + +/// Texel indices and the weight of the second, for texture coordinate `u` over `size` +/// texels, with linear filtering and clamping at the edges as the GPU does. +fn texels(u: f64, size: usize) -> (usize, usize, f64) { + let x = (u * size as f64 - 0.5).clamp(0.0, (size - 1) as f64); + let i = (x.floor() as usize).min(size - 1); + (i, (i + 1).min(size - 1), x - i as f64) +} + +/// The four texels around texture coordinates (u, v) in a table of `size` (width, +/// height), with their weights for linear filtering. +fn corners(u: f64, v: f64, size: [usize; 2]) -> [(usize, f64); 4] { + let (x0, x1, fx) = texels(u, size[0]); + let (y0, y1, fy) = texels(v, size[1]); + [ + (y0 * size[0] + x0, (1.0 - fx) * (1.0 - fy)), + (y0 * size[0] + x1, fx * (1.0 - fy)), + (y1 * size[0] + x0, (1.0 - fx) * fy), + (y1 * size[0] + x1, fx * fy), + ] +} + +/// A 2D table of `size` (width, height) read at (u, v) with linear filtering. +fn bilinear(table: &[[f32; 4]], size: [usize; 2], u: f64, v: f64) -> [f64; 4] { + let around = corners(u, v, size); + std::array::from_fn(|k| { + around + .iter() + .map(|&(i, w)| f64::from(table[i][k]) * w) + .sum() + }) +} + +/// A scattering table read as the shader does: trilinear in the 3D texture, and linear +/// between the two nu slices around `nu`. +fn lookup(table: &[[f32; 4]], r: f64, mu: f64, mu_s: f64, nu: f64, ground: bool) -> [f64; 4] { + let [size_nu, size_mu_s, size_mu, size_r] = SCATTERING_SIZE; + let uvwz = scattering_uvwz(r, mu, mu_s, nu, ground); + let x = uvwz[0] * (size_nu - 1) as f64; + let slice = x.floor().min((size_nu - 2) as f64); + let lerp = x - slice; + let width = size_nu * size_mu_s; + let read = |s: f64| { + let u = (s + uvwz[1]) / size_nu as f64; + let (x0, x1, fx) = texels(u, width); + let (y0, y1, fy) = texels(uvwz[2], size_mu); + let (z0, z1, fz) = texels(uvwz[3], size_r); + let at = |x: usize, y: usize, z: usize, k: usize| { + f64::from(table[(z * size_mu + y) * width + x][k]) + }; + std::array::from_fn::(|k| { + let plane = |z: usize| { + let row = |y: usize| at(x0, y, z, k) * (1.0 - fx) + at(x1, y, z, k) * fx; + row(y0) * (1.0 - fy) + row(y1) * fy + }; + plane(z0) * (1.0 - fz) + plane(z1) * fz + }) + }; + let (a, b) = (read(slice), read(slice + 1.0)); + std::array::from_fn(|k| a[k] * (1.0 - lerp) + b[k] * lerp) +} + +/// Rayleigh's phase function, per steradian. +pub fn rayleigh_phase(nu: f64) -> f64 { + 3.0 / (16.0 * PI) * (1.0 + nu * nu) +} + +/// The aerosols' Cornette-Shanks phase function with asymmetry [`MIE_G`], per steradian. +pub fn mie_phase(nu: f64) -> f64 { + let g = MIE_G; + let k = 3.0 / (8.0 * PI) * (1.0 - g * g) / (2.0 + g * g); + k * (1.0 + nu * nu) / (1.0 + g * g - 2.0 * g * nu).powf(1.5) +} + +/// Column densities from points in the atmosphere to its top, over the transmittance +/// table's texels: meters of air, aerosols and ozone at ground (or peak) density. They +/// do not depend on wavelength, so one table gives the transmittance at every one. +struct Columns { + values: Vec<[f64; 3]>, +} + +impl Columns { + fn new() -> Self { + let [width, height] = TRANSMITTANCE_SIZE; + let rows = parallel(height, |y| { + (0..width) + .map(|x| { + let u = (x as f64 + 0.5) / width as f64; + let v = (y as f64 + 0.5) / height as f64; + let (r, mu) = transmittance_r_mu(u, v); + let dx = distance_to_top(r, mu) / TRANSMITTANCE_STEPS as f64; + let rho: Vec<[f64; 3]> = (0..=TRANSMITTANCE_STEPS) + .map(|step| { + let d = step as f64 * dx; + densities((d * d + 2.0 * r * mu * d + r * r).sqrt() - BOTTOM_RADIUS) + }) + .collect(); + let mut column = [0.0; 3]; + // Air and aerosols fall exponentially; ozone's profile is straight + // lines, for the trapezoid rule. + for s in 0..2 { + let profile: Vec = rho.iter().map(|p| p[s]).collect(); + column[s] = exponential_weights(&profile, dx).iter().sum(); + } + column[2] = rho.windows(2).map(|w| 0.5 * (w[0][2] + w[1][2]) * dx).sum(); + column + }) + .collect::>() + }); + Self { + values: rows.concat(), + } + } + + /// The column from radius `r` to the top along cosine `mu`, linearly interpolated. + fn at(&self, r: f64, mu: f64) -> [f64; 3] { + let (u, v) = transmittance_uv(r, mu); + let around = corners(u, v, TRANSMITTANCE_SIZE); + std::array::from_fn(|s| around.iter().map(|&(i, w)| self.values[i][s] * w).sum()) + } + + /// The column between radius `r` and the point `d` along cosine `mu`, as Bruneton's + /// `GetTransmittance` takes it from two lookups. + fn between(&self, r: f64, mu: f64, d: f64, ground: bool) -> [f64; 3] { + let r_d = clamp_radius((d * d + 2.0 * r * mu * d + r * r).sqrt()); + let mu_d = clamp_cosine((r * mu + d) / r_d); + let (near, far) = if ground { + (self.at(r_d, -mu_d), self.at(r, -mu)) + } else { + (self.at(r, mu), self.at(r_d, mu_d)) + }; + std::array::from_fn(|s| (near[s] - far[s]).max(0.0)) + } + + /// The column toward the Sun from radius `r` with the Sun at cosine `mu_s` from the + /// zenith, and the fraction of the Sun's disk above the horizon there. + fn to_sun(&self, r: f64, mu_s: f64) -> ([f64; 3], f64) { + (self.at(r, mu_s), sun_fraction(r, mu_s)) + } +} + +/// The fraction of the Sun's disk above the horizon at radius `r`, with the Sun at +/// cosine `mu_s` from the zenith: a smoothstep across the disk's diameter (Bruneton's +/// approximation of a circular segment's area). +fn sun_fraction(r: f64, mu_s: f64) -> f64 { + let sin_h = BOTTOM_RADIUS / r; + let cos_h = -safe_sqrt(1.0 - sin_h * sin_h); + let edge = sin_h * SUN_ANGULAR_RADIUS; + smoothstep(-edge, edge, mu_s - cos_h) +} + +fn smoothstep(a: f64, b: f64, x: f64) -> f64 { + let t = ((x - a) / (b - a)).clamp(0.0, 1.0); + t * t * (3.0 - 2.0 * t) +} + +/// Hillaire's (2020) table of multiple scattering: per wavelength, at each radius and +/// Sun angle, the light that reaches a point after scattering more than once, per unit +/// of a scattering coefficient there and of sunlight. It assumes the light scattered +/// once arrives the same from every direction and is the same at points nearby, sums +/// second-order scattering over directions, and the higher orders as a geometric series +/// of the fraction of light each order passes on. +struct Multiple { + /// [(mu_s index + r index * MULTIPLE_SIZE) * wavelengths + wavelength]. + values: Vec, + wavelengths: usize, +} + +impl Multiple { + fn new(spectrum: &Spectral, columns: &Columns) -> Self { + let n = spectrum.lambda.len(); + let rows = parallel(MULTIPLE_SIZE, |j| { + let mut row = vec![0.0; MULTIPLE_SIZE * n]; + for i in 0..MULTIPLE_SIZE { + let mu_s = unit_from_coord((i as f64 + 0.5) / MULTIPLE_SIZE as f64, MULTIPLE_SIZE) + * 2.0 + - 1.0; + let r = BOTTOM_RADIUS + + unit_from_coord((j as f64 + 0.5) / MULTIPLE_SIZE as f64, MULTIPLE_SIZE) + * (TOP_RADIUS - BOTTOM_RADIUS); + let value = Self::texel(spectrum, columns, r, mu_s); + row[i * n..(i + 1) * n].copy_from_slice(&value); + } + row + }); + Self { + values: rows.concat(), + wavelengths: n, + } + } + + fn texel(spectrum: &Spectral, columns: &Columns, r: f64, mu_s: f64) -> Vec { + let n = spectrum.lambda.len(); + let sun = [(1.0 - mu_s * mu_s).max(0.0).sqrt(), 0.0, mu_s]; + let (mut second, mut passed) = (vec![0.0; n], vec![0.0; n]); + let directions = MULTIPLE_DIRECTIONS * MULTIPLE_DIRECTIONS; + for a in 0..MULTIPLE_DIRECTIONS { + for b in 0..MULTIPLE_DIRECTIONS { + // Uniform over the sphere. + let theta = 2.0 * PI * (a as f64 + 0.5) / MULTIPLE_DIRECTIONS as f64; + let cos_phi = 1.0 - 2.0 * (b as f64 + 0.5) / MULTIPLE_DIRECTIONS as f64; + let sin_phi = (1.0 - cos_phi * cos_phi).sqrt(); + let dir = [theta.cos() * sin_phi, theta.sin() * sin_phi, cos_phi]; + let mu = dir[2]; + let ground = intersects_ground(r, mu); + let length = if ground { + distance_to_bottom(r, mu) + } else { + distance_to_top(r, mu) + }; + let dt = length / MULTIPLE_STEPS as f64; + let mut depth = vec![0.0; n]; + for step in 0..MULTIPLE_STEPS { + let t = (step as f64 + 0.5) * dt; + let p = [t * dir[0], t * dir[1], r + t * dir[2]]; + let r_p = (p[0] * p[0] + p[1] * p[1] + p[2] * p[2]).sqrt(); + let mu_s_p = (p[0] * sun[0] + p[2] * sun[2]) / r_p; + let rho = densities(r_p - BOTTOM_RADIUS); + let (column, fraction) = columns.to_sun(clamp_radius(r_p), mu_s_p); + for l in 0..n { + let scattering = + spectrum.rayleigh[l] * rho[0] + spectrum.mie[l] * MIE_ALBEDO * rho[1]; + let extinction = spectrum.rayleigh[l] * rho[0] + + spectrum.mie[l] * rho[1] + + spectrum.ozone[l] * rho[2]; + let segment = extinction * dt; + let through = (-(depth[l] + 0.5 * segment)).exp(); + let sunlight = fraction * spectrum.transmittance(l, &column); + second[l] += through * scattering * sunlight / (4.0 * PI) * dt; + passed[l] += through * scattering * dt; + depth[l] += segment; + } + } + if ground { + let p = [length * dir[0], length * dir[1], r + length * dir[2]]; + let r_g = (p[0] * p[0] + p[1] * p[1] + p[2] * p[2]).sqrt(); + let mu_s_g = (p[0] * sun[0] + p[2] * sun[2]) / r_g; + let (column, fraction) = columns.to_sun(BOTTOM_RADIUS, mu_s_g); + for l in 0..n { + let sunlight = fraction * spectrum.transmittance(l, &column); + second[l] += + (-depth[l]).exp() * sunlight * mu_s_g.max(0.0) * GROUND_ALBEDO / PI; + } + } + } + } + // Integrated over the sphere with the isotropic phase function: the mean. + (0..n) + .map(|l| { + let second = second[l] / directions as f64; + let passed = passed[l] / directions as f64; + second / (1.0 - passed) + }) + .collect() + } + + /// Per wavelength, the table at radius `r` and Sun cosine `mu_s`, linearly + /// interpolated, added into `out` times `scale` per wavelength. + fn add(&self, r: f64, mu_s: f64, scale: &[f64], out: &mut [f64]) { + let u = coord_from_unit(mu_s * 0.5 + 0.5, MULTIPLE_SIZE); + let v = coord_from_unit( + (r - BOTTOM_RADIUS) / (TOP_RADIUS - BOTTOM_RADIUS), + MULTIPLE_SIZE, + ); + let (x0, x1, fx) = texels(u, MULTIPLE_SIZE); + let (y0, y1, fy) = texels(v, MULTIPLE_SIZE); + let n = self.wavelengths; + let base = |x: usize, y: usize| (y * MULTIPLE_SIZE + x) * n; + let corners = [ + (base(x0, y0), (1.0 - fx) * (1.0 - fy)), + (base(x1, y0), fx * (1.0 - fy)), + (base(x0, y1), (1.0 - fx) * fy), + (base(x1, y1), fx * fy), + ]; + for l in 0..n { + let value: f64 = corners.iter().map(|(b, w)| self.values[b + l] * w).sum(); + out[l] += value * scale[l]; + } + } +} + +/// The scattering tables, folded into the channels: Rayleigh (with multiple scattering +/// over Rayleigh's phase function) and Mie. +/// +/// A row of the 3D texture shares r and mu, so everything along its ray but the Sun's +/// angle is computed once per row: the samples' places, their weights in the integral +/// ([`exponential_weights`]) and the transmittance back to the ray's start. Per texel, +/// transmittance toward the Sun is read per wavelength from a table rather than +/// computed from the columns, which saves an exponential per sample and wavelength. +fn scattering( + spectrum: &Spectral, + columns: &Columns, + multiple: &Multiple, +) -> (Vec<[f32; 4]>, Vec<[f32; 4]>) { + let [size_nu, size_mu_s, size_mu, size_r] = SCATTERING_SIZE; + let width = size_nu * size_mu_s; + let n = spectrum.lambda.len(); + // Transmittance to the top over the transmittance table's texels, per wavelength. + let to_top: Vec = columns + .values + .iter() + .flat_map(|c| (0..n).map(|l| spectrum.transmittance(l, c))) + .collect(); + let samples = SCATTERING_STEPS + 1; + let rows = parallel(size_mu * size_r, |row| { + let (y, z) = (row % size_mu, row / size_mu); + let (r, mu, _, _, ground) = scattering_texel(0, y, z); + let length = if ground { + distance_to_bottom(r, mu) + } else { + distance_to_top(r, mu) + }; + let dx = length / SCATTERING_STEPS as f64; + // Along the ray: each sample's distance and radius, its weights for air and + // aerosols, the transmittance back to the start per wavelength, and that times + // the scattering coefficient and weight, for multiple scattering. + let mut along = Vec::with_capacity(samples); + let mut rho = Vec::with_capacity(samples); + let mut through = vec![0.0; samples * n]; + for step in 0..samples { + let d = step as f64 * dx; + let r_d = clamp_radius((d * d + 2.0 * r * mu * d + r * r).sqrt()); + rho.push(densities(r_d - BOTTOM_RADIUS)); + let view = columns.between(r, mu, d, ground); + for l in 0..n { + through[step * n + l] = spectrum.transmittance(l, &view); + } + along.push((d, r_d)); + } + let air = exponential_weights(&rho.iter().map(|p| p[0]).collect::>(), dx); + let aerosols = exponential_weights(&rho.iter().map(|p| p[1]).collect::>(), dx); + let mut scattered = vec![0.0; samples * n]; + for step in 0..samples { + for l in 0..n { + scattered[step * n + l] = through[step * n + l] + * (spectrum.rayleigh[l] * air[step] + + spectrum.mie[l] * MIE_ALBEDO * aerosols[step]); + } + } + let mut out = Vec::with_capacity(width); + let (mut rayleigh, mut mie, mut more) = (vec![0.0; n], vec![0.0; n], vec![0.0; n]); + let mut to_sun = vec![0.0; n]; + for x in 0..width { + let (_, _, mu_s, nu, _) = scattering_texel(x, y, z); + rayleigh.fill(0.0); + mie.fill(0.0); + more.fill(0.0); + for (step, &(d, r_d)) in along.iter().enumerate() { + let mu_s_d = clamp_cosine((r * mu_s + d * nu) / r_d); + let fraction = sun_fraction(r_d, mu_s_d); + let at = step * n; + if fraction > 0.0 { + let (u, v) = transmittance_uv(r_d, mu_s_d); + to_sun.fill(0.0); + for (c, w) in corners(u, v, TRANSMITTANCE_SIZE) { + let texel = &to_top[c * n..(c + 1) * n]; + for l in 0..n { + to_sun[l] += w * texel[l]; + } + } + for l in 0..n { + let lit = through[at + l] * fraction * to_sun[l]; + rayleigh[l] += air[step] * lit; + mie[l] += aerosols[step] * lit; + } + } + multiple.add(r_d, mu_s_d, &scattered[at..at + n], &mut more); + } + let phase = rayleigh_phase(nu); + let r_value = spectrum.fold(|l| rayleigh[l] * spectrum.rayleigh[l] + more[l] / phase); + let m_value = spectrum.fold(|l| mie[l] * spectrum.mie[l] * MIE_ALBEDO); + out.push((r_value, m_value)); + } + out + }); + rows.into_iter().flatten().unzip() +} + +/// Weights for the integral of g(x) rho(x) over samples `dx` apart, where `rho` holds +/// a density at the samples and g is smooth: rho is taken as exponential between +/// samples and g as linear, and each step's integral is exact for those. A density +/// with a scale height shorter than a step (aerosols, 1.2 km, near the ground) is then +/// integrated correctly; the trapezoid rule was 9 % off there in Mie scattering. +fn exponential_weights(rho: &[f64], dx: f64) -> Vec { + let mut weights = vec![0.0; rho.len()]; + for i in 0..rho.len().saturating_sub(1) { + let (a, b) = (rho[i], rho[i + 1]); + if a > 0.0 && b > 0.0 { + // rho = a e^(q s / dx): the integrals of rho (1 - s / dx) and rho s / dx. + let q = (b / a).ln(); + let (start, end) = if q.abs() < 1e-3 { + (0.5 + q / 6.0 + q * q / 24.0, 0.5 + q / 3.0 + q * q / 8.0) + } else { + let e = q.exp(); + ((e - 1.0 - q) / (q * q), (q * e - e + 1.0) / (q * q)) + }; + weights[i] += a * start * dx; + weights[i + 1] += a * end * dx; + } else { + // A density that reaches zero is not exponential: the trapezoid rule. + weights[i] += 0.5 * a * dx; + weights[i + 1] += 0.5 * b * dx; + } + } + weights +} + +/// The sky's irradiance on a horizontal surface, from the folded scattering tables. +fn sky_irradiance(tables: &AtmosphereTables) -> Vec<[f32; 4]> { + const STEPS: usize = 32; + let [width, height] = IRRADIANCE_SIZE; + let rows = parallel(height, |j| { + (0..width) + .map(|i| { + let mu_s = unit_from_coord((i as f64 + 0.5) / width as f64, width) * 2.0 - 1.0; + let r = BOTTOM_RADIUS + + unit_from_coord((j as f64 + 0.5) / height as f64, height) + * (TOP_RADIUS - BOTTOM_RADIUS); + let sun = [(1.0 - mu_s * mu_s).max(0.0).sqrt(), 0.0, mu_s]; + let (d_theta, d_phi) = (PI / (2.0 * STEPS as f64), PI / STEPS as f64); + let mut sum = [0.0; 4]; + for a in 0..STEPS { + let theta = (a as f64 + 0.5) * d_theta; + for b in 0..2 * STEPS { + let phi = (b as f64 + 0.5) * d_phi; + let dir = [ + phi.cos() * theta.sin(), + phi.sin() * theta.sin(), + theta.cos(), + ]; + let nu = dir[0] * sun[0] + dir[2] * sun[2]; + let light = tables.scattering(r, dir[2], mu_s, nu, false); + let weight = dir[2] * theta.sin() * d_theta * d_phi; + for k in 0..4 { + sum[k] += light[k] * weight; + } + } + } + sum.map(|x| x as f32) + }) + .collect::>() + }); + rows.concat() +} + +/// `row(i)` for i in 0..count, on every core, in order. +fn parallel(count: usize, row: impl Fn(usize) -> T + Sync) -> Vec { + use std::sync::atomic::{AtomicUsize, Ordering}; + let next = AtomicUsize::new(0); + let threads = std::thread::available_parallelism().map_or(1, |n| n.get()); + let mut done: Vec<(usize, T)> = std::thread::scope(|scope| { + let handles: Vec<_> = (0..threads) + .map(|_| { + scope.spawn(|| { + let mut mine = Vec::new(); + loop { + let i = next.fetch_add(1, Ordering::Relaxed); + if i >= count { + return mine; + } + mine.push((i, row(i))); + } + }) + }) + .collect(); + handles + .into_iter() + .flat_map(|h| h.join().expect("worker panicked")) + .collect() + }); + done.sort_by_key(|(i, _)| *i); + done.into_iter().map(|(_, value)| value).collect() +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::camera::CameraConfig; + use approx::assert_relative_eq; + + /// The column from radius `r` along `mu` to the top, by a fine direct integral. + fn column(r: f64, mu: f64) -> [f64; 3] { + let steps = 20_000; + let dx = distance_to_top(r, mu) / steps as f64; + let mut sum = [0.0; 3]; + for i in 0..steps { + let d = (i as f64 + 0.5) * dx; + let rho = densities((d * d + 2.0 * r * mu * d + r * r).sqrt() - BOTTOM_RADIUS); + for s in 0..3 { + sum[s] += rho[s] * dx; + } + } + sum + } + + #[test] + fn texture_mappings_invert() { + for (r, mu) in [ + (BOTTOM_RADIUS, 0.3), + (BOTTOM_RADIUS + 20e3, -0.1), + (TOP_RADIUS, 1.0), + ] { + let (u, v) = transmittance_uv(r, mu); + let (r2, mu2) = transmittance_r_mu(u, v); + assert_relative_eq!(r2, r, max_relative = 1e-12); + assert!((mu2 - mu).abs() < 1e-9, "{mu} {mu2}"); + } + // Each scattering texel maps to coordinates at its own center. + let [size_nu, size_mu_s, size_mu, size_r] = SCATTERING_SIZE; + for (x, y, z) in [(37, 10, 3), (200, 100, 15), (5, 64, 0), (130, 63, 8)] { + let (r, mu, mu_s, _, ground) = scattering_texel(x, y, z); + let uvwz = scattering_uvwz(r, mu, mu_s, 0.0, ground); + let center = |i: usize, size: usize| (i as f64 + 0.5) / size as f64; + assert!( + (uvwz[1] - center(x % size_mu_s, size_mu_s)).abs() < 1e-9, + "{x}" + ); + assert!((uvwz[2] - center(y, size_mu)).abs() < 1e-9, "{y}"); + assert!((uvwz[3] - center(z, size_r)).abs() < 1e-9, "{z}"); + assert_eq!(ground, y < size_mu / 2); + let _ = size_nu; + } + } + + #[test] + fn columns_match_direct_integrals() { + let columns = Columns::new(); + // The zenith from the ground: air's column is its scale height. + let zenith = columns.at(BOTTOM_RADIUS, 1.0); + assert_relative_eq!(zenith[0], RAYLEIGH_SCALE_HEIGHT, max_relative = 1e-4); + assert_relative_eq!(zenith[1], MIE_SCALE_HEIGHT, max_relative = 1e-3); + assert_relative_eq!(zenith[2], 15e3, max_relative = 1e-4); + // Rayleigh optical depth at 550 nm: about 0.1. + let depth = RAYLEIGH * 0.55f64.powi(-4) * zenith[0]; + assert_relative_eq!(depth, 0.1084, max_relative = 1e-3); + // Lower and grazing rays, between texels, against a fine direct integral. + for (r, mu) in [ + (BOTTOM_RADIUS, 0.05), + (BOTTOM_RADIUS + 3e3, 0.0), + (BOTTOM_RADIUS + 30e3, -0.05), + (BOTTOM_RADIUS + 60e3, 0.4), + ] { + let (got, want) = (columns.at(r, mu), column(r, mu)); + assert_relative_eq!(got[0], want[0], max_relative = 0.01); + // Aerosols vanish with height, so their columns from high up are tiny: within + // 1 % of the column, or of the column from the ground up. Near the ground + // and the horizon it is 0.34 % off, 1.4e-4 in transmittance. + let tolerance = 0.01 * want[1].max(MIE_SCALE_HEIGHT); + assert!((got[1] - want[1]).abs() < tolerance, "{got:?} {want:?}"); + } + } + + /// A camera with one narrow channel at `nm`. + fn narrow(nm: f64) -> Radiometry { + let mut camera = CameraConfig::example(); + camera.channels.truncate(1); + camera.channels[0].quantum_efficiency = crate::Spectrum::gaussian(0.5, nm * 1e-9, 2e-9); + Radiometry::new(&camera) + } + + #[test] + fn a_narrow_channel_folds_to_its_wavelength() { + let c = [8000.0, 1200.0, 15e3]; + let exact = |nm: f64| { + (-(RAYLEIGH * (nm * 1e-3).powi(-4) * c[0] + + mie_extinction(nm) * c[1] + + OZONE_PEAK * ozone_cross_section(nm) * c[2])) + .exp() + }; + // On a computed wavelength, and halfway between two, where the tables are + // interpolated: 1.1e-3 off for a channel this narrow. + for (nm, tolerance) in [(560.0, 3e-4), (550.0, 2e-3)] { + let spectrum = Spectral::new(&narrow(nm)); + let sum: f64 = spectrum.weights.iter().map(|w| w[0]).sum(); + assert_relative_eq!(sum, 1.0, max_relative = 1e-9); + let folded = spectrum.fold(|l| spectrum.transmittance(l, &c))[0]; + assert_relative_eq!(f64::from(folded), exact(nm), max_relative = tolerance); + } + } + + #[test] + fn single_scattering_matches_a_direct_integral() { + // Single scattering toward the ground from the top of the atmosphere at one + // texel, without multiple scattering, against a fine march with exact optical + // depths, for a narrow channel on a computed wavelength. + let radiometry = narrow(560.0); + let spectrum = Spectral::new(&radiometry); + let columns = Columns::new(); + let (x, y, z) = (4 * SCATTERING_SIZE[1] + 20, 40, SCATTERING_SIZE[3] - 1); + let (r, mu, mu_s, nu, ground) = scattering_texel(x, y, z); + assert!(ground); + let empty = Multiple { + values: vec![0.0; MULTIPLE_SIZE * MULTIPLE_SIZE * spectrum.lambda.len()], + wavelengths: spectrum.lambda.len(), + }; + let (rayleigh, mie) = scattering(&spectrum, &columns, &empty); + let width = SCATTERING_SIZE[0] * SCATTERING_SIZE[1]; + let index = (z * SCATTERING_SIZE[2] + y) * width + x; + + let l = spectrum.lambda.iter().position(|&l| l == 560.0).unwrap(); + let beta = |rho: [f64; 3]| { + spectrum.rayleigh[l] * rho[0] + spectrum.mie[l] * rho[1] + spectrum.ozone[l] * rho[2] + }; + let steps = 4000; + let length = distance_to_bottom(r, mu); + let dx = length / steps as f64; + let (mut depth, mut sum_r, mut sum_m) = (0.0, 0.0, 0.0); + for i in 0..steps { + let d = (i as f64 + 0.5) * dx; + let r_d = (d * d + 2.0 * r * mu * d + r * r).sqrt(); + let rho = densities(r_d - BOTTOM_RADIUS); + let segment = beta(rho) * dx; + let mu_s_d = (r * mu_s + d * nu) / r_d; + let sun = column(r_d, mu_s_d); + let to_sun = (-(spectrum.rayleigh[l] * sun[0] + + spectrum.mie[l] * sun[1] + + spectrum.ozone[l] * sun[2])) + .exp(); + let t = (-(depth + 0.5 * segment)).exp() * to_sun; + sum_r += t * rho[0] * dx; + sum_m += t * rho[1] * dx; + depth += segment; + } + // Measured 1.9e-3 and 8.5e-4: mostly the aerosols' columns, interpolated in + // the transmittance table, near the ground. The channel spans a few nm around + // 560. + assert_relative_eq!( + f64::from(rayleigh[index][0]), + sum_r * spectrum.rayleigh[l], + max_relative = 3e-3 + ); + assert_relative_eq!( + f64::from(mie[index][0]), + sum_m * spectrum.mie[l] * MIE_ALBEDO, + max_relative = 5e-3 + ); + } + + #[test] + fn multiple_scattering_brightens_the_sky_by_a_fraction() { + let tables = AtmosphereTables::new(&Radiometry::new(&CameraConfig::example())); + // Looking up from the ground with the Sun high, away from it: the zenith sky, + // where multiple scattering adds a fraction of single scattering. + let r = BOTTOM_RADIUS; + let (mu, mu_s) = (1.0, 0.5); + let nu = mu_s; + let total = tables.scattering(r, mu, mu_s, nu, false); + // Blue is brighter than red, by about the Rayleigh ratio (620/460)^4 = 3.3, + // less ozone and extinction along the way. + assert!(total[2] > 2.0 * total[0], "{total:?}"); + assert!(total.iter().take(3).all(|&x| x > 0.0 && x.is_finite())); + // Sky irradiance at the ground with the Sun overhead: about a tenth to a fifth + // of the Sun's in blue, less in red. + let sky = tables.sky_irradiance(r, 1.0); + assert!(sky[2] > 0.08 && sky[2] < 0.3, "{sky:?}"); + assert!(sky[0] < sky[2], "{sky:?}"); + // Transmittance down the zenith: green about exp(-(0.11 air + 0.03 ozone + 0.12 + // aerosols)) = 0.77. + let t = tables.transmittance(r, 1.0); + assert!(t[1] > 0.74 && t[1] < 0.8, "{t:?}"); + } +} diff --git a/crates/resolved/src/bodies.rs b/crates/resolved/src/bodies.rs index 08aea23..fef8756 100644 --- a/crates/resolved/src/bodies.rs +++ b/crates/resolved/src/bodies.rs @@ -14,16 +14,18 @@ pub const MOON_RADIUS: f64 = 1_737_400.0; /// Without their maps (lighting L0) both are Lambertian with one reflectance spectrum /// over the whole body. That gives phase-dependent shading, but no oceans, clouds or /// maria, and not the Moon's steep brightening toward full Moon. At lighting L1 both -/// are read from maps in the colors set here, and the Moon is scaled to its measured -/// phase curve ([`crate::moon::phase_factor`]). +/// are read from maps in the colors set here, the Moon is scaled to its measured phase +/// curve ([`crate::moon::phase_factor`]), and Earth is seen through its atmosphere +/// ([`crate::atmosphere`]). #[derive(Clone, Debug, PartialEq)] pub struct Appearance { /// Earth's reflectance, used when no Earth texture is set. pub earth_albedo: Spectrum, - /// Added to an Earth texture's reflectance: a stand-in for the light Earth's - /// atmosphere scatters toward the camera, until the atmosphere is modeled. - /// Not used with `earth_albedo`, which already includes the atmosphere's effect. - pub earth_haze: Spectrum, + /// Whether Earth is seen through its atmosphere at lighting L1: the light it + /// scatters, the light it takes out, and the glow above the limb. On by default; + /// off leaves the bare surface, for comparisons and for tests of the surface alone. + /// `earth_albedo` (below L1) already includes the atmosphere's effect. + pub earth_atmosphere: bool, /// The Moon's reflectance, used when the Moon has no map (below lighting L1). pub moon_albedo: Spectrum, /// The shape of the Moon's reflectance spectrum, 1 at 643 nm, which a Moon map's @@ -50,25 +52,13 @@ impl Default for Appearance { Self { earth_albedo: Spectrum::sampled(nm.map(|x| x * 1e-9).to_vec(), earth.to_vec()) .expect("fixed samples"), - earth_haze: rayleigh_haze(), + earth_atmosphere: true, moon_albedo: Spectrum::constant(0.12), moon_color: crate::moon::moon_spectrum(), } } } -/// Reflectance shaped like Rayleigh scattering, falling as 1/lambda^4: 0.07 at 450 nm, -/// 0.017 at 650 nm. Sampled every 10 nm from 300 to 1100 nm, zero outside. -/// -/// The size is a round figure of the right order: clear air's Rayleigh optical depth is -/// about 0.2 at 450 nm, and the light it scatters toward a camera above it amounts to a -/// reflectance of a few hundredths to a tenth, depending on the angles. -fn rayleigh_haze() -> Spectrum { - let nm: Vec = (30..=110).map(|i| i as f64 * 10.0).collect(); - let value = nm.iter().map(|l| 0.07 * (450.0 / l).powi(4)).collect(); - Spectrum::sampled(nm.iter().map(|l| l * 1e-9).collect(), value).expect("increasing samples") -} - #[cfg(test)] mod tests { use super::*; diff --git a/crates/resolved/src/fidelity.rs b/crates/resolved/src/fidelity.rs index 7be1243..35833d8 100644 --- a/crates/resolved/src/fidelity.rs +++ b/crates/resolved/src/fidelity.rs @@ -19,7 +19,9 @@ 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 Moon's surface from + /// Blue Marble) seen through Earth's atmosphere, which scatters sunlight toward the + /// camera, dims what lies behind it and glows over the limb + /// ([`crate::atmosphere`]; a camera inside it is refused); the Moon's surface from /// the LRO Camera's map (`RendererOptions::moon_texture`), turned by /// `Scene::moon_attitude`; shadows that targets cast on themselves and on each /// other, with the penumbra the Sun's disk gives them; sunlight reflected onto diff --git a/crates/resolved/src/gpu/atmosphere.rs b/crates/resolved/src/gpu/atmosphere.rs new file mode 100644 index 0000000..baafe6c --- /dev/null +++ b/crates/resolved/src/gpu/atmosphere.rs @@ -0,0 +1,162 @@ +//! Earth's atmosphere's tables on the GPU ([`crate::atmosphere`]), which the bodies pass +//! and the Earthshine pass read through `shaders/atmosphere.wgsl`. + +use half::f16; + +use super::GpuContext; +use crate::atmosphere::{AtmosphereTables, IRRADIANCE_SIZE, SCATTERING_SIZE, TRANSMITTANCE_SIZE}; + +/// The tables as f16 textures, and the sampler that reads them: linear, clamped to the +/// edges, as the tables' texture coordinates assume. +#[derive(Debug)] +pub struct AtmosphereTextures { + transmittance: wgpu::TextureView, + rayleigh: wgpu::TextureView, + mie: wgpu::TextureView, + irradiance: wgpu::TextureView, + sampler: wgpu::Sampler, +} + +impl AtmosphereTextures { + /// Uploads `tables`. + pub fn new(ctx: &GpuContext, tables: &AtmosphereTables) -> Self { + let [tw, th] = TRANSMITTANCE_SIZE; + let [nu, mu_s, mu, r] = SCATTERING_SIZE; + let [iw, ih] = IRRADIANCE_SIZE; + let scattering = [(nu * mu_s) as u32, mu as u32, r as u32]; + Self::from_texels( + ctx, + [ + (&tables.transmittance, [tw as u32, th as u32, 1]), + (&tables.rayleigh, scattering), + (&tables.mie, scattering), + (&tables.irradiance, [iw as u32, ih as u32, 1]), + ], + ) + } + + /// One-texel tables of zeros, which passes bind until the atmosphere is set; the + /// shaders do not read them then. + pub fn placeholder(ctx: &GpuContext) -> Self { + let zero = vec![[0.0f32; 4]]; + Self::from_texels(ctx, [(&zero, [1, 1, 1]); 4]) + } + + fn from_texels(ctx: &GpuContext, tables: [(&Vec<[f32; 4]>, [u32; 3]); 4]) -> Self { + let upload = |label: &str, (texels, [w, h, d]): (&Vec<[f32; 4]>, [u32; 3]), three: bool| { + let half: Vec<[f16; 4]> = texels.iter().map(|t| t.map(f16::from_f32)).collect(); + let size = wgpu::Extent3d { + width: w, + height: h, + depth_or_array_layers: d, + }; + let texture = ctx.device.create_texture(&wgpu::TextureDescriptor { + label: Some(label), + size, + mip_level_count: 1, + sample_count: 1, + dimension: if three { + wgpu::TextureDimension::D3 + } else { + wgpu::TextureDimension::D2 + }, + format: wgpu::TextureFormat::Rgba16Float, + usage: wgpu::TextureUsages::TEXTURE_BINDING | wgpu::TextureUsages::COPY_DST, + view_formats: &[], + }); + ctx.queue.write_texture( + texture.as_image_copy(), + bytemuck::cast_slice(&half), + wgpu::TexelCopyBufferLayout { + offset: 0, + bytes_per_row: Some(8 * w), + rows_per_image: Some(h), + }, + size, + ); + texture.create_view(&Default::default()) + }; + let [transmittance, rayleigh, mie, irradiance] = tables; + let sampler = ctx.device.create_sampler(&wgpu::SamplerDescriptor { + label: Some("atmosphere"), + mag_filter: wgpu::FilterMode::Linear, + min_filter: wgpu::FilterMode::Linear, + ..Default::default() + }); + Self { + transmittance: upload("atmosphere transmittance", transmittance, false), + rayleigh: upload("atmosphere rayleigh", rayleigh, true), + mie: upload("atmosphere mie", mie, true), + irradiance: upload("atmosphere irradiance", irradiance, false), + sampler, + } + } + + /// The bind group entries for the tables and sampler, at bindings `first` to + /// `first + 4` in the order the shaders declare them: transmittance, Rayleigh, + /// Mie, irradiance, sampler. + pub fn entries(&self, first: u32) -> [wgpu::BindGroupEntry<'_>; 5] { + let view = |binding, view| wgpu::BindGroupEntry { + binding, + resource: wgpu::BindingResource::TextureView(view), + }; + [ + view(first, &self.transmittance), + view(first + 1, &self.rayleigh), + view(first + 2, &self.mie), + view(first + 3, &self.irradiance), + wgpu::BindGroupEntry { + binding: first + 4, + resource: wgpu::BindingResource::Sampler(&self.sampler), + }, + ] + } + + /// The matching layout entries, visible to `stages`. + pub fn layout_entries( + first: u32, + stages: wgpu::ShaderStages, + ) -> [wgpu::BindGroupLayoutEntry; 5] { + let texture = |binding, view_dimension| wgpu::BindGroupLayoutEntry { + binding, + visibility: stages, + ty: wgpu::BindingType::Texture { + sample_type: wgpu::TextureSampleType::Float { filterable: true }, + view_dimension, + multisampled: false, + }, + count: None, + }; + [ + texture(first, wgpu::TextureViewDimension::D2), + texture(first + 1, wgpu::TextureViewDimension::D3), + texture(first + 2, wgpu::TextureViewDimension::D3), + texture(first + 3, wgpu::TextureViewDimension::D2), + wgpu::BindGroupLayoutEntry { + binding: first + 4, + visibility: stages, + ty: wgpu::BindingType::Sampler(wgpu::SamplerBindingType::Filtering), + count: None, + }, + ] + } +} + +/// The shader declarations of the tables and sampler at bindings `first` to `first + 4` +/// of group 0, the constants they need, and `atmosphere.wgsl`. +pub(crate) fn wgsl(first: u32) -> String { + format!( + "{}@group(0) @binding({}) var atmosphere_transmittance_table: texture_2d;\n\ + @group(0) @binding({}) var atmosphere_rayleigh_table: texture_3d;\n\ + @group(0) @binding({}) var atmosphere_mie_table: texture_3d;\n\ + @group(0) @binding({}) var atmosphere_irradiance_table: texture_2d;\n\ + @group(0) @binding({}) var atmosphere_sampler: sampler;\n{}", + crate::atmosphere::wgsl_constants(), + first, + first + 1, + first + 2, + first + 3, + first + 4, + include_str!("../shaders/atmosphere.wgsl"), + ) +} diff --git a/crates/resolved/src/gpu/bodies.rs b/crates/resolved/src/gpu/bodies.rs index f6accf1..3cbe7e0 100644 --- a/crates/resolved/src/gpu/bodies.rs +++ b/crates/resolved/src/gpu/bodies.rs @@ -8,9 +8,11 @@ use nalgebra::Vector3; +use super::atmosphere::{self, AtmosphereTextures}; use super::profile::PassTimestamps; use super::targets::{FrameTargets, DEPTH_FORMAT, RADIANCE_FORMAT}; use super::GpuContext; +use crate::atmosphere::AtmosphereTables; use crate::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS, MOON_RADIUS}; use crate::camera::Intrinsics; use crate::earth::EarthTexture; @@ -77,8 +79,9 @@ pub struct BodyRadiometry { /// Per channel, the albedo of each RGB basis spectrum (see [`rgb_spectrum`]): column /// j maps a texture's channel j (red, green, blue) to the camera's channels. pub earth_rgb: [[f32; 4]; 3], - /// Per channel, the albedo of `Appearance::earth_haze`, added to a texture's. - pub earth_haze: [f32; 4], + /// Whether Earth is seen through its atmosphere: lighting L1 with + /// `Appearance::earth_atmosphere`, once the tables are set. + pub atmosphere: bool, /// Whether Earth is drawn from a texture. pub earth_textured: bool, /// Whether Earth and the Moon shadow each other: lunar and solar eclipses. @@ -138,7 +141,7 @@ impl BodyRadiometry { earth_lambert_1au: lambert(radiometry.effective_albedo(&appearance.earth_albedo)), earth_irradiance_1au: lambert(vec![1.0; n]), earth_rgb: [basis(0), basis(1), basis(2)], - earth_haze: albedo(&appearance.earth_haze), + atmosphere: false, earth_textured: false, bodies_eclipsed: false, earth_albedo: albedo(&appearance.earth_albedo), @@ -167,7 +170,9 @@ pub struct BodiesUniforms { earth_irradiance: [f32; 4], // mat3x4 columns: texture red, green, blue to channel albedo. earth_rgb: [[f32; 4]; 3], - earth_haze: [f32; 4], + // xyz: unit vector toward the Sun in Earth-fixed axes stretched along the polar + // axis, the atmosphere's frame. + earth_sun_stretched: [f32; 4], moon_center: [f32; 4], moon_sun: [f32; 4], moon_lambert: [f32; 4], @@ -221,7 +226,7 @@ impl BodiesUniforms { near as f32, f32::from(u8::from(radiometry.earth_textured)), f32::from(u8::from(radiometry.bodies_eclipsed)), - 0.0, + f32::from(u8::from(radiometry.atmosphere)), ], earth_center: vec4(rel.earth.translation, EARTH_EQUATORIAL_RADIUS), earth_axes: columns(rel.earth.rotation.to_rotation_matrix()), @@ -235,7 +240,7 @@ impl BodiesUniforms { earth_lambert: scaled(radiometry.earth_lambert_1au, scene.sun.norm()), earth_irradiance: scaled(radiometry.earth_irradiance_1au, scene.sun.norm()), earth_rgb: radiometry.earth_rgb, - earth_haze: radiometry.earth_haze, + earth_sun_stretched: vec4(sun_stretched(scene).cast::(), 0.0), moon_center: vec4(rel.moon, MOON_RADIUS), moon_sun: vec4(sun_from(scene.moon), 0.0), moon_lambert: scaled(radiometry.moon_lambert_1au, moon_distance), @@ -257,6 +262,18 @@ impl BodiesUniforms { } } +/// The unit vector toward the Sun from Earth's center in Earth-fixed axes stretched +/// along the polar axis, where Earth is a sphere: the frame of `atmosphere.wgsl`. +pub(crate) fn sun_stretched(scene: &Scene) -> Vector3 { + let fixed = scene.earth_attitude.inverse() * scene.sun.normalize(); + Vector3::new( + fixed.x, + fixed.y, + fixed.z * EARTH_EQUATORIAL_RADIUS / EARTH_POLAR_RADIUS, + ) + .normalize() +} + /// The compiled bodies pipeline and its uniform buffer. #[derive(Debug)] pub struct BodiesPass { @@ -268,6 +285,7 @@ pub struct BodiesPass { earth_height: u32, moon: wgpu::TextureView, moon_height: u32, + atmosphere: AtmosphereTextures, bind_group: wgpu::BindGroup, } @@ -317,6 +335,7 @@ impl BodiesPass { source: wgpu::ShaderSource::Wgsl( (prelude + include_str!("../shaders/common.wgsl") + + &atmosphere::wgsl(5) + include_str!("../shaders/bodies.wgsl")) .into(), ), @@ -374,6 +393,7 @@ impl BodiesPass { let sun_disk = SunDisk::buffer(ctx, &LimbDarkening::uniform()); let earth_view = upload_earth(ctx, &earth); let moon_view = upload_moon(ctx, &moon); + let atmosphere = AtmosphereTextures::placeholder(ctx); let bind_group = bind( ctx, &pipeline, @@ -382,6 +402,7 @@ impl BodiesPass { &sampler, &earth_view, &moon_view, + &atmosphere, ); Self { pipeline, @@ -392,6 +413,7 @@ impl BodiesPass { earth_height: earth.height(), moon: moon_view, moon_height: moon.height(), + atmosphere, bind_group, } } @@ -410,6 +432,18 @@ impl BodiesPass { self.rebind(ctx); } + /// Uploads Earth's atmosphere's tables and binds them for later frames; the + /// uniforms' atmosphere flag turns them on. + pub fn set_atmosphere(&mut self, ctx: &GpuContext, tables: &AtmosphereTables) { + self.atmosphere = AtmosphereTextures::new(ctx, tables); + self.rebind(ctx); + } + + /// The atmosphere's tables (placeholders until set), for the Earthshine pass. + pub fn atmosphere(&self) -> &AtmosphereTextures { + &self.atmosphere + } + fn rebind(&mut self, ctx: &GpuContext) { self.bind_group = bind( ctx, @@ -419,6 +453,7 @@ impl BodiesPass { &self.sampler, &self.earth, &self.moon, + &self.atmosphere, ); } @@ -580,6 +615,7 @@ fn upload<'a>( } /// The bodies pass's bind group. +#[allow(clippy::too_many_arguments)] fn bind( ctx: &GpuContext, pipeline: &wgpu::RenderPipeline, @@ -588,31 +624,34 @@ fn bind( sampler: &wgpu::Sampler, earth: &wgpu::TextureView, moon: &wgpu::TextureView, + atmosphere: &AtmosphereTextures, ) -> wgpu::BindGroup { + let entries = [ + wgpu::BindGroupEntry { + binding: 0, + resource: uniforms.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 1, + resource: wgpu::BindingResource::TextureView(earth), + }, + wgpu::BindGroupEntry { + binding: 2, + resource: wgpu::BindingResource::Sampler(sampler), + }, + wgpu::BindGroupEntry { + binding: 3, + resource: sun_disk.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 4, + resource: wgpu::BindingResource::TextureView(moon), + }, + ]; + let entries: Vec<_> = entries.into_iter().chain(atmosphere.entries(5)).collect(); ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { label: Some("bodies"), layout: &pipeline.get_bind_group_layout(0), - entries: &[ - wgpu::BindGroupEntry { - binding: 0, - resource: uniforms.as_entire_binding(), - }, - wgpu::BindGroupEntry { - binding: 1, - resource: wgpu::BindingResource::TextureView(earth), - }, - wgpu::BindGroupEntry { - binding: 2, - resource: wgpu::BindingResource::Sampler(sampler), - }, - wgpu::BindGroupEntry { - binding: 3, - resource: sun_disk.as_entire_binding(), - }, - wgpu::BindGroupEntry { - binding: 4, - resource: wgpu::BindingResource::TextureView(moon), - }, - ], + entries: &entries, }) } diff --git a/crates/resolved/src/gpu/earthshine.rs b/crates/resolved/src/gpu/earthshine.rs index 6f31d61..c33f2f9 100644 --- a/crates/resolved/src/gpu/earthshine.rs +++ b/crates/resolved/src/gpu/earthshine.rs @@ -7,6 +7,7 @@ use nalgebra::{UnitQuaternion, Vector3}; +use super::atmosphere::{self, AtmosphereTextures}; use super::bodies::{BodiesPass, BodyRadiometry, SunDisk}; use super::profile::PassTimestamps; use super::GpuContext; @@ -104,7 +105,6 @@ impl TargetLight { struct Globals { earth_albedo: [f32; 4], earth_rgb: [[f32; 4]; 3], - earth_haze: [f32; 4], moon_albedo: [f32; 4], params: [f32; 4], grid: [u32; 4], @@ -141,48 +141,56 @@ impl EarthshinePass { has_dynamic_offset: false, min_binding_size: None, }; + let entries = [ + entry(0, buffer(wgpu::BufferBindingType::Uniform)), + entry( + 1, + buffer(wgpu::BufferBindingType::Storage { read_only: true }), + ), + entry( + 2, + buffer(wgpu::BufferBindingType::Storage { read_only: false }), + ), + entry( + 3, + wgpu::BindingType::Texture { + sample_type: wgpu::TextureSampleType::Float { filterable: true }, + view_dimension: wgpu::TextureViewDimension::D2, + multisampled: false, + }, + ), + entry( + 4, + wgpu::BindingType::Sampler(wgpu::SamplerBindingType::Filtering), + ), + entry( + 5, + wgpu::BindingType::StorageTexture { + access: wgpu::StorageTextureAccess::WriteOnly, + format: CUBE_FORMAT, + view_dimension: wgpu::TextureViewDimension::D2Array, + }, + ), + entry(6, buffer(wgpu::BufferBindingType::Uniform)), + entry( + 7, + wgpu::BindingType::Texture { + sample_type: wgpu::TextureSampleType::Float { filterable: true }, + view_dimension: wgpu::TextureViewDimension::D2, + multisampled: false, + }, + ), + ]; + let entries: Vec<_> = entries + .into_iter() + .chain(AtmosphereTextures::layout_entries( + 8, + wgpu::ShaderStages::COMPUTE, + )) + .collect(); let layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { label: Some("earthshine"), - entries: &[ - entry(0, buffer(wgpu::BufferBindingType::Uniform)), - entry( - 1, - buffer(wgpu::BufferBindingType::Storage { read_only: true }), - ), - entry( - 2, - buffer(wgpu::BufferBindingType::Storage { read_only: false }), - ), - entry( - 3, - wgpu::BindingType::Texture { - sample_type: wgpu::TextureSampleType::Float { filterable: true }, - view_dimension: wgpu::TextureViewDimension::D2, - multisampled: false, - }, - ), - entry( - 4, - wgpu::BindingType::Sampler(wgpu::SamplerBindingType::Filtering), - ), - entry( - 5, - wgpu::BindingType::StorageTexture { - access: wgpu::StorageTextureAccess::WriteOnly, - format: CUBE_FORMAT, - view_dimension: wgpu::TextureViewDimension::D2Array, - }, - ), - entry(6, buffer(wgpu::BufferBindingType::Uniform)), - entry( - 7, - wgpu::BindingType::Texture { - sample_type: wgpu::TextureSampleType::Float { filterable: true }, - view_dimension: wgpu::TextureViewDimension::D2, - multisampled: false, - }, - ), - ], + entries: &entries, }); let pipeline_layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor { label: Some("earthshine"), @@ -192,10 +200,9 @@ impl EarthshinePass { let module = device.create_shader_module(wgpu::ShaderModuleDescriptor { label: Some("earthshine"), source: wgpu::ShaderSource::Wgsl( - concat!( - include_str!("../shaders/common.wgsl"), - include_str!("../shaders/earthshine.wgsl") - ) + (include_str!("../shaders/common.wgsl").to_owned() + + &atmosphere::wgsl(8) + + include_str!("../shaders/earthshine.wgsl")) .into(), ), }); @@ -316,13 +323,12 @@ impl EarthshinePass { let globals = Globals { earth_albedo: radiometry.earth_albedo, earth_rgb: radiometry.earth_rgb, - earth_haze: radiometry.earth_haze, moon_albedo: radiometry.moon_albedo, params: [ f32::from(u8::from(radiometry.earth_textured)), (std::f64::consts::PI * EARTH_EQUATORIAL_RADIUS / height as f64) as f32, (EARTH_EQUATORIAL_RADIUS / EARTH_POLAR_RADIUS) as f32, - 0.0, + f32::from(u8::from(radiometry.atmosphere)), ], grid: [RINGS, SECTORS, CUBE_SIZE, lights.len() as u32], moon_color: radiometry.moon_color, @@ -341,43 +347,48 @@ impl EarthshinePass { dimension: Some(wgpu::TextureViewDimension::D2Array), ..Default::default() }); + let entries = [ + wgpu::BindGroupEntry { + binding: 0, + resource: self.globals.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 1, + resource: self.lights.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 2, + resource: self.samples.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 3, + resource: wgpu::BindingResource::TextureView(texture), + }, + wgpu::BindGroupEntry { + binding: 4, + resource: wgpu::BindingResource::Sampler(sampler), + }, + wgpu::BindGroupEntry { + binding: 5, + resource: wgpu::BindingResource::TextureView(&cube), + }, + wgpu::BindGroupEntry { + binding: 6, + resource: self.sun_disk.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 7, + resource: wgpu::BindingResource::TextureView(moon.0), + }, + ]; + let entries: Vec<_> = entries + .into_iter() + .chain(bodies.atmosphere().entries(8)) + .collect(); let group = ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { label: Some("earthshine"), layout: &self.layout, - entries: &[ - wgpu::BindGroupEntry { - binding: 0, - resource: self.globals.as_entire_binding(), - }, - wgpu::BindGroupEntry { - binding: 1, - resource: self.lights.as_entire_binding(), - }, - wgpu::BindGroupEntry { - binding: 2, - resource: self.samples.as_entire_binding(), - }, - wgpu::BindGroupEntry { - binding: 3, - resource: wgpu::BindingResource::TextureView(texture), - }, - wgpu::BindGroupEntry { - binding: 4, - resource: wgpu::BindingResource::Sampler(sampler), - }, - wgpu::BindGroupEntry { - binding: 5, - resource: wgpu::BindingResource::TextureView(&cube), - }, - wgpu::BindGroupEntry { - binding: 6, - resource: self.sun_disk.as_entire_binding(), - }, - wgpu::BindGroupEntry { - binding: 7, - resource: wgpu::BindingResource::TextureView(moon.0), - }, - ], + entries: &entries, }); let n = lights.len() as u32; pass.set_bind_group(0, &group, &[]); diff --git a/crates/resolved/src/gpu/mod.rs b/crates/resolved/src/gpu/mod.rs index 044d5f3..d7ee200 100644 --- a/crates/resolved/src/gpu/mod.rs +++ b/crates/resolved/src/gpu/mod.rs @@ -1,5 +1,6 @@ //! GPU plumbing: adapter selection, device setup, render targets and passes. +pub mod atmosphere; pub mod blur; pub mod bodies; pub mod context; diff --git a/crates/resolved/src/lib.rs b/crates/resolved/src/lib.rs index 613d55b..93592ea 100644 --- a/crates/resolved/src/lib.rs +++ b/crates/resolved/src/lib.rs @@ -17,6 +17,7 @@ pub use nalgebra; /// [`gpu::GpuContext`] hands out wgpu types. pub use wgpu; +pub mod atmosphere; pub mod bodies; pub mod camera; pub mod earth; diff --git a/crates/resolved/src/render.rs b/crates/resolved/src/render.rs index 30e3fb6..7dba0ba 100644 --- a/crates/resolved/src/render.rs +++ b/crates/resolved/src/render.rs @@ -15,6 +15,7 @@ use crate::bodies::{Appearance, MOON_RADIUS}; use crate::camera::{Blooming, CameraConfig, Intrinsics, Psf}; use std::sync::Arc; +use crate::atmosphere::AtmosphereTables; use crate::earth::EarthTexture; use crate::eclipse::{sun_visible_fraction, sun_visible_light}; use crate::fidelity::{CameraTier, Fidelity, LightingTier}; @@ -431,6 +432,11 @@ impl Renderer { renderer .earthshine .set_limb_darkening(&renderer.ctx, &renderer.limb); + if options.appearance.earth_atmosphere { + let tables = AtmosphereTables::new(&renderer.radiometry); + renderer.bodies.set_atmosphere(&renderer.ctx, &tables); + renderer.body_radiometry.atmosphere = true; + } } if let Some(texture) = &options.earth_texture { renderer.set_earth_texture(texture); @@ -541,6 +547,18 @@ impl Renderer { ))); } let camera = observer.then(&self.camera.mount); + if self.body_radiometry.atmosphere { + let altitude = crate::atmosphere::altitude(scene, camera.position); + if altitude < crate::atmosphere::TOP_ALTITUDE { + return Err(Error::Unsupported(format!( + "the camera is {:.0} km up, inside Earth's atmosphere, which lighting {} \ + draws only from above its top at {:.0} km", + altitude / 1e3, + self.fidelity.lighting, + crate::atmosphere::TOP_ALTITUDE / 1e3 + ))); + } + } let intrinsics = self.render_intrinsics; let near = self.camera.near_clip; let occulters = [ diff --git a/crates/resolved/src/shaders/atmosphere.wgsl b/crates/resolved/src/shaders/atmosphere.wgsl new file mode 100644 index 0000000..bf33c21 --- /dev/null +++ b/crates/resolved/src/shaders/atmosphere.wgsl @@ -0,0 +1,235 @@ +// Earth's atmosphere, from the tables atmosphere.rs computes: the sunlight it scatters +// toward a viewer in space, the light it lets through, and the skylight on the ground. +// Bruneton's lookups (functions.glsl of his 2017 implementation), in kilometers. +// +// Positions are relative to Earth's center in Earth-fixed axes stretched along the +// polar axis, where Earth is a sphere of radius ATMOSPHERE_BOTTOM (common.wgsl, +// `ellipsoid_hit`). Scattered light is radiance per unit of solar irradiance, per +// channel; transmittance is a fraction. +// +// The code that includes this file first declares the ATMOSPHERE_ constants +// (atmosphere.rs, `wgsl_constants`), the textures `atmosphere_transmittance_table` and +// `atmosphere_irradiance_table` (texture_2d), `atmosphere_rayleigh_table` and +// `atmosphere_mie_table` (texture_3d), and the sampler `atmosphere_sampler`, +// linear and clamped to the edges (gpu/atmosphere.rs, `wgsl`). + +fn atmosphere_coord(x: f32, size: f32) -> f32 { + return 0.5 / size + x * (1.0 - 1.0 / size); +} + +// The distance to the top of the atmosphere for a horizontal ray at the ground. +fn atmosphere_span() -> f32 { + return sqrt(ATMOSPHERE_TOP * ATMOSPHERE_TOP - ATMOSPHERE_BOTTOM * ATMOSPHERE_BOTTOM); +} + +fn atmosphere_to_top(r: f32, mu: f32) -> f32 { + let discriminant = r * r * (mu * mu - 1.0) + ATMOSPHERE_TOP * ATMOSPHERE_TOP; + return max(-r * mu + sqrt(max(discriminant, 0.0)), 0.0); +} + +fn atmosphere_transmittance_uv(r: f32, mu: f32) -> vec2 { + let h = atmosphere_span(); + let rho = sqrt(max(r * r - ATMOSPHERE_BOTTOM * ATMOSPHERE_BOTTOM, 0.0)); + let d = atmosphere_to_top(r, mu); + let d_min = ATMOSPHERE_TOP - r; + let d_max = rho + h; + return vec2( + atmosphere_coord((d - d_min) / (d_max - d_min), ATMOSPHERE_TRANSMITTANCE_SIZE.x), + atmosphere_coord(rho / h, ATMOSPHERE_TRANSMITTANCE_SIZE.y), + ); +} + +// Per channel, the transmittance from radius `r` to the top along cosine `mu`, for a +// ray that does not meet the ground. +fn atmosphere_transmittance_to_top(r: f32, mu: f32) -> vec4 { + let uv = atmosphere_transmittance_uv(r, mu); + return textureSampleLevel(atmosphere_transmittance_table, atmosphere_sampler, uv, 0.0); +} + +// Per channel, the transmittance between radius `r` and the point `d` along cosine +// `mu`, from two lookups (`ground`: the ray meets the ground). +fn atmosphere_transmittance(r: f32, mu: f32, d: f32, ground: bool) -> vec4 { + let r_d = clamp(sqrt(d * d + 2.0 * r * mu * d + r * r), ATMOSPHERE_BOTTOM, ATMOSPHERE_TOP); + let mu_d = clamp((r * mu + d) / r_d, -1.0, 1.0); + if ground { + return min( + atmosphere_transmittance_to_top(r_d, -mu_d) / atmosphere_transmittance_to_top(r, -mu), + vec4(1.0), + ); + } + return min( + atmosphere_transmittance_to_top(r, mu) / atmosphere_transmittance_to_top(r_d, mu_d), + vec4(1.0), + ); +} + +// Per channel, the transmittance toward the Sun from radius `r`, with the Sun at cosine +// `mu_s` from the zenith, times the fraction of its disk above the horizon. +fn atmosphere_transmittance_to_sun(r: f32, mu_s: f32) -> vec4 { + let sin_h = ATMOSPHERE_BOTTOM / r; + let cos_h = -sqrt(max(1.0 - sin_h * sin_h, 0.0)); + let edge = sin_h * ATMOSPHERE_SUN_ANGULAR_RADIUS; + return atmosphere_transmittance_to_top(r, mu_s) * smoothstep(-edge, edge, mu_s - cos_h); +} + +fn atmosphere_scattering_uvwz(r: f32, mu: f32, mu_s: f32, nu: f32, ground: bool) -> vec4 { + let size = ATMOSPHERE_SCATTERING_SIZE; + let h = atmosphere_span(); + let rho = sqrt(max(r * r - ATMOSPHERE_BOTTOM * ATMOSPHERE_BOTTOM, 0.0)); + let u_r = atmosphere_coord(rho / h, size.w); + let r_mu = r * mu; + let discriminant = r_mu * r_mu - r * r + ATMOSPHERE_BOTTOM * ATMOSPHERE_BOTTOM; + var u_mu: f32; + if ground { + let d = -r_mu - sqrt(max(discriminant, 0.0)); + let d_min = r - ATMOSPHERE_BOTTOM; + let d_max = rho; + let x = select((d - d_min) / (d_max - d_min), 0.0, d_max == d_min); + u_mu = 0.5 - 0.5 * atmosphere_coord(x, size.z / 2.0); + } else { + let d = -r_mu + sqrt(max(discriminant + h * h, 0.0)); + let d_min = ATMOSPHERE_TOP - r; + let d_max = rho + h; + u_mu = 0.5 + 0.5 * atmosphere_coord((d - d_min) / (d_max - d_min), size.z / 2.0); + } + let d = atmosphere_to_top(ATMOSPHERE_BOTTOM, mu_s); + let d_min = ATMOSPHERE_TOP - ATMOSPHERE_BOTTOM; + let d_max = h; + let a = (d - d_min) / (d_max - d_min); + let big_a = (atmosphere_to_top(ATMOSPHERE_BOTTOM, ATMOSPHERE_MU_S_MIN) - d_min) / (d_max - d_min); + let u_mu_s = atmosphere_coord(max(1.0 - a / big_a, 0.0) / (1.0 + a), size.y); + return vec4((nu + 1.0) / 2.0, u_mu_s, u_mu, u_r); +} + +// A scattering table at `uvwz`: trilinear in the 3D texture, and linear between the +// two slices of nu around it. +fn atmosphere_lookup(table: texture_3d, uvwz: vec4) -> vec4 { + let size = ATMOSPHERE_SCATTERING_SIZE; + let x = uvwz.x * (size.x - 1.0); + let slice = min(floor(x), size.x - 2.0); + let lerp = x - slice; + let a = vec3((slice + uvwz.y) / size.x, uvwz.z, uvwz.w); + let b = vec3((slice + 1.0 + uvwz.y) / size.x, uvwz.z, uvwz.w); + return mix( + textureSampleLevel(table, atmosphere_sampler, a, 0.0), + textureSampleLevel(table, atmosphere_sampler, b, 0.0), + lerp, + ); +} + +fn atmosphere_rayleigh_phase(nu: f32) -> f32 { + return 3.0 / (16.0 * PI) * (1.0 + nu * nu); +} + +fn atmosphere_mie_phase(nu: f32) -> f32 { + let g = ATMOSPHERE_MIE_G; + let k = 3.0 / (8.0 * PI) * (1.0 - g * g) / (2.0 + g * g); + return k * (1.0 + nu * nu) / pow(1.0 + g * g - 2.0 * g * nu, 1.5); +} + +// The sky's irradiance on a horizontal surface at radius `r`, with the Sun at cosine +// `mu_s` from the zenith, as a fraction of the Sun's irradiance. +fn atmosphere_sky_irradiance(r: f32, mu_s: f32) -> vec4 { + let uv = vec2( + atmosphere_coord(mu_s * 0.5 + 0.5, ATMOSPHERE_IRRADIANCE_SIZE.x), + atmosphere_coord( + (r - ATMOSPHERE_BOTTOM) / (ATMOSPHERE_TOP - ATMOSPHERE_BOTTOM), + ATMOSPHERE_IRRADIANCE_SIZE.y, + ), + ); + return textureSampleLevel(atmosphere_irradiance_table, atmosphere_sampler, uv, 0.0); +} + +// What the atmosphere does to a ray from a viewer in space: the light it scatters into +// the ray, per unit of sunlight, and the transmittance along it. +struct AtmosphereLight { + scattered: vec4, + transmittance: vec4, +} + +// Where the ray from `camera` along unit `view` enters the atmosphere: the distance to +// it, or MISS if the ray passes outside. The miss distance comes from the point of +// closest approach, as `sphere_hit` takes it, so a camera far from Earth keeps its +// precision. A camera inside the atmosphere enters at distance 0. +fn atmosphere_entry(camera: vec3, view: vec3) -> f32 { + let b = dot(camera, view); + let q = length(camera - b * view); + if q > ATMOSPHERE_TOP { + return MISS; + } + if length(camera) <= ATMOSPHERE_TOP { + return 0.0; + } + let t = -b - sqrt((ATMOSPHERE_TOP - q) * (ATMOSPHERE_TOP + q)); + return select(MISS, t, t > 0.0); +} + +// Along a ray from `camera` in unit direction `view` that crosses the atmosphere without +// meeting the ground, with the Sun along unit `sun`. Zero light and full transmittance +// for a ray that misses it. +fn atmosphere_sky(camera: vec3, view: vec3, sun: vec3) -> AtmosphereLight { + let entry = atmosphere_entry(camera, view); + if entry < 0.0 { + return AtmosphereLight(vec4(0.0), vec4(1.0)); + } + let x = camera + entry * view; + let r = max(length(x), ATMOSPHERE_BOTTOM); + 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 uvwz = atmosphere_scattering_uvwz(r, mu, mu_s, nu, false); + let rayleigh = atmosphere_lookup(atmosphere_rayleigh_table, uvwz); + let mie = atmosphere_lookup(atmosphere_mie_table, uvwz); + var light: AtmosphereLight; + light.scattered = rayleigh * atmosphere_rayleigh_phase(nu) + mie * atmosphere_mie_phase(nu); + light.transmittance = atmosphere_transmittance_to_top(r, mu); + 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 { + 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); + 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)); + // 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)) + * smoothstep(0.0, 0.01, mu_s); + light.scattered = rayleigh_between * atmosphere_rayleigh_phase(nu) + + mie_between * atmosphere_mie_phase(nu); + return light; +} + +// 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 +// surface there. +fn atmosphere_ground_light(point: vec3, sun: vec3) -> mat2x4 { + let r = clamp(length(point), ATMOSPHERE_BOTTOM, ATMOSPHERE_TOP); + let mu_s = clamp(dot(point, sun) / length(point), -1.0, 1.0); + return mat2x4(atmosphere_transmittance_to_sun(r, mu_s), atmosphere_sky_irradiance(r, mu_s)); +} diff --git a/crates/resolved/src/shaders/bodies.wgsl b/crates/resolved/src/shaders/bodies.wgsl index 0dfc647..2299067 100644 --- a/crates/resolved/src/shaders/bodies.wgsl +++ b/crates/resolved/src/shaders/bodies.wgsl @@ -16,7 +16,7 @@ struct Uniforms { // fx, fy, cx, cy in pixels. intrinsics: vec4, // x: near clip distance. y: 1 if Earth is drawn from the texture. z: 1 if Earth - // and the Moon shadow each other. + // and the Moon shadow each other. w: 1 if Earth is seen through its atmosphere. params: vec4, // xyz: Earth center. w: equatorial radius. earth_center: vec4, @@ -33,8 +33,9 @@ struct Uniforms { earth_irradiance: vec4, // Columns: per-channel albedo of texture red, green and blue reflectance. earth_rgb: mat3x4, - // Per channel: albedo added to the texture's (stand-in for the atmosphere). - earth_haze: vec4, + // xyz: unit vector toward the Sun in Earth-fixed axes stretched along the polar + // axis, the atmosphere's frame. + earth_sun_stretched: vec4, // xyz: Moon center. w: radius. moon_center: vec4, moon_sun: vec4, @@ -138,23 +139,57 @@ fn read_map(map: texture_2d, at: MapLookup) -> vec4 { return textureSampleGrad(map, earth_sampler, at.uv, at.dx, at.dy); } +// The camera's position relative to Earth's center in the atmosphere's frame, in km. +fn camera_stretched() -> vec3 { + return earth_stretched_km(transpose(u.earth_axes) * (-u.earth_center.xyz), u.earth_shape.x); +} + // Earth's channel radiance at camera-frame point `p`, where its Earth-fixed normal is -// `n_fixed`, with the texture read at `at`. +// `n_fixed`, with the texture read at `at`. Through the atmosphere, the ground is lit +// by the sunlight that reaches it and by the sky, and seen through the air between it +// and the camera, which adds the sunlight it scatters along the way. fn earth_radiance(p: vec3, n_fixed: vec3, at: MapLookup) -> vec4 { let n = u.earth_axes * n_fixed; let lit = max(dot(n, u.earth_sun.xyz), 0.0); - var radiance: vec4; - if u.params.y > 0.5 { - let rgb = read_map(earth_texture, at).rgb; - radiance = u.earth_irradiance * textured_albedo(rgb, u.earth_rgb, u.earth_haze) * lit; - } else { - radiance = u.earth_lambert * lit; + // The Moon's shadow: a solar eclipse. Through the atmosphere it also dims the sky + // and the air between, taken as at the ground point. + var eclipse = vec4(1.0); + if u.params.z > 0.5 && (lit > 0.0 || u.params.w > 0.5) { + eclipse = sun_visible(p, u.sun_center, u.moon_center); } - // The Moon's shadow: a solar eclipse. - if u.params.z > 0.5 && lit > 0.0 { - radiance *= sun_visible(p, u.sun_center, u.moon_center); + if u.params.y < 0.5 { + return u.earth_lambert * lit * eclipse; } - return radiance; + let albedo = textured_albedo(read_map(earth_texture, at).rgb, u.earth_rgb); + if u.params.w < 0.5 { + return u.earth_irradiance * albedo * lit * eclipse; + } + 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; +} + +// The light Earth's atmosphere scatters toward the camera along `d` where the ray does +// not meet the ground: the glow over the limb. Zero where the nearest body along the +// ray, at distance `nearest`, lies in front of the atmosphere, or without it. +fn earth_glow(d: vec3, nearest: f32) -> vec4 { + if u.params.w < 0.5 { + return vec4(0.0); + } + let stretched = transpose(u.earth_axes) * seen(d, 0u) * vec3(1.0, 1.0, u.earth_shape.x); + let length_s = length(stretched); + let view = stretched / length_s; + let camera = camera_stretched(); + let entry = atmosphere_entry(camera, view); + // The entry's distance in meters along the unstretched ray. + if entry < 0.0 || nearest < entry * 1000.0 / length_s { + return vec4(0.0); + } + return PI * u.earth_irradiance * atmosphere_sky(camera, view, u.earth_sun_stretched.xyz).scattered; } // The Moon's channel radiance at camera-frame point `p`, seen along `d`, where its @@ -263,6 +298,9 @@ fn trace(d: vec3, spacing: f32) -> Hit { best = ts; radiance = sun_radiance(ds); } + if !(te > 0.0 && te == best) { + radiance += earth_glow(d, best); + } return Hit(radiance, best); } @@ -381,6 +419,11 @@ fn fs_main(@builtin(position) pos: vec4) -> Out { best = ts; radiance = sun_radiance(ds); } + var glow = vec4(0.0); + if !(te > 0.0 && te == best) { + glow = earth_glow(d, best); + radiance += glow; + } let body_depth = select(0.0, u.params.x / (best * d.z), best < FAR); let pixel_index = vec2(pos.xy); @@ -444,7 +487,7 @@ fn fs_main(@builtin(position) pos: vec4) -> Out { } } } - if best == FAR && mesh_count == 0u { + if best == FAR && mesh_count == 0u && all(glow == vec4(0.0)) { discard; } let covered = f32(mesh_count) / f32(SAMPLES); diff --git a/crates/resolved/src/shaders/common.wgsl b/crates/resolved/src/shaders/common.wgsl index fe8f3d4..cd24534 100644 --- a/crates/resolved/src/shaders/common.wgsl +++ b/crates/resolved/src/shaders/common.wgsl @@ -67,9 +67,15 @@ fn moon_phase_curve(to_sun: vec3, to_viewer: vec3) -> f32 { // Per-channel albedo of Earth's surface from a texel's linear red, green and blue // reflectance: `rgb_to_channels` folds the RGB basis spectra into the camera's -// channels, and `haze` stands in for the atmosphere. -fn textured_albedo(rgb: vec3, rgb_to_channels: mat3x4, haze: vec4) -> vec4 { - return rgb_to_channels * rgb + haze; +// channels. +fn textured_albedo(rgb: vec3, rgb_to_channels: mat3x4) -> vec4 { + return rgb_to_channels * rgb; +} + +// Earth-fixed coordinates stretched along the polar axis by `flattening` (equatorial +// over polar radius), where Earth is a sphere: the frame of atmosphere.wgsl, in km. +fn earth_stretched_km(fixed: vec3, flattening: f32) -> vec3 { + return fixed * vec3(1.0, 1.0, flattening) * 0.001; } // How the solar disk's brightness is spread, per channel (sun.rs, `LimbDarkening`). A diff --git a/crates/resolved/src/shaders/earthshine.wgsl b/crates/resolved/src/shaders/earthshine.wgsl index 701fcb1..ae69c77 100644 --- a/crates/resolved/src/shaders/earthshine.wgsl +++ b/crates/resolved/src/shaders/earthshine.wgsl @@ -49,11 +49,10 @@ struct Globals { earth_albedo: vec4, // Columns: per-channel albedo of texture red, green and blue reflectance. earth_rgb: mat3x4, - // Per channel: albedo added to the texture's (stand-in for the atmosphere). - earth_haze: vec4, moon_albedo: vec4, // x: 1 if Earth is textured. y: meters per texel of the texture in latitude, at - // its full size. z: Earth's equatorial / polar radius. + // its full size. z: Earth's equatorial / polar radius. w: 1 if Earth is seen + // through its atmosphere. params: vec4, // Rings and sectors per body, cube face size in texels, number of targets. grid: vec4, @@ -107,7 +106,7 @@ fn earth_surface(t: TargetLight, d: vec3, te: f32) -> mat2x4 { let footprint = te * sqrt(t.earth.cone.y / slant); let level = max(log2(footprint / g.params.y), 0.0); let rgb = textureSampleLevel(earth_texture, earth_sampler, equirectangular(n_fixed), level).rgb; - albedo = textured_albedo(rgb, g.earth_rgb, g.earth_haze); + albedo = textured_albedo(rgb, g.earth_rgb); } return mat2x4(albedo, vec4(n, 0.0)); } @@ -127,6 +126,21 @@ fn moon_albedo(t: TargetLight, d: vec3, tm: f32, n: vec3) -> vec4 return g.moon_color * (reflectance * moon_phase_curve(t.moon.sun.xyz, -d)); } +// 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 { + 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)); + 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; +} + @compute @workgroup_size(64) fn sources(@builtin(global_invocation_id) id: vec3) { let cells = g.grid.x * g.grid.y; @@ -169,8 +183,12 @@ fn sources(@builtin(global_invocation_id) id: vec3) { visible = sun_visible(tm * d, t.sun, t.earth.center); } 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; + } var out: Sample; - out.light = albedo * lit * (source.sun.w * source.cone.y); + out.light = light * (source.sun.w * source.cone.y); out.direction = vec4(d, 0.0); samples[target_index * 2u * cells + id.x] = out; } diff --git a/crates/resolved/tests/gpu/atmosphere.rs b/crates/resolved/tests/gpu/atmosphere.rs new file mode 100644 index 0000000..9fa1d70 --- /dev/null +++ b/crates/resolved/tests/gpu/atmosphere.rs @@ -0,0 +1,227 @@ +//! Earth's atmosphere at lighting L1: what the shaders draw against the tables they +//! read (evaluated in f64 on the CPU), and the limb's glow against the physics, on the +//! GPU and on the CPU (lavapipe). + +use crate::common; + +use resolved::atmosphere::{AtmosphereTables, BOTTOM_RADIUS, TOP_RADIUS}; +use resolved::gpu::AdapterChoice; +use resolved::nalgebra::{UnitQuaternion, Vector3}; +use resolved::radiometry::{Radiometry, AU}; +use std::f64::consts::PI; +use std::sync::Arc; + +use resolved::{ + CameraConfig, EarthTexture, Fidelity, LightingTier, Pose, Renderer, RendererOptions, Scene, +}; + +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 { + let pixels = [GRAY, GRAY, GRAY, 255].repeat(2); + let options = RendererOptions { + adapter: adapter.clone(), + fidelity: Fidelity { + lighting: LightingTier::L1, + ..Default::default() + }, + earth_texture: Some(Arc::new(EarthTexture::from_rgba8(2, 1, pixels).unwrap())), + moon_texture: Some(common::uniform_moon()), + read_radiance: true, + ..Default::default() + }; + Renderer::new(CameraConfig::example(), options).unwrap() +} + +/// 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) { + let scene = Scene { + earth_attitude: UnitQuaternion::identity(), + sun: sun.normalize() * AU, + moon: Vector3::new(0.0, 0.0, -3.844e8), + moon_attitude: UnitQuaternion::identity(), + targets: vec![], + velocities: Default::default(), + }; + (scene, Pose::looking_at(camera, at, Vector3::z()).unwrap()) +} + +/// The mean of the four pixels around the boresight, whose lines of sight are +/// 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] { + 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)] { + for (k, v) in radiance[y * 1920 + x].iter().enumerate() { + sum[k] += v / 4.0; + } + } + sum +} + +/// Per channel: the Sun's irradiance at 1 AU over pi, and the gray surface's albedo. +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); + let albedo = r.effective_albedo(&resolved::model::rgb_spectrum([linear; 3])); + let sun = [0, 1, 2].map(|k| r.lambert_radiance(k, 1.0, 1.0, AU)); + (r, sun, [albedo[0], albedo[1], albedo[2]]) +} + +/// Where the ray from `c` along unit `v` meets a sphere of radius `radius` first. +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() +} + +/// What atmosphere.wgsl computes for the ground seen from `camera` along unit `view`, +/// with the Sun along unit `sun`, from the tables in f64. +fn ground_reference( + tables: &AtmosphereTables, + camera: Vector3, + view: Vector3, + sun: Vector3, +) -> [f64; 3] { + let (_, irradiance, albedo) = constants(); + let x = camera + hit(camera, view, TOP_RADIUS) * view; + let ground = camera + hit(camera, view, BOTTOM_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 mu_d = (r * mu + d) / r_d; + let mu_s_d = (r * mu_s + d * nu) / r_d; + let (near, far) = ( + tables.transmittance(r_d, -mu_d), + tables.transmittance(r, -mu), + ); + 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); + 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) + }) +} + +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!( + error.abs() < tolerance, + "{context} channel {k}: {} vs {} ({:+.3} %)", + got[k], + want[k], + 100.0 * error + ); + } +} + +#[test] +fn 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, with the Sun + // 30 degrees from its zenith on the other side. + 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 want = ground_reference(&tables, camera, -toward, sun); + // The air adds blue and takes red: more scattered than transmitted in blue. + let bare = constants(); + let surface = [0, 1, 2].map(|k| bare.1[k] * bare.2[k] * 30f64.to_radians().cos()); + assert!( + want[2] > surface[2] && want[0] < 1.3 * surface[0], + "{want:?} {surface:?}" + ); + for adapter in &ADAPTERS { + let mut renderer = renderer(adapter); + let (scene, pose) = view(sun, camera, ground); + let frame = renderer.submit(&scene, &pose).unwrap().wait().unwrap(); + // Measured within 2.5e-4 on both devices, with the tables in f16. + assert_close(center(&frame), want, 1e-3, &renderer.adapter_info().name); + } +} + +/// The light the atmosphere scatters toward a camera looking horizontally past the +/// tangent point `height` above the ground at 0 N 0 E, with the Sun at its zenith: +/// from the renderer's center pixel, and from the tables. +fn limb(renderer: &mut Renderer, tables: &AtmosphereTables, height: f64) -> ([f32; 4], [f64; 3]) { + let tangent = Vector3::new(BOTTOM_RADIUS + height, 0.0, 0.0); + let along = Vector3::y(); + let camera = tangent - 2.5e6 * along; + let sun = Vector3::x(); + let (scene, pose) = view(sun, camera, tangent); + let frame = renderer.submit(&scene, &pose).unwrap().wait().unwrap(); + let x = camera + hit(camera, along, TOP_RADIUS) * along; + let r = TOP_RADIUS; + let s = tables.scattering( + r, + x.dot(&along) / r, + x.dot(&sun) / r, + along.dot(&sun), + false, + ); + let irradiance = constants().1; + ( + center(&frame), + std::array::from_fn(|k| PI * irradiance[k] * s[k]), + ) +} + +#[test] +fn the_limb_glows_and_fades_with_the_scale_height() { + let _devices = common::device_lock(); + let tables = AtmosphereTables::new(&constants().0); + for adapter in &ADAPTERS { + let mut renderer = renderer(adapter); + let name = renderer.adapter_info().name.clone(); + let (low, low_want) = limb(&mut renderer, &tables, 40e3); + let (high, high_want) = limb(&mut renderer, &tables, 50e3); + assert_close(low, low_want, 1e-3, &format!("{name} 40 km")); + assert_close(high, high_want, 1e-3, &format!("{name} 50 km")); + // Thin air seen edge on: the glow is in proportion to the density at the + // tangent point, which falls by e over Rayleigh's 8 km scale height, so 10 km + // higher is e^1.25 = 3.49 times fainter. + let ratio = f64::from(low[1]) / f64::from(high[1]); + assert!((ratio - 3.49).abs() < 0.2, "{name}: {ratio}"); + // And blue: air scatters blue light more. + assert!(low[2] > 2.0 * low[0], "{name}: {low:?}"); + } +} + +#[test] +fn space_beside_the_limb_stays_dark() { + // 150 km above the ground, over the top of the atmosphere: nothing. + let _devices = common::device_lock(); + let tables = AtmosphereTables::new(&constants().0); + let mut renderer = renderer(&AdapterChoice::Gpu); + let (glow, _) = limb(&mut renderer, &tables, 150e3); + assert_eq!(glow, [0.0; 4]); +} + +#[test] +fn a_camera_inside_the_atmosphere_is_refused() { + let _devices = common::device_lock(); + let mut renderer = renderer(&AdapterChoice::Gpu); + let ground = Vector3::new(BOTTOM_RADIUS, 0.0, 0.0); + let (scene, pose) = view(Vector3::x(), ground * (1.0 + 80e3 / BOTTOM_RADIUS), ground); + let error = renderer.submit(&scene, &pose).err().unwrap().to_string(); + assert!( + error.contains("80 km up, inside Earth's atmosphere"), + "{error}" + ); +} diff --git a/crates/resolved/tests/gpu/earth_texture.rs b/crates/resolved/tests/gpu/earth_texture.rs index c110c2a..8094245 100644 --- a/crates/resolved/tests/gpu/earth_texture.rs +++ b/crates/resolved/tests/gpu/earth_texture.rs @@ -12,7 +12,6 @@ use std::sync::Arc; use resolved::{ CameraConfig, EarthTexture, Fidelity, LightingTier, Pose, Renderer, RendererOptions, Scene, - Spectrum, }; const ADAPTERS: [AdapterChoice; 2] = [AdapterChoice::Gpu, AdapterChoice::Cpu]; @@ -20,12 +19,13 @@ const ADAPTERS: [AdapterChoice; 2] = [AdapterChoice::Gpu, AdapterChoice::Cpu]; /// The camera's distance from Earth's center: Earth is about 580 px across. const RANGE: f64 = 2.0e8; -/// A lighting tier L1 renderer drawing Earth from `texture`. -fn renderer(adapter: &AdapterChoice, haze: bool, texture: &EarthTexture) -> Renderer { - let mut appearance = Appearance::default(); - if !haze { - appearance.earth_haze = Spectrum::constant(0.0); - } +/// A lighting tier L1 renderer drawing Earth's bare surface from `texture`, without its +/// atmosphere (tests/gpu/atmosphere.rs tests that). +fn renderer(adapter: &AdapterChoice, texture: &EarthTexture) -> Renderer { + let appearance = Appearance { + earth_atmosphere: false, + ..Default::default() + }; let options = RendererOptions { adapter: adapter.clone(), fidelity: Fidelity { @@ -78,7 +78,7 @@ fn linear(byte: u8) -> f64 { } /// Channel radiance of Earth's surface with sRGB-encoded reflectance `rgb`, facing the -/// Sun at 1 AU, without haze. +/// Sun at 1 AU. fn expected(rgb: [u8; 3]) -> [f64; 3] { let r = Radiometry::new(&CameraConfig::example()); let albedo = r.effective_albedo(&rgb_spectrum([rgb[0], rgb[1], rgb[2]].map(linear))); @@ -133,7 +133,7 @@ fn texel_under_the_camera_sets_the_radiance() { // Over 45 N 45 E the camera sees the middle of region 2 (longitude 0..90, north). let _devices = common::device_lock(); for adapter in &ADAPTERS { - let mut renderer = renderer(adapter, false, ®ions()); + let mut renderer = renderer(adapter, ®ions()); let (scene, camera) = view(45.0, 45.0, UnitQuaternion::identity()); let frame = renderer.submit(&scene, &camera).unwrap().wait().unwrap(); assert_close( @@ -151,7 +151,7 @@ fn earth_rotation_moves_the_surface() { let _devices = common::device_lock(); let turn = UnitQuaternion::from_axis_angle(&Vector3::z_axis(), 90f64.to_radians()); for adapter in &ADAPTERS { - let mut renderer = renderer(adapter, false, ®ions()); + let mut renderer = renderer(adapter, ®ions()); for (lat, region) in [(45.0, 1), (-45.0, 5)] { let (scene, camera) = view(lat, 45.0, turn); let frame = renderer.submit(&scene, &camera).unwrap().wait().unwrap(); @@ -165,23 +165,6 @@ fn earth_rotation_moves_the_surface() { } } -#[test] -fn haze_adds_its_spectrum() { - let _devices = common::device_lock(); - let r = Radiometry::new(&CameraConfig::example()); - let haze = r.effective_albedo(&Appearance::default().earth_haze); - let surface = r.effective_albedo(&rgb_spectrum(REGIONS[2].map(linear))); - let want = [0, 1, 2].map(|k| r.lambert_radiance(k, surface[k] + haze[k], 1.0, AU)); - // Bluer haze: it raises blue the most. - assert!(haze[2] > haze[1] && haze[1] > haze[0], "{haze:?}"); - for adapter in &ADAPTERS { - let mut renderer = renderer(adapter, true, ®ions()); - let (scene, camera) = view(45.0, 45.0, UnitQuaternion::identity()); - let frame = renderer.submit(&scene, &camera).unwrap().wait().unwrap(); - assert_close(center(&frame), want, &renderer.adapter_info().name); - } -} - #[test] fn no_seam_at_180_degrees() { // The four texel columns either side of 180 degrees (22.5 degrees each way) are one @@ -206,7 +189,7 @@ fn no_seam_at_180_degrees() { let texture = EarthTexture::from_rgba8(w, h, pixels).unwrap(); let want = expected(edge); for adapter in &ADAPTERS { - let mut renderer = renderer(adapter, false, &texture); + let mut renderer = renderer(adapter, &texture); // Over 179.8 E the seam is about one pixel right of center, inside a 2 x 2 pixel // quad (screen-space derivatives are taken per quad). Over exactly 180 it would // fall on a quad boundary, where even unhandled it shows nothing. @@ -251,7 +234,7 @@ fn far_views_average_rather_than_alias() { // 50 % linear reflectance is sRGB 188. let gray = EarthTexture::from_rgba8(1, 1, vec![188, 188, 188, 255]).unwrap(); for adapter in &ADAPTERS { - let mut renderer = renderer(adapter, false, &gray); + let mut renderer = renderer(adapter, &gray); let (scene, mut camera) = view(10.0, 30.0, UnitQuaternion::identity()); camera.position *= 2.0; let mut block = |texture: &EarthTexture| { diff --git a/crates/resolved/tests/gpu/earthshine.rs b/crates/resolved/tests/gpu/earthshine.rs index b10aa1e..96196f9 100644 --- a/crates/resolved/tests/gpu/earthshine.rs +++ b/crates/resolved/tests/gpu/earthshine.rs @@ -53,8 +53,8 @@ fn plate() -> Mesh { Mesh::from_stl(&mut Cursor::new(s), "plate", &options).unwrap() } -/// What a test varies about the bodies: Earth's texture (one sRGB gray level, no haze) -/// and the Moon's albedo. +/// What a test varies about the bodies: Earth's texture (one sRGB gray level, without +/// the atmosphere) and the Moon's albedo. #[derive(Clone, Copy)] struct Surfaces { earth_level: u8, @@ -69,7 +69,7 @@ impl Surfaces { fn appearance(&self) -> Appearance { Appearance { - earth_haze: Spectrum::constant(0.0), + earth_atmosphere: false, moon_albedo: Spectrum::constant(self.moon_albedo), moon_color: Spectrum::constant(1.0), ..Default::default() diff --git a/crates/resolved/tests/gpu/main.rs b/crates/resolved/tests/gpu/main.rs index 0808125..2dd6146 100644 --- a/crates/resolved/tests/gpu/main.rs +++ b/crates/resolved/tests/gpu/main.rs @@ -6,6 +6,7 @@ mod common; mod analytic; mod apparent; +mod atmosphere; mod blooming; mod bodies; mod convolve; diff --git a/docs/optics.md b/docs/optics.md index 02cdd8d..7428dad 100644 --- a/docs/optics.md +++ b/docs/optics.md @@ -1012,8 +1012,9 @@ lit) right. It misses some things that are visible in real images: old basalt lava, the "seas" seen with the naked eye) and its bright highlands also differ in albedo by about a factor of two, which from lighting L1 the Moon's map shows. -- Earth's atmosphere softens the limb and adds a thin bright band around it - (not modeled yet). Here Earth's edge is a sharp line. +- Earth's atmosphere softens the limb and adds a thin glowing band around it. + At lighting L0 Earth's edge is a sharp line; from L1 the atmosphere is drawn + ("Earth's atmosphere"). ### Earth's surface @@ -1057,12 +1058,8 @@ Two details matter for brightness: those values. A coastline seen up close blends correctly from sea to land. - Blue Marble is surface reflectance with no atmosphere, and open ocean is dark, a few percent. From orbit the ocean looks bright blue because the air - above it scatters sunlight toward the camera: Rayleigh scattering, which - grows as 1 / wavelength^4, which is also why the sky is blue. Until the - atmosphere is modeled, a stand-in adds a reflectance of that shape, - 0.07 at 450 nm falling to 0.017 at 650 nm (`Appearance::earth_haze`). It - brightens and blues the whole lit disk evenly, where the real effect grows - toward the limb, where the line of sight crosses more air. + above it scatters sunlight toward the camera. The renderer adds that air + itself, next. The texture's red, green and blue are reflectances in the same three bands a glTF base color uses (see the target model section), folded into the camera's @@ -1073,6 +1070,143 @@ calibrated reflectance product itself, so its values are close but not exact reflectances. A calibrated alternative (MODIS's albedo product) is a possible later refinement. +### Earth's atmosphere + +Code: `atmosphere.rs`, `shaders/atmosphere.wgsl`, `shaders/bodies.wgsl`. + +At lighting L1 Earth is seen through its atmosphere +(`Appearance::earth_atmosphere`, on by default). Air does three things to +light passing through it: + +- **Molecules scatter it**, much more at short wavelengths: Rayleigh + scattering grows as 1 / wavelength^4, so blue light at 450 nm scatters 3.6 + times as much as red at 620 nm. Straight up from the ground, air scatters + about 10 % of green light (an optical depth of 0.11 at 550 nm) and 21 % of + blue (0.24 at 450 nm). That is why the sky is blue and the setting Sun red. +- **Aerosols** (dust, sea salt, smoke droplets) scatter colors far more + evenly than air does, so they make haze white rather than blue, and mostly + forward, within a few tens of degrees of the light's direction. They sit + low: half of them in the lowest 0.8 km. Straight up they scatter or absorb + about as much light as the air does, an optical depth of 0.12 at 550 nm on + average over Earth, and far more over deserts and fires. +- **Ozone absorbs** light in a broad band centered near 600 nm (the Chappuis + band). From space it matters mainly at the limb, where the line of sight + crosses the ozone layer edge on, and tints the glow there. + +What leaves the ground is reduced by the light scattered or absorbed out of +it on the way to the camera: the **transmittance**, between 0 and 1. And the +air along the line of sight sends sunlight it scatters toward the camera: the +**path radiance**, which adds to whatever is behind it. Both grow with the +amount of air the line of sight crosses, so the ground looks hazier and bluer +toward Earth's edge, where the line of sight is long and slanting. Past the +edge the line of sight crosses only air, which glows: the **limb**. There the +air is thin enough that the glow is in proportion to the density at the +lowest point of the line of sight, so it fades by a factor of e (2.7) for +every 8 km of height. + +The ground itself is lit by the sunlight that reaches it, which the air has +also dimmed (more at a low Sun, and more in blue, so ground near sunset is lit +redder), and by **skylight**: the light the whole sky scatters down, about a +tenth of the Sun's in blue with the Sun overhead and less in red. Past the +terminator the sky still lights the ground for a while: twilight. + +Light can also scatter more than once, from air to air or off the ground and +back. With the Sun high that adds a few tens of percent to the sky's +brightness. + +Adding up the scattered light along every line of sight for every pixel would +be far too slow. Instead the renderer precomputes tables once, at setup, with +the method of Bruneton and Neyret (2008), in the form of Bruneton's 2017 +implementation and with its atmosphere: what a point at a given height sees +looking in a given direction, with the Sun at a given angle, is a function of +four numbers (height, view angle from the zenith, Sun angle from the zenith, +and the angle between the view and the Sun), so it fits a table. Three tables: + +- the scattered light (4D, stored as a 3D texture of 256 x 128 x 16 texels, + Rayleigh and aerosols separately); +- the transmittance from any point to the top of the atmosphere (256 x 64), + from which the transmittance between any two points follows; +- the skylight on the ground (64 x 16). + +Each pixel that sees Earth then reads them about a dozen times: the +transmittance to the ground and the light scattered along the way (the light +scattered along the whole line of sight, less what lies beyond the ground), +and the sunlight and skylight on the ground there. + +Light scattered more than once is approximated as Hillaire (2020) does: +assume the light arriving at a point after one scattering comes equally from +every direction and is the same at nearby points. Then all the higher orders +together follow from a small table per height and Sun angle, instead of being +computed order by order. + +The air is Bruneton's, with more aerosols: Rayleigh scattering falling with +an 8 km scale height; aerosols with a 1.2 km scale height, an optical depth of +0.12 at 550 nm (the global annual mean of Kinne et al.'s 2013 climatology) +falling as wavelength^-1.3, 90 % of it scattering, with Bruneton's forward +phase function (asymmetry 0.8); and 300 Dobson units of ozone between 10 and +40 km. Bruneton's aerosol depth of 0.0053 is air cleaner than almost anywhere +on Earth: with it the ground looked bluer and sharper toward the limb than it +should. The top of the atmosphere is set +100 km up. Earth is the sphere of its equatorial radius, in axes stretched +along the polar axis, where its ellipsoid is that sphere (heights are off by +at most 0.3 % near the poles). + +The tables are computed every 20 nm and folded into the camera's channels +with sunlight and each channel's QE, like every other spectrum +("Folding spectra into per-channel constants"). The scattered light folds +exactly. Transmittance does not: a channel whose band spans very different +transmittances gets their average weighted by sunlight, which is wrong when +the light passing through is a different color from sunlight. It matters most +for blue light along long, grazing paths. Setup takes 0.25 s on the 24-core +desktop CPU these numbers come from, and 1.1 s on four of its cores. + +Tests: the tables against direct integrals, in `atmosphere.rs` (the amount of +air from the ground straight up within 1e-4 of the scale height, Rayleigh +scattered light within 3.4e-4, aerosols within 0.31 %); the shaders against +the tables, on the GPU and lavapipe (`tests/gpu/atmosphere.rs`): the ground +seen at 60 degrees from its zenith and the limb at 40 and 50 km, within +2.5e-4. The limb fades by 3.38 over those 10 km, against 3.49 for the density +alone; the air between takes the rest. + +Checked against a real satellite: GOES-16's full disk of 2024-07-01 at +15:00 UTC, in its blue (470 nm) and red (640 nm) bands, over 640,000 cloud-free +ocean pixels (NOAA's cloud mask), away from the Sun's glint. Over ocean, +which is dark, most of the light the satellite measures comes from the air, +and GOES sees the sea from straight above out to its limb: the same slanting +paths through the air as a camera in low orbit looking toward the horizon. +The tables, evaluated at each pixel's angles with Blue Marble's ocean color +and channels shaped like GOES's bands, against the measured reflectance +(the median in each band of viewing angle from the zenith): + +| Viewing angle | GOES blue | Model | Clean air | GOES red | Model | Clean air | +|---|---|---|---|---|---|---| +| 20-40 degrees | 0.097 | 0.090 | 0.077 | 0.029 | 0.030 | 0.023 | +| 55-65 | 0.146 | 0.145 | 0.123 | 0.044 | 0.051 | 0.037 | +| 70-75 | 0.205 | 0.191 | 0.169 | 0.076 | 0.071 | 0.056 | +| 78-81 | 0.249 | 0.239 | 0.219 | 0.093 | 0.095 | 0.078 | +| 81-84 | 0.267 | 0.266 | 0.246 | 0.104 | 0.110 | 0.092 | + +The model is within 7 % in blue (low: Blue Marble's ocean is darker in blue +than real open ocean) and 15 % in red out to 84 degrees, and the haze grows +toward the limb as measured. With Bruneton's clean air ("Clean air") it is 8 +to 29 % too dark. NOAA's own retrieval put the aerosols' optical depth over +these pixels at 0.13 (median), close to the 0.12 the model assumes. + +On the RTX 3060 the atmosphere adds 0.8 to 1.0 ms to the bodies pass at +1920 x 1080 with camera C1's 2 x 2 supersampling, and 14 to 27 ms on lavapipe. + +Limitations: + +- **The camera must be above the atmosphere**, over 100 km up; below that a + frame is refused. +- **Sunlight reaching the target** through the atmosphere, and the Sun and + Moon seen through it at the limb, are not dimmed or reddened yet. +- **Clouds** are not drawn: Blue Marble is a cloud-free mosaic. +- **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. + ### The Moon's surface Code: `moon.rs`, `shaders/bodies.wgsl`; data from `scripts/fetch-moon.py`. From @@ -1531,7 +1665,8 @@ arrives from Earth and the Moon on a surface facing any direction: through each cell's center finds the point it meets on the body, and that point is shaded exactly as the camera's view of the body is shaded: reflectance from the Blue Marble image (or the stand-in spectrum), lit by - the cosine of the Sun's angle there, dark on the night side. Each texture + the cosine of the Sun's angle there, dark on the night side, and at + lighting L1 seen through the atmosphere. Each texture read averages the patch of ground the cell covers, about 1000 square kilometers straight down from 500 km and more toward Earth's edge, so a cell sees the mean color of that patch rather than one point. A line @@ -1588,8 +1723,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 or atmosphere yet - beyond the stand-in haze. + phase curve as seen from the target. Earth has no clouds yet. ## Where the Sun and Moon appear @@ -1762,5 +1896,6 @@ is a run of the script with `--gaia 16`, at twice the size. - Non-Lambertian materials (BRDFs), glints, and the full Moon's edge as bright as its center - Clouds on Earth -- Earth's atmosphere: scattering, the bright limb +- Sunlight through Earth's atmosphere onto targets, and the Sun and Moon seen + through it; the atmosphere seen from inside it (a camera below 100 km) - Motion blur, rolling shutter, stray light diff --git a/examples/scenes/moonrise.toml b/examples/scenes/moonrise.toml index 682c12f..c2a4d95 100644 --- a/examples/scenes/moonrise.toml +++ b/examples/scenes/moonrise.toml @@ -6,8 +6,8 @@ # # resolved render --scene examples/scenes/moonrise.toml --lighting L1 --preview moonrise.png # -# At lighting L1 Earth's surface and the Moon's come from their maps. There is no -# atmosphere yet, so the limb is a sharp edge rather than a thin glowing band. +# At lighting L1 Earth's surface and the Moon's come from their maps, and Earth's +# atmosphere hazes the ground toward the limb and glows over it. sun = [129555556378, 74798935350, 0] moon = [-378712487.94, 64656759.20, -12637838.88] diff --git a/scripts/make-example-scenes.py b/scripts/make-example-scenes.py index c7f02f3..7c8d1b2 100755 --- a/scripts/make-example-scenes.py +++ b/scripts/make-example-scenes.py @@ -421,8 +421,8 @@ def moonrise(): # # resolved render --scene examples/scenes/moonrise.toml --lighting L1 --preview moonrise.png # -# At lighting L1 Earth's surface and the Moon's come from their maps. There is no -# atmosphere yet, so the limb is a sharp edge rather than a thin glowing band. +# At lighting L1 Earth's surface and the Moon's come from their maps, and Earth's +# atmosphere hazes the ground toward the limb and glows over it. sun = {fmt([AU * x for x in sun], 0)} moon = {fmt(moon)}