From 17406173f253a78c61312a187dfbed75a7e17bcc Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 7 Jul 2026 19:37:18 +0000 Subject: [PATCH 1/5] geo/garmin: Rust IMG decoder (FAT + TRE + RGN + LBL), byte-exact vs the prototype MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Ports scripts/garmin_proto.py — the validated Grand Canyon decoder — into a dep-free, pure-std `geo/src/garmin/` module. Four stages, one submodule each: - container.rs — the IMG container: XOR de-obfuscation, the 512-byte FAT from 0x600, and multi-part subfile reassembly (blocks concatenated in part*240+i order) → named `name.TYPE` byte slices. - tre.rs — the TRE: N/E/S/W int24 bbox, the LOD level pyramid, and the subdivision tree (16-byte records, 14-byte on the most-detailed level). - rgn.rs — the RGN: points / indexed-points / polylines / polygons, including the LSB-first delta bitstream with QMapShack CShiftReg continuation + sign semantics (the marker naive decoders drop, which turns fine levels into random-walk mush — the plan's documented level-4 bug). - lbl.rs — the LBL: 6-bit-packed and 8-bit-latin1 label text (road names, elevation-labelled contours). Public API in mod.rs: `Img::{parse,read,decode}`, `Feature`, `Kind`, `Decoded`, plus shared bounds-checked little-endian byte readers. Parity is proven, not asserted: an in-module test folds every decoded coordinate through FNV-1a-64 and matches the Python golden (0xadb368a3b063c74d) over all 120,174 features of the village tile 47505316, with exact per-kind (line 114616 / poly 3185 / point 2373) and per-level counts, a concrete first-feature check, and a decoded road name ("US Hwy 180"). Three more banked tiles decode to their exact prototype feature counts. This replaces raster colour-guessing with typed lookups for the geo bake pipeline (tasks #13-#16): every polygon/polyline carries a Garmin type code (building / water / forest / street / path / contour) plus contour elevations. geo is its own cargo workspace (isolated from the q2 workspace); clippy (-D warnings) and rustfmt clean on the new files. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_012jEwwaT5JZ5x8qWvcnaMYC --- .../2026-07-07-garmin-img-typed-rebuild.md | 13 +- geo/src/garmin/container.rs | 128 ++++++ geo/src/garmin/lbl.rs | 84 ++++ geo/src/garmin/mod.rs | 382 ++++++++++++++++++ geo/src/garmin/rgn.rs | 311 ++++++++++++++ geo/src/garmin/tre.rs | 125 ++++++ geo/src/lib.rs | 5 + 7 files changed, 1045 insertions(+), 3 deletions(-) create mode 100644 geo/src/garmin/container.rs create mode 100644 geo/src/garmin/lbl.rs create mode 100644 geo/src/garmin/mod.rs create mode 100644 geo/src/garmin/rgn.rs create mode 100644 geo/src/garmin/tre.rs diff --git a/claude-notes/plans/2026-07-07-garmin-img-typed-rebuild.md b/claude-notes/plans/2026-07-07-garmin-img-typed-rebuild.md index a4eaafffa..b713d543e 100644 --- a/claude-notes/plans/2026-07-07-garmin-img-typed-rebuild.md +++ b/claude-notes/plans/2026-07-07-garmin-img-typed-rebuild.md @@ -1,7 +1,10 @@ # Garmin IMG typed rebuild — ground-up geo pipeline (canyon / iceland / berlin) **Status:** format decoder VALIDATED (Python prototype, Grand Canyon renders -recognizably). Rust port + bake pipeline pending. +recognizably) **and PORTED to Rust** (`geo/src/garmin/`, byte-for-byte parity +with the prototype — FNV-1a-64 fold `0xadb368a3b063c74d` over all 120,174 +features of the village tile, plus exact counts on 3 more tiles). Bake pipeline +(tasks #13–#16) pending. ## Why (operator direction, 2026-07-07) @@ -121,8 +124,12 @@ or a second draw), since extruded footprints are genuinely non-deterministic. ## Remaining work (tasks #12-#16) -1. **Rust port** — `geo/src/garmin/` (mirror scripts/garmin_proto.py; the - prototype IS the executable spec, diff outputs against it). +1. **Rust port** — ✅ DONE. `geo/src/garmin/` — `container.rs` (FAT), `tre.rs` + (levels + subdivs), `rgn.rs` (points/lines/polys + the CShiftReg delta + bitstream), `lbl.rs` (6-bit + 8-bit labels), `mod.rs` (`Img` / `Feature` / + `Kind` / `Decoded` API). Dep-free pure `std`, clippy + fmt clean. Parity + asserted in-module against the prototype's golden (FNV-1a-64 over every + coord + per-kind/per-level counts + a decoded road name). 2. **Type codes → KIND** — typed lookup table per kind; sweet palette like `/ice` `Kind::color()`. One typed polygon = one building/house. 3. **Contours → heightfield** — grid interpolation (or TIN) from contour diff --git a/geo/src/garmin/container.rs b/geo/src/garmin/container.rs new file mode 100644 index 000000000..9c427d3bf --- /dev/null +++ b/geo/src/garmin/container.rs @@ -0,0 +1,128 @@ +//! The IMG container: XOR de-obfuscation, the FAT, and multi-part subfile +//! reassembly. A direct port of `read_img` in `scripts/garmin_proto.py`. +//! +//! Layout facts (hard-earned, see the plan): byte 0 = XOR key (0x00 on our +//! files); `DSKIMG` at 0x10; `blocksize = 1 << (b[0x61] + b[0x62])`; the FAT is +//! 512-byte entries from 0x600 (flag 0x01, 8-byte name + 3-byte type, size u32 +//! at 0x0C **valid in the part-0 entry**, part u16 at 0x10, 240 u16 block +//! pointers at 0x20). A subfile split across parts concatenates its blocks in +//! `part*240 + i` order. + +use std::collections::BTreeMap; + +use super::{u16le, u32le, GarminError}; + +/// One accumulating FAT entry: the declared size (from the part-0 record) plus +/// the block-index → physical-block map, keyed so a `BTreeMap` iterates them in +/// `part*240 + i` order. +struct Entry { + size: u32, + blocks: BTreeMap, +} + +/// Parse the container and return `name.TYPE → bytes` for every subfile. +/// +/// De-obfuscates in place when the XOR key is non-zero, verifies the `DSKIMG` +/// signature, then walks the FAT reassembling each (possibly multi-part) +/// subfile from its block list. +pub fn read_img(raw: &[u8]) -> Result>, GarminError> { + if raw.len() < 0x600 { + return Err(GarminError::Truncated("container header")); + } + + // Byte 0 is the XOR obfuscation key; de-obfuscate the whole image if set. + let xor = raw[0]; + let owned: Vec; + let b: &[u8] = if xor != 0 { + owned = raw.iter().map(|x| x ^ xor).collect(); + &owned + } else { + raw + }; + + if b.get(0x10..0x16) != Some(&b"DSKIMG"[..]) { + return Err(GarminError::NotImg); + } + + let blocksize: usize = 1usize << (u32::from(b[0x61]) + u32::from(b[0x62])); + if blocksize == 0 { + return Err(GarminError::Truncated("zero blocksize")); + } + + // FAT: 512-byte entries from 0x600. flag 0x00 ends the table; only flag + // 0x01 entries are live; the header pseudo-entry (blank / space-led name) + // is skipped. + let mut subs: BTreeMap<(String, String), Entry> = BTreeMap::new(); + let mut o = 0x600usize; + while o + 512 <= b.len() { + let flag = b[o]; + if flag != 0x01 { + if flag == 0x00 { + break; + } + o += 512; + continue; + } + let first = b[o + 1]; + let name = ascii_field(&b[o + 1..o + 9]); + let typ = ascii_field(&b[o + 9..o + 12]); + if name.is_empty() || first == b' ' || first == 0 { + o += 512; + continue; + } + + let size = u32le(b, o + 12); + let part = u32::from(u16le(b, o + 16)); + let ent = subs.entry((name, typ)).or_insert(Entry { + size: 0, + blocks: BTreeMap::new(), + }); + // The declared size lives in the part-0 record; keep it when non-zero. + if part == 0 && size != 0 { + ent.size = size; + } + // 240 block pointers per FAT entry; part N owns block indices [N*240..). + for i in 0..240u32 { + let blk = u16le(b, o + 0x20 + (i as usize) * 2); + if blk != 0xFFFF { + ent.blocks.insert(part * 240 + i, blk); + } + } + o += 512; + } + + // Reassemble each subfile by concatenating its blocks in index order, then + // truncating to the declared size (falling back to the block-span length). + let mut out = BTreeMap::new(); + for ((name, typ), ent) in subs { + let mut data = Vec::new(); + // BTreeMap::values() iterates in key order (part*240 + i), i.e. the + // correct block concatenation order. + for &blk in ent.blocks.values() { + let start = blk as usize * blocksize; + let end = (start + blocksize).min(b.len()); + if start < end { + data.extend_from_slice(&b[start..end]); + } + } + let size = if ent.size != 0 { + ent.size as usize + } else { + data.len() + }; + data.truncate(size); + out.insert(format!("{name}.{typ}"), data); + } + Ok(out) +} + +/// Decode a fixed-width ASCII name/type field (space-padded) and trim it, the +/// Rust equivalent of `b[...].decode('ascii','replace').strip()`. +fn ascii_field(b: &[u8]) -> String { + let text: String = b + .iter() + .map(|&c| if c.is_ascii() { c as char } else { '\u{FFFD}' }) + .collect(); + text.trim_matches(|c: char| c.is_ascii_whitespace()) + .to_string() +} diff --git a/geo/src/garmin/lbl.rs b/geo/src/garmin/lbl.rs new file mode 100644 index 000000000..0bc61f3b6 --- /dev/null +++ b/geo/src/garmin/lbl.rs @@ -0,0 +1,84 @@ +//! The LBL subfile: label text. A direct port of `parse_lbl` / `label_text` in +//! `scripts/garmin_proto.py`. +//! +//! Data section at `u32@0x15`; label offsets are multiplied by `1 << b[0x1D]`; +//! encoding at 0x1E (`6` = 6-bit packed, `9` = 8-bit latin1, nul-terminated). +//! Contour labels are elevations — feet in US topo maps, metres in OTM. + +use super::u32le; + +/// The 6-bit label alphabet: space, A–Z, five unused, 0–9, six unused. A code +/// above `0x2F` terminates the string. +const T6: &[u8] = b" ABCDEFGHIJKLMNOPQRSTUVWXYZ~~~~~0123456789~~~~~~"; + +/// A parsed LBL subfile: the data-section geometry plus a borrow of the bytes, +/// from which [`Lbl::text`] resolves a label offset to a string. +#[derive(Debug, Clone)] +pub struct Lbl<'a> { + pub off: usize, + pub len: usize, + pub mult: usize, + pub enc: u8, + pub data: &'a [u8], +} + +/// Parse an LBL subfile header. +#[must_use] +pub fn parse(lbl: &[u8]) -> Lbl<'_> { + let shift = u32::from(lbl.get(0x1D).copied().unwrap_or(0)); + Lbl { + off: u32le(lbl, 0x15) as usize, + len: u32le(lbl, 0x19) as usize, + mult: 1usize.checked_shl(shift).unwrap_or(1), + enc: lbl.get(0x1E).copied().unwrap_or(0), + data: lbl, + } +} + +impl Lbl<'_> { + /// Resolve a label offset (as carried on a feature) to its text. Offset `0` + /// means "no label" and yields the empty string. + #[must_use] + pub fn text(&self, lbloff: u32) -> String { + if lbloff == 0 { + return String::new(); + } + let start = self.off + lbloff as usize * self.mult; + let d = self.data; + + if self.enc == 9 { + // 8-bit latin1, nul-terminated. + match d.get(start..) { + Some(rest) => { + let end = rest.iter().position(|&c| c == 0).unwrap_or(rest.len()); + rest[..end].iter().map(|&c| c as char).collect() + } + None => String::new(), + } + } else { + // 6-bit packed: 3 bytes → 4 chars, high-first. + let mut out = String::new(); + let mut o = start; + while o + 3 <= d.len() { + let (b0, b1, b2) = (u32::from(d[o]), u32::from(d[o + 1]), u32::from(d[o + 2])); + for c in [ + b0 >> 2, + ((b0 & 3) << 4) | (b1 >> 4), + ((b1 & 0xF) << 2) | (b2 >> 6), + b2 & 0x3F, + ] { + if c > 0x2F { + return out; + } + out.push(if (c as usize) < T6.len() { + T6[c as usize] as char + } else { + '?' + }); + } + o += 3; + } + out + } + } +} diff --git a/geo/src/garmin/mod.rs b/geo/src/garmin/mod.rs new file mode 100644 index 000000000..82463fad8 --- /dev/null +++ b/geo/src/garmin/mod.rs @@ -0,0 +1,382 @@ +//! Garmin IMG map decoder — a dep-free, pure-`std` port of the validated +//! `scripts/garmin_proto.py` prototype (which renders the Grand Canyon tile +//! recognizably). Four stages, one per submodule: +//! +//! - [`container`] — the IMG container: XOR de-obfuscation, the FAT, and +//! multi-part subfile reassembly → named subfile byte slices. +//! - [`tre`] — the TRE subfile: bounding box, the LOD level pyramid, and the +//! subdivision tree. +//! - [`rgn`] — the RGN subfile: points / polylines / polygons, including the +//! LSB-first delta bitstream with QMapShack `CShiftReg` continuation +//! semantics. +//! - [`lbl`] — the LBL subfile: 6-bit-packed and 8-bit-latin1 label text. +//! +//! Garmin IMG gives the geo pipeline **typed** features — every polygon and +//! polyline carries a type code (building / water / forest / street / path / +//! contour) instead of raster colour-guessing — plus contour polylines with +//! elevation labels. Parity with the prototype is asserted byte-for-byte in the +//! tests via an FNV-1a-64 fold over every decoded coordinate. +//! +//! ```no_run +//! use geo_hhtl::garmin::{Img, Kind}; +//! +//! let img = Img::read("tile.img").expect("valid IMG"); +//! let decoded = img.decode().expect("TRE + RGN present"); +//! let lbl = img.lbl().map(geo_hhtl::garmin::lbl::parse); +//! for f in &decoded.features { +//! if f.kind == Kind::Line { +//! let name = lbl.as_ref().map(|l| l.text(f.lbl)).unwrap_or_default(); +//! let _ = (f.type_code, &f.coords, name); +//! } +//! } +//! ``` + +use std::collections::BTreeMap; +use std::fmt; + +pub mod container; +pub mod lbl; +pub mod rgn; +pub mod tre; + +pub use tre::{Bbox, Level, Subdiv, Tre}; + +// ── shared little-endian byte readers ─────────────────────────────────────── +// All return 0 on an out-of-bounds read; every hot-path caller guards its +// offset with an explicit length check first (mirroring the prototype's +// implicit Python slice bounds), so a 0 is only ever produced on genuinely +// malformed input, never on valid data. + +pub(crate) fn u16le(b: &[u8], o: usize) -> u16 { + match b.get(o..o + 2) { + Some(s) => u16::from_le_bytes([s[0], s[1]]), + None => 0, + } +} + +pub(crate) fn u32le(b: &[u8], o: usize) -> u32 { + match b.get(o..o + 4) { + Some(s) => u32::from_le_bytes([s[0], s[1], s[2], s[3]]), + None => 0, + } +} + +pub(crate) fn i16le(b: &[u8], o: usize) -> i16 { + match b.get(o..o + 2) { + Some(s) => i16::from_le_bytes([s[0], s[1]]), + None => 0, + } +} + +pub(crate) fn u24le(b: &[u8], o: usize) -> u32 { + match b.get(o..o + 3) { + Some(s) => u32::from(s[0]) | u32::from(s[1]) << 8 | u32::from(s[2]) << 16, + None => 0, + } +} + +pub(crate) fn i24le(b: &[u8], o: usize) -> i32 { + let v = u24le(b, o) as i32; + if v >= 1 << 23 { + v - (1 << 24) + } else { + v + } +} + +/// Convert signed int24 mapunits to degrees (`deg = mu * 360 / 2^24`). +#[must_use] +pub fn mu2deg(mu: i32) -> f64 { + f64::from(mu) * 360.0 / f64::from(1i32 << 24) +} + +// ── feature model ─────────────────────────────────────────────────────────── + +/// The object kind of a decoded feature. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Kind { + Point, + IPoint, + Line, + Poly, +} + +impl Kind { + /// A stable one-byte tag (used in the parity fold and for compact keys). + #[must_use] + pub fn tag(self) -> u8 { + match self { + Kind::Point => 0, + Kind::IPoint => 1, + Kind::Line => 2, + Kind::Poly => 3, + } + } + + #[must_use] + pub fn as_str(self) -> &'static str { + match self { + Kind::Point => "point", + Kind::IPoint => "ipoint", + Kind::Line => "line", + Kind::Poly => "poly", + } + } +} + +/// One decoded map feature: its kind, Garmin type code, LOD level, coordinate +/// ring in mapunits, and label offset (`0` = no label; resolve via [`lbl::Lbl`]). +#[derive(Debug, Clone)] +pub struct Feature { + pub kind: Kind, + pub type_code: u8, + pub level: u8, + /// `(lon, lat)` mapunits; a single point for [`Kind::Point`], a ring + /// otherwise. + pub coords: Vec<(i64, i64)>, + pub lbl: u32, +} + +/// A fully decoded IMG: the TRE tree plus every RGN feature. +#[derive(Debug, Clone)] +pub struct Decoded { + pub tre: Tre, + pub features: Vec, +} + +/// Error decoding a Garmin IMG. +#[derive(Debug)] +pub enum GarminError { + /// The `DSKIMG` container signature is missing. + NotImg, + /// A required subfile (`TRE` / `RGN`) is absent. + MissingSubfile(&'static str), + /// The image ended before a required structure could be read. + Truncated(&'static str), + /// The backing file could not be read. + Io(std::io::Error), +} + +impl fmt::Display for GarminError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + GarminError::NotImg => write!(f, "not a Garmin IMG (missing DSKIMG signature)"), + GarminError::MissingSubfile(s) => write!(f, "missing {s} subfile"), + GarminError::Truncated(s) => write!(f, "truncated: {s}"), + GarminError::Io(e) => write!(f, "io error: {e}"), + } + } +} + +impl std::error::Error for GarminError {} + +impl From for GarminError { + fn from(e: std::io::Error) -> Self { + GarminError::Io(e) + } +} + +/// A parsed IMG container: its subfiles keyed by `name.TYPE`. +#[derive(Debug, Clone)] +pub struct Img { + pub subfiles: BTreeMap>, +} + +impl Img { + /// Parse an in-memory IMG image. + pub fn parse(bytes: &[u8]) -> Result { + Ok(Img { + subfiles: container::read_img(bytes)?, + }) + } + + /// Read and parse an IMG file from disk. + pub fn read(path: impl AsRef) -> Result { + let bytes = std::fs::read(path)?; + Self::parse(&bytes) + } + + /// The first subfile whose name ends with `suffix` (e.g. `".TRE"`). + #[must_use] + pub fn subfile(&self, suffix: &str) -> Option<&[u8]> { + self.subfiles + .iter() + .find(|(k, _)| k.ends_with(suffix)) + .map(|(_, v)| v.as_slice()) + } + + #[must_use] + pub fn tre(&self) -> Option<&[u8]> { + self.subfile(".TRE") + } + + #[must_use] + pub fn rgn(&self) -> Option<&[u8]> { + self.subfile(".RGN") + } + + #[must_use] + pub fn lbl(&self) -> Option<&[u8]> { + self.subfile(".LBL") + } + + /// Decode the TRE tree and every RGN feature. + pub fn decode(&self) -> Result { + let tre = tre::parse(self.tre().ok_or(GarminError::MissingSubfile("TRE"))?)?; + let rgn_bytes = self.rgn().ok_or(GarminError::MissingSubfile("RGN"))?; + let features = rgn::parse(rgn_bytes, &tre, None); + Ok(Decoded { tre, features }) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Byte-for-byte parity fold over every decoded coordinate — identical to + /// the Python golden generator (kind tag, type byte, then each `(lon, lat)` + /// as little-endian `i64`). + fn fnv1a64_coords(feats: &[Feature]) -> u64 { + let mut h: u64 = 0xcbf2_9ce4_8422_2325; + let byte = |b: u8, h: &mut u64| { + *h ^= u64::from(b); + *h = h.wrapping_mul(0x0000_0100_0000_01b3); + }; + for f in feats { + byte(f.kind.tag(), &mut h); + byte(f.type_code, &mut h); + for &(lo, la) in &f.coords { + for b in lo.to_le_bytes() { + byte(b, &mut h); + } + for b in la.to_le_bytes() { + byte(b, &mut h); + } + } + } + h + } + + fn tile(name: &str) -> Img { + let path = format!( + "{}/../.claude/maps/garmin-grand-canyon/{name}.img", + env!("CARGO_MANIFEST_DIR") + ); + Img::read(&path).unwrap_or_else(|e| panic!("read {path}: {e}")) + } + + #[test] + fn base_bits_matches_imgformat() { + // base 0-9 → base+2; base >9 → 2*base-9. + assert_eq!(super::rgn::testonly_base_bits(0), 2); + assert_eq!(super::rgn::testonly_base_bits(9), 11); + assert_eq!(super::rgn::testonly_base_bits(10), 11); + assert_eq!(super::rgn::testonly_base_bits(15), 21); + } + + #[test] + fn village_tile_matches_prototype_golden() { + let img = tile("47505316"); + + // Container: exact subfile sizes. + assert_eq!(img.rgn().map(<[u8]>::len), Some(7_705_413)); + assert_eq!(img.tre().map(<[u8]>::len), Some(28_349)); + assert_eq!(img.lbl().map(<[u8]>::len), Some(42_669)); + + let dec = img.decode().expect("decode"); + + // TRE: bounding box (degrees), level pyramid, subdivision count. + assert!((mu2deg(dec.tre.bbox.north) - 36.3428).abs() < 1e-3); + assert!((mu2deg(dec.tre.bbox.east) - -111.7941).abs() < 1e-3); + assert!((mu2deg(dec.tre.bbox.south) - 35.3265).abs() < 1e-3); + assert!((mu2deg(dec.tre.bbox.west) - -112.8268).abs() < 1e-3); + + let expect_levels = [ + Level { + zoom: 4, + inherited: true, + bits: 17, + nsubdiv: 1, + }, + Level { + zoom: 3, + inherited: false, + bits: 18, + nsubdiv: 2, + }, + Level { + zoom: 2, + inherited: false, + bits: 20, + nsubdiv: 30, + }, + Level { + zoom: 1, + inherited: false, + bits: 22, + nsubdiv: 220, + }, + Level { + zoom: 0, + inherited: false, + bits: 23, + nsubdiv: 763, + }, + ]; + assert_eq!(dec.tre.levels, expect_levels); + assert_eq!(dec.tre.subdivs.len(), 1016); + + // RGN: feature counts by kind (bitstream-sensitive). + assert_eq!(dec.features.len(), 120_174); + let count = |k: Kind| dec.features.iter().filter(|f| f.kind == k).count(); + assert_eq!(count(Kind::Line), 114_616); + assert_eq!(count(Kind::Poly), 3_185); + assert_eq!(count(Kind::Point), 2_373); + + // Per-level feature counts (level 0's lone subdiv carries no data). + let per_level = |lv: u8| dec.features.iter().filter(|f| f.level == lv).count(); + assert_eq!(per_level(1), 371); + assert_eq!(per_level(2), 6_461); + assert_eq!(per_level(3), 29_829); + assert_eq!(per_level(4), 83_513); + + // First feature, concrete. + let f0 = &dec.features[0]; + assert_eq!(f0.kind, Kind::Line); + assert_eq!(f0.type_code, 2); + assert_eq!(f0.level, 1); + assert_eq!(f0.coords.len(), 3); + assert_eq!(f0.coords[0], (-5_223_552, 1_659_840)); + + // The killer assertion: every coordinate, byte-for-byte vs the Python + // golden. If the CShiftReg continuation/sign logic drifts at all, this + // diverges immediately. + assert_eq!(fnv1a64_coords(&dec.features), 0xadb3_68a3_b063_c74d); + + // LBL: the first labelled highway resolves to a real road name. + let lbl = lbl::parse(img.lbl().unwrap()); + assert_eq!(lbl.enc, 9); + assert_eq!(lbl.mult, 1); + assert_eq!(lbl.text(19_437), "US Hwy 180"); + } + + #[test] + fn other_canyon_tiles_decode_to_expected_counts() { + for (name, rgn_len, feats) in [ + ("47505310", 9_024_634usize, 116_845usize), + ("47505311", 3_871_943, 74_948), + ("47505317", 4_150_505, 63_149), + ] { + let img = tile(name); + assert_eq!( + img.rgn().map(<[u8]>::len), + Some(rgn_len), + "tile {name} RGN size" + ); + let dec = img + .decode() + .unwrap_or_else(|e| panic!("decode {name}: {e}")); + assert_eq!(dec.features.len(), feats, "tile {name} feature count"); + } + } +} diff --git a/geo/src/garmin/rgn.rs b/geo/src/garmin/rgn.rs new file mode 100644 index 000000000..43037d13e --- /dev/null +++ b/geo/src/garmin/rgn.rs @@ -0,0 +1,311 @@ +//! The RGN subfile: points, polylines, and polygons — including the LSB-first +//! delta **bitstream**. A direct port of `decode_bitstream` / `parse_rgn` in +//! `scripts/garmin_proto.py`, following QMapShack `CShiftReg` semantics exactly +//! (the continuation marker is the part naive decoders drop, turning fine +//! levels into random-walk mush). + +use super::tre::{Subdiv, Tre}; +use super::{i16le, u16le, u24le, u32le, Feature, Kind}; + +/// LSB-first bit reader over a byte slice. +struct BitReader<'a> { + d: &'a [u8], + pos: usize, +} + +impl<'a> BitReader<'a> { + fn new(d: &'a [u8]) -> Self { + Self { d, pos: 0 } + } + + fn remaining(&self) -> usize { + self.d.len() * 8 - self.pos + } + + /// Read `n` bits LSB-first. Caller must ensure `remaining() >= n`. + fn take(&mut self, n: usize) -> u32 { + let mut v = 0u32; + for i in 0..n { + let byte = self.d[self.pos >> 3]; + let bit = (byte >> (self.pos & 7)) & 1; + v |= u32::from(bit) << i; + self.pos += 1; + } + v + } +} + +/// Bit-width of a delta from its 4-bit base: `base+2` for `base <= 9`, else +/// `2*base - 9` (imgformat). +fn base_bits(base: u32) -> u32 { + if base <= 9 { + base + 2 + } else { + 2 * base - 9 + } +} + +#[cfg(test)] +pub(crate) fn testonly_base_bits(base: u32) -> u32 { + base_bits(base) +} + +/// How one axis of a polyline is encoded, read from the two leading info bits. +struct Axis { + /// `true` → per-delta two's complement (an extra width bit, with the + /// continuation marker). `false` → constant-sign unsigned magnitude. + signed: bool, + /// Fixed sign multiplier in constant-sign mode (unused when `signed`). + konst: i64, +} + +fn axis_info(br: &mut BitReader) -> Option { + if br.remaining() < 1 { + return None; + } + if br.take(1) != 0 { + // Constant-sign mode: one more bit picks the shared sign. + if br.remaining() < 1 { + return None; + } + let konst = if br.take(1) != 0 { -1 } else { 1 }; + Some(Axis { + signed: false, + konst, + }) + } else { + Some(Axis { + signed: true, + konst: 0, + }) + } +} + +/// Decode one delta on a signed axis, resolving the continuation marker: a raw +/// value equal to the sign bit alone (`1 << (n-1)`) accumulates `2^(n-1) - 1` +/// and reads again; the closing value ends the run. Returns `None` on underrun. +fn signed_delta(br: &mut BitReader, n: usize, sign: u32) -> Option { + let mut acc: i64 = 0; + loop { + if br.remaining() < n { + return None; + } + let tmp = br.take(n); + if tmp != sign { + return Some(if tmp < sign { + acc + i64::from(tmp) + } else { + (i64::from(tmp) - (i64::from(sign) << 1)) - acc + }); + } + acc += i64::from(tmp) - 1; + } +} + +/// Decode a polyline's delta bitstream into `(dlon, dlat)` pairs in level units. +fn decode_bitstream(data: &[u8], lon_base: u32, lat_base: u32, extra_bit: bool) -> Vec<(i64, i64)> { + let mut br = BitReader::new(data); + let mut out = Vec::new(); + + let Some(x_axis) = axis_info(&mut br) else { + return out; + }; + let Some(y_axis) = axis_info(&mut br) else { + return out; + }; + + let nx = (base_bits(lon_base) + u32::from(x_axis.signed)) as usize; + let ny = (base_bits(lat_base) + u32::from(y_axis.signed)) as usize; + let xsign = 1u32 << (nx - 1); + let ysign = 1u32 << (ny - 1); + + loop { + if extra_bit { + if br.remaining() < 1 { + break; + } + br.take(1); // routing-node flag, discarded + } + + let dlon = if x_axis.signed { + match signed_delta(&mut br, nx, xsign) { + Some(v) => v, + None => break, + } + } else { + if br.remaining() < nx { + break; + } + i64::from(br.take(nx)) * x_axis.konst + }; + + let dlat = if y_axis.signed { + match signed_delta(&mut br, ny, ysign) { + Some(v) => v, + None => break, + } + } else { + if br.remaining() < ny { + break; + } + i64::from(br.take(ny)) * y_axis.konst + }; + + out.push((dlon, dlat)); + } + out +} + +/// Decode every feature in an RGN subfile, driven by the TRE subdivision tree. +/// +/// `want_levels`, when `Some`, restricts decoding to those LOD levels; `None` +/// decodes all of them. +pub fn parse(rgn: &[u8], tre: &Tre, want_levels: Option<&[u8]>) -> Vec { + let data_off = u32le(rgn, 0x15) as usize; + let data_len = u32le(rgn, 0x19) as usize; + + // Subdivisions that own RGN data, in ascending rgn_off order; each one's + // data ends where the next live one begins (or at the section end). + let mut live: Vec = tre + .subdivs + .iter() + .filter(|s| s.objtypes != 0) + .copied() + .collect(); + live.sort_by_key(|s| s.rgn_off); + + let mut feats = Vec::new(); + for i in 0..live.len() { + let s = live[i]; + if let Some(wl) = want_levels { + if !wl.contains(&s.level) { + continue; + } + } + let rgn_end = if i + 1 < live.len() { + live[i + 1].rgn_off as usize + } else { + data_len + }; + let bits = tre.levels.get(s.level as usize).map_or(24, |l| l.bits); + let shift = 24u32.saturating_sub(u32::from(bits)); + + let block_start = data_off + s.rgn_off as usize; + let block_end = (data_off + rgn_end).min(rgn.len()); + if block_start >= block_end { + continue; + } + let block = &rgn[block_start..block_end]; + + // Object kinds present, in canonical order; if K kinds, (K-1) u16 + // pointers prefix the block and split it into per-kind spans. + let mut kinds: Vec = Vec::new(); + for (kind, mask) in [ + (Kind::Point, 0x10), + (Kind::IPoint, 0x20), + (Kind::Line, 0x40), + (Kind::Poly, 0x80), + ] { + if s.objtypes & mask != 0 { + kinds.push(kind); + } + } + if kinds.is_empty() { + continue; + } + let nptr = kinds.len() - 1; + let ptrs: Vec = (0..nptr).map(|j| u16le(block, 2 * j) as usize).collect(); + + let mut starts = Vec::with_capacity(kinds.len()); + starts.push(2 * nptr); + starts.extend(ptrs.iter().copied()); + let mut ends = ptrs; + ends.push(block.len()); + + for (ki, &kind) in kinds.iter().enumerate() { + let s0 = starts[ki].min(block.len()); + let e0 = ends[ki].min(block.len()); + if s0 >= e0 { + continue; + } + decode_segment(&block[s0..e0], kind, &s, shift, &mut feats); + } + } + feats +} + +/// Decode one per-kind span of a subdivision block. +fn decode_segment(seg: &[u8], kind: Kind, s: &Subdiv, shift: u32, feats: &mut Vec) { + let mut o = 0usize; + match kind { + // Indexed points decode identically to points and are reported as points + // (matching the prototype's `'kind': 'point'` for both). + Kind::Point | Kind::IPoint => { + while o + 9 <= seg.len() { + let t = seg[o]; + let lbl = u24le(seg, o + 1); + let has_sub = lbl & 0x80_0000 != 0; + let dlon = i64::from(i16le(seg, o + 4)); + let dlat = i64::from(i16le(seg, o + 6)); + o += if has_sub { 9 } else { 8 }; + let lon = i64::from(s.clon) + (dlon << shift); + let lat = i64::from(s.clat) + (dlat << shift); + feats.push(Feature { + kind: Kind::Point, + type_code: t, + level: s.level, + coords: vec![(lon, lat)], + lbl: lbl & 0x3F_FFFF, + }); + } + } + Kind::Line | Kind::Poly => { + let type_mask: u8 = if kind == Kind::Line { 0x3F } else { 0x7F }; + while o + 10 <= seg.len() { + let b0 = seg[o]; + let t = b0 & type_mask; + let two = b0 & 0x80 != 0; + let lblraw = u24le(seg, o + 1); + let extra = kind == Kind::Line && lblraw & 0x40_0000 != 0; + let lbl = lblraw & 0x3F_FFFF; + let dlon = i64::from(i16le(seg, o + 4)); + let dlat = i64::from(i16le(seg, o + 6)); + + // Bitstream byte length: two-byte when b0 bit7 is set. The info + // byte (lon/lat base nibbles) is separate from the length count. + let (blen, o2) = if two { + (u16le(seg, o + 8) as usize, o + 10) + } else { + (seg[o + 8] as usize, o + 9) + }; + let info = seg.get(o2).copied().unwrap_or(0); + let lon_base = u32::from(info & 0x0F); + let lat_base = u32::from(info >> 4); + let bs_start = o2 + 1; + let bs_end = bs_start.saturating_add(blen).min(seg.len()); + let bs: &[u8] = if bs_start <= seg.len() { + &seg[bs_start..bs_end] + } else { + &[] + }; + o = o2 + 1 + blen; + + let mut lon = i64::from(s.clon) + (dlon << shift); + let mut lat = i64::from(s.clat) + (dlat << shift); + let mut coords = vec![(lon, lat)]; + for (dlo, dla) in decode_bitstream(bs, lon_base, lat_base, extra) { + lon += dlo << shift; + lat += dla << shift; + coords.push((lon, lat)); + } + feats.push(Feature { + kind, + type_code: t, + level: s.level, + coords, + lbl, + }); + } + } + } +} diff --git a/geo/src/garmin/tre.rs b/geo/src/garmin/tre.rs new file mode 100644 index 000000000..5b9b725c4 --- /dev/null +++ b/geo/src/garmin/tre.rs @@ -0,0 +1,125 @@ +//! The TRE subfile: map bounding box, the LOD level pyramid, and the +//! subdivision tree. A direct port of `parse_tre` in `scripts/garmin_proto.py`. +//! +//! Layout: bbox N/E/S/W as signed int24 mapunits at 0x15/0x18/0x1B/0x1E; +//! the level table at `u32@0x21` (4 bytes each: `zoom | bit7 inherited`, bits, +//! `nsubdiv u16`); the subdivision table at `u32@0x29` — 16-byte records, with +//! the most-detailed (last) level using 14-byte records that omit `next`. + +use super::{i24le, u16le, u24le, u32le, GarminError}; + +/// One LOD level of the map (coarsest first). `bits` is the coordinate +/// resolution used to shift deltas back to full mapunits in the RGN. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct Level { + pub zoom: u8, + pub inherited: bool, + pub bits: u8, + pub nsubdiv: u16, +} + +/// One subdivision: a spatial cell owning a span of RGN data, tagged with which +/// object kinds it contains (`objtypes`: 0x10 point / 0x20 indexed-point / +/// 0x40 line / 0x80 polygon). +#[derive(Debug, Clone, Copy)] +pub struct Subdiv { + /// 1-based index in file order (subdivision references are 1-based). + pub n: u32, + /// Byte offset of this cell's data within the RGN data section. + pub rgn_off: u32, + pub objtypes: u8, + /// Cell centre longitude / latitude in mapunits. + pub clon: i32, + pub clat: i32, + pub width: u16, + pub height: u16, + pub terminate: bool, + pub next: u16, + /// Index into [`Tre::levels`]. + pub level: u8, +} + +/// Map bounding box in signed int24 mapunits (`deg = mu * 360 / 2^24`). +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct Bbox { + pub north: i32, + pub east: i32, + pub south: i32, + pub west: i32, +} + +/// The decoded TRE: bounding box, level pyramid, and every subdivision. +#[derive(Debug, Clone)] +pub struct Tre { + pub bbox: Bbox, + pub levels: Vec, + pub subdivs: Vec, +} + +/// Parse a TRE subfile. +pub fn parse(tre: &[u8]) -> Result { + if tre.len() < 0x31 { + return Err(GarminError::Truncated("TRE header")); + } + let bbox = Bbox { + north: i24le(tre, 0x15), + east: i24le(tre, 0x18), + south: i24le(tre, 0x1B), + west: i24le(tre, 0x1E), + }; + let lvl_off = u32le(tre, 0x21) as usize; + let lvl_size = u32le(tre, 0x25) as usize; + let sub_off = u32le(tre, 0x29) as usize; + let sub_size = u32le(tre, 0x2D) as usize; + + // Level table: 4-byte records over [lvl_off, lvl_off + lvl_size). + let mut levels = Vec::new(); + let lvl_end = lvl_off.saturating_add(lvl_size); + let mut o = lvl_off; + while o + 4 <= lvl_end && o + 4 <= tre.len() { + levels.push(Level { + zoom: tre[o] & 0x7F, + inherited: tre[o] & 0x80 != 0, + bits: tre[o + 1], + nsubdiv: u16le(tre, o + 2), + }); + o += 4; + } + + // Subdivision table: 16-byte records, last (most-detailed) level 14-byte. + let mut subdivs = Vec::new(); + let sub_end = sub_off.saturating_add(sub_size); + let mut o = sub_off; + let mut n = 1u32; + let nlev = levels.len(); + for (li, level) in levels.iter().enumerate() { + let last = li + 1 == nlev; + let rec = if last { 14 } else { 16 }; + for _ in 0..level.nsubdiv { + if o + rec > sub_end || o + rec > tre.len() { + break; + } + let wraw = u16le(tre, o + 10); + subdivs.push(Subdiv { + n, + rgn_off: u24le(tre, o) & 0x0FFF_FFFF, + objtypes: tre[o + 3], + clon: i24le(tre, o + 4), + clat: i24le(tre, o + 7), + width: wraw & 0x7FFF, + height: u16le(tre, o + 12) & 0x7FFF, + terminate: wraw & 0x8000 != 0, + next: if last { 0 } else { u16le(tre, o + 14) }, + level: li as u8, + }); + o += rec; + n += 1; + } + } + + Ok(Tre { + bbox, + levels, + subdivs, + }) +} diff --git a/geo/src/lib.rs b/geo/src/lib.rs index 454cdf2d9..0f6854915 100644 --- a/geo/src/lib.rs +++ b/geo/src/lib.rs @@ -19,6 +19,11 @@ pub mod extrude; pub mod hhtl; pub mod spm1; +/// Garmin IMG map decoder (dep-free, pure `std`): typed features + contours +/// from the banked `.img` sources — the ground-up replacement for raster +/// colour-guessing in the geo bake pipeline. +pub mod garmin; + /// Cesium-substrate seed (dep-free): the ratified TMS-quadkey addressing for the /// `cesium-osm-substrate-v1` path — the "version 2, for later" the `/helix` bake /// is a throwaway proof of. From 16943df7bb0fe815c24e9c2db3e981146c5cab99 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 7 Jul 2026 21:00:52 +0000 Subject: [PATCH 2/5] geo/kurvenlineal: continuous inter-family CurveRuler residue (fixes the needle bug) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The prior terrain bake seeded a fresh helix::CurveRuler per lattice cell (iceland_dem::ruler_phase → from_place(mix(floor(pos·detail)))). Adjacent cells drew uncorrelated phases, so the golden-spiral residue STEPPED discontinuously at every cell boundary. Harmless while it only tinted colour; pushed into geometry it IS the "needle field" — every cell its own isolated spike, never connected to its neighbour. That is the intra-family failure. This lifts the residue into a shared `kurvenlineal` module that samples the exact stride-4-over-17 CurveRuler value at the integer lattice CORNERS and smoothstep-interpolates between them (trilinear value-noise). The field is now C1-continuous ACROSS cell boundaries — inter-family: neighbouring cells share corner values, so relief flows smoothly to the next spike. It is NOT a gaussian blur — the field passes through the exact CurveRuler value at every lattice corner; only the cell interior is blended. Deterministic (phase is convention, not data — OGAR D-QUANTGATE). Test `inter_family_is_continuous_where_intra_family_jumps` proves the fix by reproducing the old per-cell construction alongside: the old residue's max adjacent step is > 0.3 (a seam jump); the new field's is < 0.05 (continuous). Plus determinism/range and pass-through-lattice-values tests. helix-gated (uses the real helix::CurveRuler, the same crate the sculpt tool + anatomy body use). clippy clean. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_012jEwwaT5JZ5x8qWvcnaMYC --- geo/src/kurvenlineal.rs | 159 ++++++++++++++++++++++++++++++++++++++++ geo/src/lib.rs | 5 ++ 2 files changed, 164 insertions(+) create mode 100644 geo/src/kurvenlineal.rs diff --git a/geo/src/kurvenlineal.rs b/geo/src/kurvenlineal.rs new file mode 100644 index 000000000..1185e86db --- /dev/null +++ b/geo/src/kurvenlineal.rs @@ -0,0 +1,159 @@ +//! The Kurvenlineal — the real `helix::CurveRuler` golden-spiral residue +//! (stride-4-over-17), as a **continuous inter-family** field. +//! +//! ## Why this exists — the intra-family needle bug +//! +//! The prior terrain bake seeded a *fresh* `CurveRuler` per lattice cell +//! (`from_place(mix(floor(pos·detail)))`, see the old `ruler_phase` in +//! `bin/iceland_dem.rs`). Adjacent cells drew *uncorrelated* phases, so the +//! residue **stepped discontinuously at every cell boundary**. That was +//! harmless while it only tinted colour — but pushed into *geometry* it is the +//! "needle field": every cell its own isolated spike, never connected to its +//! neighbour. That is the **intra-family** failure. +//! +//! This module samples the exact stride-4-over-17 residue at the integer +//! lattice **corners** and smoothstep-interpolates between them (value noise), +//! so the field is C1-continuous *across* cell boundaries — **inter-family**. +//! Neighbouring cells share corner values, so the relief flows smoothly to the +//! next spike instead of jumping. It is **not** a Gaussian blur: the field +//! passes through the exact CurveRuler value at every lattice corner; only the +//! cell interior is blended. Same `pos` + `detail` → same value on every bake +//! (phase is convention, not data — OGAR D-QUANTGATE). + +use helix::CurveRuler; + +// Odd 64-bit mixing constants (fractional golden / √2 / √3), one per axis, so a +// lattice corner decorrelates into a place anchor. +const MIX_X: u64 = 0x9E37_79B9_7F4A_7C15; +const MIX_Y: u64 = 0xC2B2_AE3D_27D4_EB4F; +const MIX_Z: u64 = 0x1656_67B1_9E37_79F9; + +/// Decorrelate an integer lattice corner into a `u64` place anchor. +fn mix(c: [i64; 3]) -> u64 { + (c[0] as u64) + .wrapping_mul(MIX_X) + .wrapping_add((c[1] as u64).wrapping_mul(MIX_Y)) + .wrapping_add((c[2] as u64).wrapping_mul(MIX_Z)) +} + +/// The exact stride-4-over-17 CurveRuler residue at one integer lattice corner, +/// mapped to bipolar `[-1, 1]` (index 8 → 0). Deterministic per corner. +fn corner_residue(corner: [i64; 3]) -> f32 { + let place = mix(corner); + // A stable per-corner arc index, taken from the high bits of the same hash + // so it is decorrelated from `start = place % 17`. + let k = ((place >> 40) % 17) as u32; + let idx = CurveRuler::from_place(place).index(k); + (f32::from(idx) / 16.0) * 2.0 - 1.0 +} + +/// Smoothstep (cubic Hermite): `smooth(0) = 0`, `smooth(1) = 1`, zero derivative +/// at both ends — so the interpolated field has no facet kinks at the lattice. +fn smooth(t: f32) -> f32 { + t * t * (3.0 - 2.0 * t) +} + +/// Continuous inter-family kurvenlineal relief in `[-1, 1]` at `pos`, sampled at +/// spatial frequency `detail` (cells per unit). Trilinear value-noise over the 8 +/// surrounding lattice corners, smoothstep-weighted per axis. +/// +/// Feed the return value into a vertex's height (bipolar relief) or into a +/// per-vertex colour brightness — it is C1-continuous across cell boundaries +/// either way, so it connects rather than spikes. +#[must_use] +pub fn relief(pos: [f32; 3], detail: f32) -> f32 { + let p = [pos[0] * detail, pos[1] * detail, pos[2] * detail]; + let base = [p[0].floor(), p[1].floor(), p[2].floor()]; + // `floor` rounds toward -inf, so the fraction is in [0, 1) even for negatives. + let w = [ + smooth(p[0] - base[0]), + smooth(p[1] - base[1]), + smooth(p[2] - base[2]), + ]; + let b = [base[0] as i64, base[1] as i64, base[2] as i64]; + + let mut acc = 0.0f32; + for (dz, wz) in [(0i64, 1.0 - w[2]), (1, w[2])] { + for (dy, wy) in [(0i64, 1.0 - w[1]), (1, w[1])] { + for (dx, wx) in [(0i64, 1.0 - w[0]), (1, w[0])] { + let v = corner_residue([b[0] + dx, b[1] + dy, b[2] + dz]); + acc += v * wx * wy * wz; + } + } + } + acc +} + +#[cfg(test)] +mod tests { + use super::*; + + /// The OLD per-cell construction (verbatim shape of `iceland_dem::ruler_phase`) + /// — reproduced here only to PROVE it is discontinuous where `relief` is not. + fn old_per_cell(pos: [f32; 3], detail: f32) -> f32 { + let cell = [ + (pos[0] * detail).floor() as i64, + (pos[1] * detail).floor() as i64, + (pos[2] * detail).floor() as i64, + ]; + let sub = [ + (pos[0] * detail * 3.0).floor() as i64, + (pos[1] * detail * 3.0).floor() as i64, + (pos[2] * detail * 3.0).floor() as i64, + ]; + let k = (mix(sub) % 17) as u32; + let idx = CurveRuler::from_place(mix(cell)).index(k); + (f32::from(idx) / 16.0) * 2.0 - 1.0 + } + + fn max_adjacent_delta(f: impl Fn([f32; 3], f32) -> f32) -> f32 { + // Sweep a line crossing ~10 cell boundaries at fine resolution. + let detail = 1.0; + let mut prev = f([0.0, 0.37, 0.11], detail); + let mut worst = 0.0f32; + let mut x = 0.0f32; + while x <= 10.0 { + let cur = f([x, 0.37, 0.11], detail); + worst = worst.max((cur - prev).abs()); + prev = cur; + x += 0.002; + } + worst + } + + #[test] + fn relief_is_deterministic_and_bounded() { + for i in 0..1000 { + let p = [i as f32 * 0.013, i as f32 * -0.021, i as f32 * 0.007]; + let v = relief(p, 2.5); + assert_eq!(v, relief(p, 2.5), "deterministic"); + assert!((-1.0..=1.0).contains(&v), "in range: {v}"); + } + } + + #[test] + fn relief_passes_through_the_lattice_values() { + // At an exact integer lattice point the fractions are 0, so `relief` + // reduces to that corner's exact CurveRuler residue (value-noise, not blur). + for c in [[0i64, 0, 0], [3, -2, 5], [-7, 11, -4]] { + let pos = [c[0] as f32, c[1] as f32, c[2] as f32]; + assert!((relief(pos, 1.0) - corner_residue(c)).abs() < 1e-6); + } + } + + #[test] + fn inter_family_is_continuous_where_intra_family_jumps() { + // THE fix, proven: the old per-cell residue jumps hard at cell seams + // (the needle field); the new inter-family field flows smoothly. + let old_jump = max_adjacent_delta(old_per_cell); + let new_jump = max_adjacent_delta(relief); + assert!( + old_jump > 0.3, + "old per-cell should step discontinuously, got {old_jump}" + ); + assert!( + new_jump < 0.05, + "new inter-family field should be continuous, got {new_jump}" + ); + } +} diff --git a/geo/src/lib.rs b/geo/src/lib.rs index 0f6854915..1cdf00496 100644 --- a/geo/src/lib.rs +++ b/geo/src/lib.rs @@ -37,6 +37,11 @@ pub mod osm_read; #[cfg(feature = "helix")] pub mod bso2; +/// The Kurvenlineal — the real `helix::CurveRuler` golden-spiral residue as a +/// CONTINUOUS INTER-FAMILY field (fixes the per-cell "needle" discontinuity). +#[cfg(feature = "helix")] +pub mod kurvenlineal; + pub use extrude::{extrude_into, Footprint, Layer}; pub use hhtl::{point_to_hhtl, tile_to_hhtl, Hhtl, GEO_DOMAIN, HHTL_DEPTH}; pub use spm1::{Mesh, Vertex}; From 8d2ab53d211e0b911817be93b2a32f481644a659 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 7 Jul 2026 21:14:17 +0000 Subject: [PATCH 3/5] =?UTF-8?q?geo/garmin:=20typed=20(kind,type=5Fcode)=20?= =?UTF-8?q?=E2=86=92=20GeoKind=20classifier=20+=20palette=20(heuristics=20?= =?UTF-8?q?stage)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The "heuristics" stage of the bake pipeline: replaces raster colour-guessing with a typed lookup. A Garmin polygon IS a building/water/forest and a polyline IS a contour/street/stream, by its type code — no inference, no colour thresholds. GeoKind { Building, Water, Stream, Woods, Park, Street, Path, Contour, Other } with a sweet /helix-look palette (water reads as water, vegetation as vegetation, no muddy browns), a stable PALETTE index order for the ver-8 wire palette block, is_area() (building/water/woods/park fill; rest are lines), and classify(kind, type_code) over the OpenTopoMap/US-topo code family. Feature gains geo_kind() and is_contour() (the contour lines are the elevation source for the terrain heightfield). Validated against the real village tile: the full-tile KIND histogram is asserted exactly (Contour 62574, Street 32154, Stream 15161, Path 3600, Other 3664, Water 2255, Park 766; Building/Woods 0 — wilderness tile), plus a partition check that every one of the 120174 features classifies to exactly one KIND. Dep-free, clippy clean. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_012jEwwaT5JZ5x8qWvcnaMYC --- geo/src/garmin/classify.rs | 169 +++++++++++++++++++++++++++++++++++++ geo/src/garmin/mod.rs | 2 + 2 files changed, 171 insertions(+) create mode 100644 geo/src/garmin/classify.rs diff --git a/geo/src/garmin/classify.rs b/geo/src/garmin/classify.rs new file mode 100644 index 000000000..c5c627532 --- /dev/null +++ b/geo/src/garmin/classify.rs @@ -0,0 +1,169 @@ +//! Typed classification: Garmin `(kind, type_code)` → a semantic [`GeoKind`] +//! with a clean palette. This is the "heuristics" stage — it replaces raster +//! colour-guessing with a **typed lookup**: a Garmin polygon *is* a building / +//! water / forest, and a polyline *is* a contour / street / stream, by its type +//! code. No inference, no colour thresholds. +//! +//! Type codes are for the OpenTopoMap / gpsfiledepot US-topo family the banked +//! `.img` sources come from (see the plan's type-code table). The classifier is +//! validated against the real village tile: the KIND histogram it produces is +//! asserted exactly in the tests. + +use super::{Feature, Kind}; + +/// A semantic surface class for a decoded Garmin feature. The discriminant +/// order is the palette index order (see [`GeoKind::PALETTE`]). +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum GeoKind { + Building, + Water, + Stream, + Woods, + Park, + Street, + Path, + Contour, + Other, +} + +impl GeoKind { + /// Canonical palette order — `GeoKind as u8` indexes this, and it is the + /// order a bake writes into the ver-8 wire's palette block. + pub const PALETTE: [GeoKind; 9] = [ + GeoKind::Building, + GeoKind::Water, + GeoKind::Stream, + GeoKind::Woods, + GeoKind::Park, + GeoKind::Street, + GeoKind::Path, + GeoKind::Contour, + GeoKind::Other, + ]; + + /// Palette index (`0..9`), stable across bakes. + #[must_use] + pub fn tag(self) -> u8 { + self as u8 + } + + /// A sweet, clean base colour (sRGB 0..255) — the `/helix` look: no muddy + /// browns, water reads as water, vegetation as vegetation. + #[must_use] + pub fn color(self) -> [u8; 3] { + match self { + GeoKind::Building => [200, 170, 130], // warm sandstone + GeoKind::Water => [64, 120, 180], // clean river blue + GeoKind::Stream => [96, 156, 206], // lighter running blue + GeoKind::Woods => [72, 116, 68], // forest green + GeoKind::Park => [132, 176, 96], // fresh grass green + GeoKind::Street => [96, 96, 102], // asphalt grey + GeoKind::Path => [190, 150, 104], // trail tan + GeoKind::Contour => [168, 138, 104], // topo brown + GeoKind::Other => [150, 150, 150], // neutral grey + } + } + + /// Whether this class fills/extrudes as an area (polygon) versus a line. + /// Buildings extrude; water / woods / parks are flat fills; everything else + /// is a line or point. + #[must_use] + pub fn is_area(self) -> bool { + matches!( + self, + GeoKind::Building | GeoKind::Water | GeoKind::Woods | GeoKind::Park + ) + } + + /// Classify a Garmin `(kind, type_code)` into a semantic class. + #[must_use] + pub fn classify(kind: Kind, type_code: u8) -> GeoKind { + match kind { + Kind::Poly => match type_code { + 0x13 | 0x60..=0x6f => GeoKind::Building, // building + urban variants + 0x3c..=0x4f => GeoKind::Water, // sea / lake / river-fill + 0x50 => GeoKind::Woods, // woods / forest + 0x14..=0x1a => GeoKind::Park, // park / reserve / grass + _ => GeoKind::Other, + }, + Kind::Line => match type_code { + 0x00..=0x07 => GeoKind::Street, // roads + 0x0a | 0x0b | 0x16 => GeoKind::Path, // trails / walking paths + 0x18 | 0x1f | 0x26 => GeoKind::Stream, // streams / rivers + 0x20..=0x25 => GeoKind::Contour, // land + depth contours + _ => GeoKind::Other, + }, + Kind::Point | Kind::IPoint => GeoKind::Other, // POIs — not a surface + } + } +} + +impl Feature { + /// The semantic surface class for this feature. + #[must_use] + pub fn geo_kind(&self) -> GeoKind { + GeoKind::classify(self.kind, self.type_code) + } + + /// `true` when this feature is a land/depth contour line — the elevation + /// source for the terrain heightfield (its [`Feature::lbl`] resolves to the + /// contour's elevation label). + #[must_use] + pub fn is_contour(&self) -> bool { + self.geo_kind() == GeoKind::Contour + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn palette_tag_round_trips() { + for (i, k) in GeoKind::PALETTE.iter().enumerate() { + assert_eq!(k.tag() as usize, i); + } + } + + #[test] + fn spot_check_type_mappings() { + // Representative codes from the real tile's histograms. + assert_eq!(GeoKind::classify(Kind::Line, 0x20), GeoKind::Contour); // minor land contour + assert_eq!(GeoKind::classify(Kind::Line, 0x22), GeoKind::Contour); // major land contour + assert_eq!(GeoKind::classify(Kind::Line, 0x06), GeoKind::Street); // road + assert_eq!(GeoKind::classify(Kind::Line, 0x26), GeoKind::Stream); // stream + assert_eq!(GeoKind::classify(Kind::Line, 0x0a), GeoKind::Path); // trail + assert_eq!(GeoKind::classify(Kind::Poly, 0x41), GeoKind::Water); // lake + assert_eq!(GeoKind::classify(Kind::Poly, 0x14), GeoKind::Park); // park + assert_eq!(GeoKind::classify(Kind::Poly, 0x13), GeoKind::Building); // building + assert_eq!(GeoKind::classify(Kind::Poly, 0x50), GeoKind::Woods); // woods + assert_eq!(GeoKind::classify(Kind::Point, 0x66), GeoKind::Other); // POI + } + + #[test] + fn village_tile_kind_histogram_matches_golden() { + // The whole tile, classified — asserted against the Python golden. + let path = format!( + "{}/../.claude/maps/garmin-grand-canyon/47505316.img", + env!("CARGO_MANIFEST_DIR") + ); + let img = super::super::Img::read(&path).expect("read tile"); + let dec = img.decode().expect("decode"); + + let count = |g: GeoKind| dec.features.iter().filter(|f| f.geo_kind() == g).count(); + assert_eq!(count(GeoKind::Contour), 62_574); + assert_eq!(count(GeoKind::Street), 32_154); + assert_eq!(count(GeoKind::Stream), 15_161); + assert_eq!(count(GeoKind::Path), 3_600); + assert_eq!(count(GeoKind::Other), 3_664); + assert_eq!(count(GeoKind::Water), 2_255); + assert_eq!(count(GeoKind::Park), 766); + assert_eq!(count(GeoKind::Building), 0); // wilderness tile — buildings are urban + assert_eq!(count(GeoKind::Woods), 0); + + // Every feature classifies to exactly one KIND (partition check). + let total: usize = GeoKind::PALETTE.iter().map(|&g| count(g)).sum(); + assert_eq!(total, dec.features.len()); + assert_eq!(total, 120_174); + } +} diff --git a/geo/src/garmin/mod.rs b/geo/src/garmin/mod.rs index 82463fad8..4a6312cb5 100644 --- a/geo/src/garmin/mod.rs +++ b/geo/src/garmin/mod.rs @@ -34,11 +34,13 @@ use std::collections::BTreeMap; use std::fmt; +pub mod classify; pub mod container; pub mod lbl; pub mod rgn; pub mod tre; +pub use classify::GeoKind; pub use tre::{Bbox, Level, Subdiv, Tre}; // ── shared little-endian byte readers ─────────────────────────────────────── From 1a2fe40a739a741c81f4dd5d38b08a1b47791b80 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 7 Jul 2026 21:19:32 +0000 Subject: [PATCH 4/5] =?UTF-8?q?geo/garmin:=20contours=20=E2=86=92=20connec?= =?UTF-8?q?ted=20heightfield=20(terrain=20stage,=20method=20B)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The "terrain" stage: labeled contour polylines → a connected grid heightfield. Method B — rasterize each contour polyline onto the W×H grid at its elevation, then fill the gaps by harmonic relaxation (SOR) with the contour cells held FIXED (Dirichlet BC). The contours stay razor-sharp; only the interior fills, so the surface is one connected sheet (grid neighbours = structural connectivity), never a field of isolated spikes. This is NOT gaussian smoothing — a blur would smear the contours; here they are immovable and relaxation only interpolates the unconstrained cells. The banked canyon data is ideal: 100% of contours are labeled (62573/62574), uniform 40-ft interval, distributed across TRE levels 2/3/4 (3690/8068/50816) — so LOD is the level pyramid for free: heightfield_for_level bakes one surface per level, the renderer swaps tier by camera distance. API: HeightField, grid_from_contours (domain-agnostic gridder), contours_for_level + heightfield_for_level (Garmin adapter: is_contour filter + LBL elevation, feet→m). Dep-free, clippy clean. Tests: (1) label parse feet/metres; (2) synthetic two-contour ramp interpolates linearly (mid ≈ 50, monotone); (3) real village-tile level-4 reconstruction — 50816 contours → a connected canyon: elevation stays in the real envelope (560..3158 m), >1500 m relief, and isolated needle spikes are <0.1% of cells (the discrete maximum principle the harmonic fill guarantees — the exact property the old per-cell needle bake violated). All finite, no NaN. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_012jEwwaT5JZ5x8qWvcnaMYC --- geo/src/garmin/mod.rs | 2 + geo/src/garmin/terrain.rs | 314 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 316 insertions(+) create mode 100644 geo/src/garmin/terrain.rs diff --git a/geo/src/garmin/mod.rs b/geo/src/garmin/mod.rs index 4a6312cb5..a70929dd1 100644 --- a/geo/src/garmin/mod.rs +++ b/geo/src/garmin/mod.rs @@ -38,9 +38,11 @@ pub mod classify; pub mod container; pub mod lbl; pub mod rgn; +pub mod terrain; pub mod tre; pub use classify::GeoKind; +pub use terrain::HeightField; pub use tre::{Bbox, Level, Subdiv, Tre}; // ── shared little-endian byte readers ─────────────────────────────────────── diff --git a/geo/src/garmin/terrain.rs b/geo/src/garmin/terrain.rs new file mode 100644 index 000000000..b463e7de9 --- /dev/null +++ b/geo/src/garmin/terrain.rs @@ -0,0 +1,314 @@ +//! The terrain stage: labeled contour polylines → a connected grid heightfield. +//! +//! Method **B** (grid, not a point TIN, not Gaussian splat): rasterize each +//! contour polyline onto the `W×H` grid at its elevation, then fill the gaps +//! between them by **harmonic relaxation with the contour cells held fixed** +//! (Dirichlet boundary conditions). The contours stay razor-sharp — only the +//! interior fills — so the surface is one connected sheet (grid neighbours = +//! structural connectivity), never a field of isolated spikes. This is *not* +//! Gaussian smoothing: a Gaussian blur would smear the contours; here they are +//! immovable and the relaxation only interpolates the unconstrained cells. +//! +//! LOD is the TRE level pyramid, for free: each level carries a generalized +//! contour set (this canyon tile: level 2 = 3690, level 3 = 8068, level 4 = +//! 50816), so [`heightfield_for_level`] bakes one surface per level and the +//! renderer swaps tier by camera distance. +//! +//! The sub-contour relief (the `kurvenlineal` residue) and the palette are +//! applied downstream by the bake; this module produces the base geometry. + +use super::tre::Bbox; +use super::{Decoded, Feature}; + +/// A row-major heightfield over a lon/lat bbox. `z[r*w + c]` is the surface +/// elevation in metres; row 0 is the north edge (equirectangular, matching the +/// OSM / Iceland bakes). +#[derive(Debug, Clone)] +pub struct HeightField { + pub w: usize, + pub h: usize, + pub z: Vec, +} + +impl HeightField { + #[must_use] + pub fn at(&self, col: usize, row: usize) -> f32 { + self.z[row * self.w + col] + } + + /// `(min, max)` elevation over the field. + #[must_use] + pub fn range(&self) -> (f32, f32) { + let mut lo = f32::INFINITY; + let mut hi = f32::NEG_INFINITY; + for &v in &self.z { + lo = lo.min(v); + hi = hi.max(v); + } + (lo, hi) + } +} + +/// Project a `(lon, lat)` mapunit onto fractional grid `(col, row)` over `bbox`. +/// Row 0 = north. Degenerate spans collapse to 0. +fn project(bbox: Bbox, lon: i64, lat: i64, w: usize, h: usize) -> (f32, f32) { + let ew = f64::from(bbox.east) - f64::from(bbox.west); + let ns = f64::from(bbox.north) - f64::from(bbox.south); + let fx = if ew != 0.0 { + (lon as f64 - f64::from(bbox.west)) / ew + } else { + 0.0 + }; + let fy = if ns != 0.0 { + (f64::from(bbox.north) - lat as f64) / ns + } else { + 0.0 + }; + ( + (fx * (w - 1) as f64) as f32, + (fy * (h - 1) as f64) as f32, + ) +} + +/// Stamp a straight segment `(c0,r0)→(c1,r1)` into the grid at `elev` (DDA walk, +/// one sample per max-axis step), marking each touched cell `known`. +#[allow(clippy::too_many_arguments)] +fn stamp_segment( + z: &mut [f32], + known: &mut [bool], + w: usize, + h: usize, + p0: (f32, f32), + p1: (f32, f32), + elev: f32, +) { + let steps = ((p1.0 - p0.0).abs().max((p1.1 - p0.1).abs())).ceil() as i64; + let n = steps.max(1); + for i in 0..=n { + let t = i as f32 / n as f32; + let c = (p0.0 + (p1.0 - p0.0) * t).round() as i64; + let r = (p0.1 + (p1.1 - p0.1) * t).round() as i64; + if c >= 0 && r >= 0 && (c as usize) < w && (r as usize) < h { + let idx = r as usize * w + c as usize; + z[idx] = elev; + known[idx] = true; + } + } +} + +/// Fill the unknown cells by successive over-relaxation (SOR) of Laplace's +/// equation, holding `known` (contour) cells fixed. Converges when the largest +/// update drops below `tol`, or after `max_iters`. +fn relax(z: &mut [f32], known: &[bool], w: usize, h: usize, max_iters: usize, tol: f32) { + const OMEGA: f32 = 1.8; // over-relaxation factor (1 < ω < 2) + for _ in 0..max_iters { + let mut max_delta = 0.0f32; + for r in 0..h { + for c in 0..w { + let idx = r * w + c; + if known[idx] { + continue; + } + // Average of the 4-neighbourhood, clamping at the grid edge. + let left = z[if c > 0 { idx - 1 } else { idx }]; + let right = z[if c + 1 < w { idx + 1 } else { idx }]; + let up = z[if r > 0 { idx - w } else { idx }]; + let down = z[if r + 1 < h { idx + w } else { idx }]; + let avg = 0.25 * (left + right + up + down); + let new = z[idx] + OMEGA * (avg - z[idx]); + max_delta = max_delta.max((new - z[idx]).abs()); + z[idx] = new; + } + } + if max_delta < tol { + break; + } + } +} + +/// Build a heightfield from contour polylines. Each entry is `(polyline in +/// mapunits, elevation in metres)`. Unconstrained cells are harmonically +/// interpolated between the contours. +#[must_use] +pub fn grid_from_contours( + contours: &[(&[(i64, i64)], f32)], + bbox: Bbox, + w: usize, + h: usize, +) -> HeightField { + let mut z = vec![0.0f32; w * h]; + let mut known = vec![false; w * h]; + + // Seed the unknown cells at the mean contour elevation so relaxation starts + // from a sensible level (fewer iterations, no dependence on the 0.0 default). + let mean = if contours.is_empty() { + 0.0 + } else { + contours.iter().map(|(_, e)| *e).sum::() / contours.len() as f32 + }; + z.fill(mean); + + for (line, elev) in contours { + if line.len() < 2 { + if let Some(&(lon, lat)) = line.first() { + let (c, r) = project(bbox, lon, lat, w, h); + stamp_segment(&mut z, &mut known, w, h, (c, r), (c, r), *elev); + } + continue; + } + for pair in line.windows(2) { + let (a, b) = (pair[0], pair[1]); + let p0 = project(bbox, a.0, a.1, w, h); + let p1 = project(bbox, b.0, b.1, w, h); + stamp_segment(&mut z, &mut known, w, h, p0, p1, *elev); + } + } + + // Relaxation tolerance: ~1 cm of the elevation span is plenty for terrain. + relax(&mut z, &known, w, h, 400, 0.01); + HeightField { w, h, z } +} + +/// Parse a contour label ("6000") to an elevation. `feet` converts US-topo feet +/// to metres (OTM labels are already metres). +fn label_elevation(text: &str, feet: bool) -> Option { + let digits: String = text + .chars() + .skip_while(|c| !c.is_ascii_digit() && *c != '-') + .take_while(|c| c.is_ascii_digit() || *c == '-') + .collect(); + let raw: f32 = digits.parse().ok()?; + Some(if feet { raw * 0.3048 } else { raw }) +} + +/// The labeled contour polylines of one LOD `level`, paired with their elevation +/// in metres. Borrows the decoded features. +#[must_use] +pub fn contours_for_level<'a>( + dec: &'a Decoded, + lbl: &super::lbl::Lbl<'_>, + level: u8, + feet: bool, +) -> Vec<(&'a [(i64, i64)], f32)> { + dec.features + .iter() + .filter(|f: &&Feature| f.is_contour() && f.level == level && f.lbl != 0) + .filter_map(|f| label_elevation(&lbl.text(f.lbl), feet).map(|e| (f.coords.as_slice(), e))) + .collect() +} + +/// Bake one LOD level's contours into a `W×H` heightfield over the TRE bbox. +#[must_use] +pub fn heightfield_for_level( + dec: &Decoded, + lbl: &super::lbl::Lbl<'_>, + level: u8, + w: usize, + h: usize, + feet: bool, +) -> HeightField { + let contours = contours_for_level(dec, lbl, level, feet); + grid_from_contours(&contours, dec.tre.bbox, w, h) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn label_elevation_feet_and_metres() { + assert!((label_elevation("6000", true).unwrap() - 1828.8).abs() < 0.1); + assert_eq!(label_elevation("500", false), Some(500.0)); + assert_eq!(label_elevation("", true), None); + } + + #[test] + fn synthetic_two_contours_ramp_linearly() { + // Two horizontal contours (bottom row elev 0, top row elev 100) over a + // unit bbox → the fill must ramp ~linearly with row. + let bbox = Bbox { + north: 1_000_000, + east: 1_000_000, + south: 0, + west: 0, + }; + let w = 32; + let h = 32; + // North edge (row 0) = lat=north → elev 100; south edge (row h-1) = elev 0. + let top: Vec<(i64, i64)> = (0..w).map(|c| (c as i64 * 32_258, 1_000_000)).collect(); + let bot: Vec<(i64, i64)> = (0..w).map(|c| (c as i64 * 32_258, 0)).collect(); + let contours: Vec<(&[(i64, i64)], f32)> = vec![(&top, 100.0), (&bot, 0.0)]; + let hf = grid_from_contours(&contours, bbox, w, h); + + // Row 0 ≈ 100, last row ≈ 0, middle ≈ 50, monotone decreasing down. + assert!((hf.at(w / 2, 0) - 100.0).abs() < 1.0); + assert!((hf.at(w / 2, h - 1) - 0.0).abs() < 1.0); + assert!((hf.at(w / 2, h / 2) - 50.0).abs() < 6.0); + for r in 1..h { + assert!( + hf.at(w / 2, r) <= hf.at(w / 2, r - 1) + 0.5, + "monotone down at row {r}" + ); + } + } + + fn village() -> (Decoded, Vec) { + let path = format!( + "{}/../.claude/maps/garmin-grand-canyon/47505316.img", + env!("CARGO_MANIFEST_DIR") + ); + let img = super::super::Img::read(&path).expect("read tile"); + let lbl_bytes = img.lbl().expect("lbl").to_vec(); + (img.decode().expect("decode"), lbl_bytes) + } + + #[test] + fn village_tile_level4_is_connected_terrain() { + let (dec, lbl_bytes) = village(); + let lbl = super::super::lbl::parse(&lbl_bytes); + + // Level 4 = the full-detail contour set (50816 contours, 1840..10360 ft). + let contours = contours_for_level(&dec, &lbl, 4, true); + assert!(contours.len() > 40_000, "level 4 contour count: {}", contours.len()); + + let hf = heightfield_for_level(&dec, &lbl, 4, 256, 256, true); + + // Elevation stays within the real canyon envelope (1840 ft = 560.8 m, + // 10360 ft = 3157.7 m), with a little slack for the ramp. + let (lo, hi) = hf.range(); + assert!(lo >= 540.0 && hi <= 3180.0, "range {lo}..{hi} m"); + assert!(hi - lo > 1500.0, "canyon has real relief: span {} m", hi - lo); + + // CONNECTED, not a needle field. A *cliff* (a large monotone step — the + // canyon has genuine ~900 m/cell walls) is fine; a *needle* is an + // isolated spike that pokes above (or below) ALL four neighbours by a + // wide margin. The harmonic fill obeys the discrete maximum principle, + // so a connected surface has essentially none; the old per-cell bake + // produced thousands. Assert isolated spikes are < 0.1% of interior cells. + let mut spikes = 0usize; + for r in 1..hf.h - 1 { + for c in 1..hf.w - 1 { + let v = hf.at(c, r); + let n = [ + hf.at(c - 1, r), + hf.at(c + 1, r), + hf.at(c, r - 1), + hf.at(c, r + 1), + ]; + let nmin = n.iter().copied().fold(f32::INFINITY, f32::min); + let nmax = n.iter().copied().fold(f32::NEG_INFINITY, f32::max); + if v > nmax + 300.0 || v < nmin - 300.0 { + spikes += 1; + } + } + } + let interior = (hf.h - 2) * (hf.w - 2); + assert!( + spikes * 1000 < interior, + "isolated needle spikes {spikes} of {interior} interior cells (>0.1%)" + ); + + // Fidelity: the field is finite everywhere (no NaN/inf leaked from relax). + assert!(hf.z.iter().all(|v| v.is_finite()), "heightfield has non-finite cells"); + } +} From 69ceda72b9f385bed276682ec3fbb792cdb3fbb6 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 7 Jul 2026 22:53:41 +0000 Subject: [PATCH 5/5] =?UTF-8?q?geo:=20garmin=5Fbake=20bin=20=E2=80=94=20Ga?= =?UTF-8?q?rmin=20IMG=20=E2=86=92=20BSO2=20ver-8=20terrain=20wire=20(the?= =?UTF-8?q?=20bake=20pipeline)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Ties the typed pipeline together into a bakeable scene: IMG → decode → contours → heightfield (terrain::heightfield_for_level) → landcover KIND overlay → ver-8 radix-grid wire. terrain.rs additions: - kind_grid(dec, bbox, w, h, overlay): rasterize selected features' GeoKind onto the grid (areas first, lines on top); bare terrain stays Other for the shader's hypsometric tint. LANDCOVER = [Water, Stream, Woods, Park] — roads/paths/ buildings are excluded so a wilderness scene reads as TERRAIN, not a street grid (the "cluttered lines" look, avoided by construction). project() made pub. garmin_bake bin (helix-gated, mirrors iceland_dem's projection/concept/encode): level (TRE LOD tier, default 4), aspect-correct grid, equirectangular metric projection about the tile centre, true-scale elevation, HHTL-keyed row concepts, GeoKind palette, encode_grid_bso2. The rich look (hypsometric LUT, inter-family kurvenlineal brightness, Gouraud, sunset, Ice/Ocean specular) is client-derived in the shader from the stored height+kind — the ver-8 contract. Verified end-to-end: baking the Grand Canyon village tile (47505316) → 844x1024 grid, elevation 561..3121 m (real river-to-rim envelope), 864256 verts / 1.72M tris / 2.64 MB ver-8 soa. Header inspected: BSO2 ver8, nC 1024, W 844, H 1024. New kind_grid test asserts streams overlay while roads/paths do not and terrain stays the majority. clippy + fmt clean. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_012jEwwaT5JZ5x8qWvcnaMYC --- geo/Cargo.toml | 6 ++ geo/src/bin/garmin_bake.rs | 212 +++++++++++++++++++++++++++++++++++++ geo/src/garmin/terrain.rs | 106 +++++++++++++++++-- 3 files changed, 315 insertions(+), 9 deletions(-) create mode 100644 geo/src/bin/garmin_bake.rs diff --git a/geo/Cargo.toml b/geo/Cargo.toml index 4d0bd83a4..42494d992 100644 --- a/geo/Cargo.toml +++ b/geo/Cargo.toml @@ -64,6 +64,12 @@ required-features = ["osm", "helix"] name = "iceland_dem" required-features = ["helix"] +# garmin_bake → BSO2 ver-8 wire from a Garmin IMG tile: contours → heightfield +# + typed KIND overlay. Needs `helix` (bso2 + ndarray + lance-graph-contract). +[[bin]] +name = "garmin_bake" +required-features = ["helix"] + [profile.release] opt-level = 3 codegen-units = 1 diff --git a/geo/src/bin/garmin_bake.rs b/geo/src/bin/garmin_bake.rs new file mode 100644 index 000000000..65166fa79 --- /dev/null +++ b/geo/src/bin/garmin_bake.rs @@ -0,0 +1,212 @@ +//! `garmin_bake` — bake a Garmin IMG tile into the BSO2 ver-8 radix-grid wire +//! the cockpit `/garmin/` route decodes, via the typed pipeline: +//! +//! ```text +//! IMG → decode → contours → heightfield (terrain::heightfield_for_level) +//! → KIND overlay grid (rivers/roads/forest on bare terrain) +//! → ver-8 wire (height F16 + kind u8 + palette) +//! ``` +//! +//! LOD is the TRE level pyramid — pass `--level N` (default 4 = full detail; the +//! canyon carries generalized sets at levels 2/3/4). The rich look (hypsometric +//! elevation tint, the inter-family kurvenlineal brightness, Gouraud normals, +//! sunset light, Ice/Ocean specular) is CLIENT-derived in the shader from the +//! stored height + kind — the ver-8 contract: store only what the address can't +//! reconstruct. +//! +//! Usage: `garmin_bake [--level N] [--dim WxH] [--metres]` + +use geo_hhtl::bso2::{encode_grid_bso2, MeshConcept, CLASSID_GEO_V3}; +use geo_hhtl::garmin::{mu2deg, terrain, Img}; +use geo_hhtl::hhtl::point_to_hhtl4; +use geo_hhtl::osm_read::M_PER_DEG; +use lance_graph_contract::canonical_node::{NodeGuid, TailVariant}; + +const KEY_ZOOM: u32 = 32; +const FAMILY_TERRAIN: u32 = 4; + +fn main() { + let mut args = std::env::args().skip(1); + let (Some(input), Some(output)) = (args.next(), args.next()) else { + eprintln!("usage: garmin_bake [--level N] [--dim WxH] [--metres]"); + std::process::exit(2); + }; + let mut level: u8 = 4; + let mut dim: Option<(usize, usize)> = None; + let mut feet = true; // US-topo default; --metres for OTM + let rest: Vec = args.collect(); + let mut it = rest.iter(); + while let Some(a) = it.next() { + match a.as_str() { + "--level" => level = it.next().and_then(|s| s.parse().ok()).unwrap_or(4), + "--metres" | "--meters" => feet = false, + "--dim" => { + if let Some(d) = it.next() { + if let Some((w, h)) = d.split_once('x') { + dim = w.parse().ok().zip(h.parse().ok()); + } + } + } + _ => {} + } + } + + let img = match Img::read(&input) { + Ok(i) => i, + Err(e) => { + eprintln!("garmin_bake: {e}"); + std::process::exit(1); + } + }; + let dec = match img.decode() { + Ok(d) => d, + Err(e) => { + eprintln!("garmin_bake: decode: {e}"); + std::process::exit(1); + } + }; + let lbl_bytes = img.lbl().map(<[u8]>::to_vec).unwrap_or_default(); + let lbl = geo_hhtl::garmin::lbl::parse(&lbl_bytes); + let bbox = dec.tre.bbox; + + // Grid dims: keep the tile's aspect at a ~1024-cell long axis unless overridden. + let (w, h) = dim.unwrap_or_else(|| { + let deg_w = (mu2deg(bbox.east) - mu2deg(bbox.west)).abs(); + let deg_h = (mu2deg(bbox.north) - mu2deg(bbox.south)).abs(); + let cos_lat = ((mu2deg(bbox.north) + mu2deg(bbox.south)) * 0.5) + .to_radians() + .cos(); + let aspect = (deg_w * cos_lat / deg_h.max(1e-9)).clamp(0.25, 4.0); + if aspect >= 1.0 { + (1024, (1024.0 / aspect).round() as usize) + } else { + ((1024.0 * aspect).round() as usize, 1024) + } + }); + eprintln!( + "tile bbox N{:.4} E{:.4} S{:.4} W{:.4} · level {level} · grid {w}x{h}", + mu2deg(bbox.north), + mu2deg(bbox.east), + mu2deg(bbox.south), + mu2deg(bbox.west), + ); + + // ── Terrain heightfield (metres) from the labeled contours of this LOD level. ── + let hf = terrain::heightfield_for_level(&dec, &lbl, level, w, h, feet); + let (elev_lo, elev_hi) = hf.range(); + eprintln!("heightfield {w}x{h} · elevation {elev_lo:.0}..{elev_hi:.0} m"); + + // ── KIND overlay: natural landcover (rivers / lakes / forest) on bare + // terrain — roads/paths are excluded so the mountain scene reads as + // terrain, not a street grid. ── + let kinds = terrain::kind_grid(&dec, bbox, w, h, &terrain::LANDCOVER); + + // ── Equirectangular metric projection about the tile centre (matches osm_read / + // iceland_dem so the cockpit decoder is shared). ── + let lon0 = (mu2deg(bbox.west) + mu2deg(bbox.east)) * 0.5; + let lat0 = (mu2deg(bbox.north) + mu2deg(bbox.south)) * 0.5; + let cos_lat0 = lat0.to_radians().cos(); + let lon_at = |c: usize| { + mu2deg(bbox.west) + + (mu2deg(bbox.east) - mu2deg(bbox.west)) * c as f64 / (w - 1).max(1) as f64 + }; + let lat_at = |r: usize| { + mu2deg(bbox.north) + - (mu2deg(bbox.north) - mu2deg(bbox.south)) * r as f64 / (h - 1).max(1) as f64 + }; + + let mut mx = vec![0.0f32; w * h]; + let mut mz = vec![0.0f32; w * h]; + for r in 0..h { + let z = ((lat_at(r) - lat0) * M_PER_DEG) as f32; + for c in 0..w { + let x = ((lon_at(c) - lon0) * cos_lat0 * M_PER_DEG) as f32; + mx[r * w + c] = x; + mz[r * w + c] = z; + } + } + // Normalize horizontal extent to [-1,1]; elevation is TRUE-SCALE (÷ same half) + // — the shader raises it by uExag, exactly like the Iceland terrain. + let (mut lox, mut hix, mut loz, mut hiz) = (f32::MAX, f32::MIN, f32::MAX, f32::MIN); + for i in 0..w * h { + lox = lox.min(mx[i]); + hix = hix.max(mx[i]); + loz = loz.min(mz[i]); + hiz = hiz.max(mz[i]); + } + let (cx, cz) = ((lox + hix) * 0.5, (loz + hiz) * 0.5); + let half = ((hix - lox).max(hiz - loz) * 0.5).max(1.0); + let inv = 1.0 / half; + + let mut pos = vec![[0.0f32; 3]; w * h]; + for (i, p) in pos.iter_mut().enumerate() { + *p = [(mx[i] - cx) * inv, hf.z[i] * inv, (mz[i] - cz) * inv]; + } + + // ── Concepts = grid rows (contiguous vertex strips), keyed by HHTL address. ── + let mid = w / 2; + let mut concepts = Vec::with_capacity(h); + for r in 0..h { + let v_start = (r * w) as u32; + let mut ctr = [0.0f32; 3]; + for p in &pos[r * w..(r + 1) * w] { + ctr[0] += p[0]; + ctr[1] += p[1]; + ctr[2] += p[2]; + } + let invn = 1.0 / w as f32; + let centroid = [ctr[0] * invn, ctr[1] * invn, ctr[2] * invn]; + let key4 = point_to_hhtl4(lon_at(mid), lat_at(r), KEY_ZOOM); + let key = NodeGuid::mint_for( + TailVariant::V3, + CLASSID_GEO_V3, + key4.heel, + key4.hip, + key4.twig, + key4.leaf, + FAMILY_TERRAIN, + r as u32, + ); + concepts.push(MeshConcept { + key, + layer: 4, + label: 0, + centroid, + v_start, + v_count: w as u32, + }); + } + + // ── ver-8 radix-grid encode: height F16 + kind u8 + GeoKind palette. ── + let palette: Vec<[u8; 3]> = geo_hhtl::garmin::GeoKind::PALETTE + .iter() + .map(|k| k.color()) + .collect(); + let ymax = pos.iter().map(|p| p[1]).fold(1e-9f32, f32::max); + let heights: Vec = pos.iter().map(|p| p[1] / ymax).collect(); + let x0 = pos[0][0]; + let dx = if w > 1 { pos[1][0] - pos[0][0] } else { 0.0 }; + let zrow: Vec = (0..h).map(|r| pos[r * w][2]).collect(); + let labels = br#"{"names":["garmin-terrain"]}"#; + + let (soa, blocks) = encode_grid_bso2( + w as u32, h as u32, x0, dx, ymax, &zrow, &heights, &kinds, &palette, &concepts, labels, + ); + eprintln!( + "BSO2 ver8: {} concepts · {} verts · {} tris · {} B soa · {} B blocks", + concepts.len(), + w * h, + (w - 1) * (h - 1) * 2, + soa.len(), + blocks.len() + ); + + let blocks_path = format!("{}.blocks", output.strip_suffix(".soa").unwrap_or(&output)); + if let Err(e) = + std::fs::write(&output, &soa).and_then(|()| std::fs::write(&blocks_path, &blocks)) + { + eprintln!("garmin_bake: write: {e}"); + std::process::exit(1); + } + eprintln!("wrote {output} + {blocks_path}"); +} diff --git a/geo/src/garmin/terrain.rs b/geo/src/garmin/terrain.rs index b463e7de9..ebbecda0d 100644 --- a/geo/src/garmin/terrain.rs +++ b/geo/src/garmin/terrain.rs @@ -18,7 +18,7 @@ //! applied downstream by the bake; this module produces the base geometry. use super::tre::Bbox; -use super::{Decoded, Feature}; +use super::{Decoded, Feature, GeoKind, Kind}; /// A row-major heightfield over a lon/lat bbox. `z[r*w + c]` is the surface /// elevation in metres; row 0 is the north edge (equirectangular, matching the @@ -51,7 +51,7 @@ impl HeightField { /// Project a `(lon, lat)` mapunit onto fractional grid `(col, row)` over `bbox`. /// Row 0 = north. Degenerate spans collapse to 0. -fn project(bbox: Bbox, lon: i64, lat: i64, w: usize, h: usize) -> (f32, f32) { +pub fn project(bbox: Bbox, lon: i64, lat: i64, w: usize, h: usize) -> (f32, f32) { let ew = f64::from(bbox.east) - f64::from(bbox.west); let ns = f64::from(bbox.north) - f64::from(bbox.south); let fx = if ew != 0.0 { @@ -64,10 +64,7 @@ fn project(bbox: Bbox, lon: i64, lat: i64, w: usize, h: usize) -> (f32, f32) { } else { 0.0 }; - ( - (fx * (w - 1) as f64) as f32, - (fy * (h - 1) as f64) as f32, - ) + ((fx * (w - 1) as f64) as f32, (fy * (h - 1) as f64) as f32) } /// Stamp a straight segment `(c0,r0)→(c1,r1)` into the grid at `elev` (DDA walk, @@ -169,6 +166,64 @@ pub fn grid_from_contours( HeightField { w, h, z } } +/// Stamp a straight segment into a KIND grid with `tag` (DDA walk). +fn stamp_kind_line(kinds: &mut [u8], w: usize, h: usize, p0: (f32, f32), p1: (f32, f32), tag: u8) { + let steps = ((p1.0 - p0.0).abs().max((p1.1 - p0.1).abs())).ceil() as i64; + let n = steps.max(1); + for i in 0..=n { + let t = i as f32 / n as f32; + let c = (p0.0 + (p1.0 - p0.0) * t).round() as i64; + let r = (p0.1 + (p1.1 - p0.1) * t).round() as i64; + if c >= 0 && r >= 0 && (c as usize) < w && (r as usize) < h { + kinds[r as usize * w + c as usize] = tag; + } + } +} + +/// The natural-landcover overlay for a wilderness terrain scene: blue water + +/// rivers, green forest + park. Roads/paths/buildings are deliberately excluded +/// — on a mountain scene they clutter; the hypsometric terrain tint carries it. +pub const LANDCOVER: [GeoKind; 4] = [ + GeoKind::Water, + GeoKind::Stream, + GeoKind::Woods, + GeoKind::Park, +]; + +/// Rasterize selected map features' semantic KIND onto a `W×H` grid, so the +/// `overlay` classes (e.g. [`LANDCOVER`]) drape onto the terrain heightfield. +/// Bare terrain stays [`GeoKind::Other`] (the shader tints it hypsometrically by +/// height). Draw order: areas (polygons) first, then lines on top, so a river on +/// a lake reads as river. Contours are the height source and are not stamped. +#[must_use] +pub fn kind_grid(dec: &Decoded, bbox: Bbox, w: usize, h: usize, overlay: &[GeoKind]) -> Vec { + let mut kinds = vec![GeoKind::Other.tag(); w * h]; + let mut ordered: Vec<&Feature> = dec + .features + .iter() + .filter(|f| !f.is_contour() && overlay.contains(&f.geo_kind())) + .collect(); + ordered.sort_by_key(|f| match f.kind { + Kind::Poly => 0u8, + Kind::Line => 1, + _ => 2, + }); + for f in ordered { + let tag = f.geo_kind().tag(); + if f.coords.len() == 1 { + let p = project(bbox, f.coords[0].0, f.coords[0].1, w, h); + stamp_kind_line(&mut kinds, w, h, p, p, tag); + continue; + } + for pair in f.coords.windows(2) { + let p0 = project(bbox, pair[0].0, pair[0].1, w, h); + let p1 = project(bbox, pair[1].0, pair[1].1, w, h); + stamp_kind_line(&mut kinds, w, h, p0, p1, tag); + } + } + kinds +} + /// Parse a contour label ("6000") to an elevation. `feet` converts US-topo feet /// to metres (OTM labels are already metres). fn label_elevation(text: &str, feet: bool) -> Option { @@ -269,7 +324,11 @@ mod tests { // Level 4 = the full-detail contour set (50816 contours, 1840..10360 ft). let contours = contours_for_level(&dec, &lbl, 4, true); - assert!(contours.len() > 40_000, "level 4 contour count: {}", contours.len()); + assert!( + contours.len() > 40_000, + "level 4 contour count: {}", + contours.len() + ); let hf = heightfield_for_level(&dec, &lbl, 4, 256, 256, true); @@ -277,7 +336,11 @@ mod tests { // 10360 ft = 3157.7 m), with a little slack for the ramp. let (lo, hi) = hf.range(); assert!(lo >= 540.0 && hi <= 3180.0, "range {lo}..{hi} m"); - assert!(hi - lo > 1500.0, "canyon has real relief: span {} m", hi - lo); + assert!( + hi - lo > 1500.0, + "canyon has real relief: span {} m", + hi - lo + ); // CONNECTED, not a needle field. A *cliff* (a large monotone step — the // canyon has genuine ~900 m/cell walls) is fine; a *needle* is an @@ -309,6 +372,31 @@ mod tests { ); // Fidelity: the field is finite everywhere (no NaN/inf leaked from relax). - assert!(hf.z.iter().all(|v| v.is_finite()), "heightfield has non-finite cells"); + assert!( + hf.z.iter().all(|v| v.is_finite()), + "heightfield has non-finite cells" + ); + } + + #[test] + fn village_tile_kind_grid_overlays_landcover() { + let (dec, _lbl) = village(); + let kinds = kind_grid(&dec, dec.tre.bbox, 512, 512, &LANDCOVER); + assert_eq!(kinds.len(), 512 * 512); + + let n = |g: GeoKind| kinds.iter().filter(|&&k| k == g.tag()).count(); + // Rivers (streams) drape onto the terrain… + assert!(n(GeoKind::Stream) > 0, "streams stamped"); + // …roads/paths are NOT overlaid (landcover only — no urban clutter). + assert_eq!(n(GeoKind::Street), 0, "roads excluded from terrain overlay"); + assert_eq!(n(GeoKind::Path), 0, "paths excluded from terrain overlay"); + // The terrain surface still dominates (hypsometric tint carries it). + assert!( + n(GeoKind::Other) > kinds.len() / 2, + "most of the surface is bare terrain: {} of {}", + n(GeoKind::Other), + kinds.len() + ); + assert!(kinds.iter().all(|&k| (k as usize) < GeoKind::PALETTE.len())); } }