perf(dem): the planner terrain fills its tiles on every core and decodes two rows of GLO-30 squares at a time

host/obc-dem/Cargo.toml

@@ -18,7 +18,7 @@ description = "Copernicus GLO-30 GeoTIFF -> OBCT terrain cells (obc-dem fetch /
 [features]
 default = ["dem"]
 dem = ["geotiff", "fetch"]
-geotiff = ["dep:tiff", "dep:serde_json"]
+geotiff = ["dep:tiff", "dep:serde_json", "dep:rayon"]
 fetch = ["dep:obc-data"]
 # The terrain of the planner maps: Terrarium WebP tiles in MBTiles.
 terrarium = ["dem", "dep:image", "dep:rusqlite"]
@@ -59,6 +59,9 @@ serde_json = { workspace = true, optional = true }
 # `fetch` only: the store and its fetcher for the tiles. The binary also prints the Copernicus
 # credit from data/sources.toml.
 obc-data = { workspace = true, optional = true }
+# `geotiff` only: a planner terrain bake decodes the GLO-30 squares of a row and fills a row of
+# tiles on every core.
+rayon = { version = "1", optional = true }
 image = { workspace = true, features = ["webp"], optional = true }
 rusqlite = { workspace = true, features = ["bundled"], optional = true }
 

host/obc-dem/src/geotiff.rs

@@ -22,6 +22,7 @@
 //! Sampling the mosaic is `f64` throughout and bilinear over the four surrounding source posts; see
 //! [`DemMosaic::height`] for why that is a point sample and not an area aggregate.
 
+use std::cell::Cell;
 use std::fs::File;
 use std::io::BufReader;
 use std::path::{Path, PathBuf};
