sharper earth map

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]

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<Vec<u8>, 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<u8>,
+    east: Vec<u8>,
+    north: Vec<u8>,
+    south: Vec<u8>,
+}
+
+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<u8> {
+    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<f32> = (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(())
+}

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)",