From ff6a4ecf2cb92d8da7ac8cebf7d85f9f95977e03 Mon Sep 17 00:00:00 2001 From: Nikita Zavartsev Date: Fri, 12 Jun 2026 20:42:09 +0200 Subject: [PATCH 1/9] fix: [dsm-rgb-landcover-tirol] catch GDAL_NODATA and float-max NoData sentinels in DEM readers The BEV ALS DSM declares NoData = +3.4028235e38 via the GDAL_NODATA ASCII tag (42113). The old predicate (NaN / < -1000 / == 0) let those cells through as valid 3.4e38 m elevations, poisoning max_terrain_h and the raymarch sky exit. extract_window, parse_geotiff_auto and the overview cache builder now share is_nodata_value(), which also honours the file's declared GDAL_NODATA value for finite custom sentinels. Co-Authored-By: Claude Fable 5 --- crates/dem_io/src/geotiff.rs | 44 ++++++++++++--- crates/dem_io/src/overview.rs | 3 +- crates/dem_io/tests/synthetic.rs | 92 ++++++++++++++++++++++++++++++++ 3 files changed, 132 insertions(+), 7 deletions(-) diff --git a/crates/dem_io/src/geotiff.rs b/crates/dem_io/src/geotiff.rs index cbe2fa3..eb21bb3 100644 --- a/crates/dem_io/src/geotiff.rs +++ b/crates/dem_io/src/geotiff.rs @@ -8,6 +8,27 @@ use crate::heightmap::fill_nodata; use crate::lzw_lenient; use crate::{DemError, Heightmap}; +/// Read the GDAL_NODATA ASCII tag (42113) from the decoder's current IFD. +/// GDAL writes the NoData sentinel as decimal text (e.g. "3.4028235e+38" for the +/// BEV ALS DSM); absent or unparseable tags yield `None`. +pub(crate) fn read_gdal_nodata( + decoder: &mut Decoder, +) -> Option { + let value = decoder.find_tag(Tag::Unknown(42113)).ok().flatten()?; + let text = value.into_string().ok()?; + text.trim().trim_end_matches('\0').parse::().ok() +} + +/// Shared NoData predicate for float DEM samples. +/// Catches IEEE NaN, the SRTM-style large-negative sentinel, the large-positive +/// float-max sentinel family (BEV ALS DSM uses +3.4028235e38 — without this check +/// it parses as a valid 3.4e38 m elevation and poisons `max_terrain_h` and the +/// raymarch sky exit), and the file's declared GDAL_NODATA value if any. +#[inline] +pub(crate) fn is_nodata_value(v: f32, gdal_nodata: Option) -> bool { + v.is_nan() || !(-1000.0..1.0e38).contains(&v) || gdal_nodata.is_some_and(|nd| v == nd) +} + /// Returns ModelPixelScaleTag[0] from a GeoTIFF, or 0.0 on failure. pub fn geotiff_pixel_scale(path: &Path) -> f64 { let Ok(file) = File::open(path) else { @@ -76,6 +97,8 @@ pub fn parse_geotiff_auto(path: &Path) -> Result { ) }; + let gdal_nodata = read_gdal_nodata(&mut decoder); + let img = decoder.read_image()?; let raw: Vec = match img { DecodingResult::F32(v) => v, @@ -86,7 +109,13 @@ pub fn parse_geotiff_auto(path: &Path) -> Result { const NODATA: f32 = -9999.0; let mut data: Vec = raw .iter() - .map(|&v| if v.is_nan() || v < -1000.0 { NODATA } else { v }) + .map(|&v| { + if is_nodata_value(v, gdal_nodata) { + NODATA + } else { + v + } + }) .collect(); let before = data.iter().filter(|&&v| v == NODATA).count(); @@ -234,6 +263,9 @@ pub fn extract_window( let crs_origin_x = tiepoint[3]; // easting of top-left corner in EPSG:31287 metres let crs_origin_y = tiepoint[4]; + // GDAL_NODATA lives on IFD 0 like the other geo-tags. + let gdal_nodata = read_gdal_nodata(&mut decoder); + // Seek to the requested IFD level for pixel data and actual dimensions. decoder.seek_to_image(ifd_level)?; let (cols, rows): (u32, u32) = decoder.dimensions()?; @@ -345,11 +377,11 @@ pub fn extract_window( // BEV DGM uses 0.0 as NoData sentinel instead of the SRTM convention // (-32 768) or IEEE NaN. This is safe because the minimum valid elevation // in Austria is ~115 m (Hungarian border lowlands) — no real terrain pixel - // can be zero. The three conditions cover: - // v == 0.0 — BEV DGM NoData sentinel - // v.is_nan() — IEEE NaN from corrupt or partial tiles - // v < -1000.0 — large-negative sentinel (SRTM-style, defensive guard) - data[dst + i] = if v == 0.0 || v.is_nan() || v < -1000.0 { + // can be zero. The conditions cover: + // v == 0.0 — BEV DGM NoData sentinel + // is_nodata_value() — NaN, large-negative sentinel, float-max sentinel + // family, and the file's declared GDAL_NODATA value + data[dst + i] = if v == 0.0 || is_nodata_value(v, gdal_nodata) { NODATA } else { v diff --git a/crates/dem_io/src/overview.rs b/crates/dem_io/src/overview.rs index 4eb923e..03c28ff 100644 --- a/crates/dem_io/src/overview.rs +++ b/crates/dem_io/src/overview.rs @@ -48,6 +48,7 @@ pub(crate) fn build_downsampled( let src_dy = scale_tag[1]; let crs_origin_x = tiepoint[3]; let crs_origin_y = tiepoint[4]; + let gdal_nodata = crate::geotiff::read_gdal_nodata(&mut decoder); let is_geo = crs::is_geographic(&proj4); let (origin_lat, origin_lon) = if is_geo { @@ -137,7 +138,7 @@ pub(crate) fn build_downsampled( for local_c in 0..(tile_col1 - tile_col0) { let gc = tile_col0 + local_c; let v = tile_data[local_r * actual_tw + local_c]; - if v.is_nan() || v < -1000.0 { + if crate::geotiff::is_nodata_value(v, gdal_nodata) { continue; } // Scatter this pixel into every downsample level simultaneously. diff --git a/crates/dem_io/tests/synthetic.rs b/crates/dem_io/tests/synthetic.rs index c46c73f..88c25bb 100644 --- a/crates/dem_io/tests/synthetic.rs +++ b/crates/dem_io/tests/synthetic.rs @@ -42,6 +42,40 @@ fn write_min_geotiff( img.write_data(data).unwrap(); } +/// Like `write_min_geotiff` but also writes the GDAL_NODATA ASCII tag (42113), +/// the way GDAL/opals declare per-file sentinels (e.g. the BEV ALS DSM's +3.4e38). +#[allow(clippy::too_many_arguments)] +fn write_min_geotiff_with_nodata( + path: &Path, + cols: u32, + rows: u32, + scale: (f64, f64), + origin: (f64, f64), + geo_key_dir: &[u16], + data: &[f32], + nodata: &str, +) { + let file = File::create(path).unwrap(); + let mut enc = TiffEncoder::new(BufWriter::new(file)).unwrap(); + let mut img = enc.new_image::(cols, rows).unwrap(); + img.encoder() + .write_tag(Tag::Unknown(33550), &[scale.0, scale.1, 0.0_f64][..]) + .unwrap(); + img.encoder() + .write_tag( + Tag::Unknown(33922), + &[0.0_f64, 0.0, 0.0, origin.0, origin.1, 0.0][..], + ) + .unwrap(); + img.encoder() + .write_tag(Tag::Unknown(34735), geo_key_dir) + .unwrap(); + img.encoder() + .write_tag(Tag::Unknown(42113), nodata) + .unwrap(); + img.write_data(data).unwrap(); +} + /// GeoKeyDirectory for a projected CRS given by EPSG (3072). fn projected_dir(epsg: u16) -> Vec { // header: version 1.1.0, 2 keys; GTModelType=1 (projected); ProjectedCSType=epsg @@ -153,6 +187,64 @@ fn overview_cache_round_trip_preserves_crs() { assert_eq!(src_tags.geo_ascii_params, cache_tags.geo_ascii_params); } +#[test] +fn float_max_sentinel_maps_to_nodata() { + // The BEV ALS DSM declares NoData = +3.4028235e38 via GDAL_NODATA. Before the + // generalized predicate, those cells parsed as valid 3.4e38 m elevations. + let dir = tempfile::tempdir().unwrap(); + let src = dir.path().join("dsm.tif"); + let mut data = ramp(8, 8); + data[0] = 3.402_823_5e38; + data[27] = 3.402_823_5e38; + write_min_geotiff_with_nodata( + &src, + 8, + 8, + (1.0, 1.0), + (4_450_000.0, 2_700_000.0), + &projected_dir(3035), + &data, + "3.4028235e+38", + ); + + // extract_window does not infill, so sentinels must surface as exact -9999. + let win = dem_io::extract_window(&src, (4_450_004.0, 2_699_996.0), 10.0, 0).unwrap(); + assert_eq!((win.cols, win.rows), (8, 8)); + assert_eq!(win.data[0], -9999.0, "float-max sentinel → NODATA"); + assert_eq!(win.data[27], -9999.0); + assert!( + win.data.iter().all(|&v| v < 1.0e38), + "no float-max value may survive into the heightmap" + ); + assert_eq!(win.data[1], data[1], "valid neighbours pass through"); +} + +#[test] +fn declared_gdal_nodata_value_maps_to_nodata() { + // A finite sentinel inside the "plausible elevation" range is only catchable + // via the GDAL_NODATA tag — neither the NaN nor the ±magnitude checks fire. + let dir = tempfile::tempdir().unwrap(); + let src = dir.path().join("custom_nodata.tif"); + let mut data = ramp(8, 8); + data[5] = -500.0; + data[42] = -500.0; + write_min_geotiff_with_nodata( + &src, + 8, + 8, + (1.0, 1.0), + (4_450_000.0, 2_700_000.0), + &projected_dir(3035), + &data, + "-500", + ); + + let win = dem_io::extract_window(&src, (4_450_004.0, 2_699_996.0), 10.0, 0).unwrap(); + assert_eq!(win.data[5], -9999.0, "tag-declared sentinel → NODATA"); + assert_eq!(win.data[42], -9999.0); + assert_eq!(win.data[6], data[6], "valid neighbours pass through"); +} + #[test] fn no_cache_built_for_coarse_single_ifd_source() { // A ≥20 m/px source can be served directly by base-tier select_ifd, so From e2c6f4caa332442e5c142e0618d273e594ba174a Mon Sep 17 00:00:00 2001 From: Nikita Zavartsev Date: Fri, 12 Jun 2026 20:43:00 +0200 Subject: [PATCH 2/9] feat: [dsm-rgb-landcover-tirol] add DSM-over-DTM surface compositing to dem_io composite_surface_over paints the valid cells of a DSM window onto the fine-tier DTM canvas with a short BFS feather at the ALS coverage and window boundaries. The DTM stays the canvas so fine-tier coverage never shrinks where the DSM window is clipped, and real canopy/building heights are never sanded down the way fill_nodata_from_base's 30 px outward blend would. Canvas sentinel cells take the overlay at full weight so nothing ever blends toward -9999. Co-Authored-By: Claude Fable 5 --- crates/dem_io/src/heightmap.rs | 192 +++++++++++++++++++++++++++++++++ crates/dem_io/src/lib.rs | 4 +- 2 files changed, 195 insertions(+), 1 deletion(-) diff --git a/crates/dem_io/src/heightmap.rs b/crates/dem_io/src/heightmap.rs index ff64fea..123a0a0 100644 --- a/crates/dem_io/src/heightmap.rs +++ b/crates/dem_io/src/heightmap.rs @@ -439,6 +439,126 @@ pub fn fill_nodata_from_base(hm: &mut Heightmap, base: &Heightmap) { } } +/// Paint the valid cells of `overlay` (a DSM window) onto `canvas` (the fine +/// DTM window), feathering `feather_px` pixels at the overlay's validity and +/// window boundaries so canopy/building walls don't produce a hard step where +/// the ALS coverage ends. +/// +/// Direction matters: the canvas keeps its full extent, so fine-tier coverage +/// never shrinks where the DSM window is smaller (it is clipped to its tile +/// bounds near coverage edges). This is deliberately NOT `fill_nodata_from_base` +/// with swapped roles — that function's 30 px outward blend would sand real +/// canopy detail toward bare earth along every coverage edge. +/// +/// Both rasters must share the same projected CRS and pixel scale (the DSM and +/// the 1 m DTM tiles are both EPSG:3035 on the same metre lattice); anything +/// else is a no-op with a warning since cross-CRS resampling is not this +/// function's job. +pub fn composite_surface_over(canvas: &mut Heightmap, overlay: &Heightmap, feather_px: usize) { + use std::collections::VecDeque; + + if canvas.rows == 0 || canvas.cols == 0 || overlay.rows == 0 || overlay.cols == 0 { + return; + } + if canvas.crs_proj4 != overlay.crs_proj4 { + eprintln!("[surface] overlay CRS differs from canvas — skipping composite"); + return; + } + let dx = canvas.dx_meters; + let dy = canvas.dy_meters; + if (overlay.dx_meters - dx).abs() > 0.01 * dx || (overlay.dy_meters - dy).abs() > 0.01 * dy { + eprintln!( + "[surface] overlay scale {}×{} differs from canvas {}×{} — skipping composite", + overlay.dx_meters, overlay.dy_meters, dx, dy + ); + return; + } + + // Integer pixel offset of the overlay's top-left corner within the canvas. + // Both windows come from extract_window on tiles sharing one metre lattice, + // so the offset is integral up to float noise. + let off_c = ((overlay.crs_origin_x - canvas.crs_origin_x) / dx).round() as isize; + let off_r = ((canvas.crs_origin_y - overlay.crs_origin_y) / dy).round() as isize; + + let n = overlay.rows * overlay.cols; + let valid: Vec = overlay.data.iter().map(|&v| v > -1000.0).collect(); + + // Multi-source BFS: distance (in pixels) from each valid overlay cell to the + // nearest invalid cell or window edge. dist saturates at feather_px — beyond + // that the overlay applies at full weight. + let mut dist = vec![u32::MAX; n]; + let mut queue: VecDeque = VecDeque::new(); + for r in 0..overlay.rows { + for c in 0..overlay.cols { + let i = r * overlay.cols + c; + if !valid[i] { + continue; + } + let at_edge = r == 0 || c == 0 || r == overlay.rows - 1 || c == overlay.cols - 1; + let near_invalid = [(-1i32, 0i32), (1, 0), (0, -1), (0, 1)].iter().any(|&(dr, dc)| { + let nr = r as i32 + dr; + let nc = c as i32 + dc; + nr >= 0 + && nr < overlay.rows as i32 + && nc >= 0 + && nc < overlay.cols as i32 + && !valid[nr as usize * overlay.cols + nc as usize] + }); + if at_edge || near_invalid { + dist[i] = 1; + queue.push_back(i); + } + } + } + while let Some(i) = queue.pop_front() { + let d = dist[i]; + if d as usize >= feather_px { + continue; + } + let r = i / overlay.cols; + let c = i % overlay.cols; + for &(dr, dc) in &[(-1i32, 0i32), (1, 0), (0, -1), (0, 1)] { + let nr = r as i32 + dr; + let nc = c as i32 + dc; + if nr < 0 || nr >= overlay.rows as i32 || nc < 0 || nc >= overlay.cols as i32 { + continue; + } + let ni = nr as usize * overlay.cols + nc as usize; + if !valid[ni] || dist[ni] != u32::MAX { + continue; + } + dist[ni] = d + 1; + queue.push_back(ni); + } + } + + for r in 0..overlay.rows { + let cr = r as isize + off_r; + if cr < 0 || cr >= canvas.rows as isize { + continue; + } + for c in 0..overlay.cols { + let i = r * overlay.cols + c; + if !valid[i] { + continue; + } + let cc = c as isize + off_c; + if cc < 0 || cc >= canvas.cols as isize { + continue; + } + let ci = cr as usize * canvas.cols + cc as usize; + let under = canvas.data[ci]; + // A canvas sentinel must never be blended toward — take the surface. + let t = if under <= -1000.0 || feather_px == 0 { + 1.0 + } else { + (dist[i].min(feather_px as u32) as f32) / feather_px as f32 + }; + canvas.data[ci] = under * (1.0 - t) + overlay.data[i] * t; + } + } +} + fn build_grayscale_png(heightmap: &Heightmap, cols: usize, rows: usize) { let min = heightmap.data.iter().cloned().fold(f32::INFINITY, f32::min); let max = heightmap @@ -736,6 +856,78 @@ mod tests { } } + // composite_surface_over + + /// Projected 1 m/px heightmap with a CRS string so the same-CRS check passes. + fn mk_proj(rows: usize, cols: usize, ox: f64, oy: f64, data: Vec) -> Heightmap { + let mut hm = mk_hm(rows, cols, data); + hm.crs_origin_x = ox; + hm.crs_origin_y = oy; + hm.crs_proj4 = "+proj=laea +lat_0=52 +lon_0=10".to_string(); + hm + } + + #[test] + fn composite_paints_overlay_with_full_weight_past_feather() { + // 9×9 canvas at 100 m, fully-valid 9×9 overlay at 130 m, feather 2. + // The centre cell is ≥2 px from every overlay edge → full overlay value. + // An edge cell is at distance 1 → 50/50 blend. + let mut canvas = mk_proj(9, 9, 0.0, 9.0, vec![100.0; 81]); + let overlay = mk_proj(9, 9, 0.0, 9.0, vec![130.0; 81]); + composite_surface_over(&mut canvas, &overlay, 2); + assert_eq!(canvas.data[4 * 9 + 4], 130.0, "interior takes full surface"); + assert_eq!(canvas.data[0], 115.0, "window-edge cell blends 50/50"); + } + + #[test] + fn composite_respects_pixel_offset() { + // Overlay's top-left sits 2 px east / 1 px south of the canvas origin. + let mut canvas = mk_proj(6, 6, 0.0, 6.0, vec![100.0; 36]); + let overlay = mk_proj(3, 3, 2.0, 5.0, vec![200.0; 9]); + composite_surface_over(&mut canvas, &overlay, 0); + assert_eq!(canvas.data[2 * 6 + 3], 200.0, "overlay centre lands at (r2,c3)"); + assert_eq!(canvas.data[0], 100.0, "outside overlay untouched"); + assert_eq!(canvas.data[5 * 6 + 5], 100.0); + } + + #[test] + fn composite_skips_overlay_nodata_holes_and_feathers_around_them() { + // Fully valid except a centre sentinel; feather 1 ⇒ all valid cells are + // boundary cells (dist 1 == feather) so they still get full weight, while + // the hole itself keeps the canvas value. + let mut canvas = mk_proj(3, 3, 0.0, 3.0, vec![100.0; 9]); + let mut over = vec![150.0; 9]; + over[4] = -9999.0; + let overlay = mk_proj(3, 3, 0.0, 3.0, over); + composite_surface_over(&mut canvas, &overlay, 1); + assert_eq!(canvas.data[4], 100.0, "hole keeps the DTM value"); + assert_eq!(canvas.data[0], 150.0, "valid cells at dist=feather get full weight"); + } + + #[test] + fn composite_never_blends_toward_canvas_sentinel() { + // Canvas cell is the NODATA sentinel: even at the feather boundary the + // overlay must apply fully instead of mixing in -9999 (a crater). + let mut canvas = mk_proj(3, 3, 0.0, 3.0, vec![-9999.0; 9]); + let overlay = mk_proj(3, 3, 0.0, 3.0, vec![140.0; 9]); + composite_surface_over(&mut canvas, &overlay, 4); + assert!(canvas.data.iter().all(|&v| v == 140.0), "{:?}", canvas.data); + } + + #[test] + fn composite_requires_matching_crs_and_scale() { + let mut canvas = mk_proj(3, 3, 0.0, 3.0, vec![100.0; 9]); + let mut overlay = mk_proj(3, 3, 0.0, 3.0, vec![150.0; 9]); + overlay.crs_proj4 = "+proj=utm +zone=33".to_string(); + composite_surface_over(&mut canvas, &overlay, 1); + assert!(canvas.data.iter().all(|&v| v == 100.0), "CRS mismatch → no-op"); + + let mut overlay2 = mk_proj(3, 3, 0.0, 3.0, vec![150.0; 9]); + overlay2.dx_meters = 2.0; + composite_surface_over(&mut canvas, &overlay2, 1); + assert!(canvas.data.iter().all(|&v| v == 100.0), "scale mismatch → no-op"); + } + #[test] fn fill_from_base_empty_inputs_early_return() { let base = mk_geo(2, 2, 0.0, 2.0, vec![1.0; 4]); diff --git a/crates/dem_io/src/lib.rs b/crates/dem_io/src/lib.rs index 3c75bd9..57099bf 100644 --- a/crates/dem_io/src/lib.rs +++ b/crates/dem_io/src/lib.rs @@ -13,7 +13,9 @@ pub use geotiff::{ pub use grid::{ assemble_grid, crop, load_grid_from_paths, stitch_windows, stitch_windows_geographic, }; -pub use heightmap::{Heightmap, clamp_nodata_to_sea, fill_nodata_from_base, parse_bil}; +pub use heightmap::{ + Heightmap, clamp_nodata_to_sea, composite_surface_over, fill_nodata_from_base, parse_bil, +}; pub use overview::{BASE_OVERVIEW_TARGET_M, CLOSE_OVERVIEW_TARGET_M, ensure_overview_cache}; pub(crate) type DemError = Box; From e573c9b8c9200e1ade61e661b8fc907c9c8fba40 Mon Sep 17 00:00:00 2001 From: Nikita Zavartsev Date: Fri, 12 Jun 2026 20:43:24 +0200 Subject: [PATCH 3/9] feat: [dsm-rgb-landcover-tirol] add ortho RGB + land-cover color window reader to dem_io MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit extract_color_window reads a camera-centred window from the BEV 20 cm orthophoto mosaic (ModernJPEG YCbCr chunks decoded via the tiff crate's zune path, BT.601-converted on the caller's thread) and nearest-resamples the land-cover raster onto the same grid, packing one RGBA8 buffer: RGB = albedo, A = material code ladder (water 255 / high veg 192 / med veg 128 / building 64 / other 0 — codes spaced >= 64 apart so bilinear boundary mixes degrade gracefully). The LC mosaic is Gray(4): 4-bit nibble-packed chunks are unpacked in extract_u8_window. ifd_overview_levels walks the IFD chain filtered by NewSubfileType so the RGB mosaic's interleaved transparency-mask IFDs don't shift the overview indices the way raw ifd_scales would. LC class codes pinned empirically with the new inspect_color example against single-class reference spots (Achensee = 100% class 5 water, Hintertux glacier = 100% class 2 ground/ice, forest 96% class 1, meadow 88% class 6, town 44% class 4 buildings). Gated integration tests run against the local multi-GB tiles and skip cleanly elsewhere. Co-Authored-By: Claude Fable 5 --- crates/dem_io/examples/inspect_color.rs | 128 +++++++++ crates/dem_io/src/color.rs | 364 ++++++++++++++++++++++++ crates/dem_io/src/geotiff.rs | 36 +++ crates/dem_io/src/lib.rs | 9 +- crates/dem_io/tests/local_big_tiles.rs | 118 ++++++++ 5 files changed, 653 insertions(+), 2 deletions(-) create mode 100644 crates/dem_io/examples/inspect_color.rs create mode 100644 crates/dem_io/src/color.rs create mode 100644 crates/dem_io/tests/local_big_tiles.rs diff --git a/crates/dem_io/examples/inspect_color.rs b/crates/dem_io/examples/inspect_color.rs new file mode 100644 index 0000000..b4b2f5b --- /dev/null +++ b/crates/dem_io/examples/inspect_color.rs @@ -0,0 +1,128 @@ +//! Inspect a land-cover (and optionally orthophoto) mosaic: dump the mask-aware +//! overview levels and class-value histograms at known Tirol reference spots +//! (reservoir / forest slope / town / glacier). Used to pin the BEV land-cover +//! class → material mapping in `dem_io::color::lc_code_to_material` — the class +//! codes are not documented in the GeoTIFF itself. +//! +//! Usage: +//! cargo run --release -p dem_io --example inspect_color -- tiles/color/2022470_Mosaik_LC.tif [tiles/color/2019470_Mosaik_RGB.tif] + +use std::env; +use std::path::Path; + +fn main() { + let lc_arg = env::args() + .nth(1) + .expect("usage: inspect_color [ortho_rgb.tif]"); + let lc_path = Path::new(&lc_arg).to_path_buf(); + let rgb_arg = env::args().nth(2); + + let levels = dem_io::ifd_overview_levels(&lc_path).expect("walk IFDs"); + println!("land cover overview levels (ifd, scale):"); + for (ifd, scale) in &levels { + println!(" ifd {ifd:>2} {scale:.2} m/px"); + } + + let proj4 = dem_io::crs::tile_proj4(&lc_path).expect("resolve CRS"); + println!("proj4: {proj4}"); + + // (label, lat, lon, radius_m) — picked so each window is dominated by one class. + let spots: &[(&str, f64, f64, f64)] = &[ + ("Achensee centre (deep lake)", 47.4450, 11.7080, 250.0), + ("Speicher Schlegeis (reservoir)", 47.0285, 11.7090, 250.0), + ("forest slope near Mayrhofen", 47.1500, 11.8400, 150.0), + ("valley meadow near Uderns", 47.3120, 11.8740, 100.0), + ("Mayrhofen town centre", 47.1660, 11.8630, 150.0), + ("Hintertux glacier / rock", 47.0680, 11.6730, 250.0), + ]; + + // Use a mid overview (~1.6 m) so each window is a few hundred px — enough for + // a stable histogram without decoding thousands of full-res chunks. + let (ifd, scale) = levels + .iter() + .copied() + .find(|&(_, s)| s >= 1.5) + .unwrap_or(levels[0]); + println!("\nhistograms at ifd {ifd} ({scale:.2} m/px):"); + + for &(label, lat, lon, radius) in spots { + let Ok(centre) = dem_io::crs::from_wgs84(lat, lon, &proj4) else { + println!(" {label}: from_wgs84 failed"); + continue; + }; + match dem_io::landcover_histogram(&lc_path, centre, radius, ifd) { + Ok(hist) => { + let total: u64 = hist.iter().sum(); + print!(" {label}: "); + for (v, &n) in hist.iter().enumerate() { + if n > 0 { + print!("{v}:{:.0}% ", n as f64 / total as f64 * 100.0); + } + } + println!(); + } + Err(e) => println!(" {label}: {e}"), + } + } + + if let Some(rgb_arg) = rgb_arg { + let rgb_path = Path::new(&rgb_arg).to_path_buf(); + let rgb_levels = dem_io::ifd_overview_levels(&rgb_path).expect("walk RGB IFDs"); + println!("\northo overview levels (ifd, scale):"); + for (ifd, scale) in &rgb_levels { + println!(" ifd {ifd:>2} {scale:.2} m/px"); + } + let (rgb_ifd, rgb_scale) = rgb_levels + .iter() + .copied() + .find(|&(_, s)| s >= 1.5) + .unwrap_or(rgb_levels[0]); + let lc_ifd = levels + .iter() + .copied() + .find(|&(_, s)| (s - rgb_scale).abs() < 0.01) + .map(|(i, _)| i); + println!("sampling ortho at ifd {rgb_ifd} ({rgb_scale:.2} m/px), lc ifd {lc_ifd:?}"); + + for &(label, lat, lon, radius) in spots { + let proj4_rgb = dem_io::crs::tile_proj4(&rgb_path).expect("resolve RGB CRS"); + let Ok(centre) = dem_io::crs::from_wgs84(lat, lon, &proj4_rgb) else { + continue; + }; + match dem_io::extract_color_window( + &rgb_path, + Some(&lc_path), + centre, + radius, + rgb_ifd, + lc_ifd, + ) { + Ok(win) => { + let n = win.georef.rows * win.georef.cols; + let c = n / 2; // centre-ish pixel + let px = &win.rgba[c * 4..c * 4 + 4]; + let avg: [u64; 4] = win.rgba.chunks_exact(4).fold([0; 4], |mut a, p| { + for i in 0..4 { + a[i] += p[i] as u64; + } + a + }); + println!( + " {label}: {}×{} centre rgba=({},{},{},{}) avg rgba=({},{},{},{})", + win.georef.cols, + win.georef.rows, + px[0], + px[1], + px[2], + px[3], + avg[0] / n as u64, + avg[1] / n as u64, + avg[2] / n as u64, + avg[3] / n as u64, + ); + } + Err(e) => println!(" {label}: {e}"), + } + } + } +} diff --git a/crates/dem_io/src/color.rs b/crates/dem_io/src/color.rs new file mode 100644 index 0000000..c471f62 --- /dev/null +++ b/crates/dem_io/src/color.rs @@ -0,0 +1,364 @@ +//! Camera-centred windows from orthophoto + land-cover mosaics. +//! +//! Reads a window from a 3-band JPEG-in-TIFF orthophoto (BEV mosaics store +//! YCbCr; the tiff crate's zune-jpeg path returns the samples unconverted, so +//! the BT.601 → RGB step happens here) and optionally drapes a single-band +//! land-cover raster over the same CRS rect, nearest-sampled onto the ortho +//! grid. Output is one RGBA8 buffer: RGB = albedo, A = material code — the +//! shader keys water/vegetation shading off A while a single texture binding +//! carries everything. + +use std::fs::File; +use std::path::Path; + +use tiff::ColorType; +use tiff::decoder::{Decoder, DecodingResult, Limits}; +use tiff::tags::Tag; + +use crate::crs; +use crate::{DemError, Heightmap}; + +/// Material codes packed into the alpha channel. Spaced apart so bilinear +/// filtering at class boundaries degrades into the neighbouring band instead +/// of aliasing to an unrelated material. +pub const MATERIAL_NONE: u8 = 0; +pub const MATERIAL_BUILDING: u8 = 64; +pub const MATERIAL_MED_VEG: u8 = 128; +pub const MATERIAL_HIGH_VEG: u8 = 192; +pub const MATERIAL_WATER: u8 = 255; + +/// BEV Land Cover class → material code. +/// +/// Class values pinned empirically with `cargo run -p dem_io --example +/// inspect_color` against `2022470_Mosaik_LC.tif` at single-class reference +/// spots: Achensee deep centre = 100% class 5 (water), Hintertux glacier = +/// 100% class 2 (ground/ice), Mayrhofen forest slope = 96% class 1 (high +/// vegetation), Uderns valley meadow = 88% class 6 (low vegetation), Mayrhofen +/// town = 44% class 4 (buildings). Class 3 is the remainder (medium +/// vegetation); 15 = NoData. +fn lc_code_to_material(code: u8) -> u8 { + match code { + 1 => MATERIAL_HIGH_VEG, + 3 => MATERIAL_MED_VEG, + 4 => MATERIAL_BUILDING, + 5 => MATERIAL_WATER, + _ => MATERIAL_NONE, // 2 = ground/ice, 6 = low veg/grass, 15 = NoData + } +} + +/// Result of `extract_color_window`. +pub struct ColorWindow { + /// RGBA8, row-major: RGB = orthophoto albedo, A = material code. + pub rgba: Vec, + /// Georeferencing carrier for the rgba grid. `data` is intentionally empty — + /// rows/cols/dx/dy/crs_* describe the window so the viewer can reuse the + /// existing `Heightmap`-based placement code (`cross_crs_world_origin_and_extent` + /// never touches `.data`). + pub georef: Heightmap, +} + +/// Raw single- or multi-sample u8 window read straight from one IFD level. +struct U8Window { + data: Vec, // interleaved, `samples` per pixel + samples: usize, + cols: usize, + rows: usize, + dx: f64, + dy: f64, + origin_x: f64, // CRS coordinate of the window's top-left corner + origin_y: f64, + proj4: String, + epsg: u32, + color: ColorType, +} + +/// Window-read core shared by the ortho and land-cover paths. Same geometry +/// contract as `geotiff::extract_window` (centre + radius in CRS units, clipped +/// to tile bounds, selective chunk reads at `ifd_level`) but for u8 rasters. +/// JPEG/LZW/Deflate chunks all go through the tiff crate's `read_chunk`. +fn extract_u8_window( + path: &Path, + centre_crs: (f64, f64), + radius_m: f64, + ifd_level: usize, +) -> Result { + let crs_data = crs::read_geo_key_data(path)?; + let proj4 = crs::proj4_from_keys(&crs_data)?; + let epsg = crs_data + .projected_epsg + .or(crs_data.geographic_epsg) + .unwrap_or(0); + + let file = File::open(path)?; + let mut decoder = Decoder::new(std::io::BufReader::new(file))?.with_limits(Limits::unlimited()); + + // Geo-tags live on IFD 0 only. + decoder.seek_to_image(0)?; + let (full_cols, full_rows): (u32, u32) = decoder.dimensions()?; + let scale = decoder.get_tag(Tag::Unknown(33550))?.into_f64_vec()?; + let tiepoint = decoder.get_tag(Tag::Unknown(33922))?.into_f64_vec()?; + let crs_origin_x = tiepoint[3]; + let crs_origin_y = tiepoint[4]; + + decoder.seek_to_image(ifd_level)?; + let (cols, rows): (u32, u32) = decoder.dimensions()?; + let dx = scale[0] * (full_cols as f64 / cols as f64); + let dy = scale[1] * (full_rows as f64 / rows as f64); + + let color = decoder.colortype()?; + // (samples per pixel, bits per sample). 4-bit grayscale is what BEV's + // land-cover mosaic uses (6 classes + NoData fit a nibble); the tiff crate + // returns those chunks bit-packed two pixels per byte, MSB nibble first. + let (samples, bits) = match color { + ColorType::Gray(8) => (1, 8), + ColorType::Gray(4) => (1, 4), + ColorType::RGB(8) | ColorType::YCbCr(8) => (3, 8), + ColorType::RGBA(8) => (4, 8), + other => return Err(format!("unsupported color raster type {other:?}").into()), + }; + + let cx = (centre_crs.0 - crs_origin_x) / dx; + let cy = (crs_origin_y - centre_crs.1) / dy; + let radius_px_x = (radius_m / dx) as isize; + let radius_px_y = (radius_m / dy) as isize; + let px0 = (cx as isize - radius_px_x).max(0); + let px1 = (cx as isize + radius_px_x).min(cols as isize); + let py0 = (cy as isize - radius_px_y).max(0); + let py1 = (cy as isize + radius_px_y).min(rows as isize); + if px1 <= px0 || py1 <= py0 { + return Err("centre is outside tile bounds".into()); + } + let (px0, px1, py0, py1) = (px0 as usize, px1 as usize, py0 as usize, py1 as usize); + let out_w = px1 - px0; + let out_h = py1 - py0; + + let mut data = vec![0u8; out_w * out_h * samples]; + + let (tw, th) = decoder.chunk_dimensions(); + let tiles_across = (cols as usize).div_ceil(tw as usize); + let tc0 = px0 / tw as usize; + let tc1 = px1.div_ceil(tw as usize); + let tr0 = py0 / th as usize; + let tr1 = py1.div_ceil(th as usize); + + for tr in tr0..tr1 { + for tc in tc0..tc1 { + let index = (tr * tiles_across + tc) as u32; + let tile_col0 = tc * tw as usize; + let tile_row0 = tr * th as usize; + let tile_col1 = (tile_col0 + tw as usize).min(cols as usize); + let tile_row1 = (tile_row0 + th as usize).min(rows as usize); + // read_chunk trims edge tiles to their actual data dimensions. + let actual_tw = tile_col1 - tile_col0; + let actual_th = tile_row1 - tile_row0; + + let tile_data: Vec = match decoder.read_chunk(index)? { + DecodingResult::U8(v) => v, + _ => return Err("expected U8 chunk in color raster".into()), + }; + // Unpack 4-bit chunks to one sample per byte so the copy loop below + // can stay sample-addressed. Rows are byte-aligned in the packed form. + let tile_data: Vec = if bits == 4 { + let row_bytes = actual_tw.div_ceil(2); + let mut unpacked = vec![0u8; actual_tw * actual_th]; + for r in 0..actual_th { + for c in 0..actual_tw { + let byte = tile_data[r * row_bytes + c / 2]; + unpacked[r * actual_tw + c] = + if c % 2 == 0 { byte >> 4 } else { byte & 0x0F }; + } + } + unpacked + } else { + tile_data + }; + + let col_start = tile_col0.max(px0); + let col_end = tile_col1.min(px1); + let row_start = tile_row0.max(py0); + let row_end = tile_row1.min(py1); + + for row in row_start..row_end { + let src = ((row - tile_row0) * actual_tw + (col_start - tile_col0)) * samples; + let dst = ((row - py0) * out_w + (col_start - px0)) * samples; + let len = (col_end - col_start) * samples; + data[dst..dst + len].copy_from_slice(&tile_data[src..src + len]); + } + } + } + + Ok(U8Window { + data, + samples, + cols: out_w, + rows: out_h, + dx, + dy, + origin_x: crs_origin_x + px0 as f64 * dx, + origin_y: crs_origin_y - py0 as f64 * dy, + proj4, + epsg, + color, + }) +} + +/// BT.601 full-range YCbCr → RGB (the JPEG convention). +#[inline] +pub(crate) fn ycbcr_to_rgb(y: u8, cb: u8, cr: u8) -> [u8; 3] { + let y = y as f32; + let cb = cb as f32 - 128.0; + let cr = cr as f32 - 128.0; + let r = y + 1.402 * cr; + let g = y - 0.344_136 * cb - 0.714_136 * cr; + let b = y + 1.772 * cb; + [ + r.clamp(0.0, 255.0) as u8, + g.clamp(0.0, 255.0) as u8, + b.clamp(0.0, 255.0) as u8, + ] +} + +/// Extract a camera-centred RGBA8 window: RGB albedo from `rgb_path` at +/// `rgb_ifd`, material codes in A from `lc_path` at `lc_ifd` (nearest-sampled +/// onto the ortho grid by CRS coordinate; A = 0 everywhere when no land cover +/// is supplied or a pixel falls outside its coverage). +/// +/// `centre_crs` is in the ortho file's CRS; callers convert from WGS84 via +/// `crs::from_wgs84` exactly like the height-tier workers do. +pub fn extract_color_window( + rgb_path: &Path, + lc_path: Option<&Path>, + centre_crs: (f64, f64), + radius_m: f64, + rgb_ifd: usize, + lc_ifd: Option, +) -> Result { + let rgb = extract_u8_window(rgb_path, centre_crs, radius_m, rgb_ifd)?; + if rgb.samples != 3 { + return Err(format!("ortho must be 3-band, got {} samples", rgb.samples).into()); + } + let needs_ycbcr = matches!(rgb.color, ColorType::YCbCr(_)); + + let n_px = rgb.cols * rgb.rows; + let mut rgba = vec![0u8; n_px * 4]; + for i in 0..n_px { + let s = &rgb.data[i * 3..i * 3 + 3]; + let [r, g, b] = if needs_ycbcr { + ycbcr_to_rgb(s[0], s[1], s[2]) + } else { + [s[0], s[1], s[2]] + }; + rgba[i * 4] = r; + rgba[i * 4 + 1] = g; + rgba[i * 4 + 2] = b; + // alpha stays MATERIAL_NONE until the land-cover pass below + } + + if let Some(lc_path) = lc_path { + match extract_u8_window(lc_path, centre_crs, radius_m, lc_ifd.unwrap_or(0)) { + Ok(lc) if lc.samples == 1 => { + for r in 0..rgb.rows { + let y = rgb.origin_y - (r as f64 + 0.5) * rgb.dy; + let lc_r = ((lc.origin_y - y) / lc.dy) as isize; + if lc_r < 0 || lc_r >= lc.rows as isize { + continue; + } + for c in 0..rgb.cols { + let x = rgb.origin_x + (c as f64 + 0.5) * rgb.dx; + let lc_c = ((x - lc.origin_x) / lc.dx) as isize; + if lc_c < 0 || lc_c >= lc.cols as isize { + continue; + } + let code = lc.data[lc_r as usize * lc.cols + lc_c as usize]; + rgba[(r * rgb.cols + c) * 4 + 3] = lc_code_to_material(code); + } + } + } + Ok(lc) => eprintln!( + "[color] land cover {} has {} samples/px, expected 1 — skipping", + lc_path.display(), + lc.samples + ), + Err(e) => eprintln!( + "[color] land cover window failed ({e}) — continuing without material codes" + ), + } + } + + Ok(ColorWindow { + rgba, + georef: Heightmap { + data: Vec::new(), + rows: rgb.rows, + cols: rgb.cols, + nodata: -9999.0, + origin_lat: 0.0, + origin_lon: 0.0, + dx_deg: 0.0, + dy_deg: 0.0, + dx_meters: rgb.dx, + dy_meters: rgb.dy, + crs_origin_x: rgb.origin_x, + crs_origin_y: rgb.origin_y, + crs_epsg: rgb.epsg, + crs_proj4: rgb.proj4, + }, + }) +} + +/// Histogram of land-cover class values over a window — backing for the +/// `inspect_color` example that pins `lc_code_to_material`. +pub fn landcover_histogram( + lc_path: &Path, + centre_crs: (f64, f64), + radius_m: f64, + ifd: usize, +) -> Result<[u64; 256], DemError> { + let win = extract_u8_window(lc_path, centre_crs, radius_m, ifd)?; + if win.samples != 1 { + return Err("land cover must be single-band".into()); + } + let mut hist = [0u64; 256]; + for &v in &win.data { + hist[v as usize] += 1; + } + Ok(hist) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn ycbcr_known_vectors() { + // Neutral grays: Cb = Cr = 128 must pass Y through unchanged. + assert_eq!(ycbcr_to_rgb(0, 128, 128), [0, 0, 0]); + assert_eq!(ycbcr_to_rgb(255, 128, 128), [255, 255, 255]); + assert_eq!(ycbcr_to_rgb(99, 128, 128), [99, 99, 99]); + // Primary red (255,0,0) encodes to Y≈76, Cb≈84.4, Cr=255 in BT.601. + let [r, g, b] = ycbcr_to_rgb(76, 84, 255); + assert!(r >= 250, "red channel saturates, got {r}"); + assert!(g <= 10 && b <= 10, "green/blue near zero, got {g}/{b}"); + // Primary blue (0,0,255) → Y≈29, Cb=255, Cr≈107. + let [r, g, b] = ycbcr_to_rgb(29, 255, 107); + assert!(b >= 250, "blue channel saturates, got {b}"); + assert!(r <= 10 && g <= 35, "red/green low, got {r}/{g}"); + } + + #[test] + fn material_ladder_is_widely_spaced() { + // Bilinear filtering mixes adjacent texels; codes must stay ≥ 64 apart + // so a 50/50 boundary mix never lands inside another class's band. + let codes = [ + MATERIAL_NONE, + MATERIAL_BUILDING, + MATERIAL_MED_VEG, + MATERIAL_HIGH_VEG, + MATERIAL_WATER, + ]; + for pair in codes.windows(2) { + assert!(pair[1] - pair[0] >= 63, "codes too close: {pair:?}"); + } + assert_eq!(lc_code_to_material(15), MATERIAL_NONE, "LC NoData → none"); + } +} diff --git a/crates/dem_io/src/geotiff.rs b/crates/dem_io/src/geotiff.rs index eb21bb3..c6b3d13 100644 --- a/crates/dem_io/src/geotiff.rs +++ b/crates/dem_io/src/geotiff.rs @@ -231,6 +231,42 @@ pub fn ifd_scales(path: &Path) -> Result, DemError> { Ok(scales) } +/// Like `ifd_scales`, but skips non-image IFDs and reports which IFD index each +/// scale belongs to. COGs with per-dataset transparency masks (e.g. BEV ortho +/// mosaics) interleave mask IFDs (NewSubfileType bit 2) with the overview chain — +/// walking by raw index would mis-pair scales with IFDs there. Returns +/// `(ifd_index, scale)` pairs, finest first; scale is in the unit of the pixel +/// scale tag (m/px projected, deg/px geographic). +pub fn ifd_overview_levels(path: &Path) -> Result, DemError> { + const MASK_BIT: u64 = 4; // NewSubfileType (tag 254) bit 2 = transparency mask + + let file = File::open(path)?; + let mut decoder = Decoder::new(std::io::BufReader::new(file))?.with_limits(Limits::unlimited()); + decoder.seek_to_image(0)?; + let (full_cols, _) = decoder.dimensions()?; + let base_scale = decoder.get_tag(Tag::Unknown(33550))?.into_f64_vec()?[0]; + + let mut levels = Vec::new(); + let mut ifd = 0usize; + loop { + if ifd > 0 && decoder.seek_to_image(ifd).is_err() { + break; + } + let subfile_type = decoder + .find_tag(Tag::NewSubfileType) + .ok() + .flatten() + .and_then(|v| v.into_u64().ok()) + .unwrap_or(0); + if subfile_type & MASK_BIT == 0 { + let (cols, _) = decoder.dimensions()?; + levels.push((ifd, base_scale * (full_cols as f64 / cols as f64))); + } + ifd += 1; + } + Ok(levels) +} + pub fn extract_window( path: &Path, centre_crs: (f64, f64), diff --git a/crates/dem_io/src/lib.rs b/crates/dem_io/src/lib.rs index 57099bf..39ac106 100644 --- a/crates/dem_io/src/lib.rs +++ b/crates/dem_io/src/lib.rs @@ -1,3 +1,4 @@ +mod color; pub mod crs; mod geotiff; mod grid; @@ -5,10 +6,14 @@ mod heightmap; mod lzw_lenient; mod overview; +pub use color::{ + ColorWindow, MATERIAL_BUILDING, MATERIAL_HIGH_VEG, MATERIAL_MED_VEG, MATERIAL_NONE, + MATERIAL_WATER, extract_color_window, landcover_histogram, +}; pub use crs::get_tile_epsg; pub use geotiff::{ - extract_window, geotiff_pixel_scale, ifd_scales, parse_geotiff_auto, tile_bounds_wgs84, - tile_centre_crs, + extract_window, geotiff_pixel_scale, ifd_overview_levels, ifd_scales, parse_geotiff_auto, + tile_bounds_wgs84, tile_centre_crs, }; pub use grid::{ assemble_grid, crop, load_grid_from_paths, stitch_windows, stitch_windows_geographic, diff --git a/crates/dem_io/tests/local_big_tiles.rs b/crates/dem_io/tests/local_big_tiles.rs new file mode 100644 index 0000000..e3c1062 --- /dev/null +++ b/crates/dem_io/tests/local_big_tiles.rs @@ -0,0 +1,118 @@ +//! Tests gated on the multi-GB local Tirol tiles (DSM / land cover / RGB ortho). +//! These files are not committed (tiles/ is gitignored); each test skips with a +//! note when its input is absent so CI and other machines stay green. + +use std::path::PathBuf; + +fn repo_tile(rel: &str) -> Option { + let p = PathBuf::from(env!("CARGO_MANIFEST_DIR")) + .join("../../tiles") + .join(rel); + if p.exists() { + Some(p) + } else { + eprintln!("skipping — {rel} not present"); + None + } +} + +#[test] +fn rgb_ortho_overview_levels_skip_mask_ifds() { + let Some(path) = repo_tile("color/2019470_Mosaik_RGB.tif") else { + return; + }; + let levels = dem_io::ifd_overview_levels(&path).expect("walk IFDs"); + + // gdalinfo: full res + 12 overviews = 13 image IFDs; the per-dataset mask + // IFDs interleaved in the chain must not appear. + assert_eq!(levels.len(), 13, "13 image levels, masks skipped: {levels:?}"); + assert_eq!(levels[0], (0, 0.2), "IFD 0 is the 0.2 m/px full resolution"); + for pair in levels.windows(2) { + let ratio = pair[1].1 / pair[0].1; + assert!( + (ratio - 2.0).abs() < 0.05, + "overview scales must roughly double: {levels:?}" + ); + } +} + +#[test] +fn landcover_overview_levels_have_no_masks_to_skip() { + let Some(path) = repo_tile("color/2022470_Mosaik_LC.tif") else { + return; + }; + let levels = dem_io::ifd_overview_levels(&path).expect("walk IFDs"); + // gdalinfo: full res + 10 overviews, no mask bands. + assert_eq!(levels.len(), 11, "11 image levels: {levels:?}"); + assert_eq!(levels[0], (0, 0.2)); + // Without masks the IFD indices must be contiguous. + for (i, &(ifd, _)) in levels.iter().enumerate() { + assert_eq!(ifd, i); + } +} + +#[test] +fn color_window_decodes_albedo_and_water_material() { + let Some(rgb) = repo_tile("color/2019470_Mosaik_RGB.tif") else { + return; + }; + let Some(lc) = repo_tile("color/2022470_Mosaik_LC.tif") else { + return; + }; + // Achensee centre: 100% water in the land cover, blue-ish ortho pixels. + let proj4 = dem_io::crs::tile_proj4(&rgb).expect("CRS"); + let centre = dem_io::crs::from_wgs84(47.4450, 11.7080, &proj4).expect("project"); + // ifd 4 = 1.6 m/px in the RGB pyramid (mask IFDs skipped), lc ifd 3 matches. + let win = dem_io::extract_color_window(&rgb, Some(&lc), centre, 250.0, 4, Some(3)) + .expect("color window"); + + let n = win.georef.rows * win.georef.cols; + assert_eq!(win.rgba.len(), n * 4); + assert!(n > 90_000, "expected ≈312² window, got {n}"); + + let water_frac = win + .rgba + .chunks_exact(4) + .filter(|p| p[3] == dem_io::MATERIAL_WATER) + .count() as f64 + / n as f64; + assert!( + water_frac > 0.95, + "lake centre must be nearly all water material, got {water_frac:.2}" + ); + + // The YCbCr→RGB path must produce chroma, not grayscale: over a lake the + // blue channel average exceeds red by a clear margin. + let (mut r_sum, mut b_sum) = (0u64, 0u64); + for p in win.rgba.chunks_exact(4) { + r_sum += p[0] as u64; + b_sum += p[2] as u64; + } + assert!( + b_sum > r_sum + (n as u64 * 5), + "lake should be blue-tinted: avg r={} b={}", + r_sum / n as u64, + b_sum / n as u64 + ); +} + +#[test] +fn dsm_window_maps_float_max_nodata() { + let Some(path) = repo_tile("big_size/ALS_DSM_CRS3035RES50000mN2650000E4450000.tif") else { + return; + }; + // Window at the DSM's western edge (EPSG:3035) — guaranteed to include both + // valid surface pixels and +3.4e38 NoData from outside the ALS footprint. + let win = dem_io::extract_window(&path, (4_450_200.0, 2_663_000.0), 400.0, 0).expect("window"); + assert!( + win.data.iter().all(|&v| v < 1.0e38), + "no float-max sentinel may survive" + ); + let valid: Vec = win.data.iter().copied().filter(|&v| v > -9000.0).collect(); + assert!(!valid.is_empty(), "window must contain valid DSM surface"); + let max = valid.iter().cloned().fold(f32::MIN, f32::max); + assert!( + (400.0..4000.0).contains(&max), + "plausible Tirol surface height, got {max}" + ); +} From 46cc7c2286ca9b46403af6609c6672d682f1ebcb Mon Sep 17 00:00:00 2001 From: Nikita Zavartsev Date: Fri, 12 Jun 2026 20:44:14 +0200 Subject: [PATCH 4/9] feat: [dsm-rgb-landcover-tirol] add ortho albedo textures, uniforms and shader draping to render_gpu Two mipmapped Rgba8Unorm windows (fine + close) join the canonical bind group as entries 20-23, with their own origin/extent/rotation uniform blocks (CameraUniforms grows 224 -> 288 bytes; layout test pins the new offsets). upload_ortho_fine/close follow the established drop-first reload cycle (placeholder swap -> rebuild_bind_group -> poll(Wait) -> tracked alloc) and set_ortho_*_inactive eagerly frees the textures. gen_rgba_mip_bytes box-filters RGB but takes nearest alpha so the land-cover material codes survive the mip chain. Deliberately non-sRGB: the whole shader works in display-space 0-255. Shader: hit-point albedo = procedural ramp -> close ortho -> fine ortho with the same edge-smoothstep + per-window rotation as the height tiers and the existing distance-heuristic mip. Where ortho wins, the brightness floor is lifted (the photo carries baked sun shading; full diffuse on top double-shades north faces). Water material (alpha >= ~0.9) blends toward a lake tint and adds a mirror-direction sun glint. ortho_mode is threaded through dispatch_frame; the drape is a no-op while both extents are 0. GPU suites gain upload->inactive lifecycle and vram-accounting coverage for the new texture pair. Co-Authored-By: Claude Fable 5 --- crates/render_gpu/src/camera.rs | 31 ++- crates/render_gpu/src/lib.rs | 85 ++++++++ crates/render_gpu/src/scene/bind_group.rs | 16 ++ crates/render_gpu/src/scene/mod.rs | 206 +++++++++++++++++++ crates/render_gpu/src/scene/tiers.rs | 226 +++++++++++++++++++++ crates/render_gpu/src/shader_texture.wgsl | 130 +++++++++++- crates/render_gpu/tests/scene_lifecycle.rs | 44 ++++ crates/render_gpu/tests/vram_accounting.rs | 35 ++++ src/viewer/mod.rs | 6 + 9 files changed, 767 insertions(+), 12 deletions(-) diff --git a/crates/render_gpu/src/camera.rs b/crates/render_gpu/src/camera.rs index 340f705..88ac582 100644 --- a/crates/render_gpu/src/camera.rs +++ b/crates/render_gpu/src/camera.rs @@ -53,6 +53,26 @@ pub struct CameraUniforms { pub smooth_radius_m: f32, pub align_mode: u32, // 0=off, 1=tier viz (green/blue/red) pub _pad7: f32, + // fine ortho window (extent_x == 0.0 means inactive) — albedo texture draped + // over the fine tier's neighbourhood; own rect because the ortho's CRS + // (e.g. MGI GK Central) differs from both the base and fine-tier CRSes. + pub ortho_fine_origin_x: f32, + pub ortho_fine_origin_y: f32, + pub ortho_fine_extent_x: f32, + pub ortho_fine_extent_y: f32, + pub ortho_fine_cos_rot: f32, + pub ortho_fine_sin_rot: f32, + // close ortho window (extent_x == 0.0 means inactive) + pub ortho_close_origin_x: f32, + pub ortho_close_origin_y: f32, + pub ortho_close_extent_x: f32, + pub ortho_close_extent_y: f32, + pub ortho_close_cos_rot: f32, + pub ortho_close_sin_rot: f32, + pub ortho_mode: u32, // 0=off, 1=drape ortho albedo (T key) + pub _pad8: f32, + pub _pad9: f32, + pub _pad10: f32, } /// World-space camera basis `(forward, right, up)` for a +Z-up look-at camera. @@ -83,21 +103,26 @@ mod tests { /// Layout guard: the struct is mirrored byte-for-byte in WGSL (std140). Any /// field added/removed/reordered without updating the shader must trip this. - /// Total: 14 vec4-sized rows × 16 bytes = 224. + /// Total: 18 vec4-sized rows × 16 bytes = 288. /// /// The size check alone would pass a same-size field swap, so we also pin the /// offsets of the std140-sensitive `vec3` members — each must sit on a 16-byte /// boundary, which is exactly what the trailing `_padN` fields exist to enforce. /// A missing pad shifts one of these and fails here instead of silently making - /// the shader read the wrong bytes. + /// the shader read the wrong bytes. The two ortho-window blocks are pinned the + /// same way: 12 scalars starting right after `_pad7`, with `ortho_mode` + 3 pads + /// closing the final 16-byte row. #[test] fn camera_uniforms_layout_is_std140() { - assert_eq!(std::mem::size_of::(), 224); + assert_eq!(std::mem::size_of::(), 288); assert_eq!(std::mem::offset_of!(CameraUniforms, origin), 0); assert_eq!(std::mem::offset_of!(CameraUniforms, forward), 16); assert_eq!(std::mem::offset_of!(CameraUniforms, right), 32); assert_eq!(std::mem::offset_of!(CameraUniforms, up), 48); assert_eq!(std::mem::offset_of!(CameraUniforms, sun_dir), 64); + assert_eq!(std::mem::offset_of!(CameraUniforms, ortho_fine_origin_x), 224); + assert_eq!(std::mem::offset_of!(CameraUniforms, ortho_close_origin_x), 248); + assert_eq!(std::mem::offset_of!(CameraUniforms, ortho_mode), 272); } #[test] diff --git a/crates/render_gpu/src/lib.rs b/crates/render_gpu/src/lib.rs index f830639..aa210d6 100644 --- a/crates/render_gpu/src/lib.rs +++ b/crates/render_gpu/src/lib.rs @@ -76,6 +76,51 @@ pub fn gen_hm_mip_bytes( mips } +/// Generate mip level byte data for an RGBA8 ortho texture (levels 1..N-1, +/// N = `hm_mip_count(cols, rows)`). RGB channels are box-averaged; the alpha +/// channel — a discrete material-code ladder, not coverage — takes the top-left +/// sample instead, because averaging codes (0/64/128/192/255) would fabricate +/// classes that don't exist at that location. +/// Call on a background thread before `GpuScene::upload_ortho_*`. +pub fn gen_rgba_mip_bytes(base_rgba: &[u8], cols: usize, rows: usize) -> Vec<(u32, u32, Vec)> { + let extra_mips = (hm_mip_count(cols as u32, rows as u32) as usize).saturating_sub(1); + let mut mips: Vec<(u32, u32, Vec)> = Vec::with_capacity(extra_mips); + let mut prev_w = cols; + let mut prev_h = rows; + for i in 0..extra_mips { + let w = (prev_w / 2).max(1); + let h = (prev_h / 2).max(1); + let mip_bytes = { + let prev: &[u8] = if i == 0 { base_rgba } else { &mips[i - 1].2 }; + let mut out = Vec::with_capacity(w * h * 4); + for row in 0..h { + for col in 0..w { + let r0 = (row * 2).min(prev_h - 1); + let r1 = (row * 2 + 1).min(prev_h - 1); + let c0 = (col * 2).min(prev_w - 1); + let c1 = (col * 2 + 1).min(prev_w - 1); + let px = |r: usize, c: usize| -> &[u8] { + let off = (r * prev_w + c) * 4; + &prev[off..off + 4] + }; + let (p00, p01, p10, p11) = (px(r0, c0), px(r0, c1), px(r1, c0), px(r1, c1)); + for ch in 0..3 { + let sum = + p00[ch] as u16 + p01[ch] as u16 + p10[ch] as u16 + p11[ch] as u16; + out.push((sum / 4) as u8); + } + out.push(p00[3]); // material code: nearest, never averaged + } + } + out + }; + prev_w = w; + prev_h = h; + mips.push((w as u32, h as u32, mip_bytes)); + } + mips +} + /// Pack normal vectors into Rg16Snorm bytes (4 bytes/pixel) for GPU texture upload. /// Call on a background thread before `GpuScene::upload_hm5m` / `upload_hm1m`. pub fn pack_normals_rg16_bytes(nx: &[f32], ny: &[f32]) -> Vec { @@ -203,6 +248,46 @@ mod tests { } } + // gen_rgba_mip_bytes + + #[test] + fn rgba_mip_shapes_match_hm_pyramid() { + // 4×4 RGBA base → 2 generated levels (2×2, 1×1), 4 bytes/texel. + let base = vec![100u8; 4 * 4 * 4]; + let mips = gen_rgba_mip_bytes(&base, 4, 4); + assert_eq!(mips.len(), hm_mip_count(4, 4) as usize - 1); + assert_eq!((mips[0].0, mips[0].1), (2, 2)); + assert_eq!(mips[0].2.len(), 2 * 2 * 4); + assert_eq!((mips[1].0, mips[1].1), (1, 1)); + assert_eq!(mips[1].2.len(), 4); + } + + #[test] + fn rgba_mip_averages_rgb_but_takes_nearest_alpha() { + // 2×2 base where the four texels have RGB 0/40/80/120 and alphas that + // would average to a nonexistent material code (0+255+64+128)/4 = 111. + #[rustfmt::skip] + let base = vec![ + 0, 0, 0, 0, 40, 40, 40, 255, + 80, 80, 80, 64, 120, 120, 120, 128, + ]; + let mips = gen_rgba_mip_bytes(&base, 2, 2); + assert_eq!((mips[0].0, mips[0].1), (1, 1)); + let m = &mips[0].2; + assert_eq!(m[0], (0 + 40 + 80 + 120) / 4, "RGB box-averaged"); + assert_eq!(m[3], 0, "alpha = top-left sample, never an average"); + } + + #[test] + fn rgba_mip_odd_dimensions_clamp_like_hm_path() { + let base = vec![7u8; 5 * 3 * 4]; + let mips = gen_rgba_mip_bytes(&base, 5, 3); + assert_eq!((mips[0].0, mips[0].1), (2, 1)); + for (w, h, bytes) in &mips { + assert_eq!(bytes.len(), (*w * *h * 4) as usize); + } + } + // normal packing fn decode_rg16(bytes: &[u8], idx: usize) -> (f32, f32) { diff --git a/crates/render_gpu/src/scene/bind_group.rs b/crates/render_gpu/src/scene/bind_group.rs index eefa4bf..aa7b5ab 100644 --- a/crates/render_gpu/src/scene/bind_group.rs +++ b/crates/render_gpu/src/scene/bind_group.rs @@ -81,6 +81,22 @@ impl GpuScene { binding: 19, resource: self._hm1m_shadow_buf.as_entire_binding(), }, + wgpu::BindGroupEntry { + binding: 20, + resource: wgpu::BindingResource::TextureView(&self._ortho_fine_view), + }, + wgpu::BindGroupEntry { + binding: 21, + resource: wgpu::BindingResource::Sampler(&self._ortho_fine_sampler), + }, + wgpu::BindGroupEntry { + binding: 22, + resource: wgpu::BindingResource::TextureView(&self._ortho_close_view), + }, + wgpu::BindGroupEntry { + binding: 23, + resource: wgpu::BindingResource::Sampler(&self._ortho_close_sampler), + }, ], }); } diff --git a/crates/render_gpu/src/scene/mod.rs b/crates/render_gpu/src/scene/mod.rs index 9653309..a0c1b04 100644 --- a/crates/render_gpu/src/scene/mod.rs +++ b/crates/render_gpu/src/scene/mod.rs @@ -60,6 +60,31 @@ pub struct GpuScene { pub(super) hm1m_sin_rot: f32, pub(super) hm1m_buf_elems: u64, + // ortho albedo windows (placeholder until upload_ortho_*; extent_x==0 = inactive) + pub(super) _ortho_fine_tex: wgpu::Texture, + pub(super) _ortho_fine_view: wgpu::TextureView, + pub(super) _ortho_fine_sampler: wgpu::Sampler, + pub(super) ortho_fine_origin_x: f32, + pub(super) ortho_fine_origin_y: f32, + pub(super) ortho_fine_extent_x: f32, + pub(super) ortho_fine_extent_y: f32, + pub(super) ortho_fine_cos_rot: f32, + pub(super) ortho_fine_sin_rot: f32, + pub(super) ortho_fine_cols: u32, + pub(super) ortho_fine_rows: u32, + + pub(super) _ortho_close_tex: wgpu::Texture, + pub(super) _ortho_close_view: wgpu::TextureView, + pub(super) _ortho_close_sampler: wgpu::Sampler, + pub(super) ortho_close_origin_x: f32, + pub(super) ortho_close_origin_y: f32, + pub(super) ortho_close_extent_x: f32, + pub(super) ortho_close_extent_y: f32, + pub(super) ortho_close_cos_rot: f32, + pub(super) ortho_close_sin_rot: f32, + pub(super) ortho_close_cols: u32, + pub(super) ortho_close_rows: u32, + // Mutable per-frame / per-sun-update pub(super) shadow_buf: wgpu::Buffer, pub(super) cam_buf: wgpu::Buffer, @@ -309,6 +334,53 @@ pub(super) fn make_tier_size_placeholders( (texture, view, normal_tex, normal_view, shadow_buf) } +/// 1×1 Rgba8Unorm placeholder texture for an ortho albedo slot. Used both at +/// scene creation and by the drop-first reload cycle (`upload_ortho_*` / +/// `set_ortho_*_inactive`) — the existing `make_tier_size_placeholders` is +/// hardcoded to the height-tier resource trio (R16Float + Rgba8Snorm + shadow +/// buffer), so the ortho slot gets its own minimal variant. +pub(super) fn make_ortho_placeholder( + device: &wgpu::Device, + queue: &wgpu::Queue, + label: &str, +) -> (wgpu::Texture, wgpu::TextureView) { + let tex_label = format!("{}_tex", label); + let texture = vram::create_texture_tracked( + device, + &wgpu::TextureDescriptor { + label: Some(&tex_label), + 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::Rgba8Unorm, + usage: wgpu::TextureUsages::TEXTURE_BINDING | wgpu::TextureUsages::COPY_DST, + view_formats: &[], + }, + &tex_label, + ); + queue.write_texture( + texture.as_image_copy(), + &[0u8, 0, 0, 0], + wgpu::TexelCopyBufferLayout { + offset: 0, + bytes_per_row: Some(4), + rows_per_image: None, + }, + wgpu::Extent3d { + width: 1, + height: 1, + depth_or_array_layers: 1, + }, + ); + let view = texture.create_view(&wgpu::TextureViewDescriptor::default()); + (texture, view) +} + /// Generate mip levels 1..7 for a heightmap texture using a max filter. pub(super) fn write_hm_mips( queue: &wgpu::Queue, @@ -338,6 +410,35 @@ pub(super) fn write_hm_mips( } } +/// Write pre-generated RGBA8 mip levels 1..N (4 bytes/texel). +pub(super) fn write_rgba_mips( + queue: &wgpu::Queue, + texture: &wgpu::Texture, + mips: &[(u32, u32, Vec)], +) { + for (mip_idx, (w, h, data)) in mips.iter().enumerate() { + queue.write_texture( + wgpu::TexelCopyTextureInfo { + texture, + mip_level: (mip_idx + 1) as u32, + origin: wgpu::Origin3d::ZERO, + aspect: wgpu::TextureAspect::All, + }, + data, + wgpu::TexelCopyBufferLayout { + offset: 0, + bytes_per_row: Some(w * 4), + rows_per_image: None, + }, + wgpu::Extent3d { + width: *w, + height: *h, + depth_or_array_layers: 1, + }, + ); + } +} + impl GpuScene { pub fn new( gpu_ctx: GpuContext, @@ -458,6 +559,22 @@ impl GpuScene { hm1m_shadow_buf, ) = create_tier_placeholder(&gpu_ctx.device, &gpu_ctx.queue, "hm1m"); + // ortho albedo placeholders (1×1 Rgba8Unorm) — inactive until upload_ortho_* + let (ortho_fine_tex, ortho_fine_view) = + make_ortho_placeholder(&gpu_ctx.device, &gpu_ctx.queue, "ortho_fine"); + let (ortho_close_tex, ortho_close_view) = + make_ortho_placeholder(&gpu_ctx.device, &gpu_ctx.queue, "ortho_close"); + let mk_ortho_sampler = || { + gpu_ctx.device.create_sampler(&wgpu::SamplerDescriptor { + mag_filter: wgpu::FilterMode::Linear, + min_filter: wgpu::FilterMode::Linear, + mipmap_filter: wgpu::MipmapFilterMode::Linear, + ..Default::default() + }) + }; + let ortho_fine_sampler = mk_ortho_sampler(); + let ortho_close_sampler = mk_ortho_sampler(); + // normals packed buffer: bits 31–16 = nx_i16, bits 15–0 = ny_i16; nz reconstructed in shader // COPY_DST so update_heightmap can write_buffer let normals_packed: Vec = normal_map @@ -680,6 +797,39 @@ impl GpuScene { }, count: None, }, + // ortho albedo windows (fine 20/21, close 22/23) + wgpu::BindGroupLayoutEntry { + binding: 20, + visibility: wgpu::ShaderStages::COMPUTE, + ty: wgpu::BindingType::Texture { + sample_type: wgpu::TextureSampleType::Float { filterable: true }, + view_dimension: wgpu::TextureViewDimension::D2, + multisampled: false, + }, + count: None, + }, + wgpu::BindGroupLayoutEntry { + binding: 21, + visibility: wgpu::ShaderStages::COMPUTE, + ty: wgpu::BindingType::Sampler(wgpu::SamplerBindingType::Filtering), + count: None, + }, + wgpu::BindGroupLayoutEntry { + binding: 22, + visibility: wgpu::ShaderStages::COMPUTE, + ty: wgpu::BindingType::Texture { + sample_type: wgpu::TextureSampleType::Float { filterable: true }, + view_dimension: wgpu::TextureViewDimension::D2, + multisampled: false, + }, + count: None, + }, + wgpu::BindGroupLayoutEntry { + binding: 23, + visibility: wgpu::ShaderStages::COMPUTE, + ty: wgpu::BindingType::Sampler(wgpu::SamplerBindingType::Filtering), + count: None, + }, ], }); @@ -765,6 +915,23 @@ impl GpuScene { binding: 19, resource: hm1m_shadow_buf.as_entire_binding(), }, + // ortho albedo windows (placeholders) + wgpu::BindGroupEntry { + binding: 20, + resource: wgpu::BindingResource::TextureView(&ortho_fine_view), + }, + wgpu::BindGroupEntry { + binding: 21, + resource: wgpu::BindingResource::Sampler(&ortho_fine_sampler), + }, + wgpu::BindGroupEntry { + binding: 22, + resource: wgpu::BindingResource::TextureView(&ortho_close_view), + }, + wgpu::BindGroupEntry { + binding: 23, + resource: wgpu::BindingResource::Sampler(&ortho_close_sampler), + }, ], }); let render_shader = gpu_ctx @@ -834,6 +1001,28 @@ impl GpuScene { hm1m_cos_rot: 1.0, hm1m_sin_rot: 0.0, hm1m_buf_elems: 1, + _ortho_fine_tex: ortho_fine_tex, + _ortho_fine_view: ortho_fine_view, + _ortho_fine_sampler: ortho_fine_sampler, + ortho_fine_origin_x: 0.0, + ortho_fine_origin_y: 0.0, + ortho_fine_extent_x: 0.0, + ortho_fine_extent_y: 0.0, + ortho_fine_cos_rot: 1.0, + ortho_fine_sin_rot: 0.0, + ortho_fine_cols: 0, + ortho_fine_rows: 0, + _ortho_close_tex: ortho_close_tex, + _ortho_close_view: ortho_close_view, + _ortho_close_sampler: ortho_close_sampler, + ortho_close_origin_x: 0.0, + ortho_close_origin_y: 0.0, + ortho_close_extent_x: 0.0, + ortho_close_extent_y: 0.0, + ortho_close_cos_rot: 1.0, + ortho_close_sin_rot: 0.0, + ortho_close_cols: 0, + ortho_close_rows: 0, shadow_buf, cam_buf, output_buf, @@ -875,6 +1064,7 @@ impl GpuScene { lod_mode: u32, smooth_radius_m: f32, align_mode: u32, + ortho_mode: u32, ) { let (forward, right, up) = crate::camera::camera_basis(origin, look_at); let (half_w, half_h) = crate::camera::projection_half_extents(fov_deg, aspect); @@ -926,6 +1116,22 @@ impl GpuScene { smooth_radius_m, align_mode, _pad7: 0.0, + ortho_fine_origin_x: self.ortho_fine_origin_x, + ortho_fine_origin_y: self.ortho_fine_origin_y, + ortho_fine_extent_x: self.ortho_fine_extent_x, + ortho_fine_extent_y: self.ortho_fine_extent_y, + ortho_fine_cos_rot: self.ortho_fine_cos_rot, + ortho_fine_sin_rot: self.ortho_fine_sin_rot, + ortho_close_origin_x: self.ortho_close_origin_x, + ortho_close_origin_y: self.ortho_close_origin_y, + ortho_close_extent_x: self.ortho_close_extent_x, + ortho_close_extent_y: self.ortho_close_extent_y, + ortho_close_cos_rot: self.ortho_close_cos_rot, + ortho_close_sin_rot: self.ortho_close_sin_rot, + ortho_mode, + _pad8: 0.0, + _pad9: 0.0, + _pad10: 0.0, }; self.gpu_ctx diff --git a/crates/render_gpu/src/scene/tiers.rs b/crates/render_gpu/src/scene/tiers.rs index b6047c9..f1a5620 100644 --- a/crates/render_gpu/src/scene/tiers.rs +++ b/crates/render_gpu/src/scene/tiers.rs @@ -378,6 +378,232 @@ impl GpuScene { self.hm1m_sin_rot = rot_rad.sin(); } + /// Upload the fine ortho albedo window (RGBA8: RGB = orthophoto, A = material + /// code). Same drop-first reallocation strategy as `upload_hm5m`; steady-state + /// reloads with unchanged dimensions overwrite in place via `write_texture`. + // Args are the worker-packed byte slices + window geometry; kept flat by design. + #[allow(clippy::too_many_arguments)] + pub fn upload_ortho_fine( + &mut self, + origin_x: f32, + origin_y: f32, + rot_rad: f32, + extent_x: f32, + extent_y: f32, + cols: u32, + rows: u32, + rgba: &[u8], + mips: &[(u32, u32, Vec)], + ) { + let size_changed = cols != self.ortho_fine_cols || rows != self.ortho_fine_rows; + + if size_changed { + // Drop-first cycle (see upload_hm5m): placeholder swap → rebuild → + // poll(Wait) so the old texture is actually freed before the new + // allocation. extent_x = 0.0 keeps the shader off the placeholder. + self.ortho_fine_extent_x = 0.0; + + vram::track_texture_drop(&self._ortho_fine_tex, "ortho_fine_tex"); + let (ph_tex, ph_view) = super::make_ortho_placeholder( + &self.gpu_ctx.device, + &self.gpu_ctx.queue, + "ortho_fine", + ); + self._ortho_fine_tex = ph_tex; + self._ortho_fine_view = ph_view; + self.ortho_fine_cols = 0; + self.ortho_fine_rows = 0; + + self.rebuild_bind_group(); + let _ = self.gpu_ctx.device.poll(wgpu::PollType::Wait { + submission_index: None, + timeout: None, + }); + + vram::track_texture_drop(&self._ortho_fine_tex, "ortho_fine_tex"); + let texture = vram::create_texture_tracked( + &self.gpu_ctx.device, + &wgpu::TextureDescriptor { + label: Some("ortho_fine_tex"), + size: wgpu::Extent3d { + width: cols, + height: rows, + depth_or_array_layers: 1, + }, + mip_level_count: crate::hm_mip_count(cols, rows), + sample_count: 1, + dimension: wgpu::TextureDimension::D2, + format: wgpu::TextureFormat::Rgba8Unorm, + usage: wgpu::TextureUsages::TEXTURE_BINDING | wgpu::TextureUsages::COPY_DST, + view_formats: &[], + }, + "ortho_fine_tex", + ); + self._ortho_fine_view = texture.create_view(&wgpu::TextureViewDescriptor::default()); + self._ortho_fine_tex = texture; + + self.rebuild_bind_group(); + } + + self.gpu_ctx.queue.write_texture( + self._ortho_fine_tex.as_image_copy(), + rgba, + wgpu::TexelCopyBufferLayout { + offset: 0, + bytes_per_row: Some(cols * 4), + rows_per_image: None, + }, + wgpu::Extent3d { + width: cols, + height: rows, + depth_or_array_layers: 1, + }, + ); + super::write_rgba_mips(&self.gpu_ctx.queue, &self._ortho_fine_tex, mips); + + self.ortho_fine_origin_x = origin_x; + self.ortho_fine_origin_y = origin_y; + self.ortho_fine_extent_x = extent_x; + self.ortho_fine_extent_y = extent_y; + self.ortho_fine_cols = cols; + self.ortho_fine_rows = rows; + self.ortho_fine_cos_rot = rot_rad.cos(); + self.ortho_fine_sin_rot = rot_rad.sin(); + } + + /// Disable the fine ortho window and eagerly free its texture. + /// See `set_hm5m_inactive` for the mechanism. + pub fn set_ortho_fine_inactive(&mut self) { + if self.ortho_fine_cols == 0 && self.ortho_fine_rows == 0 && self.ortho_fine_extent_x == 0.0 + { + return; + } + self.ortho_fine_extent_x = 0.0; + vram::track_texture_drop(&self._ortho_fine_tex, "ortho_fine_tex"); + let (ph_tex, ph_view) = + super::make_ortho_placeholder(&self.gpu_ctx.device, &self.gpu_ctx.queue, "ortho_fine"); + self._ortho_fine_tex = ph_tex; + self._ortho_fine_view = ph_view; + self.ortho_fine_cols = 0; + self.ortho_fine_rows = 0; + self.rebuild_bind_group(); + let _ = self.gpu_ctx.device.poll(wgpu::PollType::Wait { + submission_index: None, + timeout: None, + }); + } + + /// Upload the close ortho albedo window. Same mechanism as `upload_ortho_fine`. + // Args are the worker-packed byte slices + window geometry; kept flat by design. + #[allow(clippy::too_many_arguments)] + pub fn upload_ortho_close( + &mut self, + origin_x: f32, + origin_y: f32, + rot_rad: f32, + extent_x: f32, + extent_y: f32, + cols: u32, + rows: u32, + rgba: &[u8], + mips: &[(u32, u32, Vec)], + ) { + let size_changed = cols != self.ortho_close_cols || rows != self.ortho_close_rows; + + if size_changed { + self.ortho_close_extent_x = 0.0; + + vram::track_texture_drop(&self._ortho_close_tex, "ortho_close_tex"); + let (ph_tex, ph_view) = super::make_ortho_placeholder( + &self.gpu_ctx.device, + &self.gpu_ctx.queue, + "ortho_close", + ); + self._ortho_close_tex = ph_tex; + self._ortho_close_view = ph_view; + self.ortho_close_cols = 0; + self.ortho_close_rows = 0; + + self.rebuild_bind_group(); + let _ = self.gpu_ctx.device.poll(wgpu::PollType::Wait { + submission_index: None, + timeout: None, + }); + + vram::track_texture_drop(&self._ortho_close_tex, "ortho_close_tex"); + let texture = vram::create_texture_tracked( + &self.gpu_ctx.device, + &wgpu::TextureDescriptor { + label: Some("ortho_close_tex"), + size: wgpu::Extent3d { + width: cols, + height: rows, + depth_or_array_layers: 1, + }, + mip_level_count: crate::hm_mip_count(cols, rows), + sample_count: 1, + dimension: wgpu::TextureDimension::D2, + format: wgpu::TextureFormat::Rgba8Unorm, + usage: wgpu::TextureUsages::TEXTURE_BINDING | wgpu::TextureUsages::COPY_DST, + view_formats: &[], + }, + "ortho_close_tex", + ); + self._ortho_close_view = texture.create_view(&wgpu::TextureViewDescriptor::default()); + self._ortho_close_tex = texture; + + self.rebuild_bind_group(); + } + + self.gpu_ctx.queue.write_texture( + self._ortho_close_tex.as_image_copy(), + rgba, + wgpu::TexelCopyBufferLayout { + offset: 0, + bytes_per_row: Some(cols * 4), + rows_per_image: None, + }, + wgpu::Extent3d { + width: cols, + height: rows, + depth_or_array_layers: 1, + }, + ); + super::write_rgba_mips(&self.gpu_ctx.queue, &self._ortho_close_tex, mips); + + self.ortho_close_origin_x = origin_x; + self.ortho_close_origin_y = origin_y; + self.ortho_close_extent_x = extent_x; + self.ortho_close_extent_y = extent_y; + self.ortho_close_cols = cols; + self.ortho_close_rows = rows; + self.ortho_close_cos_rot = rot_rad.cos(); + self.ortho_close_sin_rot = rot_rad.sin(); + } + + /// Disable the close ortho window and eagerly free its texture. + pub fn set_ortho_close_inactive(&mut self) { + if self.ortho_close_cols == 0 + && self.ortho_close_rows == 0 + && self.ortho_close_extent_x == 0.0 + { + return; + } + self.ortho_close_extent_x = 0.0; + vram::track_texture_drop(&self._ortho_close_tex, "ortho_close_tex"); + let (ph_tex, ph_view) = + super::make_ortho_placeholder(&self.gpu_ctx.device, &self.gpu_ctx.queue, "ortho_close"); + self._ortho_close_tex = ph_tex; + self._ortho_close_view = ph_view; + self.ortho_close_cols = 0; + self.ortho_close_rows = 0; + self.rebuild_bind_group(); + let _ = self.gpu_ctx.device.poll(wgpu::PollType::Wait { + submission_index: None, + timeout: None, + }); + } + /// Disable the 1 m fine tier and eagerly free its GPU resources. /// See `set_hm5m_inactive` for the mechanism. pub fn set_hm1m_inactive(&mut self) { diff --git a/crates/render_gpu/src/shader_texture.wgsl b/crates/render_gpu/src/shader_texture.wgsl index 5e862cc..45a2f57 100644 --- a/crates/render_gpu/src/shader_texture.wgsl +++ b/crates/render_gpu/src/shader_texture.wgsl @@ -47,6 +47,24 @@ struct CameraUniforms { smooth_radius_m: f32, align_mode: u32, _pad7: f32, + // fine ortho albedo window (extent_x == 0.0 means inactive) + ortho_fine_origin_x: f32, + ortho_fine_origin_y: f32, + ortho_fine_extent_x: f32, + ortho_fine_extent_y: f32, + ortho_fine_cos_rot: f32, + ortho_fine_sin_rot: f32, + // close ortho albedo window (extent_x == 0.0 means inactive) + ortho_close_origin_x: f32, + ortho_close_origin_y: f32, + ortho_close_extent_x: f32, + ortho_close_extent_y: f32, + ortho_close_cos_rot: f32, + ortho_close_sin_rot: f32, + ortho_mode: u32, + _pad8: f32, + _pad9: f32, + _pad10: f32, } // camera uniforms struct @@ -75,6 +93,12 @@ struct CameraUniforms { @group(0) @binding(17) var hm1m_normal_tex: texture_2d; @group(0) @binding(18) var hm1m_normal_samp: sampler; @group(0) @binding(19) var hm1m_shadow: array; +// ortho albedo windows: rgb = orthophoto, a = material code +// (water 1.0, high veg ~0.75, med veg ~0.5, building ~0.25, none 0) +@group(0) @binding(20) var ortho_fine_tex: texture_2d; +@group(0) @binding(21) var ortho_fine_samp: sampler; +@group(0) @binding(22) var ortho_close_tex: texture_2d; +@group(0) @binding(23) var ortho_close_samp: sampler; // Width of the blend zone at the edge of any tier boundary, in metres. // smoothstep from 0 (edge, fully base tier) to BLEND_MARGIN (inner, fully close tier). @@ -211,6 +235,50 @@ fn apply_convergence_1m(lx: f32, ly: f32) -> vec2 { ); } +// Ortho-window blend weight at world position (x, y): 0 outside the rectangle, +// smoothstep over BLEND_MARGIN of edge distance inside (same scheme as the +// height tiers). The returned uv is rotation-corrected and normalised. +struct OrthoSample { + weight: f32, + uv: vec2, +} + +fn ortho_fine_sample_at(x: f32, y: f32) -> OrthoSample { + var s = OrthoSample(0.0, vec2(0.0)); + if cam.ortho_fine_extent_x <= 0.0 { return s; } + let lx = x - cam.ortho_fine_origin_x; + let ly = y - cam.ortho_fine_origin_y; + if lx < 0.0 || lx >= cam.ortho_fine_extent_x || ly < 0.0 || ly >= cam.ortho_fine_extent_y { + return s; + } + let edge = min(min(lx, cam.ortho_fine_extent_x - lx), min(ly, cam.ortho_fine_extent_y - ly)); + let a = vec2( + lx * cam.ortho_fine_cos_rot + ly * cam.ortho_fine_sin_rot, + -lx * cam.ortho_fine_sin_rot + ly * cam.ortho_fine_cos_rot + ); + s.weight = smoothstep(0.0, BLEND_MARGIN, edge); + s.uv = vec2(a.x / cam.ortho_fine_extent_x, a.y / cam.ortho_fine_extent_y); + return s; +} + +fn ortho_close_sample_at(x: f32, y: f32) -> OrthoSample { + var s = OrthoSample(0.0, vec2(0.0)); + if cam.ortho_close_extent_x <= 0.0 { return s; } + let lx = x - cam.ortho_close_origin_x; + let ly = y - cam.ortho_close_origin_y; + if lx < 0.0 || lx >= cam.ortho_close_extent_x || ly < 0.0 || ly >= cam.ortho_close_extent_y { + return s; + } + let edge = min(min(lx, cam.ortho_close_extent_x - lx), min(ly, cam.ortho_close_extent_y - ly)); + let a = vec2( + lx * cam.ortho_close_cos_rot + ly * cam.ortho_close_sin_rot, + -lx * cam.ortho_close_sin_rot + ly * cam.ortho_close_cos_rot + ); + s.weight = smoothstep(0.0, BLEND_MARGIN, edge); + s.uv = vec2(a.x / cam.ortho_close_extent_x, a.y / cam.ortho_close_extent_y); + return s; +} + // Catmull-Rom 1D weight vector for fractional offset t ∈ [0, 1]. fn cr_w(t: f32) -> vec4 { let t2 = t * t; @@ -798,13 +866,57 @@ fn main(@builtin(global_invocation_id) gid: vec3) { // is narrow enough that the easing is what gives the coastline its crisp edge. let coastal = mix(sea, turquoise, clamp(pos.z / 15.5, 0.0, 1.0)); let base = mix(coastal, land, smoothstep(15.5, 17.0, pos.z)); - let base_r = base.x; - let base_g = base.y; - let base_b = base.z; - let r = u32(clamp(base_r * brightness, 0.0, 255.0)); - let g = u32(clamp(base_g * brightness, 0.0, 255.0)); - let b = u32(clamp(base_b * brightness, 0.0, 255.0)); + // ── albedo: procedural ramp, overridden by streamed ortho windows ── + // Close window first, fine window on top (same nesting as the height + // tiers). The mip level reuses the distance heuristic from the march + // loop so distant ortho texels don't shimmer; each texture clamps to + // its own level count because fine/close windows differ in size. + var albedo = base; + var ortho_w = 0.0; // strongest ortho influence at this hit + var material = 0.0; // land-cover material code from the winning sample's alpha + if cam.ortho_mode == 1u { + let lod_mip_div_o = select(select(select(8000.0, 15000.0, cam.lod_mode < 3u), 30000.0, cam.lod_mode < 2u), 1000000000.0, cam.lod_mode == 0u); + let mip_o = log2(1.0 + t / lod_mip_div_o); + + let oc = ortho_close_sample_at(pos.x, pos.y); + if oc.weight > 0.0 { + let mip_c = min(mip_o, f32(textureNumLevels(ortho_close_tex) - 1u)); + let s = textureSampleLevel(ortho_close_tex, ortho_close_samp, oc.uv, mip_c); + albedo = mix(albedo, s.rgb * 255.0, oc.weight); + material = mix(material, s.a, oc.weight); + ortho_w = max(ortho_w, oc.weight); + } + let ofs = ortho_fine_sample_at(pos.x, pos.y); + if ofs.weight > 0.0 { + let mip_f = min(mip_o, f32(textureNumLevels(ortho_fine_tex) - 1u)); + let s = textureSampleLevel(ortho_fine_tex, ortho_fine_samp, ofs.uv, mip_f); + albedo = mix(albedo, s.rgb * 255.0, ofs.weight); + material = mix(material, s.a, ofs.weight); + ortho_w = max(ortho_w, ofs.weight); + } + } + + // The orthophoto already contains the capture day's sun shading; running + // our full diffuse on top double-shades (north faces go pitch black). + // Where ortho wins, lift the brightness floor but keep enough relief + // modulation that the geometry still reads. + let lit = mix(brightness, clamp(brightness, 0.55, 1.0), ortho_w * 0.7); + var rgb = albedo * lit; + + // Water material (land-cover class baked in alpha): flatten the noisy + // ortho water pixels toward a lake tint and add a sun glint along the + // mirror direction. smoothstep over the code keeps shorelines soft. + let water = smoothstep(0.85, 0.95, material); + if water > 0.0 { + rgb = mix(rgb, vec3(45.0, 80.0, 110.0) * lit, water * 0.55); + let glint = pow(max(dot(reflect(dir, normal), normalized_sun_dir), 0.0), 64.0); + rgb += vec3(255.0, 240.0, 210.0) * glint * water * 0.5 * shadow_factor; + } + + let r = clamp(rgb.x, 0.0, 255.0); + let g = clamp(rgb.y, 0.0, 255.0); + let b = clamp(rgb.z, 0.0, 255.0); // add fog/haze in the distance — fog colour tracks the sky gradient so // looking up at a distant peak hazes into the correct overhead tone. @@ -818,9 +930,9 @@ fn main(@builtin(global_invocation_id) gid: vec3) { let fog_far = 60000.0 * alt_visibility_mul; // 60 km at sea level let fog_t = select(0.0, smoothstep(fog_near, fog_far, t), cam.fog_enabled == 1u); - let fr = mix(f32(r), sky.x, fog_t); - let fg = mix(f32(g), sky.y, fog_t); - let fb = mix(f32(b), sky.z, fog_t); + let fr = mix(r, sky.x, fog_t); + let fg = mix(g, sky.y, fog_t); + let fb = mix(b, sky.z, fog_t); output[gid.y * cam.img_width + gid.x] = u32(fb) | (u32(fg) << 8u) | (u32(fr) << 16u) | (255u << 24u); // bgra format } else { diff --git a/crates/render_gpu/tests/scene_lifecycle.rs b/crates/render_gpu/tests/scene_lifecycle.rs index fd8d6f8..4d2f456 100644 --- a/crates/render_gpu/tests/scene_lifecycle.rs +++ b/crates/render_gpu/tests/scene_lifecycle.rs @@ -34,6 +34,7 @@ fn dispatch_once(scene: &GpuScene, ctx: &render_gpu::GpuContext) { 0, // lod_mode 2000.0, // smooth_radius_m 0, // align_mode + 0, // ortho_mode ); ctx.queue.submit([enc.finish()]); drain(ctx); @@ -96,3 +97,46 @@ fn tier_upload_then_inactive_survives_validation() { scene.set_hm5m_inactive(); dispatch_once(&scene, &ctx); } + +#[test] +fn ortho_upload_then_inactive_survives_validation() { + gpu_or_skip!(ctx); + + let base = pseudo_random(64, 64, 4, 500.0, 30.0, 30.0); + let (bn, bs, bao) = derive_maps(&base); + let mut scene = GpuScene::new(ctx.clone(), &base, &bn, &bs, &bao, 96, 96); + dispatch_once(&scene, &ctx); + + // Upload both ortho windows with full mip chains — exercises the grow path + // on bindings 20–23 plus the per-mip write_texture calls. The dispatch runs + // with ortho_mode = 0 in dispatch_once, but validation covers the bindings + // regardless of the uniform flag. + let (w, h) = (96usize, 64usize); + let rgba: Vec = (0..w * h * 4).map(|i| (i % 251) as u8).collect(); + let mips = render_gpu::gen_rgba_mip_bytes(&rgba, w, h); + scene.upload_ortho_fine(0.0, 0.0, 0.01, 96.0, 64.0, w as u32, h as u32, &rgba, &mips); + dispatch_once(&scene, &ctx); + + let (cw, ch) = (40usize, 56usize); + let rgba_c: Vec = (0..cw * ch * 4).map(|i| (i % 13) as u8).collect(); + let mips_c = render_gpu::gen_rgba_mip_bytes(&rgba_c, cw, ch); + scene.upload_ortho_close( + 0.0, 0.0, -0.02, 256.0, 358.4, cw as u32, ch as u32, &rgba_c, &mips_c, + ); + dispatch_once(&scene, &ctx); + + // Same-size re-upload writes in place (no reallocation path). + scene.upload_ortho_fine(5.0, 5.0, 0.0, 96.0, 64.0, w as u32, h as u32, &rgba, &mips); + dispatch_once(&scene, &ctx); + + // Grow the fine window — drop-first reallocation on binding 20. + let (gw, gh) = (130usize, 90usize); + let rgba_g: Vec = vec![128; gw * gh * 4]; + let mips_g = render_gpu::gen_rgba_mip_bytes(&rgba_g, gw, gh); + scene.upload_ortho_fine(0.0, 0.0, 0.0, 130.0, 90.0, gw as u32, gh as u32, &rgba_g, &mips_g); + dispatch_once(&scene, &ctx); + + scene.set_ortho_fine_inactive(); + scene.set_ortho_close_inactive(); + dispatch_once(&scene, &ctx); +} diff --git a/crates/render_gpu/tests/vram_accounting.rs b/crates/render_gpu/tests/vram_accounting.rs index e3979a1..e4c2063 100644 --- a/crates/render_gpu/tests/vram_accounting.rs +++ b/crates/render_gpu/tests/vram_accounting.rs @@ -106,3 +106,38 @@ fn set_inactive_frees_tracked_textures() { "inactive tier should free its textures (residual {after} bytes)" ); } + +#[test] +#[serial] +fn ortho_inactive_frees_tracked_textures() { + gpu_or_skip!(ctx); + + let base = pseudo_random(64, 64, 6, 500.0, 30.0, 30.0); + let (bn, bs, bao) = derive_maps(&base); + let mut scene = GpuScene::new(ctx.clone(), &base, &bn, &bs, &bao, 64, 64); + + let baseline = tex_bytes(); + + let (w, h) = (128usize, 128usize); + let rgba: Vec = vec![200; w * h * 4]; + let mips = render_gpu::gen_rgba_mip_bytes(&rgba, w, h); + scene.upload_ortho_fine(0.0, 0.0, 0.0, 128.0, 128.0, w as u32, h as u32, &rgba, &mips); + let foot = tex_bytes() - baseline; + assert!(foot > 0, "ortho upload should grow tracked textures"); + + // Same-size re-upload: grow-only logic, no allocation. + scene.upload_ortho_fine(1.0, 1.0, 0.0, 128.0, 128.0, w as u32, h as u32, &rgba, &mips); + assert_eq!( + tex_bytes() - baseline, + foot, + "same-size ortho reload must not allocate" + ); + + // Deactivation drops down to the 1×1 placeholder. + scene.set_ortho_fine_inactive(); + let after = tex_bytes() - baseline; + assert!( + after.abs() < 1024, + "inactive ortho should free its texture (residual {after} bytes)" + ); +} diff --git a/src/viewer/mod.rs b/src/viewer/mod.rs index a96e7f0..62e1655 100644 --- a/src/viewer/mod.rs +++ b/src/viewer/mod.rs @@ -111,6 +111,10 @@ pub(crate) struct Viewer { /// banner silent). Drives the red "OOM prevented" warning in the HUD. fine_disabled_by_oom: bool, align_mode_viz: bool, // V key: show all 3 tiers as separate colored surfaces + /// T key: drape the streamed orthophoto albedo over the terrain. Harmless + /// when no ortho window is loaded — the shader skips draping while both + /// ortho extents are 0. + ortho_enabled: bool, } impl ApplicationHandler for Viewer { @@ -616,6 +620,7 @@ impl ApplicationHandler for Viewer { self.lod_mode, self.smooth_radius_m, self.align_mode_viz as u32, + self.ortho_enabled as u32, ); let output_buf: &wgpu::Buffer = scene.get_output_buffer(); @@ -1239,6 +1244,7 @@ impl Viewer { close_tier_disabled: false, fine_disabled_by_oom: false, align_mode_viz: false, + ortho_enabled: true, } } } From 41812ff38ea925d0103ae003472623da5477640c Mon Sep 17 00:00:00 2001 From: Nikita Zavartsev Date: Fri, 12 Jun 2026 20:44:40 +0200 Subject: [PATCH 5/9] feat: [dsm-rgb-landcover-tirol] stream DSM surface overlay and ortho albedo tiers in the viewer MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit DemoViewConfig gains surface_tile_paths / ortho_tile_paths / landcover_tile_paths (serde-defaulted to the BEV Tirol mosaics, resolved against tiles_dir; missing files degrade gracefully via build_tile_index, so old configs keep working untouched). The fine worker extracts a DSM window with the same centre/radius as its DTM window and composites it on with composite_surface_over — canopy and buildings become fine-tier geometry with normals, DDA self-shadowing and bicubic smoothing for free. StreamingTier becomes generic over its payload (TierCentre trait) so two new ortho workers (fine ~0.8 m/px, close ~6.4 m/px, radii per the new ortho_radii(VramClass) presets) reuse all drift/trigger bookkeeping. They resolve mask-aware overview levels once per file, read RGBA windows via extract_color_window, pre-pack mips, and the frame loop uploads them like any other tier. Base reloads deactivate + invalidate both windows (their origins are base-relative); the OOM ladder gains step 0 = drop ortho before any height tier; T toggles the drape. Includes a gated offscreen smoke test that renders the full DSM + ortho pipeline against the real tiles through the production placement math and asserts the frame contains ortho-textured terrain (/tmp/offscreen_dsm_ortho.png for visual inspection). Co-Authored-By: Claude Fable 5 --- src/launcher/config.rs | 47 ++++ src/launcher/mod.rs | 3 + src/viewer/mod.rs | 140 +++++++++++- src/viewer/tiers.rs | 490 ++++++++++++++++++++++++++++++++++++++++- 4 files changed, 668 insertions(+), 12 deletions(-) diff --git a/src/launcher/config.rs b/src/launcher/config.rs index b1e2100..a0ec468 100644 --- a/src/launcher/config.rs +++ b/src/launcher/config.rs @@ -18,6 +18,18 @@ pub struct DemoViewConfig { pub close_tile_paths: Vec, #[serde(default = "default_base_paths")] pub base_tile_paths: Vec, + /// DSM tiles (trees/buildings as surface height) composited over the fine + /// tier's DTM where they overlap. Missing files degrade gracefully — the + /// tile index skips paths that don't exist on disk. + #[serde(default = "default_surface_paths")] + pub surface_tile_paths: Vec, + /// Orthophoto mosaics streamed as albedo windows (fine + close ortho tiers). + #[serde(default = "default_ortho_paths")] + pub ortho_tile_paths: Vec, + /// Land-cover mosaics; material codes (water/vegetation/buildings) are baked + /// into the ortho windows' alpha channel for shader-side material shading. + #[serde(default = "default_landcover_paths")] + pub landcover_tile_paths: Vec, } fn default_demo_cam_lat() -> f64 { @@ -41,6 +53,20 @@ fn default_close_paths() -> Vec { vec![PathBuf::from("big_size/DGM_R5.tif")] } +fn default_surface_paths() -> Vec { + vec![PathBuf::from( + "big_size/ALS_DSM_CRS3035RES50000mN2650000E4450000.tif", + )] +} + +fn default_ortho_paths() -> Vec { + vec![PathBuf::from("color/2019470_Mosaik_RGB.tif")] +} + +fn default_landcover_paths() -> Vec { + vec![PathBuf::from("color/2022470_Mosaik_LC.tif")] +} + fn default_base_paths() -> Vec { // Camera tile + east + south + south-east of default camera position (N47/E011) [(47u32, 11u32), (47, 12), (46, 11), (46, 12)] @@ -90,6 +116,9 @@ impl Default for DemoViewConfig { fine_tile_paths: default_fine_paths(), close_tile_paths: default_close_paths(), base_tile_paths: default_base_paths(), + surface_tile_paths: default_surface_paths(), + ortho_tile_paths: default_ortho_paths(), + landcover_tile_paths: default_landcover_paths(), } } } @@ -310,5 +339,23 @@ mod tests { assert!(!d.base_tile_paths.is_empty()); assert!(!d.close_tile_paths.is_empty()); assert!(!d.fine_tile_paths.is_empty()); + assert!(!d.surface_tile_paths.is_empty()); + assert!(!d.ortho_tile_paths.is_empty()); + assert!(!d.landcover_tile_paths.is_empty()); + } + + #[test] + fn old_config_without_surface_or_color_fields_gets_defaults() { + // A config.toml written before the DSM/ortho fields existed must still + // parse, with the new fields falling back to their serde defaults. + let s: LauncherSettings = toml::from_str( + "[demo_view]\ncamera_lat = 47.0\nfine_tile_paths = [\"a.tif\"]\n", + ) + .expect("pre-DSM config parses"); + assert_eq!(s.demo_view.camera_lat, 47.0); + assert_eq!(s.demo_view.fine_tile_paths, vec![PathBuf::from("a.tif")]); + assert_eq!(s.demo_view.surface_tile_paths, default_surface_paths()); + assert_eq!(s.demo_view.ortho_tile_paths, default_ortho_paths()); + assert_eq!(s.demo_view.landcover_tile_paths, default_landcover_paths()); } } diff --git a/src/launcher/mod.rs b/src/launcher/mod.rs index 08f7ce7..aa25e21 100644 --- a/src/launcher/mod.rs +++ b/src/launcher/mod.rs @@ -198,6 +198,9 @@ impl LauncherApp { cfg.fine_tile_paths = resolve(cfg.fine_tile_paths); cfg.close_tile_paths = resolve(cfg.close_tile_paths); cfg.base_tile_paths = resolve(cfg.base_tile_paths); + cfg.surface_tile_paths = resolve(cfg.surface_tile_paths); + cfg.ortho_tile_paths = resolve(cfg.ortho_tile_paths); + cfg.landcover_tile_paths = resolve(cfg.landcover_tile_paths); cfg }; // Write resolved absolute paths back into settings so LauncherOutcome::Start diff --git a/src/viewer/mod.rs b/src/viewer/mod.rs index 62e1655..7cb457d 100644 --- a/src/viewer/mod.rs +++ b/src/viewer/mod.rs @@ -37,7 +37,8 @@ use crate::viewer::hud_renderer::HudRenderer; use self::geo::{latlon_to_tile_metres, sun_position}; use self::scene_init::{INIT_SIM_DAY, INIT_SIM_HOUR, compute_ao_cropped}; use self::tiers::{ - AO_DRIFT_THRESHOLD_M, BevBaseState, TierRadii, cross_crs_world_origin_and_extent, tier_radii, + AO_DRIFT_THRESHOLD_M, BevBaseState, TierRadii, cross_crs_world_origin_and_extent, ortho_radii, + tier_radii, }; use crate::consts::{GPU_SAFE_PX, M_PER_DEG}; @@ -446,9 +447,12 @@ impl ApplicationHandler for Viewer { scene.update_shadow(&data.shadow); // The fine-tier origins are offsets relative to the base heightmap origin. // After a base reload the origin shifts, so the old offsets are wrong — - // hide both fine tiers until their workers deliver fresh windows. + // hide both fine tiers (and both ortho windows, whose origins are + // base-relative too) until their workers deliver fresh windows. scene.set_hm5m_inactive(); scene.set_hm1m_inactive(); + scene.set_ortho_fine_inactive(); + scene.set_ortho_close_inactive(); } self.hm = data.hm; // Recalibrate drift threshold to match the actual loaded window. @@ -468,6 +472,12 @@ impl ApplicationHandler for Viewer { if let Some(ref mut fine) = bev_base.fine { fine.invalidate(); } + if let Some(ref mut cf) = bev_base.color_fine { + cf.invalidate(); + } + if let Some(ref mut cc) = bev_base.color_close { + cc.invalidate(); + } // Respawn shadow worker with updated heightmap. let (new_tx, new_worker_rx) = mpsc::sync_channel::<(f32, f32)>(1); let (new_worker_tx, new_rx) = mpsc::channel::(); @@ -584,6 +594,71 @@ impl ApplicationHandler for Viewer { } } } + + // ── ortho albedo windows (fine + close) ── + // Same drift/trigger discipline as the height tiers; the + // upload is a texture write (+ reallocation only on size + // change), so the main-thread cost matches the other tiers. + if let Some(ref mut cf) = bev_base.color_fine { + if let Some(data) = cf.try_recv() { + let (ox, oy, ex, ey, rot) = + cross_crs_world_origin_and_extent(&data.window.georef, &self.hm); + self.scene.as_mut().unwrap().upload_ortho_fine( + ox, + oy, + rot, + ex, + ey, + data.window.georef.cols as u32, + data.window.georef.rows as u32, + &data.window.rgba, + &data.mips, + ); + println!( + "ortho fine updated: {}×{} at {:.1}m/px", + data.window.georef.cols, + data.window.georef.rows, + data.window.georef.dx_meters + ); + } + if detail_allowed && !cf.computing { + let (lat, lon) = cam_wgs84(self.cam_pos, &self.hm); + if cf.needs_reload(lat, lon) && cf.try_trigger(lat, lon) { + println!("ortho fine reload triggered at lat={lat:.4} lon={lon:.4}"); + } + } + } + if let Some(ref mut cc) = bev_base.color_close { + if let Some(data) = cc.try_recv() { + let (ox, oy, ex, ey, rot) = + cross_crs_world_origin_and_extent(&data.window.georef, &self.hm); + self.scene.as_mut().unwrap().upload_ortho_close( + ox, + oy, + rot, + ex, + ey, + data.window.georef.cols as u32, + data.window.georef.rows as u32, + &data.window.rgba, + &data.mips, + ); + println!( + "ortho close updated: {}×{} at {:.1}m/px", + data.window.georef.cols, + data.window.georef.rows, + data.window.georef.dx_meters + ); + } + if detail_allowed && !cc.computing { + let (lat, lon) = cam_wgs84(self.cam_pos, &self.hm); + if cc.needs_reload(lat, lon) && cc.try_trigger(lat, lon) { + println!( + "ortho close reload triggered at lat={lat:.4} lon={lon:.4}" + ); + } + } + } } let surface: &wgpu::Surface = @@ -755,6 +830,16 @@ impl ApplicationHandler for Viewer { self.smooth_radius_m = presets[(cur + 1) % presets.len()]; return; } + // Toggle the streamed orthophoto albedo drape. The shader + // falls back to the procedural elevation ramp when off. + if kc == KeyCode::KeyT && event.state == winit::event::ElementState::Pressed { + self.ortho_enabled = !self.ortho_enabled; + eprintln!( + "ortho albedo: {}", + if self.ortho_enabled { "on" } else { "off" } + ); + return; + } // Debug: force close + fine tier reloads on the next frame. // Used to repro tier-swap memory peaks without flying. if kc == KeyCode::KeyR && event.state == winit::event::ElementState::Pressed { @@ -763,7 +848,13 @@ impl ApplicationHandler for Viewer { if let Some(ref mut fine) = bev_base.fine { fine.invalidate(); } - eprintln!("[vram] debug: close + fine tiers invalidated (R)"); + if let Some(ref mut cf) = bev_base.color_fine { + cf.invalidate(); + } + if let Some(ref mut cc) = bev_base.color_close { + cc.invalidate(); + } + eprintln!("[vram] debug: close + fine + ortho tiers invalidated (R)"); } return; } @@ -908,9 +999,11 @@ impl Viewer { /// the active tier set instead of letting the next allocation panic. /// /// Step-down order: + /// 0. Ortho albedo windows — cosmetic only, cheapest quality loss; both + /// textures reclaimed via `set_ortho_*_inactive`. /// 1. Fine tier — disabled, GPU memory reclaimed via `set_hm1m_inactive`. /// 2. Close tier — disabled, GPU memory reclaimed via `set_hm5m_inactive`. - /// 3. If both already disabled and we still OOM, give up — the base tier + /// 3. If all already disabled and we still OOM, give up — the base tier /// itself isn't safe to drop, so we let wgpu's default behaviour panic. fn poll_and_handle_oom(&mut self) { if !render_gpu::OOM_OBSERVED.load(std::sync::atomic::Ordering::Relaxed) { @@ -923,6 +1016,21 @@ impl Viewer { return; }; + // Step 0: drop both ortho albedo windows — terrain geometry is intact, + // the shader falls back to the procedural ramp. + if let Some(ref mut bev_base) = self.bev_base + && (bev_base.color_fine.is_some() || bev_base.color_close.is_some()) + { + eprintln!( + "[OOM #{count}] disabling ortho albedo windows — freeing both RGBA textures" + ); + scene.set_ortho_fine_inactive(); + scene.set_ortho_close_inactive(); + bev_base.color_fine = None; + bev_base.color_close = None; + return; + } + // Step 1: kill fine tier if it's still alive. if let Some(ref mut bev_base) = self.bev_base && bev_base.fine.is_some() @@ -1066,15 +1174,35 @@ impl Viewer { let fine_index = Arc::new(tile_index::build_tile_index(&dv.fine_tile_paths)); let close_index = Arc::new(tile_index::build_tile_index(&dv.close_tile_paths)); let base_index = Arc::new(tile_index::build_tile_index(&dv.base_tile_paths)); + let surface_index = Arc::new(tile_index::build_tile_index(&dv.surface_tile_paths)); + if !surface_index.is_empty() { + println!( + "DSM surface overlay active ({} tile(s))", + surface_index.len() + ); + } + let ortho_index = Arc::new(tile_index::build_tile_index(&dv.ortho_tile_paths)); + let lc_index = Arc::new(tile_index::build_tile_index(&dv.landcover_tile_paths)); + if !ortho_index.is_empty() { + println!( + "ortho albedo streaming active ({} mosaic(s), {} land-cover)", + ortho_index.len(), + lc_index.len() + ); + } let (cam_lat_d, cam_lon_d) = (dv.camera_lat, dv.camera_lon); bev_base = Some(BevBaseState::new( fine_index, close_index, base_index, + surface_index, + ortho_index, + lc_index, cam_lat_d, cam_lon_d, lat_rad, tier_radii, + ortho_radii(chosen_class), &hm, &mut scene, )); @@ -1182,10 +1310,14 @@ impl Viewer { fine_index, close_index, base_index, + Arc::new(vec![]), // single-file mode has no DSM surface overlay + Arc::new(vec![]), // …and no ortho albedo mosaics + Arc::new(vec![]), CAM_LAT, CAM_LON, lat_rad, tier_radii, + ortho_radii(chosen_class), &hm, &mut scene, )); diff --git a/src/viewer/tiers.rs b/src/viewer/tiers.rs index 1a52d54..d082f35 100644 --- a/src/viewer/tiers.rs +++ b/src/viewer/tiers.rs @@ -129,23 +129,51 @@ pub(super) struct TierData { pub(super) gpu_ao_u8: Vec, } -/// Per-tier channel state and drift-detection bookkeeping. +/// Worker-result bundle for an ortho albedo window: pre-packed RGBA bytes (in +/// `window.rgba`), CPU-generated mips, and the WGS84 centre for drift tracking. +pub(super) struct ColorData { + pub(super) window: dem_io::ColorWindow, + pub(super) mips: Vec<(u32, u32, Vec)>, + pub(super) centre_lat: f64, + pub(super) centre_lon: f64, +} + +/// Anything a streaming worker can deliver: exposes the WGS84 centre of the +/// loaded window so `StreamingTier` can do payload-agnostic drift bookkeeping. +pub(super) trait TierCentre { + fn centre(&self) -> (f64, f64); +} + +impl TierCentre for TierData { + fn centre(&self) -> (f64, f64) { + (self.centre_lat, self.centre_lon) + } +} + +impl TierCentre for ColorData { + fn centre(&self) -> (f64, f64) { + (self.centre_lat, self.centre_lon) + } +} + +/// Per-tier channel state and drift-detection bookkeeping, generic over the +/// worker's payload (`TierData` for height tiers, `ColorData` for ortho). /// /// `last_cx`/`last_cy` store WGS84 (lat, lon) in degrees. /// `drift_threshold_m` is stored in degrees (metres / M_PER_DEG). -pub(super) struct StreamingTier { +pub(super) struct StreamingTier { pub(super) tx: mpsc::SyncSender<(f64, f64)>, - rx: mpsc::Receiver, + rx: mpsc::Receiver, pub(super) computing: bool, last_cx: f64, last_cy: f64, drift_threshold_m: f64, } -impl StreamingTier { +impl StreamingTier { pub(super) fn new( tx: mpsc::SyncSender<(f64, f64)>, - rx: mpsc::Receiver, + rx: mpsc::Receiver, init_cx: f64, init_cy: f64, drift_threshold_m: f64, @@ -180,12 +208,13 @@ impl StreamingTier { /// Poll for a finished bundle. On success, clears `computing` and /// updates `last_cx`/`last_cy` from the bundle's centre coordinates. - pub(super) fn try_recv(&mut self) -> Option { + pub(super) fn try_recv(&mut self) -> Option { match self.rx.try_recv() { Ok(data) => { self.computing = false; - self.last_cx = data.centre_lat; - self.last_cy = data.centre_lon; + let (cx, cy) = data.centre(); + self.last_cx = cx; + self.last_cy = cy; Some(data) } Err(_) => None, @@ -222,11 +251,190 @@ pub(super) fn select_ifd(scales: &[f64], min_scale_m: f64, radius_m: f64, max_px scales.len().saturating_sub(1) } +/// `select_ifd` over mask-filtered `(ifd, scale)` pairs (from +/// `dem_io::ifd_overview_levels`) — needed for ortho mosaics whose IFD chain +/// interleaves transparency-mask IFDs with the overview pyramid, where a raw +/// scale index no longer equals the IFD index. +pub(super) fn select_overview_level( + levels: &[(usize, f64)], + min_scale_m: f64, + radius_m: f64, + max_px: u32, +) -> (usize, f64) { + for &(ifd, scale) in levels { + let window_px = (radius_m * 2.0 / scale) as u32; + if scale >= min_scale_m && window_px <= max_px { + return (ifd, scale); + } + } + levels.last().copied().unwrap_or((0, 1.0)) +} + +/// Ortho albedo window geometry per VRAM class. Target scales snap to the BEV +/// mosaics' 0.2·2^k overview pyramid inside `select_overview_level`; RGBA8 + +/// ⅓ mips puts the steady-state cost at ≈170 MB (Mid), ≈330 MB (High), ≈17 MB +/// (Low) for both windows together — reclaimed first by the OOM ladder. +#[derive(Clone, Copy, Debug)] +pub(super) struct OrthoRadii { + pub(super) fine_radius_m: f64, + pub(super) fine_min_scale_m: f64, + pub(super) fine_drift_m: f64, + pub(super) close_radius_m: f64, + pub(super) close_min_scale_m: f64, + pub(super) close_drift_m: f64, +} + +pub(super) fn ortho_radii(class: VramClass) -> OrthoRadii { + match class { + VramClass::High => OrthoRadii { + fine_radius_m: 3_000.0, + fine_min_scale_m: 0.75, + fine_drift_m: 800.0, + close_radius_m: 20_000.0, + close_min_scale_m: 6.0, + close_drift_m: 3_000.0, + }, + VramClass::Mid => OrthoRadii { + fine_radius_m: 2_000.0, + fine_min_scale_m: 0.75, + fine_drift_m: 600.0, + close_radius_m: 14_000.0, + close_min_scale_m: 6.0, + close_drift_m: 2_000.0, + }, + VramClass::Low => OrthoRadii { + fine_radius_m: 1_000.0, + fine_min_scale_m: 1.5, + fine_drift_m: 300.0, + close_radius_m: 8_000.0, + close_min_scale_m: 12.0, + close_drift_m: 1_500.0, + }, + } +} + +/// Spawn a background worker that streams camera-centred ortho albedo windows +/// (RGB orthophoto + land-cover material codes packed RGBA). Mirrors the height +/// tier workers: receives WGS84 `(lat, lon)`, resolves the ortho tile's CRS and +/// mask-filtered overview level, reads + converts the window, packs mips, sends. +fn spawn_color_worker( + ortho_index: Arc, + lc_index: Arc, + radius_m: f64, + min_scale_m: f64, + drift_m: f64, + label: &'static str, +) -> StreamingTier { + use super::tile_index::tiles_overlapping_wgs84; + use std::collections::HashMap; + use std::path::{Path, PathBuf}; + + let (tx, worker_rx) = mpsc::sync_channel::<(f64, f64)>(1); + let (worker_tx, rx) = mpsc::channel::(); + + std::thread::spawn(move || { + // Overview walking opens the file and touches every IFD header; cache + // per path so steady-state reloads skip it. + let mut levels_cache: HashMap> = HashMap::new(); + let levels_for = |cache: &mut HashMap>, + path: &Path| + -> Vec<(usize, f64)> { + if let Some(v) = cache.get(path) { + return v.clone(); + } + let v = dem_io::ifd_overview_levels(path).unwrap_or_else(|_| vec![(0, 1.0)]); + cache.insert(path.to_path_buf(), v.clone()); + v + }; + + while let Ok((lat, lon)) = worker_rx.recv() { + let overlapping = tiles_overlapping_wgs84(&ortho_index, lat, lon, radius_m); + let Some(&oi) = overlapping.first() else { + continue; + }; + let entry = &ortho_index[oi]; + let Ok(centre) = dem_io::crs::from_wgs84(lat, lon, &entry.crs_proj4) else { + continue; + }; + let rgb_levels = levels_for(&mut levels_cache, &entry.path); + let (rgb_ifd, rgb_scale) = + select_overview_level(&rgb_levels, min_scale_m, radius_m, GPU_SAFE_PX as u32); + + // Land cover: pick the overview whose scale sits closest to the + // chosen RGB scale so the nearest-resampling in extract_color_window + // is ~1:1. + let lc_entry = tiles_overlapping_wgs84(&lc_index, lat, lon, radius_m) + .first() + .map(|&i| &lc_index[i]); + let (lc_path, lc_ifd) = match lc_entry { + Some(e) => { + let lc_levels = levels_for(&mut levels_cache, &e.path); + let ifd = lc_levels + .iter() + .min_by(|a, b| { + (a.1 - rgb_scale) + .abs() + .partial_cmp(&(b.1 - rgb_scale).abs()) + .unwrap_or(std::cmp::Ordering::Equal) + }) + .map(|&(i, _)| i); + (Some(e.path.clone()), ifd) + } + None => (None, None), + }; + + let t0 = std::time::Instant::now(); + match dem_io::extract_color_window( + &entry.path, + lc_path.as_deref(), + centre, + radius_m, + rgb_ifd, + lc_ifd, + ) { + Ok(window) => { + let mips = render_gpu::gen_rgba_mip_bytes( + &window.rgba, + window.georef.cols, + window.georef.rows, + ); + eprintln!( + "[ortho-{label}] {}×{} at {:.1} m/px ({:.2?})", + window.georef.cols, + window.georef.rows, + window.georef.dx_meters, + t0.elapsed() + ); + if worker_tx + .send(ColorData { + window, + mips, + centre_lat: lat, + centre_lon: lon, + }) + .is_err() + { + break; + } + } + Err(e) => eprintln!("[ortho-{label}] window failed: {e}"), + } + } + }); + + // last centre (0, 0) guarantees the first drift check fires immediately. + StreamingTier::new(tx, rx, 0.0, 0.0, drift_m / M_PER_DEG) +} + /// Persistent state for BEV multi-tier streaming mode. pub(super) struct BevBaseState { pub(super) base: StreamingTier, // wide window, low resolution (IFD-2/1) pub(super) close: StreamingTier, // close window, 5 m/px (IFD-0) pub(super) fine: Option, // fine window, 1 m/px (1m tile IFD-0); None if no 1m tiles available + /// Ortho albedo streamers (fine + close windows); None when no ortho mosaic + /// is configured or its file is missing. + pub(super) color_fine: Option>, + pub(super) color_close: Option>, } impl BevBaseState { @@ -243,10 +451,14 @@ impl BevBaseState { fine_index: Arc, close_index: Arc, base_index: Arc, + surface_index: Arc, + ortho_index: Arc, + lc_index: Arc, cam_lat: f64, cam_lon: f64, lat_rad: f32, radii: TierRadii, + ortho: OrthoRadii, hm: &Arc, scene: &mut GpuScene, ) -> Self { @@ -463,6 +675,7 @@ impl BevBaseState { let (hm1m_tx, hm1m_worker_rx) = mpsc::sync_channel::<(f64, f64)>(1); let (hm1m_worker_tx, hm1m_rx) = mpsc::channel::(); let fine_idx = Arc::clone(&fine_index); + let surface_idx = Arc::clone(&surface_index); let lat_rad_1m = lat_rad; let recent_5m_w = recent_5m_fine; std::thread::spawn(move || { @@ -501,6 +714,21 @@ impl BevBaseState { { dem_io::fill_nodata_from_base(&mut raw1m, close_hm); } + // DSM overlay: composite trees/buildings over the bare-earth + // DTM wherever surface tiles cover this window. Normals, + // shadows and the bicubic march then all see the composite, + // so the canopy gets geometry and self-shadowing for free. + for &si in &tiles_overlapping_wgs84(&surface_idx, lat, lon, fine_radius_m) { + let e = &surface_idx[si]; + let Ok((set, snt)) = dem_io::crs::from_wgs84(lat, lon, &e.crs_proj4) + else { + continue; + }; + if let Ok(dsm) = extract_window(&e.path, (set, snt), fine_radius_m, e.ifd) + { + dem_io::composite_surface_over(&mut raw1m, &dsm, 6); + } + } dem_io::clamp_nodata_to_sea(&mut raw1m); let hm1m = Arc::new(raw1m); let normals = terrain::compute_normals_vector_par(&hm1m); @@ -544,6 +772,31 @@ impl BevBaseState { let effective_base_threshold = radii.base_drift_m.min(base_half_m * 0.5); let base_drift_deg = effective_base_threshold / M_PER_DEG; + // Ortho albedo streamers — only when an ortho mosaic actually resolved. + // (Land cover alone does nothing: material codes ride the ortho window.) + let (color_fine, color_close) = if ortho_index.is_empty() { + (None, None) + } else { + ( + Some(spawn_color_worker( + Arc::clone(&ortho_index), + Arc::clone(&lc_index), + ortho.fine_radius_m, + ortho.fine_min_scale_m, + ortho.fine_drift_m, + "fine", + )), + Some(spawn_color_worker( + ortho_index, + lc_index, + ortho.close_radius_m, + ortho.close_min_scale_m, + ortho.close_drift_m, + "close", + )), + ) + }; + BevBaseState { base: StreamingTier::new(base_tx, base_rx, cam_lat, cam_lon, base_drift_deg), close: StreamingTier::new( @@ -554,6 +807,8 @@ impl BevBaseState { effective_close_threshold / M_PER_DEG, ), fine, + color_fine, + color_close, } } } @@ -716,6 +971,61 @@ mod tests { assert!(hi.fine_radius_m >= mid.fine_radius_m && mid.fine_radius_m >= lo.fine_radius_m); } + // ortho_radii + + #[test] + fn ortho_radii_are_internally_ordered() { + for class in [VramClass::Low, VramClass::Mid, VramClass::High] { + let o = ortho_radii(class); + assert!( + o.fine_radius_m < o.close_radius_m, + "{class:?}: fine ortho window must sit inside the close one" + ); + assert!( + o.fine_min_scale_m < o.close_min_scale_m, + "{class:?}: fine window must target finer texels" + ); + // A drift ≥ radius would let the camera exit before a reload fires. + assert!(o.fine_drift_m < o.fine_radius_m, "{class:?} fine drift"); + assert!(o.close_drift_m < o.close_radius_m, "{class:?} close drift"); + // The selected window must fit GPU_SAFE_PX at the target scale. + assert!( + (o.fine_radius_m * 2.0 / o.fine_min_scale_m) as usize <= GPU_SAFE_PX, + "{class:?}: fine ortho window exceeds the GPU texture cap" + ); + assert!( + (o.close_radius_m * 2.0 / o.close_min_scale_m) as usize <= GPU_SAFE_PX, + "{class:?}: close ortho window exceeds the GPU texture cap" + ); + } + } + + // select_overview_level + + #[test] + fn select_overview_level_keeps_mask_filtered_ifd_indices() { + // Levels as ifd_overview_levels reports them for the BEV RGB mosaic: + // IFD 1 is a transparency mask, so the pyramid jumps 0 → 2. + let levels = [ + (0usize, 0.2), + (2, 0.4), + (3, 0.8), + (4, 1.6), + (5, 3.2), + (6, 6.4), + ]; + // Fine ortho (≥0.75 m, 2 km radius): 0.8 m level → IFD 3, not index 2. + assert_eq!(select_overview_level(&levels, 0.75, 2_000.0, 8192), (3, 0.8)); + // Close ortho (≥6 m, 14 km radius): 6.4 m level → IFD 6. + assert_eq!(select_overview_level(&levels, 6.0, 14_000.0, 8192), (6, 6.4)); + // Nothing fits → coarsest level. + assert_eq!( + select_overview_level(&levels, 0.2, 1.0e9, 8192), + (6, 6.4), + "fall through to the coarsest pair" + ); + } + // select_ifd #[test] @@ -856,4 +1166,168 @@ mod tests { let mut t = streaming_tier(0.0, 0.0, 100.0); assert!(t.try_recv().is_none()); } + + /// Offscreen end-to-end render of the DSM + ortho pipeline against the real + /// Tirol tiles: DTM window + DSM composite as geometry, RGB ortho + land + /// cover as albedo, placed with the production `cross_crs_world_origin_and_extent`, + /// rendered through the real shader, read back, and saved to + /// `/tmp/offscreen_dsm_ortho.png` for visual inspection. + /// + /// Skips when the multi-GB local tiles (gitignored) or a GPU adapter are + /// absent, so CI stays green. Run with: + /// `cargo test --release --bin dem_renderer offscreen -- --nocapture` + #[test] + fn offscreen_dsm_ortho_render_smoke() { + let tiles = std::path::Path::new(env!("CARGO_MANIFEST_DIR")).join("tiles"); + let dtm = tiles.join("big_size/CRS3035RES50000mN2650000E4450000.tif"); + let dsm = tiles.join("big_size/ALS_DSM_CRS3035RES50000mN2650000E4450000.tif"); + let rgb = tiles.join("color/2019470_Mosaik_RGB.tif"); + let lc = tiles.join("color/2022470_Mosaik_LC.tif"); + if ![&dtm, &dsm, &rgb, &lc].iter().all(|p| p.exists()) { + eprintln!("skipping — local Tirol tiles not present"); + return; + } + + // Mayrhofen valley: inside DSM, ortho and land-cover coverage. + let (lat, lon) = (47.16, 11.86); + let radius = 2_000.0; + let (width, height) = (800u32, 600u32); + + // Geometry: DTM window with the DSM surface composited over it — the + // exact sequence the fine worker runs. + let proj4 = dem_io::crs::tile_proj4(&dtm).expect("DTM CRS"); + let (e, n) = dem_io::crs::from_wgs84(lat, lon, &proj4).expect("project"); + let mut hm = extract_window(&dtm, (e, n), radius, 0).expect("DTM window"); + let dsm_win = extract_window(&dsm, (e, n), radius, 0).expect("DSM window"); + dem_io::composite_surface_over(&mut hm, &dsm_win, 6); + dem_io::clamp_nodata_to_sea(&mut hm); + let normals = terrain::compute_normals_vector_par(&hm); + let shadow = terrain::compute_shadow_vector_par_with_azimuth( + &hm, + 180.0_f32.to_radians(), + 45.0_f32.to_radians(), + 200.0, + ); + let ao = vec![1.0f32; hm.rows * hm.cols]; + + let ctx = render_gpu::GpuContext::new(); + let mut scene = GpuScene::new(ctx, &hm, &normals, &shadow, &ao, width, height); + + // Albedo: real ortho + land cover window, placed with the production + // placement function (georef stub → world rect over this heightmap). + let rgb_proj4 = dem_io::crs::tile_proj4(&rgb).expect("RGB CRS"); + let centre = dem_io::crs::from_wgs84(lat, lon, &rgb_proj4).expect("project ortho"); + let levels = dem_io::ifd_overview_levels(&rgb).expect("RGB levels"); + let (rgb_ifd, rgb_scale) = select_overview_level(&levels, 0.75, radius, 8192); + let lc_levels = dem_io::ifd_overview_levels(&lc).expect("LC levels"); + let lc_ifd = lc_levels + .iter() + .min_by(|a, b| { + (a.1 - rgb_scale) + .abs() + .partial_cmp(&(b.1 - rgb_scale).abs()) + .unwrap() + }) + .map(|&(i, _)| i); + let win = dem_io::extract_color_window(&rgb, Some(&lc), centre, radius, rgb_ifd, lc_ifd) + .expect("color window"); + let mips = render_gpu::gen_rgba_mip_bytes(&win.rgba, win.georef.cols, win.georef.rows); + let (ox, oy, ex, ey, rot) = cross_crs_world_origin_and_extent(&win.georef, &hm); + scene.upload_ortho_fine( + ox, + oy, + rot, + ex, + ey, + win.georef.cols as u32, + win.georef.rows as u32, + &win.rgba, + &mips, + ); + + // Camera 400 m above the valley floor, pitched down toward the north. + let cam_x = (e - hm.crs_origin_x) as f32; + let cam_y = (hm.crs_origin_y - n) as f32; + let gi = (cam_y as usize).min(hm.rows - 1) * hm.cols + (cam_x as usize).min(hm.cols - 1); + let ground = hm.data[gi]; + let origin = [cam_x, cam_y, ground + 400.0]; + let look_at = [cam_x, cam_y - 900.0, ground]; + + let ctx = scene.get_gpu_ctx().clone(); + let mut enc = ctx + .device + .create_command_encoder(&wgpu::CommandEncoderDescriptor { label: None }); + scene.dispatch_frame( + &mut enc, + origin, + look_at, + 70.0, + width as f32 / height as f32, + [0.3, -0.5, 0.8], + 1.0, // step_m + 20_000.0, // t_max + 0, // ao_mode off + 1, // shadows on + 0, // fog off — keep colors unmixed for the assertions + 2, // vat Mid + 0, // lod Ultra + 0.0, // no bicubic + 0, // align viz off + 1, // ortho_mode ON + ); + // Readback: copy the BGRA output into a mappable staging buffer. + let out_size = (width * height * 4) as u64; + let staging = ctx.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("offscreen_staging"), + size: out_size, + usage: wgpu::BufferUsages::MAP_READ | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }); + enc.copy_buffer_to_buffer(scene.get_output_buffer(), 0, &staging, 0, out_size); + ctx.queue.submit([enc.finish()]); + let slice = staging.slice(..); + slice.map_async(wgpu::MapMode::Read, |r| r.expect("map")); + let _ = ctx.device.poll(wgpu::PollType::Wait { + submission_index: None, + timeout: None, + }); + let bgra = slice.get_mapped_range().to_vec(); + staging.unmap(); + + // BGRA → RGBA and save for human inspection. + let mut rgba_img = vec![0u8; bgra.len()]; + for (dst, src) in rgba_img.chunks_exact_mut(4).zip(bgra.chunks_exact(4)) { + dst[0] = src[2]; + dst[1] = src[1]; + dst[2] = src[0]; + dst[3] = 255; + } + image::RgbaImage::from_raw(width, height, rgba_img.clone()) + .expect("image dims") + .save("/tmp/offscreen_dsm_ortho.png") + .expect("save png"); + eprintln!("offscreen render saved to /tmp/offscreen_dsm_ortho.png"); + + // Sanity: the frame must contain terrain (not all sky) and the terrain + // must carry ortho color variety (not the flat procedural ramp, whose + // greens/grays are far less diverse than a photo mosaic). + let n_px = (width * height) as usize; + let sky = rgba_img + .chunks_exact(4) + .filter(|p| p[2] > p[0] + 30 && p[2] > 120) + .count(); + assert!( + sky < n_px * 3 / 4, + "frame is mostly sky — camera placement or march broken ({sky}/{n_px})" + ); + let mut distinct = std::collections::HashSet::new(); + for p in rgba_img.chunks_exact(4) { + distinct.insert((p[0] >> 3, p[1] >> 3, p[2] >> 3)); + } + assert!( + distinct.len() > 300, + "terrain colors too uniform ({}) — ortho albedo likely not applied", + distinct.len() + ); + } } From 6c9ef62ef1b2366304c1ab7699ce7ff6d5c45b9f Mon Sep 17 00:00:00 2001 From: Nikita Zavartsev Date: Fri, 12 Jun 2026 20:44:51 +0200 Subject: [PATCH 6/9] docs: [dsm-rgb-landcover-tirol] document DSM/ortho layers in CLAUDE.md, add rendering-approaches learning note CLAUDE.md gains the DSM + ortho/land-cover architecture section (bind group entries 20-23, material code ladder, empirically pinned LC classes, OOM step 0), the new dem_io modules, and the updated tiers.rs inventory. The learnings note compares polygons vs SDF raymarching vs this renderer's heightfield raymarch and maps the trade-offs onto shader_texture.wgsl. Co-Authored-By: Claude Fable 5 --- CLAUDE.md | 17 +- .../polygons-vs-sdf-vs-heightfield.md | 163 ++++++++++++++++++ 2 files changed, 177 insertions(+), 3 deletions(-) create mode 100644 docs/learnings/polygons-vs-sdf-vs-heightfield.md diff --git a/CLAUDE.md b/CLAUDE.md index 96ea476..ce78e84 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -49,7 +49,7 @@ Settings persist to `dirs::config_dir() / dem_renderer / config.toml` (macOS: `~ |---|---| | `src/viewer/mod.rs` | `Viewer` (ApplicationHandler), WASD+mouse, sun animation, tile streaming dispatch, key bindings | | `src/viewer/scene_init.rs` | `prepare_scene_with_ctx` (single-tile / projected-CRS streaming), `prepare_demo_scene_with_ctx` (N×M Copernicus base, camera-centered crop to `GPU_SAFE_PX`), `compute_ao_cropped` | -| `src/viewer/tiers.rs` | `StreamingTier` (drift-detected reload), `BevBaseState` (base/close/fine workers), `TierRadii` + `tier_radii(VramClass)` preset mapping, `cross_crs_world_origin_and_extent`, `select_ifd`, `cap_to_gpu_limit` | +| `src/viewer/tiers.rs` | `StreamingTier` (drift-detected reload; payload-generic), `BevBaseState` (base/close/fine + color_fine/color_close ortho workers), `TierRadii` + `tier_radii(VramClass)`, `OrthoRadii` + `ortho_radii(VramClass)`, `cross_crs_world_origin_and_extent`, `select_ifd`, `select_overview_level`, `cap_to_gpu_limit` | | `src/viewer/tile_index.rs` | `TileEntry` + `TileIndex` — discover WGS84 bounds of multiple tiles per tier; `tiles_overlapping_wgs84` | | `src/viewer/geo.rs` | `latlon_to_tile_metres` (handles both geographic and projected CRSes), `sun_position` | | `src/viewer/hud_renderer.rs` | glyphon HUD overlay, sun indicator, settings panel | @@ -83,8 +83,9 @@ dem_renderer/ ├── crates/ │ ├── dem_io/src/ │ │ ├── lib.rs -│ │ ├── heightmap.rs # Heightmap (f32 data), parse_bil, fill_nodata, fill_nodata_from_base -│ │ ├── geotiff.rs # parse_geotiff_auto, extract_window, ifd_scales, tile_bounds_wgs84, tile_centre_crs +│ │ ├── heightmap.rs # Heightmap (f32 data), parse_bil, fill_nodata, fill_nodata_from_base, composite_surface_over (DSM-over-DTM) +│ │ ├── geotiff.rs # parse_geotiff_auto, extract_window, ifd_scales, ifd_overview_levels (mask-aware), tile_bounds_wgs84, tile_centre_crs, GDAL_NODATA(42113)-aware NoData predicate +│ │ ├── color.rs # extract_color_window: ortho RGB (JPEG YCbCr→RGB) + land-cover material codes packed RGBA8; landcover_histogram │ │ ├── grid.rs # assemble_grid (N×M), load_grid_from_paths, crop, stitch_windows[_geographic] │ │ ├── crs.rs # tile_proj4 / to_wgs84 / from_wgs84 / is_geographic / epsg_towgs84 / read_raw_crs_tags (proj4rs + proj4wkt + crs-definitions; WKT → inline-GeoKey → EPSG fallback chain) │ │ ├── overview.rs # ensure_overview_cache: build .tmp_dem_pre_calc_*.tif from large single-IFD tiles (copies source CRS tags verbatim so cache is self-describing) @@ -175,9 +176,19 @@ For projected single-file mode and demo mode, the viewer spawns three background | 8 / 9 | AO texture (R8Unorm) + sampler | | 10–14 | 5 m close tier — heightmap (R32Float) + samp + normals (Rg16Snorm) + samp + shadow buffer | | 15–19 | 1 m fine tier — same layout as 5 m | +| 20/21 | Fine ortho albedo (Rgba8Unorm, mipped; RGB = orthophoto, A = land-cover material code) + sampler | +| 22/23 | Close ortho albedo + sampler | `scene::bind_group::rebuild_bind_group()` rebuilds the BG whenever a tier's texture / buffer is recreated (size grows). Steady-state reloads do not allocate — `write_texture` / `write_buffer` overwrite in place. +### DSM + Ortho/Land-Cover Layers (Tirol demo) + +- **DSM (surface model, trees/buildings)**: `demo_view.surface_tile_paths` (default: the BEV ALS DSM, EPSG:3035 1 m). The fine worker extracts a DSM window with the same centre/radius as its DTM window and paints it on with `dem_io::composite_surface_over` (6 px feather; DTM stays the canvas so coverage never shrinks). Normals/shadows/bicubic then run on the composite — canopy and buildings get geometry and self-shadowing with zero shader changes. The DSM's `GDAL_NODATA = +3.4028235e38` is caught by the generalized NoData predicate (tag 42113 + `≥1e38` guard) — without it those cells parse as 3.4e38 m elevations and poison `max_terrain_h`. +- **Ortho albedo**: `demo_view.ortho_tile_paths` + `landcover_tile_paths` (defaults: BEV 20 cm RGB + LC mosaics, EPSG:31255). Two `StreamingTier` workers (fine ≈0.8 m/px, close ≈6.4 m/px; radii per `ortho_radii(VramClass)`) read windows via `extract_color_window` — JPEG YCbCr chunks decoded by the tiff crate's zune path, BT.601-converted on the worker, land cover nearest-resampled onto the ortho grid. One RGBA8 texture per window: RGB = albedo, **A = material code ladder** (water 255 / high veg 192 / med veg 128 / building 64 / other 0 — codes ≥64 apart so bilinear boundary mixes degrade gracefully; mip gen box-filters RGB but takes nearest A). LC classes pinned empirically (`examples/inspect_color.rs`): 1=high veg, 2=ground/ice, 3=med veg, 4=buildings, 5=water, 6=low veg, 15=NoData. Note the LC mosaic is **Gray(4)** — 4-bit nibble-packed chunks, unpacked in `extract_u8_window`; the RGB mosaic interleaves **mask IFDs** in its chain, hence `ifd_overview_levels` (NewSubfileType-filtered) instead of raw `ifd_scales`. +- **Shader**: hit-point albedo = procedural ramp → close ortho → fine ortho (edge-smoothstep over `BLEND_MARGIN`, per-window rotation like the height tiers, distance-heuristic mip). Where ortho wins, the brightness floor is lifted (`mix(brightness, clamp(brightness, 0.55, 1.0), w·0.7)`) because the photo carries baked sun shading — full diffuse on top double-shades north faces. Water (A ≥ ~0.9) blends toward a lake tint + mirror-direction sun glint. `T` toggles the drape (`ortho_mode` uniform). +- **OOM ladder step 0**: both ortho textures are dropped before the fine tier — cosmetic-first degradation. +- Missing files degrade gracefully everywhere (`build_tile_index` skips absent paths; empty index → no workers spawned). + --- ## Launcher UI — Development Rules diff --git a/docs/learnings/polygons-vs-sdf-vs-heightfield.md b/docs/learnings/polygons-vs-sdf-vs-heightfield.md new file mode 100644 index 0000000..4b20f8b --- /dev/null +++ b/docs/learnings/polygons-vs-sdf-vs-heightfield.md @@ -0,0 +1,163 @@ +# Polygons vs SDFs vs Heightfield Raymarching + +Two fundamentally different answers to the same problem: *"given a 3D scene and a +camera, what color is each pixel?"* The split is **where the geometry lives** and +**which direction you traverse** the projection. + +This renderer uses the third option — heightfield raymarching — which borrows the +*structure* of SDF raymarching but not its distance guarantee. See the bottom +section for how that maps onto `crates/render_gpu/src/shader_texture.wgsl`. + +--- + +## Polygons (rasterization) + +**Representation:** the surface is an explicit mesh — a list of vertices (x,y,z) +stitched into triangles. A sphere isn't "a sphere," it's 5,000 triangles +approximating one. The geometry is *enumerated*: you store every piece of surface +up front. + +**Algorithm — you loop over geometry, not pixels:** + +1. For each triangle, transform its 3 vertices by the model→view→projection + matrices. This lands each vertex in 2D screen space (plus a depth). +2. **Rasterize:** figure out which pixels the triangle covers. Hardware does an + edge-function test per pixel — "is this pixel inside the 3 edges?" +3. For each covered pixel, run a **fragment shader** to compute color (texture + lookup, lighting). +4. **Z-buffer:** each pixel stores the depth of the closest triangle seen so far. + A new fragment writes only if it's nearer. This resolves occlusion without + sorting. + +Mental model: rasterization is **scatter** — geometry is thrown *at* the screen, +and each triangle "lights up" the pixels it lands on. You iterate triangles, and +pixels are a side effect. + +**Why hardware loves it:** the whole pipeline is fixed-function and massively +parallel per-triangle and per-pixel. Cost scales with `triangles × pixels_covered`. +A triangle that covers 0 pixels (off-screen, back-facing) is cheap to reject. This +is what GPUs were *built* for — decades of silicon dedicated to exactly this loop. + +**Cost characteristics:** +- More detail = more triangles = more vertex work + memory. A film-quality asset + can be millions of triangles. +- Smooth/curved surfaces are *faked* — tessellate finer until the facets are + sub-pixel. +- Hard things: true reflections, refraction, soft shadows, global illumination — + rasterization only knows "this triangle covers this pixel," with no cheap way to + ask "what's along this arbitrary ray?" You bolt those on with tricks (shadow + maps, screen-space reflections) or move to ray tracing. + +--- + +## SDFs (signed distance fields, via raymarching) + +**Representation:** the surface is *implicit* — defined by a function `f(p)` that +returns, for any point `p` in space, the **distance to the nearest surface**. +Signed: negative inside the object, positive outside, zero *on* the surface. A +sphere of radius `r` at origin is literally: + +``` +f(p) = length(p) - r +``` + +That one line *is* the sphere — exactly, at infinite resolution, no triangles. You +don't store geometry; you store a *rule* for measuring distance to it. + +**Algorithm — you loop over pixels, not geometry (sphere tracing):** + +1. For each pixel, shoot a ray from the camera through it. +2. Evaluate `f(p)` at the current position. It returns `d` — distance to the + nearest surface *in any direction*. +3. That `d` is a **safety radius**: nothing is closer than `d`, so you can leap + forward exactly `d` along the ray without risk of skipping through anything. +4. Repeat. As you approach a surface, `d → 0` and steps shrink automatically. When + `d < epsilon`, you've hit. + +Mental model: raymarching is **gather** — each pixel asks "what's the first thing +along my ray?" You iterate pixels, and geometry is queried. The genius of the SDF +is that the distance value *tells you how far you're allowed to jump*, so you take +big steps in empty space and tiny careful steps near surfaces. That's why it's +called *sphere tracing* — at each point an empty sphere of radius `d` is safe to +cross. + +**Why it's powerful:** +- **Infinite precision** — no facets, curves are exact. +- **Cheap CSG** — union is `min(f1, f2)`, intersection is `max`, subtraction is + `max(f1, -f2)`. Blending/morphing is a smooth `min`. This is why Shadertoy demos + build whole worlds in a few lines. +- **Free side effects** — because you can evaluate `f(p)` anywhere, you get soft + shadows, AO, and reflections nearly for free (march another ray; AO ≈ sample `f` + at a few points along the normal). +- **Surface normal** for lighting is just the gradient of `f` (finite + differences) — no stored normals needed. + +**Cost characteristics:** +- Cost scales with `pixels × steps_per_ray × cost_of_f`. Complex scenes mean + evaluating `f` (a `min` over many primitives) dozens of times per pixel per ray. +- Grazing angles and thin/high-frequency features are the enemy — the ray creeps + in tiny steps, or overshoots and misses. +- Doesn't map to fixed-function hardware — all general compute shader work, no + dedicated silicon. + +--- + +## The core differences, side by side + +| | Polygons (raster) | SDF (raymarch) | +|---|---|---| +| Geometry | Explicit (stored triangles) | Implicit (a function) | +| Traversal | Scatter: loop geometry → pixels | Gather: loop pixels → query geometry | +| Resolution | Faceted, finite | Exact, infinite | +| Occlusion | Z-buffer | First hit along ray | +| Cost driver | triangle count × coverage | steps per ray × `f` complexity | +| Hardware | Dedicated fixed-function pipeline | General compute | +| Shadows/reflections/AO | Bolted on with tricks | Nearly free (march more rays) | +| Sweet spot | Detailed authored assets, games | Procedural/CSG scenes, smooth blends | +| Pain point | Curves, secondary rays | Thin features, grazing rays, `f` cost | + +--- + +## Where this renderer sits: heightfield raymarching + +We do the **gather** style (loop pixels, march rays) like SDFs — but **not** with a +distance function. Our `f(p)` doesn't return distance-to-nearest-surface; it +returns terrain height at an (x,y) and we test `ray.z < height`. That's a +**heightfield raymarch**. + +The difference matters precisely at the step-size problem: + +- A true SDF gives you a *guaranteed safe jump distance in all 3D directions*. +- Our heightfield only knows **vertical clearance** (`pos.z - h`) — how far the ray + is *above* the terrain straight down. That's a weaker guarantee: if terrain rises + faster horizontally than the step moves us, vertical clearance lies and the ray + punches through a ridge. + +In `shader_texture.wgsl`: +- Adaptive step (line ~586): `t += max((pos.z - h) * sphere_factor, lod_min_step);` + — step scales by vertical clearance ("sphere tracing" comment, 0.5 safety + factor). `sphere_factor` is quality-controlled (Ultra 0.1 → Low 0.8). Bigger + factor = more overshoot = thin ridges disappear on Low. +- `lod_min_step` is the floor that prevents the ray stalling near the surface. +- Once the march detects it crossed below the surface, `binary_search_hit(t_prev, + t, dir, 10)` (line ~573) refines the hit with **10 iterations** of bisection + between `t_prev` (last point above) and `t` (first below). + +So we get the best-fit structure for terrain — a heightfield query is *far* cheaper +than a general SDF (one texture sample vs. a `min` over primitives) — at the cost +of that weaker stepping guarantee, patched with the minimum-step floor and the +binary-search refinement. + +The "arc" artefacts in the foreground are the step-size error made visible: all +rays with the same step count share the same overshoot distance, creating +concentric iso-step contours. + +--- + +## The third giant: ray-traced polygons (RTX) + +Same "gather" traversal as raymarching, but instead of stepping along the ray +evaluating a function, it intersects the ray analytically against triangles using a +BVH acceleration structure. It's the convergence point: explicit geometry like +rasterization, ray-per-pixel traversal like raymarching — and the part modern GPUs +now accelerate in dedicated hardware. From efd12316f6a00b89ff8b4939b922776de7f7a0ba Mon Sep 17 00:00:00 2001 From: Nikita Zavartsev Date: Fri, 12 Jun 2026 21:03:12 +0200 Subject: [PATCH 7/9] refactor: [dsm-rgb-landcover-tirol] allow identity_op in rgba mip test for readable texel sum cargo clippy --all-targets -D warnings flags the deliberate (0 + 40 + 80 + 120) / 4 in the mip-average assertion. Keep the four contributing texels visible instead of collapsing to a bare 60, matching the existing rationale'd allow on the f16 mip shape test. Co-Authored-By: Claude Fable 5 --- crates/render_gpu/src/lib.rs | 3 +++ 1 file changed, 3 insertions(+) diff --git a/crates/render_gpu/src/lib.rs b/crates/render_gpu/src/lib.rs index aa210d6..ccbf189 100644 --- a/crates/render_gpu/src/lib.rs +++ b/crates/render_gpu/src/lib.rs @@ -262,6 +262,9 @@ mod tests { assert_eq!(mips[1].2.len(), 4); } + // The expected value is written as `(0 + 40 + 80 + 120) / 4` so the four + // contributing texels stay visible, instead of collapsing to a bare `60`. + #[allow(clippy::identity_op)] #[test] fn rgba_mip_averages_rgb_but_takes_nearest_alpha() { // 2×2 base where the four texels have RGB 0/40/80/120 and alphas that From 2d5c9a0adc136823d05ed4f722afedb68af210a8 Mon Sep 17 00:00:00 2001 From: Nikita Zavartsev Date: Sat, 13 Jun 2026 13:20:26 +0200 Subject: [PATCH 8/9] fix: [dsm-rgb-landcover-tirol] sharpen DSM buildings/trees and gate ortho black coverage Two rendering fixes for the DSM + ortho draping: 1. Buildings and trees rendered like smooth mountains. The fine tier now carries the DSM, but the Catmull-Rom bicubic normal smoothing (meant for bare-earth ridgelines within smooth_radius_m) rounded every building wall and canopy edge into a hill. Gate the smoothing by world-slope magnitude: blend back to the raw per-texel Sobel normal where the gradient is wall-steep (smoothstep 1.0->2.5 m/m) so man-made faces and canopy edges stay crisp while gentle terrain stays bicubic-smooth. 2. Big black spot where the ortho coverage ends. The BEV mosaic fills areas with no photo coverage with pure black (0,0,0), which draped onto the terrain. Gate each ortho sample by luminance (smoothstep 2->16) so near-black no-data fades out and the procedural terrain colour shows through; real ortho pixels stay well above the threshold. Co-Authored-By: Claude Fable 5 --- crates/render_gpu/src/shader_texture.wgsl | 44 ++++++++++++++++------- 1 file changed, 32 insertions(+), 12 deletions(-) diff --git a/crates/render_gpu/src/shader_texture.wgsl b/crates/render_gpu/src/shader_texture.wgsl index 45a2f57..30451bc 100644 --- a/crates/render_gpu/src/shader_texture.wgsl +++ b/crates/render_gpu/src/shader_texture.wgsl @@ -733,18 +733,29 @@ fn main(@builtin(global_invocation_id) gid: vec3) { let a1h = apply_convergence_1m(lx1, ly1); let fine_uv = vec2(a1h.x / cam.hm1m_extent_x, a1h.y / cam.hm1m_extent_y); - // Within smooth_radius: derive normal analytically from bicubic gradient - // (C1-continuous → no slope discontinuities at cell boundaries). - // Outside: fall back to the pre-computed Rg16Snorm normal texture. + // Raw per-texel Sobel normal (sharp; computed on the DSM-over-DTM + // composite so building walls and canopy edges keep their steps). + let n1_rg = textureSampleLevel(hm1m_normal_tex, hm1m_normal_samp, fine_uv, 0.0).rg; + let n1_raw = normalize(vec3(n1_rg.x, n1_rg.y, sqrt(max(0.0, 1.0 - dot(n1_rg, n1_rg))))); + + // Within smooth_radius: derive a C1-continuous normal from the + // bicubic gradient (no slope discontinuities on natural terrain). + // BUT Catmull-Rom rounds sharp steps into hills, which turns the + // DSM's buildings and trees into melted bumps. Blend back to the + // raw normal where the gradient is wall-steep so man-made / canopy + // edges stay crisp while gentle terrain stays smooth. let dist_to_cam_hit = length(pos.xy - cam.origin.xy); var n1: vec3; if dist_to_cam_hit < cam.smooth_radius_m { let hg1 = sample_h_grad_bicubic_1m(a1h.x, a1h.y); // surface normal from gradient: N = normalize(-dh/dlx, -dh/dly, 1) - n1 = normalize(vec3(-hg1.y, -hg1.z, 1.0)); + let n_smooth = normalize(vec3(-hg1.y, -hg1.z, 1.0)); + // world slope magnitude (m/m): ~0.3–0.8 alpine, ≫2 at a wall. + let slope = length(vec2(hg1.y, hg1.z)); + let sharp = smoothstep(1.0, 2.5, slope); + n1 = normalize(mix(n_smooth, n1_raw, sharp)); } else { - let n1_rg = textureSampleLevel(hm1m_normal_tex, hm1m_normal_samp, fine_uv, 0.0).rg; - n1 = normalize(vec3(n1_rg.x, n1_rg.y, sqrt(max(0.0, 1.0 - dot(n1_rg, n1_rg))))); + n1 = n1_raw; } let dx1 = cam.hm1m_extent_x / f32(cam.hm1m_cols); @@ -883,17 +894,26 @@ fn main(@builtin(global_invocation_id) gid: vec3) { if oc.weight > 0.0 { let mip_c = min(mip_o, f32(textureNumLevels(ortho_close_tex) - 1u)); let s = textureSampleLevel(ortho_close_tex, ortho_close_samp, oc.uv, mip_c); - albedo = mix(albedo, s.rgb * 255.0, oc.weight); - material = mix(material, s.a, oc.weight); - ortho_w = max(ortho_w, oc.weight); + // The BEV mosaic fills areas with no photo coverage with pure + // black; without a gate those drape as a black hole on the + // terrain. Real ortho pixels are essentially never near-black, + // so fade the drape out by luminance and let the procedural + // colour show through where the mosaic has no data. + let cover = smoothstep(2.0, 16.0, dot(s.rgb * 255.0, vec3(0.299, 0.587, 0.114))); + let w = oc.weight * cover; + albedo = mix(albedo, s.rgb * 255.0, w); + material = mix(material, s.a, w); + ortho_w = max(ortho_w, w); } let ofs = ortho_fine_sample_at(pos.x, pos.y); if ofs.weight > 0.0 { let mip_f = min(mip_o, f32(textureNumLevels(ortho_fine_tex) - 1u)); let s = textureSampleLevel(ortho_fine_tex, ortho_fine_samp, ofs.uv, mip_f); - albedo = mix(albedo, s.rgb * 255.0, ofs.weight); - material = mix(material, s.a, ofs.weight); - ortho_w = max(ortho_w, ofs.weight); + let cover = smoothstep(2.0, 16.0, dot(s.rgb * 255.0, vec3(0.299, 0.587, 0.114))); + let w = ofs.weight * cover; + albedo = mix(albedo, s.rgb * 255.0, w); + material = mix(material, s.a, w); + ortho_w = max(ortho_w, w); } } From 8c752bce104813f02b7bf529ecfbb6e7ddb0d5fb Mon Sep 17 00:00:00 2001 From: Nikita Zavartsev Date: Sat, 13 Jun 2026 13:20:37 +0200 Subject: [PATCH 9/9] test: [dsm-rgb-landcover-tirol] cover color reader, cross-CRS placement and IFD enumeration The color/ortho readers and placement math only ran against the multi-GB BEV mosaics, which CI never has, so they showed as uncovered and dropped patch coverage below the 60% gate. Add synthetic-fixture tests that exercise the same production paths without the real tiles: - color_synth.rs: tiny in-process RGB + single-band GeoTIFFs drive extract_color_window (RGB pass-through, land-cover resample, material packing, georef stub) and landcover_histogram (color.rs 18% -> 95%). - tiers.rs inline tests for cross_crs_world_origin[_and_extent]: same-CRS exact placement plus both cross-CRS branches (geographic base, projected base) through real proj4 transforms (tiers.rs 28% -> 69%). - synthetic.rs: ifd_overview_levels over a multi-IFD overview cache (geotiff.rs 40% -> 98%). Patch coverage 56% -> 81% under the CI condition (tiles absent). Co-Authored-By: Claude Fable 5 --- crates/dem_io/tests/color_synth.rs | 174 +++++++++++++++++++++++++++++ crates/dem_io/tests/synthetic.rs | 34 ++++++ src/viewer/tiers.rs | 116 +++++++++++++++++++ 3 files changed, 324 insertions(+) create mode 100644 crates/dem_io/tests/color_synth.rs diff --git a/crates/dem_io/tests/color_synth.rs b/crates/dem_io/tests/color_synth.rs new file mode 100644 index 0000000..9f21322 --- /dev/null +++ b/crates/dem_io/tests/color_synth.rs @@ -0,0 +1,174 @@ +//! Synthetic-fixture tests for the color-window reader. +//! +//! `extract_color_window` and `landcover_histogram` only run end-to-end against +//! the multi-GB BEV mosaics, which CI never has — so the worker decode loop, the +//! land-cover resample and the material mapping showed as uncovered. These build +//! tiny RGB + single-band GeoTIFFs in-process (no JPEG / 4-bit needed: the RGB +//! photometric path passes through, the Gray8 path exercises the byte reader) +//! and assert the decoded window, the packed material codes, and the histogram. + +use std::fs::File; +use std::io::BufWriter; +use std::path::Path; + +use dem_io::{MATERIAL_HIGH_VEG, MATERIAL_NONE, MATERIAL_WATER}; +use tiff::encoder::{TiffEncoder, colortype}; +use tiff::tags::Tag; + +/// GeoKeyDirectory for a projected CRS by EPSG (GeoKey 3072). +fn projected_dir(epsg: u16) -> Vec { + vec![1, 1, 0, 2, 1024, 0, 1, 1, 3072, 0, 1, epsg] +} + +/// Write a 3-band RGB byte GeoTIFF (PhotometricInterpretation = RGB, so the +/// reader takes the pass-through branch rather than YCbCr conversion). +fn write_rgb(path: &Path, cols: u32, rows: u32, origin: (f64, f64), rgb: &[u8]) { + let file = File::create(path).unwrap(); + let mut enc = TiffEncoder::new(BufWriter::new(file)).unwrap(); + let mut img = enc.new_image::(cols, rows).unwrap(); + img.encoder() + .write_tag(Tag::Unknown(33550), &[1.0_f64, 1.0, 0.0][..]) + .unwrap(); + img.encoder() + .write_tag( + Tag::Unknown(33922), + &[0.0_f64, 0.0, 0.0, origin.0, origin.1, 0.0][..], + ) + .unwrap(); + img.encoder() + .write_tag(Tag::Unknown(34735), projected_dir(3035).as_slice()) + .unwrap(); + img.write_data(rgb).unwrap(); +} + +/// Write a single-band byte GeoTIFF used as the land-cover raster. +fn write_gray(path: &Path, cols: u32, rows: u32, origin: (f64, f64), g: &[u8]) { + let file = File::create(path).unwrap(); + let mut enc = TiffEncoder::new(BufWriter::new(file)).unwrap(); + let mut img = enc.new_image::(cols, rows).unwrap(); + img.encoder() + .write_tag(Tag::Unknown(33550), &[1.0_f64, 1.0, 0.0][..]) + .unwrap(); + img.encoder() + .write_tag( + Tag::Unknown(33922), + &[0.0_f64, 0.0, 0.0, origin.0, origin.1, 0.0][..], + ) + .unwrap(); + img.encoder() + .write_tag(Tag::Unknown(34735), projected_dir(3035).as_slice()) + .unwrap(); + img.write_data(g).unwrap(); +} + +#[test] +fn color_window_passes_through_rgb_and_packs_material_codes() { + let dir = tempfile::tempdir().unwrap(); + let (cols, rows) = (64u32, 64u32); + let origin = (4_400_000.0, 2_700_000.0); // EPSG:3035 top-left (easting, northing) + + // RGB: left half pure red, right half pure blue — distinct enough to assert + // the pass-through and the column split survive the windowed read. + let mut rgb = vec![0u8; (cols * rows * 3) as usize]; + for r in 0..rows { + for c in 0..cols { + let i = ((r * cols + c) * 3) as usize; + if c < cols / 2 { + rgb[i] = 200; // red + } else { + rgb[i + 2] = 200; // blue + } + } + } + let rgb_path = dir.path().join("ortho.tif"); + write_rgb(&rgb_path, cols, rows, origin, &rgb); + + // Land cover on the identical grid: left half water (class 5), right half + // high vegetation (class 1) → alpha must be WATER on the left, HIGH_VEG on + // the right, after the nearest-resample onto the ortho grid. + let mut lc = vec![0u8; (cols * rows) as usize]; + for r in 0..rows { + for c in 0..cols { + lc[(r * cols + c) as usize] = if c < cols / 2 { 5 } else { 1 }; + } + } + let lc_path = dir.path().join("landcover.tif"); + write_gray(&lc_path, cols, rows, origin, &lc); + + // Window centred on the tile centre, radius 24 px → ~48×48 clipped window. + let centre = (origin.0 + 32.0, origin.1 - 32.0); + let win = dem_io::extract_color_window(&rgb_path, Some(&lc_path), centre, 24.0, 0, Some(0)) + .expect("color window"); + + let (wc, wr) = (win.georef.cols, win.georef.rows); + assert!(wc >= 40 && wr >= 40, "expected ~48² window, got {wc}×{wr}"); + assert_eq!(win.rgba.len(), wc * wr * 4); + assert_eq!(win.georef.dx_meters, 1.0, "1 m/px carried into georef"); + assert!(win.georef.data.is_empty(), "georef is a placement stub"); + + // Sample a pixel in the left quarter (red + water) and the right quarter + // (blue + high veg). Window origin maps to a source column ≥ 8, so quarter + // offsets stay clear of the centre split. + let left = ((wr / 2) * wc + wc / 4) * 4; + let right = ((wr / 2) * wc + (3 * wc) / 4) * 4; + assert_eq!(win.rgba[left], 200, "left half red channel passes through"); + assert_eq!(win.rgba[left + 2], 0, "left half has no blue"); + assert_eq!(win.rgba[left + 3], MATERIAL_WATER, "left half water material"); + assert_eq!(win.rgba[right + 2], 200, "right half blue channel"); + assert_eq!(win.rgba[right], 0, "right half has no red"); + assert_eq!( + win.rgba[right + 3], + MATERIAL_HIGH_VEG, + "right half high-veg material" + ); +} + +#[test] +fn color_window_without_landcover_leaves_material_zero() { + let dir = tempfile::tempdir().unwrap(); + let (cols, rows) = (32u32, 32u32); + let origin = (4_400_000.0, 2_700_000.0); + let rgb = vec![150u8; (cols * rows * 3) as usize]; + let rgb_path = dir.path().join("ortho_only.tif"); + write_rgb(&rgb_path, cols, rows, origin, &rgb); + + let centre = (origin.0 + 16.0, origin.1 - 16.0); + let win = + dem_io::extract_color_window(&rgb_path, None, centre, 10.0, 0, None).expect("color window"); + assert!( + win.rgba.chunks_exact(4).all(|p| p[3] == MATERIAL_NONE), + "no land cover → every material code stays NONE" + ); + assert!( + win.rgba.chunks_exact(4).all(|p| p[0] == 150), + "flat ortho passes through unchanged" + ); +} + +#[test] +fn landcover_histogram_counts_class_values() { + let dir = tempfile::tempdir().unwrap(); + let (cols, rows) = (40u32, 40u32); + let origin = (4_400_000.0, 2_700_000.0); + // Three horizontal bands of classes 5 / 1 / 6. + let mut g = vec![0u8; (cols * rows) as usize]; + for r in 0..rows { + let v = if r < 13 { 5 } else if r < 26 { 1 } else { 6 }; + for c in 0..cols { + g[(r * cols + c) as usize] = v; + } + } + let lc_path = dir.path().join("lc_hist.tif"); + write_gray(&lc_path, cols, rows, origin, &g); + + // Full-tile window (radius covers the whole grid). + let centre = (origin.0 + 20.0, origin.1 - 20.0); + let hist = dem_io::landcover_histogram(&lc_path, centre, 100.0, 0).expect("histogram"); + assert!(hist[5] > 0 && hist[1] > 0 && hist[6] > 0, "all three classes present"); + assert_eq!( + hist.iter().sum::(), + (cols * rows) as u64, + "every pixel counted exactly once" + ); + assert_eq!(hist[2], 0, "absent class has zero count"); +} diff --git a/crates/dem_io/tests/synthetic.rs b/crates/dem_io/tests/synthetic.rs index 88c25bb..8a67bc7 100644 --- a/crates/dem_io/tests/synthetic.rs +++ b/crates/dem_io/tests/synthetic.rs @@ -245,6 +245,40 @@ fn declared_gdal_nodata_value_maps_to_nodata() { assert_eq!(win.data[6], data[6], "valid neighbours pass through"); } +#[test] +fn ifd_overview_levels_enumerates_a_multi_ifd_pyramid() { + // The overview cache is a multi-IFD GeoTIFF with no transparency-mask IFDs, + // so ifd_overview_levels must return one (ifd, scale) pair per level with + // contiguous indices and ascending scales. (The mask-skip branch needs a + // real BEV ortho and is covered by the gated local_big_tiles suite.) + let dir = tempfile::tempdir().unwrap(); + let src = dir.path().join("src.tif"); + let (cols, rows) = (96, 96); + write_min_geotiff( + &src, + cols, + rows, + (2.0, 2.0), + (4_400_000.0, 2_700_000.0), + &projected_dir(3035), + &ramp(cols, rows), + ); + let cache = ensure_overview_cache(&src, |_, _| {}) + .expect("cache build") + .expect("sub-5m source needs a cache"); + + let levels = dem_io::ifd_overview_levels(&cache).expect("walk IFDs"); + assert!(levels.len() >= 2, "multi-IFD pyramid expected, got {levels:?}"); + assert_eq!(levels[0].0, 0, "first level is IFD 0"); + for (i, &(ifd, scale)) in levels.iter().enumerate() { + assert_eq!(ifd, i, "no mask IFDs → contiguous indices"); + assert!(scale > 0.0, "positive pixel scale"); + } + for pair in levels.windows(2) { + assert!(pair[1].1 > pair[0].1, "scales increase down the pyramid"); + } +} + #[test] fn no_cache_built_for_coarse_single_ifd_source() { // A ≥20 m/px source can be served directly by base-tier select_ifd, so diff --git a/src/viewer/tiers.rs b/src/viewer/tiers.rs index d082f35..a0ff58b 100644 --- a/src/viewer/tiers.rs +++ b/src/viewer/tiers.rs @@ -1076,6 +1076,122 @@ mod tests { assert_eq!((out.rows, out.cols), (100, 100)); } + // cross_crs_world_origin / cross_crs_world_origin_and_extent + + /// Projected heightmap with an explicit CRS + origin (1 m/px square). + fn crs_proj(proj4: &str, ox: f64, oy: f64, cols: usize, rows: usize) -> Heightmap { + let mut hm = proj_hm(rows, cols); + hm.crs_proj4 = proj4.to_string(); + hm.crs_origin_x = ox; + hm.crs_origin_y = oy; + hm + } + + /// Geographic heightmap (deg/px) with an explicit lon/lat origin. + fn crs_geo(lon0: f64, lat0: f64, dscale: f64, cols: usize, rows: usize) -> Heightmap { + let mut hm = proj_hm(rows, cols); + hm.crs_proj4 = "+proj=longlat +datum=WGS84 +no_defs".to_string(); + hm.crs_origin_x = lon0; + hm.crs_origin_y = lat0; + hm.origin_lon = lon0; + hm.origin_lat = lat0; + hm.dx_deg = dscale; + hm.dy_deg = dscale; + hm.dx_meters = dscale * M_PER_DEG * lat0.to_radians().cos(); + hm.dy_meters = dscale * M_PER_DEG; + hm + } + + #[test] + fn cross_crs_same_crs_is_pure_offset() { + // Identical CRS → no projection: origin is the metre delta of the + // top-left corners (X right, Y down), extent is cols·dx / rows·dy, rot 0. + let utm = "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"; + let base = crs_proj(utm, 620_000.0, 5_240_000.0, 2000, 2000); + let mut hm = crs_proj(utm, 630_000.0, 5_235_000.0, 1000, 800); + hm.dx_meters = 2.0; + hm.dy_meters = 2.0; + + let (ox, oy) = cross_crs_world_origin(&hm, &base); + assert_eq!((ox, oy), (10_000.0, 5_000.0), "metre offset of TL corners"); + + let (ox2, oy2, ex, ey, rot) = cross_crs_world_origin_and_extent(&hm, &base); + assert_eq!((ox2, oy2), (10_000.0, 5_000.0)); + assert_eq!(ex, 2000.0, "extent_x = cols·dx"); + assert_eq!(ey, 1600.0, "extent_y = rows·dy"); + assert_eq!(rot, 0.0, "no meridian convergence within one CRS"); + } + + #[test] + fn cross_crs_projected_over_geographic_base_is_finite_and_sized() { + // 1 km UTM-32N window placed over a geographic (lon/lat) base — the + // geographic-base branch (cos-scaled metres + meridian rotation). + let base = crs_geo(11.4, 47.4, 0.000_27, 3000, 3000); + let hm = crs_proj( + "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs", + 630_000.0, + 5_235_000.0, + 1000, + 1000, + ); + let (ox, oy, ex, ey, rot) = cross_crs_world_origin_and_extent(&hm, &base); + for v in [ox, oy, ex, ey, rot] { + assert!(v.is_finite(), "all outputs finite, got {v}"); + } + // 1000 px @ 1 m projected into the base stays ~1 km within ±40 %. + assert!((600.0..1400.0).contains(&ex), "extent_x ≈ 1 km, got {ex}"); + assert!((600.0..1400.0).contains(&ey), "extent_y ≈ 1 km, got {ey}"); + // Meridian convergence between UTM-32N and lon/lat near 11.4°E is small. + assert!(rot.abs() < 0.3, "rotation small, got {rot} rad"); + } + + #[test] + fn cross_crs_projected_over_projected_base_is_finite_and_sized() { + // 1 km EPSG:3035 (LAEA) window over a UTM-32N base — the projected-base + // branch (from_wgs84 back into the base CRS). + let base = crs_proj( + "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs", + 620_000.0, + 5_240_000.0, + 2000, + 2000, + ); + let hm = crs_proj( + "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 \ + +ellps=GRS80 +towgs84=0,0,0 +units=m +no_defs", + 4_430_000.0, + 2_695_000.0, + 1000, + 1000, + ); + let (ox, oy, ex, ey, rot) = cross_crs_world_origin_and_extent(&hm, &base); + for v in [ox, oy, ex, ey, rot] { + assert!(v.is_finite(), "all outputs finite, got {v}"); + } + assert!((600.0..1400.0).contains(&ex), "extent_x ≈ 1 km, got {ex}"); + assert!((600.0..1400.0).contains(&ey), "extent_y ≈ 1 km, got {ey}"); + assert!(rot.abs() < 0.3, "rotation small, got {rot} rad"); + } + + #[test] + fn cross_crs_world_origin_geographic_base_branch() { + // Same geographic-base path for the origin-only helper. + let base = crs_geo(11.4, 47.4, 0.000_27, 3000, 3000); + let hm = crs_proj( + "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs", + 630_000.0, + 5_235_000.0, + 1000, + 1000, + ); + let (ox, oy) = cross_crs_world_origin(&hm, &base); + // The metre offset is finite and on the order of the inter-origin + // distance (tens of km), i.e. the projection actually ran rather than + // returning the (0, 0) parse-failure fallback. + assert!(ox.is_finite() && oy.is_finite()); + assert!(ox.abs() > 1.0 || oy.abs() > 1.0, "non-degenerate offset"); + } + #[test] fn cap_crops_oversized_axis_around_camera() { // 8200 px wide (over the 8192 limit), 4 px tall (under it).