@@ -141,6 +142,36 @@ impl DemTile {
     /// free of them: a tile that survives this is a regular geographic post grid whose every sample
     /// position is one multiply away.
     pub fn open(path: &Path) -> Result<DemTile, String> {
+        let (mut dec, tile) = Self::read_header(path)?;
+        let (rows, cols, name) = (tile.rows, tile.cols, || path.display().to_string());
+        // The samples, flipped once.
+        let image = dec.read_image().map_err(|e| format!("{}: decode failed ({e})", name()))?;
+        let north_up: Vec<f32> = match image {
+            DecodingResult::F32(v) => v,
+            DecodingResult::F64(v) => v.into_iter().map(|x| x as f32).collect(),
+            DecodingResult::I16(v) => v.into_iter().map(f32::from).collect(),
+            DecodingResult::U16(v) => v.into_iter().map(f32::from).collect(),
+            _ => return Err(format!("{}: unsupported sample type — a DEM must be float or 16-bit integer", name())),
+        };
+        if north_up.len() != rows * cols {
+            return Err(format!("{}: decoded {} samples for a {cols}×{rows} raster", name(), north_up.len()));
+        }
+        let mut data = vec![0f32; rows * cols];
+        for row in 0..rows {
+            let src = &north_up[(rows - 1 - row) * cols..(rows - row) * cols];
+            data[row * cols..(row + 1) * cols].copy_from_slice(src);
+        }
+        Ok(DemTile { data, ..tile })
+    }
+
+    /// The checks of [`DemTile::open`] that need no sample: a caller that decodes the tile later
+    /// finds a file that is not a DEM tile before it starts.
+    pub fn check(path: &Path) -> Result<(), String> {
+        Self::read_header(path).map(|_| ())
+    }
+
+    /// The decoder of the tile at `path`, and the tile without its samples.
+    fn read_header(path: &Path) -> Result<(Decoder<BufReader<File>>, DemTile), String> {
         let file = File::open(path).map_err(|e| format!("{}: {e}", path.display()))?;
         let mut dec =
             Decoder::new(BufReader::new(file)).map_err(|e| format!("{}: not a readable TIFF ({e})", path.display()))?;
@@ -167,34 +198,17 @@ impl DemTile {
             .ok()
             .and_then(|s| s.trim().trim_end_matches('\0').parse::<f64>().ok());
 
-        // The samples, flipped once.
-        let image = dec.read_image().map_err(|e| format!("{}: decode failed ({e})", name()))?;
-        let north_up: Vec<f32> = match image {
-            DecodingResult::F32(v) => v,
-            DecodingResult::F64(v) => v.into_iter().map(|x| x as f32).collect(),
-            DecodingResult::I16(v) => v.into_iter().map(f32::from).collect(),
-            DecodingResult::U16(v) => v.into_iter().map(f32::from).collect(),
-            _ => return Err(format!("{}: unsupported sample type — a DEM must be float or 16-bit integer", name())),
-        };
-        if north_up.len() != rows * cols {
-            return Err(format!("{}: decoded {} samples for a {cols}×{rows} raster", name(), north_up.len()));
-        }
-        let mut data = vec![0f32; rows * cols];
-        for row in 0..rows {
-            let src = &north_up[(rows - 1 - row) * cols..(rows - row) * cols];
-            data[row * cols..(row + 1) * cols].copy_from_slice(src);
-        }
-
-        Ok(DemTile {
+        let tile = DemTile {
             south_lat_deg: north_lat_deg - step_lat_deg * (rows - 1) as f64,
             west_lon_deg,
             step_lat_deg,
             step_lon_deg,
             rows,
             cols,
-            data,
+            data: Vec::new(),
             nodata,
-        })
+        };
+        Ok((dec, tile))
     }
 
     /// Posts along each axis, `(rows, cols)`.
@@ -309,18 +323,22 @@ pub(crate) fn geo_keys<R: std::io::Read + std::io::Seek>(dec: &mut Decoder<R>) -
 #[derive(Debug, Default)]
 pub struct DemMosaic {
     tiles: Vec<DemTile>,
-    /// The tile that answered the previous query. A bake walks a cell sample by sample, so
-    /// essentially every query lands in the same tile as the one before it, and without this the
-    /// lookup is a linear scan of the mosaic per sample.
+}
+
+thread_local! {
+    /// The tile that answered this thread's previous query. A bake walks a cell sample by sample,
+    /// so essentially every query lands in the same tile as the one before it, and without this
+    /// the lookup is a linear scan of the mosaic per sample. One per thread, so threads that sample
+    /// different places of one mosaic keep their own.
     ///
     /// It cannot change an answer: tile coverage boxes are disjoint by construction, so at most one
     /// tile ever matches and checking a different one first only reorders the search.
-    last_hit: std::cell::Cell<usize>,
+    static LAST_HIT: Cell<usize> = const { Cell::new(0) };
     /// The same memo, for the corner resolver alone. `height` resolves an off-tile stencil corner
     /// through [`DemMosaic::nearest_post`], which is a different tile by definition; sharing one
     /// slot means every seam-adjacent sample overwrites the memo with the neighbour and the next
     /// sample overwrites it back, which is two linear scans per sample exactly along the seams.
-    last_corner_hit: std::cell::Cell<usize>,
+    static LAST_CORNER_HIT: Cell<usize> = const { Cell::new(0) };
 }
 
 impl DemMosaic {
@@ -355,8 +373,11 @@ impl DemMosaic {
     pub fn push(&mut self, tile: DemTile) {
         self.tiles.push(tile);
         self.tiles.sort_by_key(DemTile::order_key);
-        self.last_hit.set(0); // the sort invalidated both memo indices
-        self.last_corner_hit.set(0);
+    }
+
+    /// Keep only the tiles for which `keep` holds.
+    pub fn retain(&mut self, keep: impl FnMut(&DemTile) -> bool) {
+        self.tiles.retain(keep);
     }
 
     /// How many tiles the mosaic holds.
@@ -370,12 +391,12 @@ impl DemMosaic {
 
     /// The tile that answers for `(lat, lon)`, or `None` outside coverage.
     fn tile_for(&self, lat_deg: f64, lon_deg: f64) -> Option<&DemTile> {
-        self.tile_for_memo(lat_deg, lon_deg, &self.last_hit)
+        LAST_HIT.with(|memo| self.tile_for_memo(lat_deg, lon_deg, memo))
     }
 
     /// [`DemMosaic::tile_for`], against a caller-chosen memo slot. Which slot is used can only
     /// change how long the search takes, never its answer, because coverage boxes are disjoint.
-    fn tile_for_memo(&self, lat_deg: f64, lon_deg: f64, memo: &std::cell::Cell<usize>) -> Option<&DemTile> {
+    fn tile_for_memo(&self, lat_deg: f64, lon_deg: f64, memo: &Cell<usize>) -> Option<&DemTile> {
         if let Some(tile) = self.tiles.get(memo.get()) {
             if tile.owns(lat_deg, lon_deg) {
                 return Some(tile);
@@ -391,7 +412,7 @@ impl DemMosaic {
     /// is what lets two tiles with different post spacings meet without inventing a lattice that
     /// neither of them is on.
     fn nearest_post(&self, lat_deg: f64, lon_deg: f64) -> Option<f64> {
-        let tile = self.tile_for_memo(lat_deg, lon_deg, &self.last_corner_hit)?;
+        let tile = LAST_CORNER_HIT.with(|memo| self.tile_for_memo(lat_deg, lon_deg, memo))?;
         let row = ((lat_deg - tile.south_lat_deg) / tile.step_lat_deg).round();
         let col = ((lon_deg - tile.west_lon_deg) / tile.step_lon_deg).round();
         let row = (row as i64).clamp(0, tile.rows as i64 - 1) as usize;

host/obc-dem/src/planner.rs

@@ -1,10 +1,13 @@
 //! The heights of the planner: its terrain tiles and the slopes of its routes. The reference
 //! archive gives the height where it holds the ground, and GLO-30 elsewhere.
 
-use std::collections::{BTreeSet, VecDeque};
+use std::collections::{BTreeMap, BTreeSet, VecDeque};
+use std::ops::RangeInclusive;
 use std::path::{Path, PathBuf};
+use std::sync::Mutex;
 
-use obc_data::store::sha256_hex;
+use obc_data::store::hash_file;
+use rayon::prelude::*;
 
 use crate::fetch::TileId;
 use crate::geotiff::{DemMosaic, DemTile};
@@ -18,30 +21,55 @@ pub fn tiles([west, south, east, north]: [f64; 4]) -> Vec<TileId> {
 }
 
 /// The GLO-30 files of `bounds` in `dir`: a square without a file is sea.
-pub fn tiles_in(dir: &Path, bounds: [f64; 4]) -> Vec<PathBuf> {
-    tiles(bounds).iter().map(|tile| dir.join(tile.file_name())).filter(|path| path.is_file()).collect()
+pub fn tiles_in(dir: &Path, bounds: [f64; 4]) -> Vec<(TileId, PathBuf)> {
+    let files = tiles(bounds).into_iter().map(|tile| (tile, dir.join(tile.file_name())));
+    files.filter(|(_, path)| path.is_file()).collect()
 }
 
+/// The GLO-30 squares, by latitude, whose posts a height at `lat` in microdegrees can read: its
+/// bilinear corners lie within two posts, and a square holds the posts of its own degree.
+pub fn squares(lat: i32) -> RangeInclusive<i32> {
+    const SLACK: i32 = 10_000;
+    (lat - SLACK).div_euclid(1_000_000)..=(lat + SLACK).div_euclid(1_000_000)
+}
+
+/// The reference tiles that the last heights read. A caller that samples near its last sample
+/// keeps one, so a reference tile decodes once for many heights.
+#[derive(Default)]
+pub struct Recent(VecDeque<((u32, u32), Option<ReferenceTile>)>);
+
 pub struct Terrain {
+    /// The GLO-30 files by the latitude of their square.
+    glo30: BTreeMap<i32, Vec<PathBuf>>,
+    /// The decoded squares of the latitudes in `loaded`.
     fallback: DemMosaic,
+    loaded: BTreeSet<i32>,
     reference: Option<ReferenceArchive>,
-    cache: VecDeque<((u32, u32), Option<ReferenceTile>)>,
     /// The SHA-256 of each GLO-30 file, then of each reference tile that `bounds` reaches.
     pub identities: Vec<String>,
-    used: BTreeSet<String>,
+    used: Mutex<BTreeSet<String>>,
     /// The credit of GLO-30.
     credit: String,
 }
 
 impl Terrain {
     /// The GLO-30 files `glo30`, in the order of [`tiles`], with their `credit`, and the reference
-    /// archive at `reference`, for heights in `bounds`.
-    pub fn open(glo30: &[PathBuf], credit: &str, reference: Option<&Path>, bounds: [f64; 4]) -> Result<Self, String> {
-        let mut fallback = DemMosaic::default();
-        let mut identities = Vec::new();
-        for path in glo30 {
-            identities.push(sha256_hex(&std::fs::read(path).map_err(|e| format!("{}: {e}", path.display()))?));
-            fallback.push(DemTile::open(path)?);
+    /// archive at `reference`, for heights in `bounds`. Each file is hashed and its header checked;
+    /// no square is decoded until [`Terrain::load`].
+    pub fn open(
+        glo30: &[(TileId, PathBuf)],
+        credit: &str,
+        reference: Option<&Path>,
+        bounds: [f64; 4],
+    ) -> Result<Self, String> {
+        let hashes = glo30.par_iter().map(|(_, path)| {
+            DemTile::check(path)?;
+            hash_file(path).map(|(sha256, _)| sha256).map_err(|e| format!("{}: {e}", path.display()))
+        });
+        let mut identities = hashes.collect::<Result<Vec<_>, _>>()?;
+        let mut files = BTreeMap::<_, Vec<_>>::new();
+        for (tile, path) in glo30 {
+            files.entry(tile.lat).or_default().push(path.clone());
         }
         let reference = reference.map(ReferenceArchive::open).transpose()?;
         if let Some(archive) = &reference {
@@ -53,12 +81,34 @@ impl Terrain {
             };
             identities.extend(archive.held_digests(window).into_iter().map(|(_, _, hash)| hash.to_owned()));
         }
-        let (cache, used, credit) = (VecDeque::new(), BTreeSet::new(), credit.into());
-        Ok(Self { fallback, reference, cache, identities, used, credit })
+        let (fallback, loaded, used) = (DemMosaic::default(), BTreeSet::new(), Mutex::default());
+        Ok(Self { glo30: files, fallback, loaded, reference, identities, used, credit: credit.into() })
+    }
+
+    /// Decode the GLO-30 squares of the latitudes `lats` and drop every other square.
+    pub fn load(&mut self, lats: RangeInclusive<i32>) -> Result<(), String> {
+        let missing: Vec<_> = self
+            .glo30
+            .range(lats.clone())
+            .filter(|(lat, _)| !self.loaded.contains(lat))
+            .flat_map(|(_, paths)| paths)
+            .collect();
+        let decoded = missing.into_par_iter().map(|path| DemTile::open(path)).collect::<Result<Vec<_>, _>>()?;
+        // A GLO-30 square holds the posts north of its south edge, up to and with the next degree.
+        self.fallback.retain(|tile| lats.contains(&(tile.south_west_deg().0.floor() as i32)));
+        decoded.into_iter().for_each(|tile| self.fallback.push(tile));
+        self.loaded = self.glo30.range(lats).map(|(lat, _)| *lat).collect();
+        Ok(())
     }
 
     /// The height in metres at `lat`, `lon` in microdegrees, or `None` where no source has one.
-    pub fn height(&mut self, lat: i32, lon: i32) -> Result<Option<f64>, String> {
+    /// GLO-30 answers from the squares that [`Terrain::load`] decoded, which must hold
+    /// [`squares`] of `lat`.
+    pub fn height(&self, lat: i32, lon: i32, recent: &mut Recent) -> Result<Option<f64>, String> {
+        debug_assert!(
+            squares(lat).all(|lat| self.loaded.contains(&lat) || !self.glo30.contains_key(&lat)),
+            "the GLO-30 squares of {lat} µdeg are not loaded"
+        );
         let mut best = None;
         if let Some(archive) = &self.reference {
             let window = Window {
@@ -67,13 +117,14 @@ impl Terrain {
                 lon_lo: lon as i64 - 32,
                 lon_hi: lon as i64 + 32,
             };
+            let cache = &mut recent.0;
             for key in window.tiles() {
-                let tile = if let Some(index) = self.cache.iter().position(|(id, _)| id == &key) {
-                    self.cache.remove(index).unwrap().1
+                let tile = if let Some(index) = cache.iter().position(|(id, _)| id == &key) {
+                    cache.remove(index).unwrap().1
                 } else {
                     match archive.tile(key.0, key.1)? {
                         TileLookup::Held { tile, sources } => {
-                            self.used.extend(sources.iter().cloned());
+                            self.used.lock().unwrap().extend(sources.iter().cloned());
                             Some(tile)
                         }
                         TileLookup::Absent | TileLookup::Unknown => None,
@@ -82,9 +133,9 @@ impl Terrain {
                 if let Some(tile) = &tile {
                     tile.centres_in(window, |_, _, metres| best = Some(metres as f64));
                 }
-                self.cache.push_back((key, tile));
-                if self.cache.len() > 32 {
-                    self.cache.pop_front();
+                cache.push_back((key, tile));
+                if cache.len() > 32 {
+                    cache.pop_front();
                 }
             }
         }
@@ -95,15 +146,16 @@ impl Terrain {
     pub fn attribution(&self) -> Vec<String> {
         let mut credits = Vec::new();
         if let Some(archive) = &self.reference {
+            let used = self.used.lock().unwrap();
             credits.extend(
                 archive
                     .credits()
                     .iter()
-                    .filter(|c| self.used.contains(&c.key))
+                    .filter(|c| used.contains(&c.key))
                     .map(|c| format!("{} ({})", c.attribution, c.licence)),
             );
         }
-        if !self.fallback.is_empty() {
+        if !self.glo30.is_empty() {
             credits.push(self.credit.clone());
         }
         credits

host/obc-dem/src/step.rs

@@ -9,6 +9,7 @@ use obc_data::engine::Request;
 use serde_json::Value;
 
 use crate::bake::{published_cell, write_cell_file};
+use crate::fetch::TileId;
 use crate::geotiff::{DemMosaic, DemTile};
 use crate::reference::ReferenceArchive;
 
@@ -107,10 +108,10 @@ pub fn bounds(options: &Value) -> Result<[f64; 4], String> {
 }
 
 /// The GLO-30 files of the request, in the order of [`crate::planner::tiles`] of `bounds`.
-pub fn glo30(request: &Request, bounds: [f64; 4]) -> Vec<PathBuf> {
+pub fn glo30(request: &Request, bounds: [f64; 4]) -> Vec<(TileId, PathBuf)> {
     let files = request.snapshots.get(GLO30).into_iter().flatten();
     let files: BTreeMap<&str, &PathBuf> =
         files.map(|(name, path)| (name.rsplit('/').next().unwrap_or(name), path)).collect();
-    let tiles = crate::planner::tiles(bounds);
-    tiles.iter().filter_map(|tile| files.get(tile.file_name().as_str()).map(|path| path.to_path_buf())).collect()
+    let tiles = crate::planner::tiles(bounds).into_iter();
+    tiles.filter_map(|tile| files.get(tile.file_name().as_str()).map(|path| (tile, path.to_path_buf()))).collect()
 }