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/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 cbe2fa3..c6b3d13 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(); @@ -202,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), @@ -234,6 +299,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 +413,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/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..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,15 +6,21 @@ 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, }; -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; 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/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/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}" + ); +} diff --git a/crates/dem_io/tests/synthetic.rs b/crates/dem_io/tests/synthetic.rs index c46c73f..8a67bc7 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,98 @@ 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 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/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..ccbf189 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,49 @@ 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); + } + + // 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 + // 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..30451bc 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; @@ -665,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); @@ -798,13 +877,66 @@ 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); + // 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); + 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); + } + } + + // 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 +950,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/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. 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 a96e7f0..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}; @@ -111,6 +112,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 { @@ -442,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. @@ -464,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::(); @@ -580,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 = @@ -616,6 +695,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(); @@ -750,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 { @@ -758,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; } @@ -903,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) { @@ -918,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() @@ -1061,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, )); @@ -1177,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, )); @@ -1239,6 +1376,7 @@ impl Viewer { close_tier_disabled: false, fine_disabled_by_oom: false, align_mode_viz: false, + ortho_enabled: true, } } } diff --git a/src/viewer/tiers.rs b/src/viewer/tiers.rs index 1a52d54..a0ff58b 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] @@ -766,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). @@ -856,4 +1282,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() + ); + } }