feat(render): the settled tilted frame's mesh grows with the export's pixels, up to the DEM's cells (#1882)

crates/render/src/lib.rs

@@ -2590,9 +2590,10 @@ fn render_tilted(map: &mut Map, view: &View, full_res_rasters: bool, sources_rea
 /// triangles raised to the ground's height and mapped from its picture, every band keeping to
 /// its own stretch of ground, all drawn far to near so nearer slopes hide what lies behind
 /// them; fog with distance, and the sky behind. (Without terrain, [`Warp`] does it per pixel.)
-/// Reading the grids' heights and ordering the triangles (tens of milliseconds) doesn't look at
-/// Cancel; drawing them, the long part, does. False when an export cancelled meanwhile stopped
-/// it, the picture left half drawn.
+/// Cancel is looked at before each band's grid of heights is read (up to 641 × 641 points,
+/// `terrain::settled_cells`) and while the triangles are drawn, the long part; ordering them
+/// doesn't. False when an export cancelled meanwhile, or a settled frame gone out of date
+/// (`cancel`), stopped it, the picture left half drawn.
 fn drape_terrain(pm: &mut Pixmap, view: &View, map: &Map, bands: &[(View, (f64, f64), Pixmap)], ground: &terrain::Ground) -> bool {
     let (sky, fog) = tilt_sky(map);
     let cam = view.camera;
@@ -2617,16 +2618,20 @@ fn drape_terrain(pm: &mut Pixmap, view: &View, map: &Map, bands: &[(View, (f64,
         None => pm.fill(c(fog)),
     }
     // Every band's triangles: (depth, band, cell * 2 + which half), with each band's grid of
-    // vertices (screen point, picture point, depth). Indices, not corners: the settled mesh is
-    // fine (`terrain::SETTLED_CELLS`), and corners copied into each triangle took 112 bytes a
-    // triangle against 16 (13 MB against 8 with the grids, for a 1600 × 1100 export).
+    // vertices (screen point, picture point, depth) and its cells a side. Indices, not corners: the
+    // settled mesh is fine (`terrain::settled_cells`), and corners copied into each triangle took
+    // 112 bytes a triangle against 16 (13 MB against 8 with the grids, for a 1600 × 1100 export
+    // at 160 cells).
     type Vertex = ((f64, f64), (f64, f64), f64);
-    let n = terrain::SETTLED_CELLS;
     let (pw, ph) = (f64::from(pm.width()), f64::from(pm.height()));
-    let mut grids: Vec<Vec<Option<Vertex>>> = Vec::with_capacity(bands.len());
+    let mut grids: Vec<(usize, Vec<Option<Vertex>>)> = Vec::with_capacity(bands.len());
     let mut tris: Vec<(f64, u32, u32)> = Vec::new();
     for (bi, (v, (d0, d1), pic)) in bands.iter().enumerate() {
+        if rustgis_data::progress::check().is_err() || cancel::requested() {
+            return false;
+        }
         let kb = f64::from(pic.width()) / v.width.max(1e-9);
+        let n = band_cells(v, pic.width().max(pic.height()), ground);
         let grid: Vec<Option<Vertex>> = (0..=n)
             .flat_map(|j| (0..=n).map(move |i| (i, j)))
             .map(|(i, j)| {
@@ -2656,18 +2661,18 @@ fn drape_terrain(pm: &mut Pixmap, view: &View, map: &Map, bands: &[(View, (f64,
                 tris.push((depth, bi as u32, (cell * 2 + half) as u32));
             }
         }
-        grids.push(grid);
+        grids.push((n, grid));
     }
     tris.sort_by(|a, b| b.0.total_cmp(&a.0));
     let fog_paint = |t: f64| paint(fog, t as f32);
     for (drawn, (depth, bi, at)) in tris.into_iter().enumerate() {
         // An export cancelled meanwhile (see `rustgis_data::progress`) stops here: it writes
-        // nothing. The settled frame on screen has no run to cancel.
-        if drawn % 4096 == 0 && rustgis_data::progress::check().is_err() {
+        // nothing. So does a settled frame on screen gone out of date (see `cancel`).
+        if drawn % 4096 == 0 && (rustgis_data::progress::check().is_err() || cancel::requested()) {
             return false;
         }
-        let (Some((_, _, pic)), Some(grid)) = (bands.get(bi as usize), grids.get(bi as usize)) else { continue };
-        let cell = at as usize / 2;
+        let (Some((_, _, pic)), Some((n, grid))) = (bands.get(bi as usize), grids.get(bi as usize)) else { continue };
+        let (n, cell) = (*n, at as usize / 2);
         let Some([p, q, r]) = triangles(grid, n, cell % n, cell / n).and_then(|t| t.get(at as usize % 2).copied()) else { continue };
         let (dst, src) = ([p.0, q.0, r.0], [p.1, q.1, r.1]);
         let Some(t) = affine(src, dst) else { continue };
@@ -2685,6 +2690,14 @@ fn drape_terrain(pm: &mut Pixmap, view: &View, map: &Map, bands: &[(View, (f64,
     true
 }
 
+/// Cells along each side of the mesh laying the band `v` over `ground`, its picture
+/// `picture_px` on its longer side: as fine as the picture, up to the ground's own cells across
+/// the band (#1882).
+fn band_cells(v: &View, picture_px: u32, ground: &terrain::Ground) -> usize {
+    let corners = [(0.0, 0.0), (v.width, 0.0), (0.0, v.height), (v.width, v.height)].map(|(x, y)| v.camera.to_map(x, y, v.width, v.height));
+    terrain::settled_cells(picture_px, ground.cells_across(&corners))
+}
+
 /// The two triangles of the cell (`i`, `j`) of a band's grid of `n` cells a side, [a, b, d] and
 /// [a, d, c] of its corners a b / c d; none when a corner couldn't be placed.
 fn triangles<V: Copy>(grid: &[Option<V>], n: usize, i: usize, j: usize) -> Option<[[V; 3]; 2]> {
@@ -9231,9 +9244,52 @@ mod tests {
         assert!(drawn < flat_y - 0.5 * (flat_y - raised_y), "drawn from {drawn}: flat {flat_y}, raised {raised_y}");
     }
 
+    #[test]
+    fn the_settled_mesh_follows_the_picture_up_to_the_dem() {
+        // #1882: each band's mesh grows with its picture, within 160 and 640 cells a side, and
+        // never past the DEM's own cells across its ground; a small DEM keeps 160.
+        use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster};
+        let camera = Camera { center: Coord::new(0.0, 0.0), res: 10.0, bearing: 30.0, pitch: 65.0 };
+        let cells = |side: u32, cell: f64, ppp: f64| -> Vec<(usize, Option<f64>)> {
+            let half = f64::from(side) * cell / 2.0;
+            let gt = GeoTransform { origin_x: -half, origin_y: half, pixel_width: cell, pixel_height: -cell };
+            let heights = vec![0.0; (side * side) as usize];
+            let dem = Raster::new("dem", side, side, vec![heights], gt, Crs::WebMercator, None, ColorHint::Gray, "Float32").unwrap();
+            let mut m = Map::new("Map");
+            m.add_raster("dem", rustgis_map::RasterLayer::from_path(dem, "dem.tif"));
+            m.view.terrain =
+                Some(rustgis_map::view::Terrain { source: rustgis_map::view::TerrainSource::Layer("dem".into()), exaggeration: 1.0 });
+            let g = terrain::Ground::of(&m, &camera).expect("terrain");
+            let view = View { width: 1600.0, height: 1100.0, pixels_per_point: ppp, camera, sky: false };
+            camera
+                .tilt_bands(view.width, view.height)
+                .iter()
+                .map(|b| {
+                    let v = band_view(&view, b);
+                    let (w, h) = v.pixel_size();
+                    let c = [(0.0, 0.0), (v.width, 0.0), (0.0, v.height), (v.width, v.height)];
+                    let across = g.cells_across(&c.map(|(x, y)| v.camera.to_map(x, y, v.width, v.height)));
+                    (band_cells(&v, w.max(h), &g), across)
+                })
+                .collect()
+        };
+        // A DEM of 8 m cells, 16 km a side: finer than the bands' pictures.
+        let (screen, print) = (cells(2000, 8.0, 1.0), cells(2000, 8.0, 3.0));
+        assert!(print.iter().all(|(n, _)| (terrain::SETTLED_CELLS..=terrain::MAX_SETTLED_CELLS).contains(n)), "{print:?}");
+        assert!(screen.iter().zip(&print).all(|(a, b)| a.0 <= b.0), "{screen:?} vs {print:?}");
+        assert!(print.iter().any(|(n, _)| *n == terrain::MAX_SETTLED_CELLS), "{print:?}");
+        assert!(screen.iter().any(|(n, _)| *n > terrain::SETTLED_CELLS && *n < terrain::MAX_SETTLED_CELLS), "{screen:?}");
+        // A coarser DEM caps the bands at its cells across them, no lower than 160.
+        let coarse = cells(400, 40.0, 3.0);
+        assert!(coarse.iter().all(|(n, d)| *n <= d.unwrap_or(0.0).ceil().max(terrain::SETTLED_CELLS as f64) as usize), "{coarse:?}");
+        assert!(coarse.iter().zip(&print).any(|(c, p)| c.0 < p.0), "{coarse:?} vs {print:?}");
+        // A 120-cell DEM: 160 everywhere, as before.
+        assert!(cells(120, 100.0, 3.0).iter().all(|(n, _)| *n == terrain::SETTLED_CELLS));
+    }
+
     #[test]
     fn a_cancelled_export_stops_laying_the_map_over_its_terrain() {
-        // Cancelled once the terrain's drawing has looked once: halfway through its triangles.
+        // Cancelled once the terrain's drawing has looked once: before its second band's grid.
         struct CancelledMidway(std::sync::atomic::AtomicUsize);
         impl rustgis_data::progress::Probe for CancelledMidway {
             fn cancelled(&self) -> bool {
@@ -9265,6 +9321,32 @@ mod tests {
         assert_eq!(render_png(&mut m, &view).err().as_deref(), Some(rustgis_data::DataError::Cancelled.to_string().as_str()));
     }
 
+    #[test]
+    fn a_settled_frame_gone_out_of_date_stops_laying_its_terrain() {
+        // A settled frame's worker stops it (`cancel`): the drape, now up to 640 cells a band
+        // (#1882), stops too, before reading another band's grid.
+        let camera = Camera { center: lonlat_to_merc(Coord::new(0.0, -0.6)), res: 600.0, bearing: 0.0, pitch: 65.0 };
+        let view = View { width: 400.0, height: 300.0, pixels_per_point: 1.0, camera, sky: false };
+        let mut m = hill_map();
+        m.view.terrain =
+            Some(rustgis_map::view::Terrain { source: rustgis_map::view::TerrainSource::Layer("dem".into()), exaggeration: 2.0 });
+        let ground = terrain::Ground::of(&m, &camera).expect("terrain");
+        let bands: Vec<(View, (f64, f64), Pixmap)> = camera
+            .tilt_bands(view.width, view.height)
+            .iter()
+            .filter_map(|b| {
+                let v = band_view(&view, b);
+                let (w, h) = v.pixel_size();
+                Some((v, b.depth, Pixmap::new(w, h)?))
+            })
+            .collect();
+        let mut pm = Pixmap::new(400, 300).expect("a pixmap");
+        assert!(drape_terrain(&mut pm, &view, &m, &bands, &ground));
+        let flag = std::sync::Arc::new(std::sync::atomic::AtomicBool::new(true));
+        let _stop = cancel::install(flag);
+        assert!(!drape_terrain(&mut pm, &view, &m, &bands, &ground));
+    }
+
     #[test]
     fn a_tilted_view_draws_the_ground_in_perspective_under_a_sky() {
         let mut m = Map::new("Map");

crates/render/src/terrain.rs

@@ -17,11 +17,33 @@ use rustgis_tiles::template::{self, TileId};
 /// drape, `tilt_drape` in the app): few enough to lay out every frame.
 pub const CELLS: usize = 40;
 
-/// Cells along each side of a band laid over the relief in the frame that stays on screen once
-/// the map settles, and in exports (#1523): fine enough that a ridge's outline follows the
-/// ridge, where 40 drew it in long straight lines and let a massif sink into its surroundings.
+/// The fewest cells along each side of a band laid over the relief in the frame that stays on
+/// screen once the map settles, and in exports (#1523; [`settled_cells`] grows it with the
+/// picture, #1882): fine enough that a ridge's outline follows the ridge, where 40 drew it in long
+/// straight lines and let a massif sink into its surroundings.
 pub const SETTLED_CELLS: usize = 4 * CELLS;
 
+/// The most cells along each side of a band laid over the relief in the settled frame and in
+/// exports (#1882): 1280 showed no more of a ridge than 640, for about three times the memory.
+pub const MAX_SETTLED_CELLS: usize = 4 * SETTLED_CELLS;
+
+/// Pixels of a band's picture to each of its mesh's cells in the settled frame and in exports
+/// (#1882): a finer mesh outruns the picture draped over it.
+const PX_PER_SETTLED_CELL: u32 = 4;
+
+/// Cells along each side of a band laid over the relief in the settled frame and in exports
+/// (#1882): one every [`PX_PER_SETTLED_CELL`] pixels of the longer side of the band's picture,
+/// `picture_px`, and no more than `dem_cells`, the ground's own cells across the band
+/// ([`Ground::cells_across`], `None` when unknown); never fewer than [`SETTLED_CELLS`] (no
+/// coarser than before) nor more than [`MAX_SETTLED_CELLS`].
+pub fn settled_cells(picture_px: u32, dem_cells: Option<f64>) -> usize {
+    let mut n = (picture_px / PX_PER_SETTLED_CELL) as usize;
+    if let Some(d) = dem_cells.filter(|d| d.is_finite()) {
+        n = n.min(d.max(0.0).ceil() as usize);
+    }
+    n.clamp(SETTLED_CELLS, MAX_SETTLED_CELLS)
+}
+
 /// The deepest zoom of elevation tiles read.
 const MAX_TILE_ZOOM: u8 = 15;
 /// How many zooms up a missing tile is stood in for by a coarser one already read.
@@ -116,6 +138,50 @@ impl Ground {
         per_unit(m, self.body)
     }
 
+    /// How many of the ground's own cells span the stretch of ground with these `corners` (Web
+    /// Mercator), along the longer side of their box in the ground's grid: a DEM's cells, the
+    /// box clipped to the DEM (a geographic DEM's columns met again a whole turn round, as
+    /// [`raster_height`] reads them; a whole turn's not clipped at all); elevation tiles'
+    /// pixels at the zoom read. A turned stretch's box is up to √2 longer than its side.
+    /// `None` when no corner lands on the grid's plane.
+    pub fn cells_across(&self, corners: &[Coord]) -> Option<f64> {
+        // `turn`: a geographic DEM's columns in 360°, where `raster_height` meets it again.
+        let (pixels, (w, h), full_turn, turn): (Vec<(f64, f64)>, _, _, _) = match &self.heights {
+            Heights::Raster { raster, wrap_width } => {
+                let t = raster.transformer();
+                let pixels = corners.iter().filter_map(|m| t.from_merc(*m)).map(|n| raster.transform.to_pixel(n)).collect();
+                let turn = t.is_geographic().then(|| 360.0 / raster.transform.pixel_width.abs()).filter(|t| t.is_finite());
+                (pixels, (f64::from(raster.width), f64::from(raster.height)), wrap_width.is_some(), turn)
+            }
+            Heights::Tiles { zoom, .. } => {
+                let half = std::f64::consts::PI * MERCATOR_R;
+                let n = 256.0 * f64::from(1u32 << (*zoom).min(30));
+                let pixels = corners.iter().map(|m| ((m.x + half) / (2.0 * half) * n, (half - m.y) / (2.0 * half) * n)).collect();
+                (pixels, (n, n), false, None)
+            }
+        };
+        let (mut x0, mut y0, mut x1, mut y1) = (f64::INFINITY, f64::INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY);
+        for (x, y) in pixels.into_iter().filter(|(x, y)| x.is_finite() && y.is_finite()) {
+            (x0, y0, x1, y1) = (x0.min(x), y0.min(y), x1.max(x), y1.max(y));
+        }
+        if x0 > x1 {
+            return None;
+        }
+        let across = |a: f64, b: f64, size: f64| (b.min(size) - a.max(0.0)).max(0.0);
+        let x = match turn {
+            // A whole turn's columns go on past either edge: no edge to clip at, one turn at most.
+            _ if full_turn => (x1 - x0).clamp(0.0, w),
+            // Moved round by the whole turns that put most of the stretch on the DEM (a 0–360°
+            // grid, a stretch across the seam from a DEM beside it, a view a turn away).
+            Some(t) => {
+                let k = ((w - x0 - x1) / 2.0 / t).round();
+                [k - 1.0, k, k + 1.0].iter().map(|k| across(x0 + k * t, x1 + k * t, w)).fold(0.0, f64::max)
+            }
+            None => across(x0, x1, w),
+        };
+        Some(x.max(across(y0, y1, h)))
+    }
+
     /// Whether elevation tiles are still on their way (drawing again later shows more).
     pub fn pending(&self) -> bool {
         match &self.heights {
@@ -298,6 +364,72 @@ mod tests {
         assert!(ground.height(lonlat_to_merc(Coord::new(-22.5, 0.0))).is_none());
     }
 
+    #[test]
+    fn the_settled_mesh_follows_the_picture_within_its_bounds() {
+        // #1882: one cell every 4 px, within 160 and 640.
+        for (px, n) in [(0, 160), (400, 160), (640, 160), (1600, 400), (2000, 500), (2560, 640), (4800, 640), (u32::MAX, 640)] {
+            assert_eq!(settled_cells(px, None), n, "{px} px");
+        }
+        // Never past the DEM's cells across the band, nor coarser than 160.
+        assert_eq!(settled_cells(4800, Some(300.0)), 300);
+        assert_eq!(settled_cells(4800, Some(299.2)), 300);
+        assert_eq!(settled_cells(1600, Some(1924.0)), 400);
+        assert_eq!(settled_cells(4800, Some(100.0)), SETTLED_CELLS);
+        assert_eq!(settled_cells(4800, Some(0.0)), SETTLED_CELLS);
+        assert_eq!(settled_cells(1600, Some(f64::NAN)), 400);
+        assert_eq!(settled_cells(4800, Some(f64::INFINITY)), 640);
+    }
+
+    #[test]
+    fn a_stretch_spans_the_ground_cells_under_it() {
+        // A 4 × 3 grid of 90° by 60° cells, the whole Earth: a stretch over two columns spans 2
+        // cells, one past its edges only the grid's.
+        let ground = raster_ground(90.0, Crs::LonLat, None);
+        let at = |lon: f64, lat: f64| lonlat_to_merc(Coord::new(lon, lat));
+        let two = ground.cells_across(&[at(0.0, 30.0), at(180.0, 30.0), at(0.0, -30.0), at(180.0, -30.0)]).unwrap();
+        assert!((two - 2.0).abs() < 1e-6, "{two}");
+        let all = ground.cells_across(&[Coord::new(-1e8, 1e8), Coord::new(1e8, -1e8)]).unwrap();
+        assert!(all <= 4.0 + 1e-9, "{all}");
+        assert_eq!(ground.cells_across(&[]), None);
+        // Across the seam, the whole stretch: 135° to 225° (-135°) spans one cell, not the half
+        // cell short of the grid's edge; the same stretch a turn away likewise.
+        let seam = ground.cells_across(&[at(135.0, 30.0), at(225.0, 30.0)]).unwrap();
+        assert!((seam - 1.0).abs() < 1e-6, "{seam}");
+        let away = ground.cells_across(&[at(495.0, 30.0), at(585.0, 30.0)]).unwrap();
+        assert!((away - 1.0).abs() < 1e-6, "{away}");
+        // A regional DEM of 1° cells beside the seam: a stretch across the seam from it counts
+        // the cells it reaches there, as its heights are read, though its middle is off the DEM.
+        let regional = |west: f64| {
+            let gt = GeoTransform { origin_x: west, origin_y: 90.0, pixel_width: 1.0, pixel_height: -60.0 };
+            let raster =
+                Arc::new(Raster::new("heights", 10, 3, vec![vec![0.0; 30]], gt, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap());
+            Ground { heights: Heights::Raster { raster, wrap_width: None }, exaggeration: 1.0, body: 1.0 }
+        };
+        let east_of_seam = regional(-180.0).cells_across(&[at(170.0, 30.0), at(186.0, 30.0)]).unwrap();
+        assert!((east_of_seam - 6.0).abs() < 1e-6, "{east_of_seam}");
+        let a_turn_on = regional(0.0).cells_across(&[at(345.0, 30.0), at(365.0, 30.0)]).unwrap();
+        assert!((a_turn_on - 5.0).abs() < 1e-6, "{a_turn_on}");
+        // Elevation tiles: their pixels at the zoom read; the whole world at zoom 2 is 1024.
+        let tiles = Ground {
+            heights: Heights::Tiles {
+                url: String::new(),
+                encoding: TerrainEncoding::Terrarium,
+                zoom: 2,
+                read: RefCell::new(HashMap::new()),
+                pending: RefCell::new(false),
+                waiting: RefCell::new(std::collections::HashSet::new()),
+            },
+            exaggeration: 1.0,
+            body: 1.0,
+        };
+        let half = std::f64::consts::PI * MERCATOR_R;
+        let world = tiles.cells_across(&[Coord::new(-2.0 * half, half), Coord::new(2.0 * half, -half)]).unwrap();
+        assert!((world - 1024.0).abs() < 1e-6, "{world}");
+        // An eighth of the world across, a sixteenth down: the longer side.
+        let quarter = tiles.cells_across(&[Coord::new(0.0, 0.0), Coord::new(half / 2.0, -half / 4.0)]).unwrap();
+        assert!((quarter - 256.0).abs() < 1e-6, "{quarter}");
+    }
+
     #[test]
     fn heights_stand_in_the_metres_of_the_body() {
         // A metre on Mars spans more of the (Earth-sized) map than one on the Earth (#223).