From 8c0a1439d7228d733c0b41a978ec688af79294e3 Mon Sep 17 00:00:00 2001 From: Connor Johnstone Date: Sun, 18 Oct 2026 00:05:42 -0400 Subject: [PATCH] sharper earth map --- Cargo.lock | 39 +++ crates/resolved-cli/Cargo.toml | 6 +- crates/resolved-cli/src/earth_tiles.rs | 205 ++++++++++++ crates/resolved-cli/src/main.rs | 51 ++- crates/resolved-showcase/src/main.rs | 18 +- crates/resolved/Cargo.toml | 2 + crates/resolved/src/earth.rs | 384 +++++++++++++++++++++- crates/resolved/src/gpu/bodies.rs | 186 ++++++++++- crates/resolved/src/gpu/context.rs | 6 +- crates/resolved/src/lib.rs | 2 +- crates/resolved/src/render.rs | 28 +- crates/resolved/src/shaders/bodies.wgsl | 71 +++- crates/resolved/tests/gpu/earth_detail.rs | 229 +++++++++++++ crates/resolved/tests/gpu/main.rs | 1 + docs/optics.md | 46 ++- readme.md | 5 +- scripts/fetch-earth.sh | 189 +++++++++-- scripts/gallery-delta.py | 4 +- 18 files changed, 1426 insertions(+), 46 deletions(-) create mode 100644 crates/resolved-cli/src/earth_tiles.rs create mode 100644 crates/resolved/tests/gpu/earth_detail.rs diff --git a/Cargo.lock b/Cargo.lock index 4df9d63..9c80abe 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -729,6 +729,12 @@ version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" +[[package]] +name = "hermit-abi" +version = "0.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e17592d60ebacc7d5e169f4663c5f84f9161cc90328abcfe8456f41e4dfcb284" + [[package]] name = "hifitime" version = "4.3.1" @@ -777,12 +783,31 @@ version = "1.1.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a257582fdcde896fd96463bf2d40eefea0580021c0712a0e2b028b60b47a837a" +[[package]] +name = "intel_tex_2" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7980d359170e7da4bc2421aa74c49d338a3153313e006d732a3430f1af375a60" +dependencies = [ + "ispc_rt", +] + [[package]] name = "is_terminal_polyfill" version = "1.70.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a6cb138bb79a146c1bd460005623e142ef0181e3d0219cb493e02f7d08a35695" +[[package]] +name = "ispc_rt" +version = "2.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a396f76fd253fab9945bb4c6080f515d0c4036bc30219462d199b58886fd37b4" +dependencies = [ + "libc", + "num_cpus", +] + [[package]] name = "itertools" version = "0.13.0" @@ -1150,6 +1175,16 @@ dependencies = [ "libm", ] +[[package]] +name = "num_cpus" +version = "1.17.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91df4bbde75afed763b708b7eee1e8e7651e02d97f6d5dd763e89367e957b23b" +dependencies = [ + "hermit-abi", + "libc", +] + [[package]] name = "objc2" version = "0.6.5" @@ -1559,6 +1594,7 @@ dependencies = [ "half", "hifitime", "image", + "intel_tex_2", "log", "nalgebra", "pollster", @@ -1571,9 +1607,12 @@ dependencies = [ name = "resolved-cli" version = "0.1.0" dependencies = [ + "bytemuck", "clap", + "half", "hifitime", "image", + "intel_tex_2", "resolved", "serde", "toml", diff --git a/crates/resolved-cli/Cargo.toml b/crates/resolved-cli/Cargo.toml index 2ae756b..fd61fe5 100644 --- a/crates/resolved-cli/Cargo.toml +++ b/crates/resolved-cli/Cargo.toml @@ -25,7 +25,11 @@ resolved.workspace = true clap = { version = "4.6", features = ["derive"] } serde = { version = "1.0", features = ["derive"] } toml = "1.1" -image = { version = "0.25", default-features = false, features = ["png"] } +image = { version = "0.25", default-features = false, features = ["png", "jpeg"] } +# BC6H compression for `earth-tiles` (Intel's ISPC texture compressor, prebuilt). +intel_tex_2 = "0.5" +half = { version = "2.7", features = ["bytemuck"] } +bytemuck = "1.0" hifitime = { version = "4.3", optional = true } [dev-dependencies] diff --git a/crates/resolved-cli/src/earth_tiles.rs b/crates/resolved-cli/src/earth_tiles.rs new file mode 100644 index 0000000..2740fd7 --- /dev/null +++ b/crates/resolved-cli/src/earth_tiles.rs @@ -0,0 +1,205 @@ +//! `resolved earth-tiles`: Blue Marble's sharp 86400 x 43200 map, as NASA publishes it +//! (eight JPEG tiles of 21600 x 21600 per month), cut into the GPU's tiles +//! ([`DetailLayout::BLUE_MARBLE`]) and compressed to BC6H, once, ahead of rendering. + +use std::io::{Seek, SeekFrom, Write}; +use std::path::{Path, PathBuf}; +use std::time::Instant; + +use half::f16; +use resolved::earth::{srgb_to_linear, DetailFormat, DetailLayout, EarthDetail}; + +/// NASA's tiles are 21600 texels on a side: four columns west to east from -180 +/// degrees (A to D), two rows north to south (1, 2). +const SOURCE: u32 = 21600; +const COLUMNS: [char; 4] = ['A', 'B', 'C', 'D']; + +/// The name of NASA's tile at `column`, `row` for `month`. +fn source_name(month: u32, column: u32, row: u32) -> String { + format!( + "world.2004{month:02}.3x21600x21600.{}{}.jpg", + COLUMNS[column as usize], + row + 1 + ) +} + +fn decode(path: &Path) -> Result, String> { + let mut reader = + image::ImageReader::open(path).map_err(|e| format!("{}: {e}", path.display()))?; + reader.no_limits(); + let image = reader + .decode() + .map_err(|e| format!("{}: {e}", path.display()))? + .into_rgb8(); + if image.dimensions() != (SOURCE, SOURCE) { + return Err(format!( + "{}: {:?}, not {SOURCE} x {SOURCE}", + path.display(), + image.dimensions() + )); + } + Ok(image.into_raw()) +} + +/// The strips along each source tile's four edges that its neighbors' borders need: +/// `width` columns at the west and east, `width` rows at the north and south, as RGB. +struct Edges { + west: Vec, + east: Vec, + north: Vec, + south: Vec, +} + +impl Edges { + fn new(rgb: &[u8], width: u32) -> Self { + let (s, w) = (SOURCE as usize, width as usize); + let columns = |from: usize| { + (0..s) + .flat_map(|y| { + rgb[3 * (y * s + from)..3 * (y * s + from + w)] + .iter() + .copied() + }) + .collect() + }; + Self { + west: columns(0), + east: columns(s - w), + north: rgb[..3 * s * w].to_vec(), + south: rgb[3 * s * (s - w)..].to_vec(), + } + } + + /// The texel at `x`, `y` in its tile, which must lie in one of the strips. + fn at(&self, x: u32, y: u32, width: u32) -> [u8; 3] { + let (s, w) = (SOURCE, width); + let pick = |strip: &[u8], i: usize| [strip[3 * i], strip[3 * i + 1], strip[3 * i + 2]]; + if y < w { + pick(&self.north, (y * s + x) as usize) + } else if y >= s - w { + pick(&self.south, ((y - (s - w)) * s + x) as usize) + } else if x < w { + pick(&self.west, (y * w + x) as usize) + } else if x >= s - w { + pick(&self.east, (y * w + x - (s - w)) as usize) + } else { + panic!("texel {x}, {y} is not in a strip") + } + } +} + +/// BC6H blocks of a `size` x `size` RGBA f16 image, its edges repeated out to a +/// multiple of 4, compressed in bands on every core. +fn compress(texels: &[f16], size: u32) -> Vec { + let padded = size.div_ceil(4) * 4; + let mut image = vec![f16::ZERO; 4 * (padded * padded) as usize]; + for y in 0..padded { + let sy = y.min(size - 1); + for x in 0..padded { + let sx = x.min(size - 1); + let (to, from) = (4 * (y * padded + x) as usize, 4 * (sy * size + sx) as usize); + image[to..to + 4].copy_from_slice(&texels[from..from + 4]); + } + } + let bytes: &[u8] = bytemuck::cast_slice(&image); + let stride = 8 * padded as usize; + let rows = padded as usize / 4; + let threads = std::thread::available_parallelism().map_or(1, |n| n.get()); + let per = rows.div_ceil(threads); + // Measured on Blue Marble's July map around Sinai: the 99th percentile of the + // per-pixel error 1.0 to 2.5 %, the mean under 0.1 %; "very fast" was 1.4 to 2.9 %, + // and "basic" no better than this. + let settings = intel_tex_2::bc6h::fast_settings(); + std::thread::scope(|scope| { + let bands: Vec<_> = (0..rows) + .step_by(per) + .map(|first| { + let count = per.min(rows - first); + let data = &bytes[first * 4 * stride..(first + count) * 4 * stride]; + scope.spawn(move || { + let surface = intel_tex_2::RgbaSurface { + data, + width: padded, + height: 4 * count as u32, + stride: stride as u32, + }; + intel_tex_2::bc6h::compress_blocks(&settings, &surface) + }) + }) + .collect(); + bands + .into_iter() + .flat_map(|b| b.join().expect("compression thread")) + .collect() + }) +} + +/// Builds `out` from the month's eight tiles in `source`. +pub fn build(source: &Path, month: u32, out: &PathBuf) -> Result<(), String> { + let started = Instant::now(); + let layout = DetailLayout::BLUE_MARBLE; + let format = DetailFormat::Bc6h; + let per_source = SOURCE / layout.tile; + let decode_tile = |column: u32, row: u32| decode(&source.join(source_name(month, column, row))); + + // First pass: the edges of every tile, for the borders their neighbors need. + let mut edges = Vec::new(); + for row in 0..2 { + for column in 0..4 { + eprintln!("reading the edges of {}", source_name(month, column, row)); + edges.push(Edges::new(&decode_tile(column, row)?, layout.gutter)); + } + } + + let fail = |e: std::io::Error| format!("{}: {e}", out.display()); + let mut file = std::fs::File::create(out).map_err(fail)?; + EarthDetail::write_header(&mut file, layout, format).map_err(fail)?; + let tile_bytes: u64 = (0..layout.levels) + .map(|l| format.bytes(layout.padded(l))) + .sum(); + file.set_len(64 + tile_bytes * u64::from(layout.layers())) + .map_err(fail)?; + let decode_lut: Vec = (0..=255u8) + .map(|v| srgb_to_linear(f32::from(v) / 255.0)) + .collect(); + + // Second pass: each source tile's output tiles. + for row in 0..2 { + for column in 0..4 { + let name = source_name(month, column, row); + eprintln!("tiling {name}"); + let rgb = decode_tile(column, row)?; + let texel = |x: u32, y: u32| { + let (sc, sr) = (x / SOURCE, y / SOURCE); + let (lx, ly) = (x % SOURCE, y % SOURCE); + let v = if (sc, sr) == (column, row) { + let i = 3 * (ly as usize * SOURCE as usize + lx as usize); + [rgb[i], rgb[i + 1], rgb[i + 2]] + } else { + edges[(sr * 4 + sc) as usize].at(lx, ly, layout.gutter) + }; + v.map(|c| decode_lut[c as usize]) + }; + for dy in 0..per_source { + for dx in 0..per_source { + let (tx, ty) = (column * per_source + dx, row * per_source + dy); + let levels = EarthDetail::build_tile(&layout, tx, ty, texel); + let index = u64::from(ty * layout.tiles[0] + tx); + file.seek(SeekFrom::Start(64 + index * tile_bytes)) + .map_err(fail)?; + for (l, level) in levels.iter().enumerate() { + file.write_all(&compress(level, layout.padded(l as u32))) + .map_err(fail)?; + } + } + } + } + } + file.sync_all().map_err(fail)?; + eprintln!( + "wrote {} in {:.0} s", + out.display(), + started.elapsed().as_secs_f64() + ); + Ok(()) +} diff --git a/crates/resolved-cli/src/main.rs b/crates/resolved-cli/src/main.rs index 94b3ab9..976365c 100644 --- a/crates/resolved-cli/src/main.rs +++ b/crates/resolved-cli/src/main.rs @@ -1,5 +1,6 @@ //! The `resolved` command: render a scene file to PNG without writing a host program. +mod earth_tiles; mod scene_file; use std::path::{Path, PathBuf}; @@ -10,8 +11,8 @@ use std::time::Instant; use clap::{Args, Parser, Subcommand}; use resolved::gpu::AdapterChoice; use resolved::{ - CameraTier, EarthTexture, Fidelity, Frame, LightingTier, Look, MoonTexture, Preview, Renderer, - RendererOptions, StarCatalog, StarPhotometry, + CameraTier, EarthDetail, EarthTexture, Fidelity, Frame, LightingTier, Look, MoonTexture, + Preview, Renderer, RendererOptions, StarCatalog, StarPhotometry, }; use crate::scene_file::SceneFile; @@ -32,6 +33,19 @@ enum Command { Render(RenderArgs), /// List the adapters this machine offers. Adapters, + /// Make the sharp tiled Earth map for a month from Blue Marble's eight 21600 x 21600 + /// JPEG tiles (scripts/fetch-earth.sh --sharp runs this): about 4.9 GB of BC6H. + EarthTiles { + /// Directory with world.2004MM.3x21600x21600.{A,B,C,D}{1,2}.jpg. + #[arg(long)] + source: PathBuf, + /// Month, 1 to 12. + #[arg(long)] + month: u32, + /// Where to write the map. + #[arg(long)] + out: PathBuf, + }, } #[derive(Args)] @@ -114,6 +128,7 @@ fn main() -> ExitCode { Ok(()) } Command::Render(args) => render(&args), + Command::EarthTiles { source, month, out } => earth_tiles::build(&source, month, &out), }; match result { Ok(()) => ExitCode::SUCCESS, @@ -180,6 +195,28 @@ fn render(args: &RenderArgs) -> Result<(), String> { // What the tiers need, gathered before the renderer exists so it can check them. let start = Instant::now(); + // The sharp tiled map, if one has been made for the month (scripts/fetch-earth.sh + // --sharp): optional, so its absence is no error. + let earth_detail = if args.lighting >= LightingTier::L1 { + match file.earth_month()? { + Some(month) => { + let name = format!("earth/sharp-2004-{month:02}.tiles"); + match data_directories(args.data.as_deref()) + .into_iter() + .map(|dir| dir.join(&name)) + .find(|p| p.is_file()) + { + Some(path) => Some(Arc::new( + EarthDetail::open(&path).map_err(|e| e.to_string())?, + )), + None => None, + } + } + None => None, + } + } else { + None + }; let earth_texture = if args.lighting >= LightingTier::L1 { let month = file.earth_month()?.ok_or( "lighting L1 needs the month for Earth's surface: give the scene an epoch, or [appearance] earth_month", @@ -234,6 +271,7 @@ fn render(args: &RenderArgs) -> Result<(), String> { camera: args.camera, }, earth_texture, + earth_detail: earth_detail.clone(), moon_texture, appearance: file.appearance(), read_radiance: false, @@ -278,8 +316,13 @@ fn render(args: &RenderArgs) -> Result<(), String> { setup.as_secs_f64() * 1e3 ); println!( - "maps {:.1} ms (Earth and Moon image decode)", - texturing.as_secs_f64() * 1e3 + "maps {:.1} ms (Earth and Moon image decode){}", + texturing.as_secs_f64() * 1e3, + if earth_detail.is_some() { + "; Earth's sharp tiled map in use" + } else { + "" + } ); println!( "scene {:.1} ms (models and ephemeris)", diff --git a/crates/resolved-showcase/src/main.rs b/crates/resolved-showcase/src/main.rs index 1229d36..e6d4c6e 100644 --- a/crates/resolved-showcase/src/main.rs +++ b/crates/resolved-showcase/src/main.rs @@ -27,9 +27,9 @@ use resolved::ephemeris::Ephemeris; use resolved::gpu::AdapterChoice; use resolved::nalgebra::Vector3; use resolved::{ - CameraConfig, CameraTier, Detector, EarthTexture, Fidelity, LightingTier, LoadOptions, Look, - Mesh, MoonTexture, Preview, Renderer, RendererOptions, Scene, StarCatalog, StarPhotometry, - Target, + CameraConfig, CameraTier, Detector, EarthDetail, EarthTexture, Fidelity, LightingTier, + LoadOptions, Look, Mesh, MoonTexture, Preview, Renderer, RendererOptions, Scene, StarCatalog, + StarPhotometry, Target, }; type Result = std::result::Result; @@ -72,7 +72,8 @@ struct Cli { #[arg(long, default_value = "auto")] adapter: String, /// Directory with earth/bmng-2004-MM.jpg, moon/wac-hapke-643.png and - /// stars/stars.bin. Without it, the repository's data/. + /// stars/stars.bin, and optionally earth/sharp-2004-MM.tiles. Without it, the + /// repository's data/. #[arg(long, value_name = "DIR")] data: Option, /// Directory of astrokit's kernels. Without it: ASTROKIT_KERNELS, then @@ -178,6 +179,14 @@ impl Movie { .map_err(|e| format!("{e} (scripts/fetch-earth.sh downloads it)"))?; let moon = MoonTexture::load(&moon) .map_err(|e| format!("{e} (scripts/fetch-moon.py makes it)"))?; + // The sharp tiled map, if one has been made for the month (scripts/fetch-earth.sh + // --sharp); optional. + let sharp = data.join(format!("earth/sharp-2004-{month:02}.tiles")); + let earth_detail = if sharp.is_file() { + Some(Arc::new(EarthDetail::open(&sharp).map_err(error)?)) + } else { + None + }; let mut camera = CameraConfig::example(); camera.exposure = cli.exposure; @@ -213,6 +222,7 @@ impl Movie { camera: CameraTier::C1, }, earth_texture: Some(Arc::new(earth)), + earth_detail, moon_texture: Some(Arc::new(moon)), stars: Some(Arc::new(stars)), ..Default::default() diff --git a/crates/resolved/Cargo.toml b/crates/resolved/Cargo.toml index 93e5545..4222fed 100644 --- a/crates/resolved/Cargo.toml +++ b/crates/resolved/Cargo.toml @@ -38,6 +38,8 @@ stl_io = "0.11" approx = "0.5" image = { version = "0.25", default-features = false, features = ["png"] } criterion = "0.8" +# BC6H compression, to test the sharp Earth map in the format it ships in. +intel_tex_2 = "0.5" [[bench]] name = "render" diff --git a/crates/resolved/src/earth.rs b/crates/resolved/src/earth.rs index f14b919..21578e8 100644 --- a/crates/resolved/src/earth.rs +++ b/crates/resolved/src/earth.rs @@ -139,7 +139,318 @@ pub(crate) fn parallel_rows( }); } -fn srgb_to_linear(v: f32) -> f32 { +/// How a sharp map of Earth's surface is cut into square tiles, for GPUs that take no +/// texture as large as the whole map (lavapipe's limit is 16384 on a side). +/// +/// The map is equirectangular like [`EarthTexture`], `tiles[0] * tile` texels across +/// and `tiles[1] * tile` down. Each tile stores its interior and a border of `gutter` +/// texels on every side copied from its neighbors (wrapping across the 180 degree +/// meridian, repeating the edge row at the poles), so filtering near a tile's edge +/// reads the right texels without consulting another tile. Each level of the mip chain +/// halves the tile and its border. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub struct DetailLayout { + /// Tiles across (in longitude) and down (in latitude). + pub tiles: [u32; 2], + /// Texels across a tile's interior, at full size. + pub tile: u32, + /// Border texels on each side of a tile at full size. + pub gutter: u32, + /// Mip levels stored, from full size. + pub levels: u32, +} + +impl DetailLayout { + /// Blue Marble Next Generation's 86400 x 43200 map (15 arcseconds, about 460 m at + /// the equator) as 8 x 4 tiles of 10800, with levels at 460 m, 930 m and 1.85 km: + /// coarser than that, the global [`EarthTexture`] is as sharp. + pub const BLUE_MARBLE: Self = Self { + tiles: [8, 4], + tile: 10800, + gutter: 4, + levels: 3, + }; + + /// Texels across and down the whole map at full size. + pub fn size(&self) -> [u32; 2] { + [self.tiles[0] * self.tile, self.tiles[1] * self.tile] + } + + /// The number of tiles. + pub fn layers(&self) -> u32 { + self.tiles[0] * self.tiles[1] + } + + /// Texels across a stored tile, border included, at mip `level`. + pub fn padded(&self, level: u32) -> u32 { + (self.tile + 2 * self.gutter) >> level + } + + /// Checks that every level halves exactly and keeps a border of at least one texel, + /// and for BC6H that the full-size tile is a whole number of 4 x 4 blocks. + pub fn validate(&self, format: DetailFormat) -> Result<(), Error> { + let halving = 1 << (self.levels.max(1) - 1); + let ok = self.tiles[0] > 0 + && self.tiles[1] > 0 + && self.levels > 0 + && self.tile % halving == 0 + && self.gutter % halving == 0 + && self.gutter >= halving + && (format == DetailFormat::Rgba16Float || self.padded(0) % 4 == 0); + if ok { + Ok(()) + } else { + Err(Error::Texture(format!( + "unusable tile layout {self:?} for {format:?}" + ))) + } + } +} + +/// How a sharp map's texels are stored. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum DetailFormat { + /// Linear reflectance as RGBA f16, 8 bytes a texel: for tests and small maps. + Rgba16Float, + /// BC6H, unsigned: blocks of 4 x 4 texels in 16 bytes, 1 byte a texel, keeping + /// linear values with about f16's precision, decoded by the GPU as it reads them. + Bc6h, +} + +impl DetailFormat { + /// Bytes of one tile of `size` x `size` texels. + pub fn bytes(&self, size: u32) -> u64 { + let size = u64::from(size); + match self { + Self::Rgba16Float => size * size * 8, + Self::Bc6h => size.div_ceil(4) * size.div_ceil(4) * 16, + } + } + + fn code(&self) -> u32 { + match self { + Self::Rgba16Float => 0, + Self::Bc6h => 1, + } + } +} + +/// A sharp map of Earth's surface, tiled ([`DetailLayout`]): drawn where a view +/// resolves ground finer than the global [`EarthTexture`] does. Made from Blue Marble's +/// 86400 x 43200 images by `resolved earth-tiles` (`scripts/fetch-earth.sh --sharp`). +/// +/// A file holds a 64-byte header, then each tile's mip levels in turn, tiles in rows +/// from the north-west. [`EarthDetail::open`] reads only the header; the tiles are +/// read as they are uploaded, a tile at a time, so the 4.9 GB of Blue Marble's map is +/// never held whole in memory. +#[derive(Clone, Debug)] +pub struct EarthDetail { + layout: DetailLayout, + format: DetailFormat, + data: DetailData, +} + +#[derive(Clone, Debug)] +enum DetailData { + Memory(std::sync::Arc>), + File(std::path::PathBuf), +} + +const DETAIL_MAGIC: &[u8; 8] = b"RSLVTILE"; +const DETAIL_HEADER: u64 = 64; + +impl EarthDetail { + /// Opens a file written by [`EarthDetail::write_header`] and the tiles after it, + /// checking its header and length. + pub fn open(path: &Path) -> Result { + use std::io::Read; + let fail = |what: String| Error::Texture(format!("{}: {what}", path.display())); + let mut file = std::fs::File::open(path).map_err(|e| fail(e.to_string()))?; + let mut header = [0u8; DETAIL_HEADER as usize]; + file.read_exact(&mut header) + .map_err(|e| fail(format!("header: {e}")))?; + if &header[..8] != DETAIL_MAGIC { + return Err(fail("not a tiled Earth map".into())); + } + let word = |i: usize| u32::from_le_bytes(header[8 + 4 * i..12 + 4 * i].try_into().unwrap()); + if word(0) != 1 { + return Err(fail(format!("version {} of the format, not 1", word(0)))); + } + let format = match word(1) { + 0 => DetailFormat::Rgba16Float, + 1 => DetailFormat::Bc6h, + n => return Err(fail(format!("unknown texel format {n}"))), + }; + let layout = DetailLayout { + tiles: [word(2), word(3)], + tile: word(4), + gutter: word(5), + levels: word(6), + }; + layout.validate(format)?; + let detail = Self { + layout, + format, + data: DetailData::File(path.to_owned()), + }; + let want = DETAIL_HEADER + detail.tile_bytes() * u64::from(layout.layers()); + let have = file.metadata().map_err(|e| fail(e.to_string()))?.len(); + if have != want { + return Err(fail(format!( + "{have} bytes, not the {want} its header implies" + ))); + } + Ok(detail) + } + + /// A map held in memory: each tile's levels in turn, as in a file after its header. + pub fn from_bytes( + layout: DetailLayout, + format: DetailFormat, + bytes: Vec, + ) -> Result { + layout.validate(format)?; + let detail = Self { + layout, + format, + data: DetailData::Memory(std::sync::Arc::new(bytes)), + }; + let want = detail.tile_bytes() * u64::from(layout.layers()); + let have = match &detail.data { + DetailData::Memory(b) => b.len() as u64, + DetailData::File(_) => unreachable!(), + }; + if have != want { + return Err(Error::Texture(format!( + "{have} bytes of tiles, not the {want} the layout needs" + ))); + } + Ok(detail) + } + + /// Writes the header of a file holding a map of `layout` and `format`; the tiles + /// follow, each as [`EarthDetail::build_tile`]'s levels, encoded. + pub fn write_header( + out: &mut impl std::io::Write, + layout: DetailLayout, + format: DetailFormat, + ) -> std::io::Result<()> { + let mut header = [0u8; DETAIL_HEADER as usize]; + header[..8].copy_from_slice(DETAIL_MAGIC); + let words = [ + 1, + format.code(), + layout.tiles[0], + layout.tiles[1], + layout.tile, + layout.gutter, + layout.levels, + ]; + for (i, w) in words.iter().enumerate() { + header[8 + 4 * i..12 + 4 * i].copy_from_slice(&w.to_le_bytes()); + } + out.write_all(&header) + } + + /// The tile layout. + pub fn layout(&self) -> DetailLayout { + self.layout + } + + /// The texel format. + pub fn format(&self) -> DetailFormat { + self.format + } + + /// Bytes of one tile, all its levels. + pub fn tile_bytes(&self) -> u64 { + (0..self.layout.levels) + .map(|l| self.format.bytes(self.layout.padded(l))) + .sum() + } + + /// Tile `index`'s levels, full size first, each as its texels or blocks in rows. + pub(crate) fn tile_levels(&self, index: u32) -> Result>, Error> { + use std::io::{Read, Seek, SeekFrom}; + let size = self.tile_bytes(); + let start = u64::from(index) * size; + let bytes = match &self.data { + DetailData::Memory(b) => b[start as usize..(start + size) as usize].to_vec(), + DetailData::File(path) => { + let fail = |e: std::io::Error| Error::Texture(format!("{}: {e}", path.display())); + let mut file = std::fs::File::open(path).map_err(fail)?; + file.seek(SeekFrom::Start(DETAIL_HEADER + start)) + .map_err(fail)?; + let mut b = vec![0; size as usize]; + file.read_exact(&mut b).map_err(fail)?; + b + } + }; + let mut levels = Vec::new(); + let mut at = 0; + for l in 0..self.layout.levels { + let n = self.format.bytes(self.layout.padded(l)) as usize; + levels.push(bytes[at..at + n].to_vec()); + at += n; + } + Ok(levels) + } + + /// Tile (`x`, `y`) of `layout`, from the north-west, as linear RGBA f16 at each mip + /// level: `padded(level)` x `padded(level)` texels in rows, border included. Texels + /// come from `texel(column, row)` of the whole map at full size, with the column + /// wrapped and the row clamped to the map. Each level averages 2 x 2 texels of the + /// one above, in linear light. + pub fn build_tile( + layout: &DetailLayout, + x: u32, + y: u32, + texel: impl Fn(u32, u32) -> [f32; 3] + Sync, + ) -> Vec> { + let [width, height] = layout.size(); + let size = layout.padded(0) as usize; + let one = f16::from_f32(1.0); + let mut first = vec![one; 4 * size * size]; + let (x0, y0) = ( + i64::from(x * layout.tile) - i64::from(layout.gutter), + i64::from(y * layout.tile) - i64::from(layout.gutter), + ); + parallel_rows(&mut first, 4 * size, |row, out| { + let gy = (y0 + row as i64).clamp(0, i64::from(height) - 1) as u32; + for (col, t) in out.chunks_exact_mut(4).enumerate() { + let gx = (x0 + col as i64).rem_euclid(i64::from(width)) as u32; + let rgb = texel(gx, gy); + for c in 0..3 { + t[c] = f16::from_f32(rgb[c]); + } + } + }); + let mut levels = vec![first]; + for l in 1..layout.levels { + let above = levels.last().expect("at least one level"); + let (from, to) = (layout.padded(l - 1) as usize, layout.padded(l) as usize); + let mut below = vec![one; 4 * to * to]; + parallel_rows(&mut below, 4 * to, |row, out| { + for col in 0..to { + for c in 0..3 { + let sum: f32 = [(0, 0), (1, 0), (0, 1), (1, 1)] + .iter() + .map(|(dx, dy)| { + above[4 * ((2 * row + dy) * from + 2 * col + dx) + c].to_f32() + }) + .sum(); + out[4 * col + c] = f16::from_f32(sum / 4.0); + } + } + }); + levels.push(below); + } + levels + } +} + +/// Linear light from an sRGB-encoded value in [0, 1]. +pub fn srgb_to_linear(v: f32) -> f32 { if v <= 0.040_45 { v / 12.92 } else { @@ -188,6 +499,77 @@ mod tests { .all(|(_, _, p)| p.chunks(4).all(|t| t[..3] == [first; 3]))); } + #[test] + fn tiles_carry_their_neighbors_borders() { + let layout = DetailLayout { + tiles: [4, 2], + tile: 8, + gutter: 2, + levels: 2, + }; + // Each texel holds its own column and row. + let texel = |x: u32, y: u32| [x as f32, y as f32, 0.0]; + let size = layout.padded(0) as usize; + let at = |levels: &Vec>, col: usize, row: usize| { + let t = &levels[0][4 * (row * size + col)..]; + (t[0].to_f32(), t[1].to_f32()) + }; + // The first tile: its border wraps west to the map's last columns and repeats + // the top row above. + let first = EarthDetail::build_tile(&layout, 0, 0, texel); + assert_eq!(at(&first, 0, 0), (30.0, 0.0)); + assert_eq!(at(&first, 2, 2), (0.0, 0.0)); + assert_eq!(at(&first, 11, 11), (9.0, 9.0)); + // The last tile: its border wraps east to column 0 and repeats the bottom row. + let last = EarthDetail::build_tile(&layout, 3, 1, texel); + assert_eq!(at(&last, 11, 11), (1.0, 15.0)); + // The next level halves the tile and its border; it averages 2 x 2 texels. + assert_eq!(first[1].len(), 4 * 6 * 6); + assert_eq!(first[1][4 * (6 + 1)].to_f32(), 0.5); + } + + #[test] + fn tiled_maps_round_trip_through_a_file() { + let layout = DetailLayout { + tiles: [2, 1], + tile: 4, + gutter: 2, + levels: 2, + }; + let format = DetailFormat::Rgba16Float; + let mut bytes = Vec::new(); + for x in 0..2 { + for level in EarthDetail::build_tile(&layout, x, 0, |c, _| [c as f32; 3]) { + bytes.extend_from_slice(bytemuck::cast_slice(&level)); + } + } + let memory = EarthDetail::from_bytes(layout, format, bytes.clone()).unwrap(); + let path = std::env::temp_dir().join(format!("resolved-tiles-{}.bin", std::process::id())); + let mut file = std::fs::File::create(&path).unwrap(); + EarthDetail::write_header(&mut file, layout, format).unwrap(); + std::io::Write::write_all(&mut file, &bytes).unwrap(); + drop(file); + let opened = EarthDetail::open(&path).unwrap(); + assert_eq!(opened.layout(), layout); + assert_eq!( + opened.tile_levels(1).unwrap(), + memory.tile_levels(1).unwrap() + ); + // A truncated file is refused. + std::fs::OpenOptions::new() + .write(true) + .open(&path) + .unwrap() + .set_len(100) + .unwrap(); + assert!(EarthDetail::open(&path).is_err()); + std::fs::remove_file(&path).unwrap(); + assert_eq!(DetailLayout::BLUE_MARBLE.size(), [86400, 43200]); + assert!(DetailLayout::BLUE_MARBLE + .validate(DetailFormat::Bc6h) + .is_ok()); + } + #[test] fn wrong_sizes_are_rejected() { assert!(EarthTexture::from_rgba8(4, 2, vec![0; 31]).is_err()); diff --git a/crates/resolved/src/gpu/bodies.rs b/crates/resolved/src/gpu/bodies.rs index 3cbe7e0..4c734ed 100644 --- a/crates/resolved/src/gpu/bodies.rs +++ b/crates/resolved/src/gpu/bodies.rs @@ -15,7 +15,7 @@ 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; +use crate::earth::{DetailFormat, DetailLayout, EarthDetail, EarthTexture}; use crate::model::rgb_spectrum; use crate::moon::{MoonTexture, NORMALIZED_INCIDENCE, PHASE_SCALE}; use crate::radiometry::{Radiometry, AU, SUN_RADIUS}; @@ -23,6 +23,11 @@ use crate::scene::{Pose, Scene}; use crate::sun::{ring_rule, LimbDarkening, RINGS, SAMPLES}; use crate::Error; +/// The sharp map's tiles are split among this many texture arrays (`DETAIL_GROUPS` in +/// bodies.wgsl): lavapipe allocates no more than 2 GB at once, and Blue Marble's map +/// is 4.9 GB, so each array holds an eighth to a quarter of it. +pub const DETAIL_GROUPS: usize = 4; + /// Channels one radiance target holds. pub const MAX_CHANNELS: usize = 4; @@ -84,6 +89,8 @@ pub struct BodyRadiometry { pub atmosphere: bool, /// Whether Earth is drawn from a texture. pub earth_textured: bool, + /// The sharp tiled map's layout, once one is set. + pub earth_detail: Option, /// Whether Earth and the Moon shadow each other: lunar and solar eclipses. pub bodies_eclipsed: bool, /// Per channel, Earth's albedo from `Appearance::earth_albedo` (no texture). @@ -142,6 +149,7 @@ impl BodyRadiometry { earth_irradiance_1au: lambert(vec![1.0; n]), earth_rgb: [basis(0), basis(1), basis(2)], atmosphere: false, + earth_detail: None, earth_textured: false, bodies_eclipsed: false, earth_albedo: albedo(&appearance.earth_albedo), @@ -186,6 +194,11 @@ pub struct BodiesUniforms { // Per body, Earth, Moon, Sun: the camera's velocity relative to it over c, in the // camera frame. beta: [[f32; 4]; 3], + // The sharp tiled map: x 1 if set, y tiles across, z tiles down, w mip levels. + earth_detail: [f32; 4], + // x: a tile's interior over its stored size, y: its border over its stored size, + // z, w: the whole map's texels across and down. + earth_detail_size: [f32; 4], } impl BodiesUniforms { @@ -250,6 +263,19 @@ impl BodiesUniforms { moon_irradiance: scaled(radiometry.earth_irradiance_1au, moon_distance), moon_color: radiometry.moon_color, moon_params: [f32::from(u8::from(radiometry.moon_textured)), 0.0, 0.0, 0.0], + earth_detail: radiometry.earth_detail.map_or([0.0; 4], |l| { + [1.0, l.tiles[0] as f32, l.tiles[1] as f32, l.levels as f32] + }), + earth_detail_size: radiometry.earth_detail.map_or([0.0; 4], |l| { + let padded = l.padded(0) as f32; + let [w, h] = l.size(); + [ + l.tile as f32 / padded, + l.gutter as f32 / padded, + w as f32, + h as f32, + ] + }), beta: scene.velocities.relative().map(|b| { let b = if apparent { to_camera * b @@ -286,6 +312,8 @@ pub struct BodiesPass { moon: wgpu::TextureView, moon_height: u32, atmosphere: AtmosphereTextures, + detail: [wgpu::TextureView; DETAIL_GROUPS], + detail_sampler: wgpu::Sampler, bind_group: wgpu::BindGroup, } @@ -394,6 +422,15 @@ impl BodiesPass { let earth_view = upload_earth(ctx, &earth); let moon_view = upload_moon(ctx, &moon); let atmosphere = AtmosphereTextures::placeholder(ctx); + // Clamped: each tile carries its neighbors' texels in its border. + let detail_sampler = device.create_sampler(&wgpu::SamplerDescriptor { + label: Some("earth detail"), + mag_filter: wgpu::FilterMode::Linear, + min_filter: wgpu::FilterMode::Linear, + mipmap_filter: wgpu::MipmapFilterMode::Linear, + ..Default::default() + }); + let detail = std::array::from_fn(|_| placeholder_detail(ctx)); let bind_group = bind( ctx, &pipeline, @@ -403,6 +440,7 @@ impl BodiesPass { &earth_view, &moon_view, &atmosphere, + (&detail, &detail_sampler), ); Self { pipeline, @@ -414,6 +452,8 @@ impl BodiesPass { moon: moon_view, moon_height: moon.height(), atmosphere, + detail, + detail_sampler, bind_group, } } @@ -432,6 +472,106 @@ impl BodiesPass { self.rebind(ctx); } + /// Uploads the sharp tiled map of Earth's surface, a tile at a time, and binds it + /// for later frames. A BC6H map needs a device with BC compression. + pub fn set_earth_detail( + &mut self, + ctx: &GpuContext, + detail: &EarthDetail, + ) -> Result<(), Error> { + let layout = detail.layout(); + let format = match detail.format() { + DetailFormat::Rgba16Float => wgpu::TextureFormat::Rgba16Float, + DetailFormat::Bc6h => { + if !ctx + .device + .features() + .contains(wgpu::Features::TEXTURE_COMPRESSION_BC) + { + return Err(Error::Unsupported(format!( + "{} has no BC texture compression, which the sharp Earth map needs", + ctx.info.name + ))); + } + wgpu::TextureFormat::Bc6hRgbUfloat + } + }; + let size = layout.padded(0); + let per_group = layout.layers().div_ceil(DETAIL_GROUPS as u32); + let mut views = Vec::new(); + for group in 0..DETAIL_GROUPS as u32 { + let first = group * per_group; + let count = per_group.min(layout.layers().saturating_sub(first)).max(1); + let texture = ctx.device.create_texture(&wgpu::TextureDescriptor { + label: Some("earth detail"), + size: wgpu::Extent3d { + width: size, + height: size, + depth_or_array_layers: count, + }, + mip_level_count: layout.levels, + sample_count: 1, + dimension: wgpu::TextureDimension::D2, + format, + usage: wgpu::TextureUsages::TEXTURE_BINDING | wgpu::TextureUsages::COPY_DST, + view_formats: &[], + }); + for local in 0..count { + let layer = first + local; + if layer >= layout.layers() { + break; + } + for (level, bytes) in detail.tile_levels(layer)?.iter().enumerate() { + let s = layout.padded(level as u32); + // A compressed level is written in whole blocks, past its edge if + // need be. + let (width, rows, row_bytes) = match detail.format() { + DetailFormat::Rgba16Float => (s, s, 8 * s), + DetailFormat::Bc6h => { + (s.div_ceil(4) * 4, s.div_ceil(4), 16 * s.div_ceil(4)) + } + }; + ctx.queue.write_texture( + wgpu::TexelCopyTextureInfo { + texture: &texture, + mip_level: level as u32, + origin: wgpu::Origin3d { + x: 0, + y: 0, + z: local, + }, + aspect: wgpu::TextureAspect::All, + }, + bytes, + wgpu::TexelCopyBufferLayout { + offset: 0, + bytes_per_row: Some(row_bytes), + rows_per_image: Some(rows), + }, + wgpu::Extent3d { + width, + height: width, + depth_or_array_layers: 1, + }, + ); + } + // Hand each tile to the device before reading the next, so the staging + // copies do not pile up: the whole map is 4.9 GB. + ctx.queue.submit([]); + ctx.device + .poll(wgpu::PollType::wait_indefinitely()) + .map_err(|e| Error::Gpu(e.to_string()))?; + } + views.push(texture.create_view(&wgpu::TextureViewDescriptor { + dimension: Some(wgpu::TextureViewDimension::D2Array), + ..Default::default() + })); + } + self.detail = views.try_into().expect("one view per group"); + self.rebind(ctx); + Ok(()) + } + /// 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) { @@ -454,6 +594,7 @@ impl BodiesPass { &self.earth, &self.moon, &self.atmosphere, + (&self.detail, &self.detail_sampler), ); } @@ -538,6 +679,29 @@ impl BodiesPass { } } +/// A one-texel, one-layer stand-in for the sharp map until one is set; the shader does +/// not read it then. +fn placeholder_detail(ctx: &GpuContext) -> wgpu::TextureView { + let texture = ctx.device.create_texture(&wgpu::TextureDescriptor { + label: Some("earth detail placeholder"), + size: wgpu::Extent3d { + width: 1, + height: 1, + depth_or_array_layers: 1, + }, + mip_level_count: 1, + sample_count: 1, + dimension: wgpu::TextureDimension::D2, + format: wgpu::TextureFormat::Rgba16Float, + usage: wgpu::TextureUsages::TEXTURE_BINDING, + view_formats: &[], + }); + texture.create_view(&wgpu::TextureViewDescriptor { + dimension: Some(wgpu::TextureViewDimension::D2Array), + ..Default::default() + }) +} + /// Uploads Earth's texture with its mip chain as linear RGBA f16. fn upload_earth(ctx: &GpuContext, texture: &EarthTexture) -> wgpu::TextureView { let levels = texture.mip_chain(); @@ -625,6 +789,7 @@ fn bind( earth: &wgpu::TextureView, moon: &wgpu::TextureView, atmosphere: &AtmosphereTextures, + detail: (&[wgpu::TextureView; DETAIL_GROUPS], &wgpu::Sampler), ) -> wgpu::BindGroup { let entries = [ wgpu::BindGroupEntry { @@ -648,7 +813,24 @@ fn bind( resource: wgpu::BindingResource::TextureView(moon), }, ]; - let entries: Vec<_> = entries.into_iter().chain(atmosphere.entries(5)).collect(); + let entries: Vec<_> = entries + .into_iter() + .chain(atmosphere.entries(5)) + .chain( + detail + .0 + .iter() + .enumerate() + .map(|(i, view)| wgpu::BindGroupEntry { + binding: 10 + i as u32, + resource: wgpu::BindingResource::TextureView(view), + }), + ) + .chain([wgpu::BindGroupEntry { + binding: 10 + DETAIL_GROUPS as u32, + resource: wgpu::BindingResource::Sampler(detail.1), + }]) + .collect(); ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { label: Some("bodies"), layout: &pipeline.get_bind_group_layout(0), diff --git a/crates/resolved/src/gpu/context.rs b/crates/resolved/src/gpu/context.rs index febad82..2f9d0ad 100644 --- a/crates/resolved/src/gpu/context.rs +++ b/crates/resolved/src/gpu/context.rs @@ -86,11 +86,13 @@ impl GpuContext { label: Some("resolved"), // Timestamps for profiling; mappable storage so the DN image can be read // where the sensor pass writes it, without a copy (see `FrameTargets`); - // the adapter's own format support, for multisampled f32 radiance. + // the adapter's own format support, for multisampled f32 radiance; BC + // compression, for the sharp Earth map (`EarthDetail`). required_features: adapter.features() & (wgpu::Features::TIMESTAMP_QUERY | wgpu::Features::MAPPABLE_PRIMARY_BUFFERS - | wgpu::Features::TEXTURE_ADAPTER_SPECIFIC_FORMAT_FEATURES), + | wgpu::Features::TEXTURE_ADAPTER_SPECIFIC_FORMAT_FEATURES + | wgpu::Features::TEXTURE_COMPRESSION_BC), required_limits: adapter.limits(), ..Default::default() }))?; diff --git a/crates/resolved/src/lib.rs b/crates/resolved/src/lib.rs index 93592ea..aa7f1e7 100644 --- a/crates/resolved/src/lib.rs +++ b/crates/resolved/src/lib.rs @@ -44,7 +44,7 @@ pub mod sun; pub use camera::{ Blooming, BrownConrady, CameraConfig, Channel, Detector, Intrinsics, Psf, PsfImage, }; -pub use earth::EarthTexture; +pub use earth::{EarthDetail, EarthTexture}; pub use error::Error; pub use fidelity::{CameraTier, Fidelity, LightingTier}; pub use fixed_pattern::FixedPattern; diff --git a/crates/resolved/src/render.rs b/crates/resolved/src/render.rs index 7dba0ba..6222090 100644 --- a/crates/resolved/src/render.rs +++ b/crates/resolved/src/render.rs @@ -16,7 +16,7 @@ use crate::camera::{Blooming, CameraConfig, Intrinsics, Psf}; use std::sync::Arc; use crate::atmosphere::AtmosphereTables; -use crate::earth::EarthTexture; +use crate::earth::{EarthDetail, EarthTexture}; use crate::eclipse::{sun_visible_fraction, sun_visible_light}; use crate::fidelity::{CameraTier, Fidelity, LightingTier}; use crate::fixed_pattern::FixedPattern; @@ -52,6 +52,10 @@ pub struct RendererOptions { pub fidelity: Fidelity, /// Earth's surface, needed from lighting tier L1 up; ignored below it. pub earth_texture: Option>, + /// A sharp tiled map of Earth's surface, drawn from lighting tier L1 where a view + /// resolves ground finer than `earth_texture` does; optional. Uploading Blue + /// Marble's takes about 4.9 GB of device memory. + pub earth_detail: Option>, /// The Moon's surface, needed from lighting tier L1 up; ignored below it. pub moon_texture: Option>, /// How Earth and the Moon reflect light. @@ -441,6 +445,9 @@ impl Renderer { if let Some(texture) = &options.earth_texture { renderer.set_earth_texture(texture); } + if let Some(detail) = &options.earth_detail { + renderer.set_earth_detail(detail)?; + } if let Some(texture) = &options.moon_texture { renderer.set_moon_texture(texture); } @@ -475,6 +482,25 @@ impl Renderer { self.body_radiometry.earth_textured = true; } + /// Draws Earth's surface from the sharp tiled map `detail` from now on where a view + /// resolves ground finer than the Earth texture does, replacing any given in + /// [`RendererOptions::earth_detail`]: for switching months. Uploads it a tile at a + /// time from its file. A BC6H map needs a device with BC compression. + /// + /// Below lighting tier L1 there is no Earth texture, and this does nothing. + pub fn set_earth_detail(&mut self, detail: &EarthDetail) -> Result<(), Error> { + if self.fidelity.lighting < LightingTier::L1 { + log::warn!( + "lighting {} draws no Earth texture; ignoring the sharp map", + self.fidelity.lighting + ); + return Ok(()); + } + self.bodies.set_earth_detail(&self.ctx, detail)?; + self.body_radiometry.earth_detail = Some(detail.layout()); + Ok(()) + } + /// Draws the Moon's surface from `texture` from now on, replacing the one given in /// [`RendererOptions::moon_texture`]. Uploads the map with its mip chain, about /// 360 MB of GPU memory for 16384 x 8192. diff --git a/crates/resolved/src/shaders/bodies.wgsl b/crates/resolved/src/shaders/bodies.wgsl index 2299067..1acbbac 100644 --- a/crates/resolved/src/shaders/bodies.wgsl +++ b/crates/resolved/src/shaders/bodies.wgsl @@ -56,6 +56,12 @@ struct Uniforms { // Per body, Earth, Moon, Sun: the camera's velocity relative to it over c, in the // camera frame; zero below camera C1. beta: array, 3>, + // The sharp tiled map (earth.rs, `DetailLayout`): x 1 if set, y tiles across, z + // tiles down, w mip levels. + earth_detail: vec4, + // x: a tile's interior over its stored size, y: its border over its stored size, + // z, w: the whole map's texels across and down. + earth_detail_size: vec4, } @group(0) @binding(0) var u: Uniforms; @@ -66,6 +72,16 @@ struct Uniforms { // Equirectangular 643 nm reflectance of the Moon, as Earth's. @group(0) @binding(4) var moon_texture: texture_2d; +// The sharp map's tiles, one per layer, split in order among DETAIL_GROUPS arrays +// (lavapipe allocates no more than 2 GB at once), and the clamped sampler that reads +// them. +const DETAIL_GROUPS: u32 = 4u; +@group(0) @binding(10) var earth_detail_0: texture_2d_array; +@group(0) @binding(11) var earth_detail_1: texture_2d_array; +@group(0) @binding(12) var earth_detail_2: texture_2d_array; +@group(0) @binding(13) var earth_detail_3: texture_2d_array; +@group(0) @binding(14) var detail_sampler: sampler; + // The mesh pass's targets, SAMPLES per pixel. @group(1) @binding(0) var mesh_radiance: MeshRadiance; @group(1) @binding(1) var mesh_depth: MeshDepth; @@ -139,6 +155,59 @@ fn read_map(map: texture_2d, at: MapLookup) -> vec4 { return textureSampleGrad(map, earth_sampler, at.uv, at.dx, at.dy); } +// Earth's surface at `at` from the sharp map's tile there: the tile's layer, and the +// coordinates and gradients scaled from the whole map to the tile, inside its border. +fn read_detail(at: MapLookup) -> vec3 { + let tiles = u.earth_detail.yz; + let uv = vec2(fract(at.uv.x), clamp(at.uv.y, 0.0, 0.99999994)) * tiles; + let cell = min(floor(uv), tiles - 1.0); + let scale = tiles * u.earth_detail_size.x; + let st = u.earth_detail_size.y + (uv - cell) * u.earth_detail_size.x; + let layer = u32(cell.y * tiles.x + cell.x); + let per_group = (u32(tiles.x * tiles.y) + DETAIL_GROUPS - 1u) / DETAIL_GROUPS; + let local = i32(layer % per_group); + let dx = at.dx * scale; + let dy = at.dy * scale; + switch layer / per_group { + case 0u: { + return textureSampleGrad(earth_detail_0, detail_sampler, st, local, dx, dy).rgb; + } + case 1u: { + return textureSampleGrad(earth_detail_1, detail_sampler, st, local, dx, dy).rgb; + } + case 2u: { + return textureSampleGrad(earth_detail_2, detail_sampler, st, local, dx, dy).rgb; + } + default: { + return textureSampleGrad(earth_detail_3, detail_sampler, st, local, dx, dy).rgb; + } + } +} + +// Earth's surface reflectance at `at`. With the sharp map, it is read where a pixel +// covers less ground than its coarsest level's texels resolve, the global map beyond +// that, and the two are blended across half a level so the change does not show. +fn read_earth(at: MapLookup) -> vec3 { + if u.earth_detail.x < 0.5 { + return read_map(earth_texture, at).rgb; + } + // The pixel's footprint in the sharp map's full-size texels, along its longer axis. + let texels = u.earth_detail_size.zw; + let fx = at.dx * texels; + let fy = at.dy * texels; + let level = 0.5 * log2(max(max(dot(fx, fx), dot(fy, fy)), 1e-12)); + let coarsest = u.earth_detail.w - 1.0; + let global = smoothstep(coarsest - 0.5, coarsest, level); + var rgb = vec3(0.0); + if global < 1.0 { + rgb = read_detail(at) * (1.0 - global); + } + if global > 0.0 { + rgb += read_map(earth_texture, at).rgb * global; + } + return rgb; +} + // 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); @@ -160,7 +229,7 @@ fn earth_radiance(p: vec3, n_fixed: vec3, at: MapLookup) -> vec4 if u.params.y < 0.5 { return u.earth_lambert * lit * eclipse; } - let albedo = textured_albedo(read_map(earth_texture, at).rgb, u.earth_rgb); + let albedo = textured_albedo(read_earth(at), u.earth_rgb); if u.params.w < 0.5 { return u.earth_irradiance * albedo * lit * eclipse; } diff --git a/crates/resolved/tests/gpu/earth_detail.rs b/crates/resolved/tests/gpu/earth_detail.rs new file mode 100644 index 0000000..b03c857 --- /dev/null +++ b/crates/resolved/tests/gpu/earth_detail.rs @@ -0,0 +1,229 @@ +//! The sharp tiled Earth map (`EarthDetail`): read across tile edges and the 180 degree +//! meridian as one map would be, handed over to the global map in far views, and +//! compressed to BC6H within a measured error; on the GPU and on the CPU (lavapipe). + +use crate::common; + +use half::f16; +use resolved::bodies::Appearance; +use resolved::earth::{srgb_to_linear, DetailFormat, DetailLayout, EarthDetail}; +use resolved::gpu::AdapterChoice; +use resolved::nalgebra::{UnitQuaternion, Vector3}; +use resolved::radiometry::AU; +use resolved::{ + CameraConfig, EarthTexture, Fidelity, LightingTier, Pose, Renderer, RendererOptions, Scene, +}; +use std::sync::Arc; + +const ADAPTERS: [AdapterChoice; 2] = [AdapterChoice::Gpu, AdapterChoice::Cpu]; + +/// 4 x 2 tiles of 256: a 1024 x 512 map, 39 km a texel at the equator. +const LAYOUT: DetailLayout = DetailLayout { + tiles: [4, 2], + tile: 256, + gutter: 4, + levels: 3, +}; + +/// The map's sRGB value at column `x`, row `y`: dark like ocean in the west half of +/// each tile and bright like desert in the east, with smooth detail at scales from 3 +/// to 60 texels (real ground has edges too, but few as sharp as a checker's). +fn pattern(x: u32, y: u32) -> [u8; 3] { + let (fx, fy) = (f64::from(x), f64::from(y)); + let detail = (fx / 3.1).sin() * (fy / 4.3).cos() * 0.3 + + (fx / 11.0 + fy / 7.0).sin() * 0.4 + + (fy / 60.0).cos() * (fx / 45.0).sin() * 0.3; + let level = |base: f64, swing: f64| (base + swing * detail).round().clamp(0.0, 255.0) as u8; + if (x / 128) % 2 == 1 { + [level(190.0, 40.0), level(160.0, 35.0), level(120.0, 30.0)] + } else { + [level(15.0, 8.0), level(30.0, 12.0), level(60.0, 20.0)] + } +} + +fn global() -> EarthTexture { + let [w, h] = LAYOUT.size(); + let pixels = (0..h) + .flat_map(|y| (0..w).flat_map(move |x| pattern(x, y).into_iter().chain([255]))) + .collect(); + EarthTexture::from_rgba8(w, h, pixels).unwrap() +} + +/// BC6H blocks of a `size` x `size` RGBA f16 image, its edges repeated out to a +/// multiple of 4, at the encoder's fast settings (as `resolved earth-tiles` +/// compresses). +fn bc6h(texels: &[f16], size: u32) -> Vec { + let padded = size.div_ceil(4) * 4; + let mut image = vec![f16::ZERO; 4 * (padded * padded) as usize]; + for y in 0..padded { + for x in 0..padded { + let from = 4 * (y.min(size - 1) * size + x.min(size - 1)) as usize; + let to = 4 * (y * padded + x) as usize; + image[to..to + 4].copy_from_slice(&texels[from..from + 4]); + } + } + let surface = intel_tex_2::RgbaSurface { + data: bytemuck::cast_slice(&image), + width: padded, + height: padded, + stride: 8 * padded, + }; + intel_tex_2::bc6h::compress_blocks(&intel_tex_2::bc6h::fast_settings(), &surface) +} + +fn detail(format: DetailFormat) -> EarthDetail { + let texel = |x: u32, y: u32| pattern(x, y).map(|v| srgb_to_linear(f32::from(v) / 255.0)); + let mut bytes = Vec::new(); + for ty in 0..LAYOUT.tiles[1] { + for tx in 0..LAYOUT.tiles[0] { + for (level, texels) in EarthDetail::build_tile(&LAYOUT, tx, ty, texel) + .iter() + .enumerate() + { + match format { + DetailFormat::Rgba16Float => { + bytes.extend_from_slice(bytemuck::cast_slice(texels)) + } + DetailFormat::Bc6h => bytes.extend(bc6h(texels, LAYOUT.padded(level as u32))), + } + } + } + } + EarthDetail::from_bytes(LAYOUT, format, bytes).unwrap() +} + +/// Lighting L1 with the global map `global()`, and the sharp map in `format` if given, +/// with the bare surface (no atmosphere). +fn renderer(adapter: &AdapterChoice, format: Option) -> Renderer { + let options = RendererOptions { + adapter: adapter.clone(), + fidelity: Fidelity { + lighting: LightingTier::L1, + ..Default::default() + }, + earth_texture: Some(Arc::new(global())), + earth_detail: format.map(|f| Arc::new(detail(f))), + moon_texture: Some(common::uniform_moon()), + appearance: Appearance { + earth_atmosphere: false, + ..Default::default() + }, + read_radiance: true, + ..Default::default() + }; + Renderer::new(CameraConfig::example(), options).unwrap() +} + +/// Looking straight down from `altitude` over geocentric `lat`, `lon` in degrees, the +/// Sun behind the camera. +fn view(lat: f64, lon: f64, altitude: f64) -> (Scene, Pose) { + let (la, lo) = (lat.to_radians(), lon.to_radians()); + let up = Vector3::new(la.cos() * lo.cos(), la.cos() * lo.sin(), la.sin()); + let scene = Scene { + earth_attitude: UnitQuaternion::identity(), + sun: up * AU, + moon: Vector3::new(0.0, 0.0, -3.844e8), + moon_attitude: UnitQuaternion::identity(), + targets: vec![], + velocities: Default::default(), + }; + let camera = Pose::looking_at(up * (6.4e6 + altitude), Vector3::zeros(), Vector3::z()); + (scene, camera.unwrap()) +} + +fn radiance(renderer: &mut Renderer, view: &(Scene, Pose)) -> Vec<[f32; 4]> { + let frame = renderer.submit(&view.0, &view.1).unwrap().wait().unwrap(); + frame.radiance.unwrap() +} + +/// The largest and the mean difference between two images, over the brightest value +/// in them. +fn differences(a: &[[f32; 4]], b: &[[f32; 4]]) -> (f64, f64) { + let peak = a + .iter() + .flat_map(|p| &p[..3]) + .fold(0.0f32, |m, &v| m.max(v)); + let all: Vec = a + .iter() + .zip(b) + .flat_map(|(p, q)| (0..3).map(move |c| f64::from((p[c] - q[c]).abs() / peak))) + .collect(); + let worst = all.iter().fold(0.0, |m: f64, &v| m.max(v)); + (worst, all.iter().sum::() / all.len() as f64) +} + +#[test] +fn tiles_read_as_one_map() { + // 5000 km up a pixel covers 0.55 km, far finer than a texel, so the sharp map is + // read, magnified. Where four tiles meet, and across the 180 degree meridian, it + // must read as the global map of the same size does. + let _devices = common::device_lock(); + for adapter in &ADAPTERS { + let mut plain = renderer(adapter, None); + let mut tiled = renderer(adapter, Some(DetailFormat::Rgba16Float)); + let name = plain.adapter_info().name.clone(); + for (lat, lon) in [(0.0, 0.0), (0.0, 90.0), (30.0, 180.0)] { + let v = view(lat, lon, 5e6); + let (pa, pb) = (radiance(&mut plain, &v), radiance(&mut tiled, &v)); + // Texture filtering blends neighboring texels in steps of 1/256 on the + // RTX 3060, and the two textures' coordinates fall on different sides of a + // step here and there: up to 1/256 of the contrast between neighbors, 0.45 % + // of the peak in this pattern, spread over the image. A wrong border would + // be off by the whole contrast along the tiles' edges. + let (worst, mean) = differences(&pa, &pb); + assert!(worst < 0.01, "{name} at {lat}, {lon}: {worst}"); + assert!(mean < 1e-4, "{name} at {lat}, {lon}: mean {mean}"); + } + } +} + +#[test] +fn far_views_read_the_global_map() { + // From 3e9 m a pixel covers 330 km at the center of the disk, over 8 of the sharp + // map's 39 km texels: past its coarsest level, so the global map alone, to the bit. + let _devices = common::device_lock(); + for adapter in &ADAPTERS { + let mut plain = renderer(adapter, None); + let mut tiled = renderer(adapter, Some(DetailFormat::Rgba16Float)); + let v = view(20.0, 30.0, 3e9); + let (a, b) = (radiance(&mut plain, &v), radiance(&mut tiled, &v)); + let differing = a.iter().zip(&b).filter(|(p, q)| p != q).count(); + assert_eq!(differing, 0, "{}", plain.adapter_info().name); + } +} + +#[test] +fn bc6h_stays_close_to_the_f16_map() { + // Dark (ocean-like) and bright (desert-like) halves of a tile, magnified. + let _devices = common::device_lock(); + for adapter in &ADAPTERS { + let mut exact = renderer(adapter, Some(DetailFormat::Rgba16Float)); + let mut compressed = renderer(adapter, Some(DetailFormat::Bc6h)); + let name = exact.adapter_info().name.clone(); + for (lat, lon, what) in [(30.0, -160.0, "dark"), (30.0, -110.0, "bright")] { + let v = view(lat, lon, 5e6); + let (a, b) = (radiance(&mut exact, &v), radiance(&mut compressed, &v)); + // Relative error of the mean over the image, and per pixel at the 99th + // percentile. Measured on both devices: means within 0.16 %, the 99th + // percentile 2.2 % on the dark half and 0.77 % on the bright. On Blue + // Marble's own July map around Sinai, 1.0 to 2.5 %, means under 0.1 %. + let mean = |img: &[[f32; 4]], c: usize| { + img.iter().map(|p| f64::from(p[c])).sum::() / img.len() as f64 + }; + let mut errors: Vec = a + .iter() + .zip(&b) + .flat_map(|(p, q)| { + (0..3).map(move |c| (f64::from(q[c]) / f64::from(p[c]) - 1.0).abs()) + }) + .collect(); + errors.sort_by(f64::total_cmp); + let p99 = errors[errors.len() * 99 / 100]; + for c in 0..3 { + let bias = mean(&b, c) / mean(&a, c) - 1.0; + assert!(bias.abs() < 5e-3, "{name} {what} channel {c}: mean {bias}"); + } + assert!(p99 < 0.04, "{name} {what}: {p99}"); + } + } +} diff --git a/crates/resolved/tests/gpu/main.rs b/crates/resolved/tests/gpu/main.rs index 2dd6146..56db604 100644 --- a/crates/resolved/tests/gpu/main.rs +++ b/crates/resolved/tests/gpu/main.rs @@ -11,6 +11,7 @@ mod blooming; mod bodies; mod convolve; mod distortion; +mod earth_detail; mod earth_texture; mod earthshine; mod eclipses; diff --git a/docs/optics.md b/docs/optics.md index 7428dad..926615d 100644 --- a/docs/optics.md +++ b/docs/optics.md @@ -1033,21 +1033,53 @@ point's geodetic latitude and longitude from the ellipsoid's normal there, in Earth-fixed axes, and reads the image at that place. Earth's orientation from the scene turns the continents with the epoch. -The full image is 500 m per pixel at the equator, resized to 660 m to fit the -largest texture lavapipe allows (16384 x 8192). Whether that is fine or coarse -depends on the view: +The image comes in two sizes: + +- The global map, always there: NASA's 21600 x 10800 image (1 arcminute, + 1.85 km per pixel at the equator), resized to 16384 x 8192 (2.45 km), the + largest texture lavapipe allows. +- A sharp map, optional: NASA's 86400 x 43200 image (15 arcseconds, 460 m), + made once per month by `scripts/fetch-earth.sh --sharp` (below). + +Whether a map is fine or coarse depends on the view: - From low orbit looking straight down, a camera pixel of the example camera - covers 77 m of ground, so each texture pixel spans about 8 image pixels and - looks soft. + covers 77 m of ground. A texel of the global map spans 32 image pixels and + looks blurred; one of the sharp map, 6. - From geostationary orbit, a camera pixel covers about 4 km, so each image - pixel averages about 35 texture pixels (6 x 6). + pixel averages about 2 x 2 texels of the global map. For the second case the texture is stored with mipmaps: copies at half, a quarter, an eighth of the size, and so on. The GPU reads from the copy whose pixels best match the image pixel, so a far view shows the true average color of a region rather than a shimmering sample of individual texture pixels. +The sharp map is too large for one texture on either device (lavapipe's limit +is 16384 on a side) and for one allocation on lavapipe (2 GB), and 30 GB as +16-bit floating point. So it is cut into 8 x 4 square tiles of 10800 texels, +each stored with a border of 4 texels copied from its neighbors, so that +blending near a tile's edge reads the right texels without a second tile, and +compressed to BC6H: blocks of 4 x 4 texels in 16 bytes, which the GPU decodes +as it reads them. BC6H keeps linear values with about the precision of 16-bit +floating point; on Blue Marble's July map around Sinai, 99 % of pixels are +within 1 to 2.5 % of the uncompressed value and the mean within 0.1 %. Three +mipmap levels are stored (460 m, 930 m, 1.85 km), 4.9 GB in all, in four +texture arrays of 8 tiles. A pixel reads the sharp map where it covers less +ground than the coarsest level's texels, the global map beyond, and the two +blended across half a level so the change does not show. + +Its cost on the RTX 3060: 0.1 to 0.2 ms more in the bodies pass where it is +read, and 1.6 s more setup to upload it. On lavapipe, which decodes BC6H in +software, about 100 ms more in the bodies pass where it is read and 5 s more +setup, with 8 GB of memory at its peak. From geostationary orbit it is not +read, and costs nothing per frame. + +The GPU and lavapipe choose mipmap levels and blend texels slightly +differently, so a textured Earth differs between them by about 1 DN on +average at lighting L1 and camera C1 (0.8 DN with the global map alone, +1.3 DN with the sharp map; 2.7 DN before the atmosphere, which lowers the +contrast). The tests of the maps compare each device with itself. + Two details matter for brightness: - Averaging has to happen on reflectance, not on the encoded values image files @@ -1896,6 +1928,8 @@ 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 surface finer than Blue Marble's 460 m: mosaics of tens of meters + (Landsat, Sentinel-2) need tiles streamed in as a view needs them - 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/readme.md b/readme.md index 4d3fe31..a9f59bc 100644 --- a/readme.md +++ b/readme.md @@ -125,7 +125,10 @@ because most NASA models are Draco-compressed and are decoded on the way in. At lighting L1, Earth's surface comes from NASA's Blue Marble, which `scripts/fetch-earth.sh` downloads and resizes (with ImageMagick) into `data/earth/`, one image per month. The month comes from the scene's epoch, or -from `[appearance] earth_month` in a scene without one. The Moon's surface +from `[appearance] earth_month` in a scene without one. `scripts/fetch-earth.sh +--sharp 07` also makes July's sharp map (460 m per texel instead of 2.45 km; 330 MB +to download and about 3 minutes to make 4.9 GB), which the CLI and the movie use +when it is there. The Moon's surface comes from the LRO Camera's normalized-reflectance mosaic, which `scripts/fetch-moon.py` downloads (1.2 GB) and assembles (with numpy and Pillow) into `data/moon/wac-hapke-643.png`. The Moon is turned by the scene's epoch, or diff --git a/scripts/fetch-earth.sh b/scripts/fetch-earth.sh index d0a67a9..062e81d 100755 --- a/scripts/fetch-earth.sh +++ b/scripts/fetch-earth.sh @@ -2,25 +2,158 @@ # Downloads NASA's Blue Marble Next Generation monthly Earth images and resizes # them for the renderer. # +# scripts/fetch-earth.sh [--sharp MM ...] +# # Source: Blue Marble: Next Generation (2004), base map without topography or -# bathymetry, 21600 x 10800 (500 m per pixel), one image per month, from +# bathymetry, one image per month, from # https://science.nasa.gov/earth/earth-observatory/blue-marble-next-generation/. # Produced by Reto Stöckli, NASA Earth Observatory; NASA asks that users credit # "NASA Earth Observatory". # -# The originals go to data/earth/source/ and are checked by sha256. Each is -# resized to 16384 x 8192 (660 m per pixel), the largest texture lavapipe takes, -# into data/earth/bmng-2004-MM.jpg. Resizing uses ImageMagick (`magick`). -# Both directories are gitignored. +# The global maps: NASA's 21600 x 10800 images (1 arcminute, 1.85 km per pixel at +# the equator), resized to 16384 x 8192 (2.45 km), the largest texture lavapipe +# takes, into data/earth/bmng-2004-MM.jpg. Resizing uses ImageMagick (`magick`). +# +# With --sharp and month numbers (01 to 12), also the sharp maps for those months: +# NASA's 86400 x 43200 images (15 arcseconds, 460 m), published as eight JPEG tiles +# of 21600 x 21600 (about 330 MB a month), made into data/earth/sharp-2004-MM.tiles +# by `resolved earth-tiles` (about 3 minutes and 4.9 GB a month; it builds the CLI +# first). +# +# The originals go to data/earth/source/ and are checked by sha256. Both +# directories are gitignored. set -eu +sharp="" +if [ "${1:-}" = "--sharp" ]; then + shift + sharp="$*" +fi + base='https://assets.science.nasa.gov/content/dam/science/esd/eo/images/bmng/bmng-base' root="$(cd "$(dirname "$0")/.." && pwd)/data/earth" mkdir -p "$root/source" # month number, folder name, sha256 of the original -while read -r n month sum; do +months=$(cat <<'EOF' +01 january c2f2c80a352518cfea23d462dad6f6e778375ddd74bbec5786bcfc85cc49c61d +02 february 8800081f38b53bfa1e3f0c2fa17be694b7646942fe5add009ea0dcd64c03e70b +03 march af37f61a9715c59dc86763df19f588a20c9d92b3d0a8034a0b00acd6b97ccb7d +04 april 3bd13f173555eb24451c45d3a83ca177baca02e9a8bdd9fa3e4989f7b496ce6b +05 may 0de40c868680ab9e86fef42406d76386b243b32444ff50edfd406b296638c206 +06 june 282608580a290723655bfbe917e69ea2f93f4da068dfa96ff4a152b1c6ec36be +07 july dea8b4dc8a4f93f5f8bce0c8c85a508a178e7901e9ed8e6bf86e6ce7ef6d61e2 +08 august 926be7253be90a1fd57071a33afbfe7a712a2dfe0d5b5b4602bc5dfbcfdc2c89 +09 september 7cf788e13a3a7b4a926b524f8f71d635c56a18c48ba0bfe35695b05a305db3f3 +10 october bf568ceca2ad55e9c8dd7bcedfd0e419201d03bb560206a7bfa2cc57c1e04c64 +11 november 1959bbe72897d3e9526c440dae3ec04cfdbdc4080c8ff0b14632925ae31167b1 +12 december b0ffc9d12205689aeaec66116bbae2c4c263fcd98a0cbadba8d5488567b53278 +EOF +) + +# The sharp maps' tiles: file name, sha256. +sharp_sums=$(cat <<'EOF' +world.200401.3x21600x21600.A1.jpg 8953cad5e491003b5d82fb230c9338bb88bf02af91a24c1fa80ab72fd0bb55ff +world.200401.3x21600x21600.A2.jpg 8c81cc3f212c7edfd3592d4395394fec48dd65ebffb32f1f5de73c0417b87c58 +world.200401.3x21600x21600.B1.jpg 8e89d3871fbdaff5dc8ede580c9f04e2ce5498eee284505bed23427229f180ec +world.200401.3x21600x21600.B2.jpg 7221c0de3eedb184f94f7e2319ce6de73f62c82b27367fe191f5f88cb626859d +world.200401.3x21600x21600.C1.jpg 1e03fa28fd3d638e99fd9d350964d8ca2e4e0b4285148be57d7081417203082f +world.200401.3x21600x21600.C2.jpg 50f16c3af137e38e8fefb86beb732e35eb45549d8a69a7028d1deb38f6e17b1d +world.200401.3x21600x21600.D1.jpg 2a29d4ac59e60f5c8f8d94100f70c38516607c5242109b1b48c6cb7e34ec3e93 +world.200401.3x21600x21600.D2.jpg f5e2366d9a79a81313a839e702c8aa620ce0201aa7f4e75bccb8cf5f451115cb +world.200402.3x21600x21600.A1.jpg 00a1fc049c50da3e69e6ba540e2e2e8214455073c9a3f8bfb7f9bf895f7c896f +world.200402.3x21600x21600.A2.jpg 1396f92a7de89abfa226acf4b94232b1111cc8eb03675060da9d1f691f3c4071 +world.200402.3x21600x21600.B1.jpg d817a0747c6c85d80ed6e888e75bc56c0fa86c9187623e222a87472a33e5bbe9 +world.200402.3x21600x21600.B2.jpg 5b3268f2d49eae7d3f12cd6694a8dc8adca56f30830fe55a31c1780d725761b6 +world.200402.3x21600x21600.C1.jpg e937d30375aec54060b1916e788886942fdc4362e2a64bcebcc4a3b29c9e864f +world.200402.3x21600x21600.C2.jpg 8f375bd01443969c3e4d3e53fd4fc1e553e78a940ee387af73ba34621180d753 +world.200402.3x21600x21600.D1.jpg 992144f1cdb8afb496b357ccff1d049cac43838e27c5ff2ccc5c094eb0c27995 +world.200402.3x21600x21600.D2.jpg c10ac7eb90ed4dee3613e6a8ac9cfc4be68f853506166e81c1d48c719f060269 +world.200403.3x21600x21600.A1.jpg 28d54db1428452e7250cf2a75d1ab82aa26fac8831290c2fa8a9ae08ea079c10 +world.200403.3x21600x21600.A2.jpg 8de3f51e27f27639a226805a62b17f50bcbade3595435e5545f068c93375e864 +world.200403.3x21600x21600.B1.jpg 04d7d7118f1b7fb37f308da68cd9e311eb9b8c0920c79cd95ab41986547b7dd6 +world.200403.3x21600x21600.B2.jpg 6f696d87d0e227639af6671c5abd8ab9fe4574493cadefd572246583910ae224 +world.200403.3x21600x21600.C1.jpg 120238fb5f29d3685aaa1a26194423adea19bd26bf5781550c37852573e0bc2b +world.200403.3x21600x21600.C2.jpg 5939731d8f1677490dedf80cc7153c1c13936ff15546fcf647ddd9aceb73e901 +world.200403.3x21600x21600.D1.jpg c74ce51a78a7a058385bc0779fe2fb9446a5fc9fb82079a12d48260606fd4e4f +world.200403.3x21600x21600.D2.jpg 79577840e4a3fe458018dd77989c1c0a61498150bb62ec683837c47b8d062a9a +world.200404.3x21600x21600.A1.jpg 993d148cdbae790c9a4908c1d292b81ff3b85a1b2bfdaf82a6d56b35d96bf035 +world.200404.3x21600x21600.A2.jpg cb105e494c8f82b1b02d7ccbef65e06038f3bb9968e7ea4a77c6c047da059ecd +world.200404.3x21600x21600.B1.jpg df6ba4cca1403dea3696e0614cc05d9823ce84da44be561502e88c060363eff8 +world.200404.3x21600x21600.B2.jpg 8f15641e63f3520b1e188715258a866a65ecc170b5ba79a3046461eac2dfba73 +world.200404.3x21600x21600.C1.jpg c2f06c539c182fe1206a8bcdb86abd50451d0f78ec905a02542745fb0255addb +world.200404.3x21600x21600.C2.jpg e892b4f40b447cd0c691ec05cb7575976b7f013bc42e64ad60e952a48a3940a8 +world.200404.3x21600x21600.D1.jpg b4eb36215675adf5c610f0648bb8c0747364d5e65a91280878935d4a52eff57e +world.200404.3x21600x21600.D2.jpg d39bc15b5608d039b7ec4a12d657359ed4a7d790672874cde86f4f469e252973 +world.200405.3x21600x21600.A1.jpg d59031bf900a16e71ab696551fc0c9c4d76e4a0ced034647a700fca80deb11e8 +world.200405.3x21600x21600.A2.jpg a3698e369922d2b1a58796f121a7c904b74b4a63a83cb327c2cbbb352c8f3ade +world.200405.3x21600x21600.B1.jpg b5787874d6d2b7636b709ea394ea988714a85f13eac18e116f9fa7ebe26139e7 +world.200405.3x21600x21600.B2.jpg 1ad9ef53e9b00b358de21cd439e1791e0451795b1ec663aed1149df8cb1fc1a0 +world.200405.3x21600x21600.C1.jpg b3773bcc5ceb876f4a217db2145709f095f1209009c504fa0540370e3c338e5d +world.200405.3x21600x21600.C2.jpg 4d2cd09c5b1e54ff0ac4c8556750076f8fdc8e1c0513febfe996ef10312e0b22 +world.200405.3x21600x21600.D1.jpg 035ed9a398c33fee0cd84cb2ad5647ab0fbd454208d170ecdd3ee5087e008b01 +world.200405.3x21600x21600.D2.jpg a24ea5c3b4d8e007b6e6396ed9707aa67a0a395a22db394f44ce3b33c3026dd6 +world.200406.3x21600x21600.A1.jpg df597db0ea202a4ed68a27e9023077bcb8a1c70f870503aa3d8d6595433b26da +world.200406.3x21600x21600.A2.jpg 14fc4556cf8c371291c2c1f6bb1ea6c99ca78f4009c10b370721942ca91382ae +world.200406.3x21600x21600.B1.jpg 2fd24d69ab7859164fa12885cbf808b7870725c5fc33cb79e63299bda5fda185 +world.200406.3x21600x21600.B2.jpg ecec960dd81d2fe4653099eaf9e09a9ca2ce9b45a9783c473321d52c0e7503e7 +world.200406.3x21600x21600.C1.jpg e85aa5230be16689b52dc5b4b3b70fb47c2eb830acd242c54ee75ab879050de6 +world.200406.3x21600x21600.C2.jpg 37ba0e21027ec9e207d97f5d52d3768c91525efbd73969c5fb06376696aed2cf +world.200406.3x21600x21600.D1.jpg 12b98d5892bc4b3e9822dbd134702da88d85d521e25778c1ebef225112a267bd +world.200406.3x21600x21600.D2.jpg 669ef6b76fcf2009bd7eedcaa1101c7c4f30d17314fd82ec679ac4abb0d1b0c0 +world.200407.3x21600x21600.A1.jpg 045bdba424fb4ab5a385b6282fe36d1ee43ceb336630bc9d1197db8d08038456 +world.200407.3x21600x21600.A2.jpg 6d4d048ee6f3aaf28f4b8bbf727d9d635217d81171783111ae09b706e8a7d108 +world.200407.3x21600x21600.B1.jpg 4fe552a037a000ff778df08900f185e491534f2b2a8f1779a85512f3ac8fd625 +world.200407.3x21600x21600.B2.jpg ba7f8e4b9ef4858e14d74b5a7d69be142524f6838775ed0ec00296ca26c89e99 +world.200407.3x21600x21600.C1.jpg 18c4b869510f63bb2261182e26e5c072faee7bc6553842bf9d52175da539e956 +world.200407.3x21600x21600.C2.jpg 15fe50659ddb7b3f6395a76d7d880bfa1c78c7c957b9ed3b57dda60cac0866a1 +world.200407.3x21600x21600.D1.jpg d2358197c6cbc479051141a8455ee456425c6479e94498888338813daecade68 +world.200407.3x21600x21600.D2.jpg 5a5b804bacae0dd4303753a8bd35f26187270374205ab1a85ddab9ad2040cf9d +world.200408.3x21600x21600.A1.jpg 80351fd9ce7eb7f2fbd827fdc5f23526429baefe6eaf091876042fc94c80013c +world.200408.3x21600x21600.A2.jpg b996913af7858f94365ced0e189579a706c30495cde85cfcee1967a63a46fd8b +world.200408.3x21600x21600.B1.jpg 654dc2c927ba44830edfbd843d4dff0264f03f8016eff89260de5b7b766fb963 +world.200408.3x21600x21600.B2.jpg 10ba13cd1d1751d16f7f8c375900a6f3dad826a7dba95486410154fb72fa4224 +world.200408.3x21600x21600.C1.jpg 3d7fe4fa9bb2c56d233e6e01e9212348604d9b6490f472710745a3b9ea4cd5cc +world.200408.3x21600x21600.C2.jpg 725bced59b4a90ed14d0ff05bc3b0526fe12f65a891a959afe7d496860accbed +world.200408.3x21600x21600.D1.jpg 20f5b09793b0bdb949e9a3a996efa6eb6eacf8c98e40c646ea1f3be3059def5b +world.200408.3x21600x21600.D2.jpg 00a840df7a226fec4cdb82847cf6d877782bf0f1f3167ede2083b633d3361310 +world.200409.3x21600x21600.A1.jpg a8a71afa554fb5810756d42b71ad6e18c6fbfdb16176adc91654793b0211deff +world.200409.3x21600x21600.A2.jpg f158edd3ec90b78405d230482ccdbfe40d777877b5b1ef2eefbed2071c17baa7 +world.200409.3x21600x21600.B1.jpg 42167660a4f09cbbf5219682967a14510ba9c44a9e0ecffd1aaacd28b5f8fdb4 +world.200409.3x21600x21600.B2.jpg 78036a95153a5db851c568a29567518d0f61a2b5c06135780981c3f208a3950a +world.200409.3x21600x21600.C1.jpg 905c30a23836312bbb19f0087ff2e18fa54dca606b06e30bb12ddfb0340ec46e +world.200409.3x21600x21600.C2.jpg 2f2e84bdb740ea408e9709c77d78bd87d63dbfc726c9c91d0c41c3510daca46d +world.200409.3x21600x21600.D1.jpg 54fa2c1b6417e05cc2af635abe88b1254af51e035285a787ba79e85d8ba29cf6 +world.200409.3x21600x21600.D2.jpg 739fcda45577910b4c6f2d71693e133a14f4a10b14ac37456b190f6dc4505be9 +world.200410.3x21600x21600.A1.jpg 0578984d7b37979587d58eb274de1d107a61bedd453d6762aa5df1eba78a8a8c +world.200410.3x21600x21600.A2.jpg a3132f72c63f78f37020faafa871d31616e94f3673a4760d8a2b1567d4de5dd5 +world.200410.3x21600x21600.B1.jpg 7515176ed467c16049680eca30a4973155e53e9f0ab0ba209b14b1a19974211b +world.200410.3x21600x21600.B2.jpg 4bf0e010d844ac29a5602cbab96b69a351c78d97f3f5455a652cf801fd783d4a +world.200410.3x21600x21600.C1.jpg 02cc6a2f699bf16a40dd652410eb1aa24cf93dbe8395c1fd2ba38978560a278f +world.200410.3x21600x21600.C2.jpg 6964e8c97b66432ea7c917d2b144b2210d9dd9b46d55304e9fda5f9bc7db21c5 +world.200410.3x21600x21600.D1.jpg 8ae0622dc19e3d02413120ada7cecbd3ae9dab20980acf70471f13c5a6af1257 +world.200410.3x21600x21600.D2.jpg c5cb23e3d4e1b45c5019ef38984cc9a851242df0321248fc8e9370a97579b447 +world.200411.3x21600x21600.A1.jpg 0431cdb10c77ff309837b11b228d2bca5a56badcd109ea1a5bc2751da5f54ad9 +world.200411.3x21600x21600.A2.jpg f9b4b6dd2e6f104ef1a0a6d69a88cc639a8c171e4d6fef214afcf777c58b24e2 +world.200411.3x21600x21600.B1.jpg e2eb1b81c9884cb9cbb9fc583d6bcdf41bfa0fb94132073b12494990bbaeda37 +world.200411.3x21600x21600.B2.jpg 17cf75b5931e2c8b7a3a74b1f8f5d9504a92853f48afde425704487d016fbdb0 +world.200411.3x21600x21600.C1.jpg 6d19ac865ae48cb39eceea447d246bb84077c316565737088f92fb5723f6182b +world.200411.3x21600x21600.C2.jpg b5f1fb519b56d784127e44cd661ea0053c3a571db30f7a54aa128b5ab33741af +world.200411.3x21600x21600.D1.jpg f83fa6f494ec8cbfdc527befde834dbe28ee3513b6982c02cf2088a3d6101edf +world.200411.3x21600x21600.D2.jpg 9cb66f1666d2e7540faab834dd488d6a12c954d901a1cffc9238d9fca712f305 +world.200412.3x21600x21600.A1.jpg 40150990562535ce70bd9ef8271b1493ced90d351101c115a63c79e28764c474 +world.200412.3x21600x21600.A2.jpg d4feffe2584e754a3b7fe7ba22ab0ee85f99d27f8452e8b72bc5f931a25f647a +world.200412.3x21600x21600.B1.jpg 115ad94094e501556aab585bf5ecf371cd6d180f187d3d260db445b5fd8a90dc +world.200412.3x21600x21600.B2.jpg b26573ae3bb0784c60ad9684a457a2e9532dbbade15688c0bad6276d2a94e88b +world.200412.3x21600x21600.C1.jpg 103062deb2d87b5fbbe6e0fbe704c53fd2f66f676d2cc6e8de7b8e680659d18b +world.200412.3x21600x21600.C2.jpg c67b3a314d57fb80218557871f2424f13791dd6aced67f9fe83bad5058e5f7e1 +world.200412.3x21600x21600.D1.jpg 6499c86513d02077cf309b2b65344af246bb9b2006017a3aecc3e5da4cb0f906 +world.200412.3x21600x21600.D2.jpg eb545d1a728133b065c5466b668c63abb7402de1bd0bbf4b67e2a9095869f368 +EOF +) + +echo "$months" | while read -r n month sum; do original="$root/source/world.2004$n.3x21600x10800.jpg" if [ ! -f "$original" ]; then echo "fetching $month" @@ -35,19 +168,35 @@ while read -r n month sum; do magick "$original" -colorspace RGB -filter Lanczos -resize 16384x8192! \ -colorspace sRGB -quality 95 "$resized" fi -done <<'EOF' -01 january c2f2c80a352518cfea23d462dad6f6e778375ddd74bbec5786bcfc85cc49c61d -02 february 8800081f38b53bfa1e3f0c2fa17be694b7646942fe5add009ea0dcd64c03e70b -03 march af37f61a9715c59dc86763df19f588a20c9d92b3d0a8034a0b00acd6b97ccb7d -04 april 3bd13f173555eb24451c45d3a83ca177baca02e9a8bdd9fa3e4989f7b496ce6b -05 may 0de40c868680ab9e86fef42406d76386b243b32444ff50edfd406b296638c206 -06 june 282608580a290723655bfbe917e69ea2f93f4da068dfa96ff4a152b1c6ec36be -07 july dea8b4dc8a4f93f5f8bce0c8c85a508a178e7901e9ed8e6bf86e6ce7ef6d61e2 -08 august 926be7253be90a1fd57071a33afbfe7a712a2dfe0d5b5b4602bc5dfbcfdc2c89 -09 september 7cf788e13a3a7b4a926b524f8f71d635c56a18c48ba0bfe35695b05a305db3f3 -10 october bf568ceca2ad55e9c8dd7bcedfd0e419201d03bb560206a7bfa2cc57c1e04c64 -11 november 1959bbe72897d3e9526c440dae3ec04cfdbdc4080c8ff0b14632925ae31167b1 -12 december b0ffc9d12205689aeaec66116bbae2c4c263fcd98a0cbadba8d5488567b53278 -EOF +done + +for n in $sharp; do + case "$n" in + 0[1-9] | 1[0-2]) ;; + *) + echo "--sharp takes months 01 to 12, not $n" >&2 + exit 1 + ;; + esac + month=$(echo "$months" | awk -v n="$n" '$1 == n { print $2 }') + for t in A1 B1 C1 D1 A2 B2 C2 D2; do + name="world.2004$n.3x21600x21600.$t.jpg" + original="$root/source/$name" + if [ ! -f "$original" ]; then + echo "fetching $month $t" + curl -fsSL "$base/$month/$name" -o "$original.part" + mv "$original.part" "$original" + fi + sum=$(echo "$sharp_sums" | awk -v f="$name" '$1 == f { print $2 }') + echo "$sum $original" | sha256sum -c --quiet - + done + tiles="$root/sharp-2004-$n.tiles" + if [ ! -f "$tiles" ]; then + echo "tiling $month" + (cd "$root/../.." && cargo run --quiet --release -p resolved-cli -- earth-tiles \ + --source "$root/source" --month "$n" --out "$tiles.part") + mv "$tiles.part" "$tiles" + fi +done echo "Earth images in $root" diff --git a/scripts/gallery-delta.py b/scripts/gallery-delta.py index 1761e48..9c00b3b 100755 --- a/scripts/gallery-delta.py +++ b/scripts/gallery-delta.py @@ -46,8 +46,8 @@ OUT = ROOT / "gallery" / "delta" TIERS = ["L0", "L1"] CAMERA = "C1" # Renders at once. Each loads its own copy of the Earth and Moon maps; four at once ran -# the RTX 3060 (12 GB) out of memory. -WORKERS = 2 +# the RTX 3060 (12 GB) out of memory, and with a sharp Earth map (4.9 GB) two do. +WORKERS = 1 if any((ROOT / "data" / "earth").glob("sharp-*.tiles")) else 2 BRIGHTER = np.array([255, 150, 40], dtype=float) DARKER = np.array([60, 140, 255], dtype=float)