diff --git a/.claude/maps/SOURCES.md b/.claude/maps/SOURCES.md new file mode 100644 index 000000000..c772ce458 --- /dev/null +++ b/.claude/maps/SOURCES.md @@ -0,0 +1,30 @@ +# .claude/maps — source provenance + +Git-tracked map artifacts (each file < 100 MB — GitHub blob limit). The +Dockerfile copies/reassembles what the deploy serves; the bakers read the +rest at bake time. + +## Garmin IMG sources (the ground-up typed rebuild, 2026-07-07) + +Garmin IMG carries **typed vector data** (every polygon/polyline has a type +code — building / water / forest / street / walking path / contour line) plus +multi-level generalization (TRE levels ≈ LOD pyramid) and, in newer builds, +per-tile `.DEM` elevation subfiles. Typed codes replace color heuristics: +one typed polygon = one building; contours → interpolated heightfield. + +| file | source | contents | +|---|---|---| +| `otm-iceland.zip` | | `otm-iceland.img` (71.9 MB): 7 map tiles `534008xx`, each with TRE/RGN/LBL/NET/NOD **/DEM**, + `OPENTOPO.TYP` style. mkgmap-built, unencrypted (xor 0x00), blocksize 2048. | +| `otm-iceland-contours.zip` | | `otm-iceland-contours.img` (26.8 MB): 16 contour tiles `533500xx` (TRE/RGN/LBL) + `CONTOURS.TYP`. Elevation contour polylines. | +| `garmin-grand-canyon/475053{10,11,16,17}.img` | (Arizona Topo, NSIS installer, 186 MB — NOT committed; tiles extracted via `7z x`) | The 4 tiles covering the Grand Canyon. `47505316` = Grand Canyon Village (N 36.343 S 35.327 W −112.827 E −111.794); `47505317` = North Rim strip; `47505310`/`47505311` = western canyon. 2016 build: no DEM subfile — heights are **contour polylines** in RGN. | + +Other free Garmin map source (not yet used): . + +## Existing artifacts (pre-Garmin arc) + +| file | role | +|---|---| +| `iceland_dem_16m.z01` + `iceland_dem_16m.zip` | 2-part split zip of `iceland_dem.helix.soa.gz` (151 MB, BSO2 ver-7, 16.5M verts). Reassemble: `zip -s 0 iceland_dem_16m.zip --out full.zip && unzip full.zip`. The Dockerfile does this into `cockpit/dist/`. | +| `berlin-latest.osm.pbf` | Raw OSM PBF for the Berlin bakes (94 MB). | +| `iceland-dem-cop30-100m-isn93.tif` / `iceland-hillshade-cop30-100m-isn93.tif` | Copernicus 30 m DEM + hillshade (ISN93). | +| `iceland.helix.soa.gz` | OLD sparse Iceland scatter (unreferenced by the deploy). | diff --git a/cockpit/public/iceland_dem.helix.soa.gz b/.claude/maps/garmin-grand-canyon/47505310.img similarity index 52% rename from cockpit/public/iceland_dem.helix.soa.gz rename to .claude/maps/garmin-grand-canyon/47505310.img index a9848139e..da15e8cdd 100644 Binary files a/cockpit/public/iceland_dem.helix.soa.gz and b/.claude/maps/garmin-grand-canyon/47505310.img differ diff --git a/.claude/maps/garmin-grand-canyon/47505311.img b/.claude/maps/garmin-grand-canyon/47505311.img new file mode 100644 index 000000000..6f915f608 Binary files /dev/null and b/.claude/maps/garmin-grand-canyon/47505311.img differ diff --git a/.claude/maps/garmin-grand-canyon/47505316.img b/.claude/maps/garmin-grand-canyon/47505316.img new file mode 100644 index 000000000..3bc111e94 Binary files /dev/null and b/.claude/maps/garmin-grand-canyon/47505316.img differ diff --git a/.claude/maps/garmin-grand-canyon/47505317.img b/.claude/maps/garmin-grand-canyon/47505317.img new file mode 100644 index 000000000..af0539c7f Binary files /dev/null and b/.claude/maps/garmin-grand-canyon/47505317.img differ diff --git a/.claude/maps/iceland_dem_16m.z01 b/.claude/maps/iceland_dem_16m.z01 new file mode 100644 index 000000000..1c656a62e Binary files /dev/null and b/.claude/maps/iceland_dem_16m.z01 differ diff --git a/.claude/maps/iceland_dem_16m.zip b/.claude/maps/iceland_dem_16m.zip new file mode 100644 index 000000000..a5316b331 Binary files /dev/null and b/.claude/maps/iceland_dem_16m.zip differ diff --git a/.claude/maps/otm-iceland-contours.zip b/.claude/maps/otm-iceland-contours.zip new file mode 100644 index 000000000..4e0008ef6 Binary files /dev/null and b/.claude/maps/otm-iceland-contours.zip differ diff --git a/.claude/maps/otm-iceland.zip b/.claude/maps/otm-iceland.zip new file mode 100644 index 000000000..243ff1d91 Binary files /dev/null and b/.claude/maps/otm-iceland.zip differ diff --git a/.gitignore b/.gitignore index 6e5e9d97f..2c6809e7c 100644 --- a/.gitignore +++ b/.gitignore @@ -60,3 +60,9 @@ __pycache__/ cockpit/public/body.soa cockpit/public/body.soa.gz scratch-fma/ +cockpit/public/*.blocks + +# Iceland DEM bake is 151MB (full z9 + imagery) — over GitHub 100MB; ship via Release asset. +cockpit/public/iceland_dem.helix.soa.gz +# in-repo scratch (demgrids, screenshots, harness) — belongs out of git. +/scratchpad/ diff --git a/Dockerfile b/Dockerfile index c91e147e9..8e142fcd1 100644 --- a/Dockerfile +++ b/Dockerfile @@ -33,7 +33,7 @@ RUN apt-get update && apt-get install -y \ git curl build-essential cmake clang \ libssl-dev pkg-config python3 \ protobuf-compiler libprotobuf-dev \ - ca-certificates lld \ + ca-certificates lld zip unzip \ && rm -rf /var/lib/apt/lists/* # Rust 1.94.0 @@ -69,13 +69,21 @@ RUN curl -fSL https://github.com/AdaWorldAPI/q2/releases/download/fma-body-soa-v -o /build/q2/cockpit/dist/body.20260629c.v6helix.soa.gz \ && ls -lh /build/q2/cockpit/dist/body.20260629c.v6helix.soa.gz -# The Iceland OSM helix bake (BSO2 ver 6) ships IN-REPO under .claude/maps/, not from -# the release — this session type can't upload release assets, so the artifact is -# version-controlled instead. BodyHelix fetches /iceland.helix.soa.gz same-origin -# (manifest iceland_latest); include_dir! embeds cockpit/dist/, so copy it there. The -# build context already carries .claude/maps via `COPY . /build/q2` above. -RUN cp /build/q2/.claude/maps/iceland.helix.soa.gz /build/q2/cockpit/dist/iceland.helix.soa.gz \ - && ls -lh /build/q2/cockpit/dist/iceland.helix.soa.gz +# The Iceland DEM terrain bake (BSO2 ver 7 = F16 pos + Signed360 normal + per-vertex +# RGB texture): 16,515,072 verts, ESRI imagery draped as KIND × helix-residue colour. +# It ships IN-REPO under .claude/maps/, not from the release — this session type can't +# upload release assets, so the 151 MB artifact is version-controlled as a 2-part split +# zip (iceland_dem_16m.z01 + iceland_dem_16m.zip, each < GitHub's 100 MB blob limit). +# Reassemble both parts into one full zip, then extract iceland_dem.helix.soa.gz into +# dist/; include_dir! embeds cockpit/dist/, so BodyHelix resolves it SAME-ORIGIN +# (manifest iceland_latest = iceland_dem.helix.soa.gz). Same-origin is required for the +# same CORS reason as the body wires — the github releases redirect sends no CORS header. +# The build context already carries .claude/maps via `COPY . /build/q2` above. +RUN cd /build/q2/.claude/maps \ + && zip -s 0 iceland_dem_16m.zip --out /tmp/iceland_dem_full.zip \ + && unzip -o /tmp/iceland_dem_full.zip -d /build/q2/cockpit/dist/ \ + && rm -f /tmp/iceland_dem_full.zip \ + && ls -lh /build/q2/cockpit/dist/iceland_dem.helix.soa.gz # The Berlin OSM helix bake — the /geo + /helix?scene=osm scene (manifest osm_latest). # It's 92 MB, so it lives in the release, not git; download it into dist/ so BodyHelix 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 new file mode 100644 index 000000000..a4eaafffa --- /dev/null +++ b/claude-notes/plans/2026-07-07-garmin-img-typed-rebuild.md @@ -0,0 +1,146 @@ +# 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. + +## Why (operator direction, 2026-07-07) + +The `/geo` building heuristics were judged "extremely terrible" against the +working `/osm` (Cesium-style slippy map) and `/helix` (surfel-interpolated +anatomy). Operator's spec, verbatim requirements: + +- "helix uses surfel interpolation" — the geo mesh needs that look. +- "you need building heuristics · anything color > same building · buildings + become houses · green becomes grass or forest texture · blue becomes water · + walking pathways vs street texture" +- "you have neither mastered the terrain nor the kurvenlineal AFTER + heuristics, nor the goureaud shading in the end" — **pipeline order is + mandatory: heuristics → kurvenlineal → terrain → Gouraud.** +- "there are so many algorithms to fix height overlay · start from basics and + use Garmin format — that has all the data you need" + +Garmin IMG replaces raster color-guessing with **typed lookups**: every +polygon/polyline carries a type code (building / water / forest / street / +walking path / contour), multi-level generalization (TRE levels = LOD +pyramid), contour polylines with elevation labels (and per-tile `.DEM` +subfiles in newer builds). + +## Sources (banked in `.claude/maps/`, see SOURCES.md) + +- `otm-iceland.zip` + `otm-iceland-contours.zip` — garmin.opentopomap.org +- `garmin-grand-canyon/475053{10,11,16,17}.img` — gpsfiledepot Arizona Topo + (map id 1). `47505316` = Grand Canyon Village tile. +- Other free source noted by operator: alternativaslibres.org. + +## Decoder — VALIDATED (scripts/garmin_proto.py) + +Python prototype, proven by rendering tile 47505316 at full detail: the +Colorado gorge, side canyons, village road grid, Bright Angel Trail, and +elevation-labelled contours all appear correctly. + +Format facts (hard-earned, keep): + +- **Container:** byte 0 = XOR key (0x00 on all our files); `DSKIMG` @0x10; + blocksize = `1 << (b[0x61] + b[0x62])`; FAT = 512-B entries from 0x600 + (flag 0x01, name 8B + typ 3B, size u32 @0x0C **valid in part-0 entry**, + part u16 @0x10, 240 u16 block pointers @0x20). Multi-part subfiles: + concatenate blocks by `part*240 + i` order. RGN sizes read 0 in naive + parsers — slice by block list, trust part-0 size when nonzero. +- **TRE:** bbox N/E/S/W int24 mapunits @0x15/18/1B/1E (deg = mu·360/2²⁴); + levels @u32 0x21 (4 B each: zoom|bit7=inherited, bits, nsubdiv u16); + subdivisions @u32 0x29 — 16-B records, **last level 14-B** (no `next`): + rgn_off u24, objtypes (0x10 pt / 0x20 idx-pt / 0x40 line / 0x80 poly), + center lon/lat int24, width u16 (bit15 = terminate), height u16, next u16. +- **RGN:** data section @u32 0x15, len @0x19. Subdiv block = sorted-by-rgn_off + spans; if N object kinds present, (N−1) u16 pointers prefix the block, kinds + in order pt/idx-pt/line/poly. Point: type, lbl u24 (bit23 = has-subtype), + dlon/dlat i16 (<< 24−bits). Line: byte0 = type(0..5)|bit6 dir|bit7 two-byte + len; lbl u24 (bit22 = extra-bit-per-node, bit23 = NET ref); dlon/dlat i16; + len u8/u16 (**bitstream bytes only, info byte separate**); info byte = + lon_base lo-nibble, lat_base hi-nibble. Poly: same but type mask 0x7F. +- **Bitstream (exact QMapShack CShiftReg semantics — the part everyone gets + wrong):** LSB-first bits. Per axis: 1 bit "same-sign"; if set, 1 more bit = + constant sign, values unsigned; else per-delta two's complement **with one + extra bit added to the width**. Width from base: `n = base+2` (base ≤ 9) + else `2·base−9`, then +1 if signed. **Continuation:** in signed mode raw == + `1<<(n−1)` (sign bit only) accumulates `2^(n−1)−1` and reads again; final + value < sign → `acc + tmp`; ≥ sign → `(tmp − 2^n) − acc`. Missing this + turns fine levels into random-walk mush (our level-4 bug). extra_bit = 1 + skip bit per pair before x. +- **LBL:** data @u32 0x15, offset multiplier `1 << b[0x1D]`, encoding @0x1E + (6 = 6-bit packed, 9 = 8-bit latin1 0-terminated). 6-bit: 3 bytes → 4 chars + hi-first; table `" A..Z~~~~~0..9~~~~~~"`; value > 0x2F terminates. + Contour labels = elevation, **feet** in US topo maps (÷3.28084 → m), + meters in OTM. +- **Type codes (this map family):** lines: 0x20/0x21/0x22 = minor/inter/major + land contour (0x23-0x25 depth), 0x18/0x1f/0x26 streams/rivers, ≤0x07 roads, + 0x0a/0x0b/0x16 trails/walking paths. Polys: 0x3c-0x49 water, 0x50 woods, + 0x14-0x16 park/reserve. Buildings (urban maps): poly 0x13 (+0x6x variants). + +## The ver-8 radix-grid wire (operator-derived, 2026-07-07) + +Operator insight chain, applied to the wire format: + +1. *"helix vector = 'zenith is tilted this axis' — even one hemisphere would + suffice"* — terrain normals live strictly in the upper hemisphere + (`n_up ∈ [0,1]`); tilt (slope) + azimuth (aspect) is ALL Gouraud needs to + reconstruct the polygons from surfels. +2. *"HHTL provides height + location + connectedness"* — for a grid + heightfield, connectedness is the address structure itself: row-major + neighbors ARE the topology. The 12 B/tri index (396 MB of Iceland's 710 MB + wire = 56%) encodes nothing the addressing doesn't know. Don't ship it. +3. *"or simply radix-like free cartesian deterministic location"* — the final + collapse: position isn't stored either. `i → (row, col)` by radix + decomposition; `(row, col) → (x, z)` deterministically from the header + bbox. The OGAR canon applied to the wire: the key prerenders the node with + zero value decode; location is convention, not data. + +What remains stored per vertex: **height (F16, 2 B)** — the one +non-deterministic coordinate — and **KIND (u8, 1 B)** into a header palette. +Everything else is derived at decode: positions (radix), index (grid loop), +normals (one gradient pass — the `terrain_normals` port), color (palette × +CurveRuler residue, both deterministic; CurveRuler stride-4-over-17 is +bit-exact integer, ~25 lines of TS). + +Iceland at 16.5 M verts: 710 MB → ~50 MB raw → ~20-30 MB gz (smooth +heightfields compress well). Single file, no split-zip workaround, 5-7× +wire collapse. Layout sketch (concept block kept ver-6/7-shaped for the +sidebar + LOD): + +``` +BSO2 | ver=8 u16 | nC u32 | W u32 | H u32 (nV = W·H, nT derived) +concepts block (guid | material | layer | label | centroid | v_range) +grid header: x0 f32 | z0 f32 | dx f32 | dz f32 | yscale f32 +palette: nK u8 | nK × rgb u8 +height F16 [W·H] | kind u8 [W·H] +labels_json | materials_json | HXFL trailer +``` + +Buildings (Berlin) ride on top as the one non-grid layer (ver-7 mesh section +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). +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 + polylines + elevations; all features drape on it. Iceland tiles also + carry `.DEM` subfiles (format: mkgmap DEM — parse later if wanted). +4. **Bake pipeline in the MANDATED order:** typed heuristics → kurvenlineal + (`helix::CurveRuler` per-class surfel texture, AFTER heuristics) → terrain + (houses extruded from terrain height; water/grass/paths flat) → **Gouraud + smooth per-vertex normals LAST** (surfel-interpolation look, not flat + facets). Encode `encode_mesh_bso2` ver-7. +5. **Scenes:** `/canyon` (Grand Canyon — new), Iceland re-bake, Berlin + (needs a Garmin Germany/Berlin extract or mkgmap build from + `.claude/maps/berlin-latest.osm.pbf`; java available for mkgmap). +6. **Verify visually per scene** (scratchpad shot.js harness, small-first). + +## References + +- QMapShack `CGarminPolygon.cpp` (GPL, consulted for format semantics only; + our code is an original implementation of the documented format). +- John Mechalas, "The Garmin IMG File Format" (imgformat.pdf). +- mkgmap (the encoder these maps were built with). diff --git a/claude-notes/plans/2026-07-07-iceland-dem-loader-wiring.md b/claude-notes/plans/2026-07-07-iceland-dem-loader-wiring.md new file mode 100644 index 000000000..60bbec643 --- /dev/null +++ b/claude-notes/plans/2026-07-07-iceland-dem-loader-wiring.md @@ -0,0 +1,146 @@ +# /ice Iceland DEM loader wiring — 2-part split zip → same-origin serve + +**Status:** wiring landed (commit `d0a0e347`), **not yet deploy-verified** +(no Docker/Railway in the authoring session). This note is the handover for +the session that verifies the deploy. + +## Overview + +The `/ice` cockpit scene renders a full-resolution **textured** Iceland DEM. +The bake is `iceland_dem.helix.soa.gz` — **151 MB**, BSO2 **ver 7** (F16 pos + +Signed360 normal + per-vertex RGB texture), **16,515,072 verts / 32,997,378 +tris / 3,584 HHTL tiles**. Colour = KIND classification +(ocean/green/rock/scree/ice/lava, from ESRI World Imagery luminance+saturation ++ elevation) × `helix::CurveRuler` golden-spiral residue brightness × imagery +luminance. + +Two hard constraints shaped how it ships: + +1. **> GitHub's 100 MB blob limit** → cannot be a normal git-tracked file. +2. **This session type cannot upload release assets** (verified: `urllib` and + `pygithub` both return 403 "Creating, editing, or deleting releases is not + permitted for this session type") → cannot live in the release like the + body/berlin wires. + +So it is version-controlled as a **2-part split zip** under `.claude/maps/`, +each part < 100 MB: + +``` +.claude/maps/iceland_dem_16m.z01 76 MB (part 1) +.claude/maps/iceland_dem_16m.zip 69 MB (last part = central directory) +``` + +## The loader chain (how /ice loads the bake) + +``` +.claude/maps/iceland_dem_16m.{z01,zip} (committed, 2-part split zip) + │ Dockerfile builder stage: reassemble + extract + ▼ +cockpit/dist/iceland_dem.helix.soa.gz (151 MB, single gzip) + │ include_dir! embeds cockpit/dist/ into the q2-cockpit binary + ▼ +server serves it SAME-ORIGIN at /iceland_dem.helix.soa.gz + │ BodyHelix.tsx fetchSoa(): reads body.manifest.json → iceland_latest + │ → fetch(`/${iceland_latest}`) same-origin FIRST (BodyHelix.tsx:568) + │ → inflate (DecompressionStream 'gzip' on raw .gz bytes) + │ → decode BSO2 ver 7 (per-vertex aColor from the rgb block) + ▼ +/ice renders the 16.5M-vert textured terrain +``` + +The release URL (`${REL}/iceland_dem.helix.soa.gz`) is only a **fallback** that +BodyHelix never reaches when the same-origin copy is present — and it would fail +anyway (the github releases redirect sends no CORS header; that is the whole +reason every big wire is staged same-origin). + +## What landed (commit `d0a0e347`) + +**`Dockerfile`** (builder stage, `debian:bookworm`): +- Added `zip unzip` to the `apt-get install` line. +- Replaced the old `cp .../iceland.helix.soa.gz` step (which staged the sparse + ~1M-vert scatter under the wrong bit of the manifest) with: + ```dockerfile + RUN cd /build/q2/.claude/maps \ + && zip -s 0 iceland_dem_16m.zip --out /tmp/iceland_dem_full.zip \ + && unzip -o /tmp/iceland_dem_full.zip -d /build/q2/cockpit/dist/ \ + && rm -f /tmp/iceland_dem_full.zip \ + && ls -lh /build/q2/cockpit/dist/iceland_dem.helix.soa.gz + ``` + `zip -s 0 --out ` is the canonical reassembly of a split + archive — `cat`-ing the parts does NOT work for `zip -s` splits. The `.z01` + part must sit in the same dir as the `.zip` (it does — both under + `.claude/maps/`, carried into the build by `COPY . /build/q2`). + +**`cockpit/public/body.manifest.json`**: +- `iceland_latest` was already `iceland_dem.helix.soa.gz` (unchanged). +- Corrected the stale `iceland_note` (it still described the old Terrarium + ~1M-vert bake) to the current 16.5M textured one. + +## Verified in the authoring session (native, no Docker) + +- Reassembly `zip -s 0 iceland_dem_16m.zip --out full.zip && unzip full.zip` + produces `iceland_dem.helix.soa.gz` — a **valid gzip** (`gzip -t` OK), + **710,102,109 B uncompressed**, header `42 53 4f 32 07 00` = **BSO2 ver 7**, + **sha256 `cb3ef34ed6f0e153315b9881ec9dbb15b7ef24364097176effecbd6fb46378d4`**. +- Header decodes to nc=3584, nv=16,515,072, nt=32,997,378. +- BodyHelix `fetchSoa()` fetches `/${iceland_latest}` same-origin first + (`BodyHelix.tsx:568`), matching the exact filename the Dockerfile stages. +- No stale references to the old name remain anywhere the deploy reads. + +## NOT verified (the handover work) + +The **Docker build + Railway deploy** was not run (no Docker daemon / Railway +access in the authoring session). To close the loop: + +1. **Build the image** (or push and let Railway build): + ```bash + docker build -t q2-ice-test . + ``` + Watch the reassembly step's `ls -lh` line print + `... 151M ... /build/q2/cockpit/dist/iceland_dem.helix.soa.gz`. + If `zip`/`unzip` are missing, the apt line didn't take — re-check line ~36. + +2. **Confirm the binary embeds it.** `include_dir!` embeds `cockpit/dist/` at + compile time, so the 151 MB file makes the binary ~151 MB larger (on top of + berlin 92 MB + the body wires). This is the intended, established pattern — + the user explicitly wants the full 16.5M-vert bake ("20M is not a budget, + it's what works flawlessly"). If the build OOMs or the embed is rejected, + that is the thing to report back, not silently downsample. + +3. **Runtime smoke test:** + ```bash + docker run -p 8080:8080 q2-ice-test + curl -sI http://localhost:8080/iceland_dem.helix.soa.gz # expect 200, ~151 MB + ``` + Then open `/ice` in a browser (WebGL) and confirm the textured terrain + renders — green/rock/ice/ocean, smooth, not needles, not muddy. + +## Gotchas / notes + +- **Do NOT re-split or re-shrink.** The 2-part zip is exactly what the user + asked for ("2 part zip"). If a rebake changes the bytes, regenerate the split + the same way and keep two parts: + ```bash + # from a fresh iceland_dem.helix.soa.gz: + zip -s 76m .claude/maps/iceland_dem_16m.zip iceland_dem.helix.soa.gz + # → iceland_dem_16m.z01 (76 MB) + iceland_dem_16m.zip (remainder) + ``` + Keep both parts < 100 MB. +- The old sparse `.claude/maps/iceland.helix.soa.gz` (15 MB) is still tracked + but now **unreferenced** by the deploy — left on disk deliberately (the + manifest note says so). Safe to delete in a later cleanup if desired. +- A live dev bake at `cockpit/public/iceland_dem.helix.soa.gz` is gitignored + (`.gitignore:66`) so a 151 MB working copy never gets committed by accident. +- If a release upload ever becomes possible in a future session, the artifact + could move to the `fma-body-soa-v3-v1` release like berlin (92 MB) and the + Dockerfile could `curl` it instead of reassembling the split zip — but that + is optional; the same-origin serve is what matters, and the split-zip path + already delivers it with zero external dependency. + +## Provenance / re-bake recipe + +- Fetch DEM + imagery: `scripts/fetch_iceland_dem.py` → `.demgrid` (DEMG v2: + elevation + ESRI RGB on one grid). +- Bake: `geo/src/bin/iceland_dem.rs` (feature `helix`) → + `bso2::encode_mesh_bso2(pos, nrm, rows, tris, concepts, labels, colors)` + (non-empty `colors` ⇒ ver 7). diff --git a/cockpit/public/body.manifest.json b/cockpit/public/body.manifest.json index 7ae58e69e..53415057d 100644 --- a/cockpit/public/body.manifest.json +++ b/cockpit/public/body.manifest.json @@ -2,7 +2,11 @@ "helix_latest": "body.20260629c.v6helix.soa.gz", "osm_latest": "berlin.helix.soa.gz", "iceland_latest": "iceland_dem.helix.soa.gz", - "iceland_note": "iceland_dem.helix.soa.gz = real Terrarium DEM heightfield (z9 keyless terrain-RGB, 1152x896 grid = 1,032,192 verts / 2,060,290 tris, equirectangular, true-scale elevation 0..0.0067 in [-1,1], sea level -> y=0). Solid landmass, replaces the sparse iceland.helix.soa.gz scatter (kept on disk, unreferenced). Baked by geo/src/bin/iceland_dem.rs via bso2::encode_mesh_bso2; DEM fetched by scripts/fetch_iceland_dem.py.", + "garmin_scenes": { + "iceland": "iceland_dem.v8grid.soa.gz" + }, + "garmin_note": "garmin_scenes = the /garmin/:location registry (mod-rewrite style). The page route /garmin/ is SPA-handled (BodyHelix parses the path); the DATA is served by the cockpit-server param route /api/garmin/:location, which resolves through this map to an embedded dist file. Adding a scene = bake the .soa.gz into dist (Dockerfile) + one entry here — no server code change. 'iceland' aliases the existing DEM bake; 'canyon' lands with the Garmin typed-rebuild bake (claude-notes/plans/2026-07-07-garmin-img-typed-rebuild.md).", + "iceland_note": "iceland_dem.helix.soa.gz = full-resolution TEXTURED Iceland DEM (BSO2 ver 7 = F16 pos + Signed360 normal + per-vertex RGB drape): 16,515,072 verts / 32,997,378 tris / 3,584 HHTL tiles. Colour = KIND classification (ocean/green/rock/scree/ice/lava, from ESRI World Imagery luminance+saturation + elevation) x helix::CurveRuler golden-spiral residue brightness x imagery luminance. Equirectangular, true-scale elevation, sea level -> y=0. Baked by geo/src/bin/iceland_dem.rs (encode_mesh_bso2, feature 'helix'); DEM elevation + ESRI imagery fetched on one grid by scripts/fetch_iceland_dem.py (DEMG v2). 151 MB gz -> version-controlled as a 2-part split zip under .claude/maps/ (iceland_dem_16m.z01 + .zip, each < 100 MB); the Dockerfile reassembles + extracts it into cockpit/dist/ so it serves same-origin. Replaces the sparse iceland.helix.soa.gz scatter (kept on disk, unreferenced).", "note": "20260629c re-emit from soa_v2 (geometry identical to 20260629b): 39 connective structures (ligaments / tendons / interosseous membranes / fascia / retinacula / iliotibial tract) reclassified out of the ORGAN and SKIN layers into the now-live CONNECTIVE layer 7 — they were FMA-filed under /viscera/solid_organ/ligament_organ, so the is_a walk tagged them viscus->organ and they floated in the organ view as tan limb-shaped strays (interosseous membrane of leg/forearm, calcaneal tendon, long plantar ligament). Carries the 20260629b fixes (teeth->skeleton, per-vessel slicer-fill diameter). BSO2 ver 6 = F16 pos + Signed360 NORMAL + HXFL trailer; Gouraud per-vertex shading. Published to fma-body-soa-v3-v1; Dockerfile pulls same-origin.", "verts": 4283525 -} +} \ No newline at end of file diff --git a/cockpit/public/iceland_dem.v8grid.soa.gz b/cockpit/public/iceland_dem.v8grid.soa.gz new file mode 100644 index 000000000..441c279f2 Binary files /dev/null and b/cockpit/public/iceland_dem.v8grid.soa.gz differ diff --git a/cockpit/src/BodyHelix.tsx b/cockpit/src/BodyHelix.tsx index 7a660a304..6b9c06df5 100644 --- a/cockpit/src/BodyHelix.tsx +++ b/cockpit/src/BodyHelix.tsx @@ -29,6 +29,11 @@ function pathScene(): string | null { const p = window.location.pathname; if (p === '/geo') return 'osm'; if (p === '/ice') return 'iceland'; + // Mod-rewrite style: /garmin/ → scene "garmin:". The slug is + // resolved SERVER-side (/api/garmin/:location → manifest garmin_scenes → dist file), + // so a new scene is a manifest entry + a bake — no client change. + const m = p.match(/^\/garmin\/([a-z0-9-]+)$/); + if (m) return `garmin:${m[1]}`; return null; } @@ -108,11 +113,131 @@ interface Decoded { } interface ConceptMeta { row: number; name: string; layer: number; cx: number; cy: number; cz: number; } +// ── ver-8 radix-grid decode ───────────────────────────────────────────────── +// The wire stores ONLY height (F16) + kind (u8 → header palette). Everything +// else is DETERMINISTIC from the address and reconstructed here: position by +// radix (i → (row, col) → x0 + col·dx, zrow[row]), the triangle index by the +// grid loop, normals by one gradient pass (true-scale, so display slope = +// real-world slope), and colour by palette[kind] × the CurveRuler golden-spiral +// residue (stride-4-over-17, bit-exact 64-bit integer — "phase is convention, +// not data"). 710 MB (ver-7) → ~50 MB raw for the same 16.5 M-vert Iceland. +const MIX_X = 0x9E3779B97F4A7C15n, MIX_Y = 0xC2B2AE3D27D4EB4Fn, MIX_Z = 0x165667B19E3779F9n; +const U64 = (v: bigint) => BigInt.asUintN(64, v); +/// mix(cell) % 17 — the only reading of the mix the ruler needs (start / k). +function mix17(cx: number, cy: number, cz: number): number { + const m = U64( + U64(BigInt.asUintN(64, BigInt(cx)) * MIX_X) + ^ U64(BigInt.asUintN(64, BigInt(cy)) * MIX_Y) + ^ U64(BigInt.asUintN(64, BigInt(cz)) * MIX_Z)); + return Number(m % 17n); +} +function decodeGrid(buf: ArrayBuffer): Decoded { + const dv = new DataView(buf); + const nC = dv.getUint32(6, true), W = dv.getUint32(10, true), H = dv.getUint32(14, true); + const nV = W * H, nT = (W - 1) * (H - 1) * 2; + let o = 18; + o += 16 * nC; // guid + o += nC; // material (unused) + const layerOff = o; o += nC; // LAYER u8 + const labelOff = o; o += 4 * nC; // label idx + const cenOff = o; o += 12 * nC; // centroid 3f (pre-remapped like ver-6/7) + o += 8 * nC; // vrange + const x0 = dv.getFloat32(o, true), dx = dv.getFloat32(o + 4, true), yscale = dv.getFloat32(o + 8, true); + o += 12; + const zrow = new Float32Array(buf.slice(o, o + 4 * H)); o += 4 * H; + const nK = dv.getUint8(o); o += 1; + const pal = new Uint8Array(buf.slice(o, o + 3 * nK)); o += 3 * nK; + const hf = new Uint16Array(buf.slice(o, o + 2 * nV)); o += 2 * nV; + const kinds = new Uint8Array(buf.slice(o, o + nV)); o += nV; + const labLen = dv.getUint32(o, true); o += 4; + let names: string[] = []; + try { const lj = JSON.parse(new TextDecoder().decode(new Uint8Array(buf.slice(o, o + labLen)))); names = lj.names ?? lj; } catch { /* names optional */ } + const cLayer = new Uint8Array(buf.slice(layerOff, layerOff + nC)); + + // positions: radix x, tabulated z, decoded F16 height. Display frame DIRECT — + // ver-8 synthesizes display coords; no (-x, z, y) source remap round-trip. + const heights = new Float32Array(nV); + for (let i = 0; i < nV; i++) heights[i] = HALF_LUT[hf[i]] * yscale; + const positions = new Float32Array(nV * 3); + const rowArr = new Uint32Array(nV); + const layer = new Float32Array(nV); + for (let r = 0, i = 0; r < H; r++) { + const z = zrow[r], cr = Math.min(r, nC - 1), li = cLayer[cr] || 8; + for (let c = 0; c < W; c++, i++) { + positions[i * 3] = x0 + c * dx; + positions[i * 3 + 1] = heights[i]; + positions[i * 3 + 2] = z; + rowArr[i] = cr; + layer[i] = li; + } + } + // normals: one central-difference gradient pass (one-sided at edges). True-scale + // heights → true-world slopes, byte-parity with the ver-7 baker's terrain_normals. + const normals = new Int8Array(nV * 3); + for (let r = 0; r < H; r++) { + const rm = Math.max(r - 1, 0), rp = Math.min(r + 1, H - 1); + const dzr = (zrow[rp] - zrow[rm]) || 1e-6; + for (let c = 0; c < W; c++) { + const i = r * W + c; + const cm = Math.max(c - 1, 0), cp = Math.min(c + 1, W - 1); + const gx = (heights[r * W + cp] - heights[r * W + cm]) / (((cp - cm) * dx) || 1e-6); + const gz = (heights[rp * W + c] - heights[rm * W + c]) / dzr; + const il = 127 / Math.hypot(gx, 1, gz); + normals[i * 3] = Math.round(-gx * il); + normals[i * 3 + 1] = Math.round(il); + normals[i * 3 + 2] = Math.round(-gz * il); + } + } + // colour: palette[kind] × the deterministic CurveRuler residue (±18% within-kind + // texture — the /helix surfel look). Lattice mixes cached per cell (≈ tens of + // thousands of distinct cells for 16.5 M verts → the BigInt cost is amortized). + const colors = new Uint8Array(nV * 3); + const DETAIL = 48; + const cache = new Map(); + const cell17 = (px: number, py: number, pz: number, d: number): number => { + const cx = Math.floor(px * d), cy = Math.floor(py * d), cz = Math.floor(pz * d); + const key = (cx + 512) + (cy + 512) * 1024 + (cz + 512) * 1048576; + let v = cache.get(key); + if (v === undefined) { v = mix17(cx, cy, cz); cache.set(key, v); } + return v; + }; + for (let i = 0; i < nV; i++) { + const px = positions[i * 3], py = positions[i * 3 + 1], pz = positions[i * 3 + 2]; + const start = cell17(px, py, pz, DETAIL); + const k = cell17(px, py, pz, DETAIL * 3); + const phase = (((start + 4 * k) % 17) / 16) * 2 - 1; + const sweet = 0.90 + 0.18 * phase; + const kb = kinds[i] * 3; + colors[i * 3] = Math.max(0, Math.min(255, Math.round(pal[kb] * sweet))); + colors[i * 3 + 1] = Math.max(0, Math.min(255, Math.round(pal[kb + 1] * sweet))); + colors[i * 3 + 2] = Math.max(0, Math.min(255, Math.round(pal[kb + 2] * sweet))); + } + // index: the grid loop — connectedness IS the address structure (baker winding). + const index = new Uint32Array(nT * 3); + let wI = 0; + for (let r = 0; r < H - 1; r++) { + for (let c = 0; c < W - 1; c++) { + const a = r * W + c, b = a + 1, d2 = a + W, e = d2 + 1; + index[wI++] = a; index[wI++] = b; index[wI++] = d2; + index[wI++] = b; index[wI++] = e; index[wI++] = d2; + } + } + const labelIdx = new Uint32Array(buf.slice(labelOff, labelOff + 4 * nC)); + const cen = new Float32Array(buf.slice(cenOff, cenOff + 12 * nC)); + const conceptList: ConceptMeta[] = []; + for (let c = 0; c < nC; c++) { + conceptList.push({ row: c, name: names[labelIdx[c]] ?? `concept ${c}`, layer: cLayer[c] || 8, + cx: -cen[c * 3], cy: cen[c * 3 + 2], cz: cen[c * 3 + 1] }); // source → display (-x,z,y) + } + return { nVerts: nV, nTris: nT, positions, index, colors, normals, layer, vrow: rowArr, concepts: nC, conceptList }; +} + function decode(buf: ArrayBuffer): Decoded { const dv = new DataView(buf); const magic = String.fromCharCode(dv.getUint8(0), dv.getUint8(1), dv.getUint8(2), dv.getUint8(3)); if (magic !== 'BSO2') throw new Error(`bad magic "${magic}"`); const ver = dv.getUint16(4, true); + if (ver === 8) return decodeGrid(buf); // radix-grid wire: height + kind only const posBytes = ver >= 4 ? 6 : 12; const nC = dv.getUint32(6, true), nV = dv.getUint32(10, true), nT = dv.getUint32(14, true); let o = 18; @@ -125,6 +250,8 @@ function decode(buf: ArrayBuffer): Decoded { const posOff = o; o += posBytes * nV; const helixOff = o; o += 6 * nV; // pos3 | nrm3 — we read the nrm half const rowOff = o; o += 4 * nV; + const colorOff = ver >= 7 ? o : -1; // ver-7: per-vertex RGB drape (real kind/imagery colour) + if (ver >= 7) o += 3 * nV; const idxOff = o; o += 12 * nT; void matOff; @@ -160,12 +287,17 @@ function decode(buf: ArrayBuffer): Decoded { const colors = new Uint8Array(nV * 3); const normals = new Int8Array(nV * 3); // rim-decoded unit normal (display frame), cheap i8 const layer = new Float32Array(nV); + // ver-7 carries a real per-vertex colour drape (the DEM baker's kind × imagery × helix texture); + // older bakes synthesize the colour from the concept layer. When present, the wire colour wins. + const wireColors = ver >= 7 ? new Uint8Array(buf.slice(colorOff, colorOff + nV * 3)) : null; for (let i = 0; i < nV; i++) { positions[i * 3] = -srcPos[i * 3]; positions[i * 3 + 1] = srcPos[i * 3 + 2]; positions[i * 3 + 2] = srcPos[i * 3 + 1]; const r0 = rowArr[i], li = cLayer[r0] || 8; - const rgb = conceptColor(li, r0); + const rgb = wireColors + ? ([wireColors[i * 3], wireColors[i * 3 + 1], wireColors[i * 3 + 2]] as [number, number, number]) + : conceptColor(li, r0); colors[i * 3] = rgb[0]; colors[i * 3 + 1] = rgb[1]; colors[i * 3 + 2] = rgb[2]; // Signed360 → unit normal: r=sinθ from the Fisher-Z RIM (its strength; saturated cliff // falls back to the polar partition), hemisphere sign from polar, φ from azimuth. Same @@ -242,24 +374,15 @@ uniform float uYMin; // decoded height range (display.y), measured once at loa uniform float uYMax; uniform float uExag; // geo relief exaggeration: the Iceland bake is true-scale (span ~0.0074 in the // [-1,1] frame), so raise the geometry to read as terrain. 1 for anatomy (untouched). -uniform float uTime; // seconds — drives the Kurvenlineal breathing (0 for anatomy). -uniform float uRuler; // Kurvenlineal displacement amplitude (0 for anatomy). +uniform float uTime; // retained-but-0: the Kurvenlineal residue is baked into the mesh, not animated. +uniform float uRuler; // retained-but-0: no shader ruler (the golden-spiral residue is baked in). varying vec3 vColor; -// THE KURVENLINEAL (helix::CurveRuler) — a GLSL port of the sculpt crate's deterministic -// stride-4-over-17 phase (gcd(4,17)=1 → a full 17-residue permutation). Not 100% bit-exact to -// the Rust CurveRuler (that walks a u64 place with integer mixing; here we hash floats), but -// faithful in spirit: a deterministic bipolar phase in ~[-1,1] regenerated from the vertex's -// lattice address — "phase is convention, not data" (OGAR D-QUANTGATE). It synthesises relief -// detail where the sparse DEM bake has gaps, and — animated by uTime — makes the terrain BREATHE. -float rhash(vec3 p){ return fract(sin(dot(floor(p), vec3(127.1, 311.7, 74.7))) * 43758.5453); } -float kurvenlineal(vec3 pos, float detail){ - vec3 cell = floor(pos * detail); // the lattice cell = the address - vec3 sub = floor(pos * detail * 3.0); // finer sub-lattice picks the step k - float place = floor(rhash(cell) * 65536.0); - float k = mod(floor(rhash(sub) * 991.0), 17.0); // k ∈ [0,17) - float idx = mod(mod(place, 17.0) + 4.0 * k, 17.0); // stride-4-over-17 coprime walk - return idx / 8.0 - 1.0; // ~[-1, 1] bipolar -} +// THE KURVENLINEAL is now baked into the mesh, not approximated here. The real +// helix::CurveRuler golden-spiral residue (stride-4-over-17) is applied at BAKE time in +// geo/src/bin/iceland_dem.rs (::ruler_phase) as per-vertex surface displacement + recomputed +// normals — so the residue is carried by the decoded position + Signed360 normal the same way +// the anatomy body carries it. The vertex shader below therefore only lifts the terrain by uExag +// and Gouraud-shades; the earlier GLSL float-hash approximation of the ruler has been removed. // Height-profile terrain palette for the geo bakes (Iceland DEM, OSM). Elevation is display.y. // In the Iceland bake it is TRUE-SCALE (not vertically exaggerated) → a tiny span (~[0, 0.0074]) // with ~39% ocean at EXACTLY 0 and a heavily-quantized lowland plateau holding ~58% of verts, @@ -299,18 +422,16 @@ void main(){ const vec3 L = vec3(-0.401, 0.783, 0.476); float ndl = max(abs(dot(n, L)), 0.0); float shade = min(0.34 + 0.20*(abs(n.y)*0.5+0.5) + 0.12*(-n.x*0.5+0.5) + 0.92*ndl, 1.3); - // uGeo is a dynamically-uniform branch (a uniform, not a varying): coherent, no divergence. - // uGeo == 0 -> EXACTLY the pre-existing aColor*shade, so anatomy scenes are byte-identical. - vColor = (uGeo > 0.5) ? terrainColor(position.y) * shade : aColor * shade; - // Geo relief + the KURVENLINEAL: raise the true-scale island by uExag, then add the deterministic - // stride-4-over-17 phase as extra relief that BREATHES over uTime — it synthesises detail where the - // sparse DEM bake has gaps and gives the horizon a living, breathing motion (not 100% physical: a - // mental horizon that motivates the next, DEM-true, improvement). Anatomy (uGeo==0): untouched. + // Colour is ALWAYS aColor now: for the DEM terrain aColor is the baked kind × imagery × helix + // texture (real green / ice / volcano / rock), for buildings + anatomy it is the concept tint — + // the shader no longer guesses colour from height. Gouraud shade applies the 3D lighting form. + vColor = aColor * shade; + // Geo relief: raise the true-scale island by uExag. The Kurvenlineal golden-spiral residue is + // already in position + normal (baked by the DEM baker), so the shader adds no ruler of its own. + // Anatomy (uGeo==0): untouched. vec3 dpos = position; if (uGeo > 0.5) { - float rp = kurvenlineal(position, 90.0); - float breathe = 0.55 + 0.45 * sin(uTime * 0.5 + rp * 6.2831); - dpos.y = position.y * uExag + rp * uRuler * breathe; + dpos.y = position.y * uExag; } gl_Position = projectionMatrix * modelViewMatrix * vec4(dpos, 1.0); }`; @@ -357,6 +478,15 @@ function mount(container: HTMLDivElement, d: Decoded, enabled: Float32Array, dir // sky dome is added. An empty `?scene=` resolves falsy → anatomy body → all three stay off. const sceneParam = new URLSearchParams(window.location.search).get('scene'); const isGeoScene = Boolean(sceneParam ?? pathScene()); + // TERRAIN vs BUILDINGS. A terrain scene (Iceland DEM) is a dense, watertight heightfield — + // the same kind of surface as the anatomy body, so it earns the height-recolour + relief + // exaggeration + breathing that make the body read in 3D. A building scene (OSM/Berlin) is + // thin extruded footprints: exaggerating those 10x and recolouring them by height turns each + // building into a black needle (the #88 regression the operator flagged). So beautification is + // gated on TERRAIN, never on "any geo" — buildings render plain baked colours (uGeo/uExag/uRuler + // all neutral, exactly as anatomy), which is a solid city, not a needle field. + const sceneName = sceneParam ?? pathScene(); + const isTerrainScene = sceneName === 'iceland' || Boolean(sceneName?.startsWith('garmin:')); // Height range for the geo palette: measured ONCE from the decoded position buffer (display.y). // The Iceland bake is true-scale so the span is tiny — normalizing against [-1,1] would flatten @@ -397,7 +527,11 @@ function mount(container: HTMLDivElement, d: Decoded, enabled: Float32Array, dir const applyIndex = () => { geom.setDrawRange(0, rebuild()); idxAttr.needsUpdate = true; }; geom.setDrawRange(0, rebuild()); - const uniforms = { uAlpha: { value: 1 }, uGeo: { value: isGeoScene ? 1 : 0 }, uYMin: { value: yMin }, uYMax: { value: yMax }, uExag: { value: isGeoScene ? 10 : 1 }, uTime: { value: 0 }, uRuler: { value: isGeoScene ? 0.03 : 0 } }; + // uGeo (height-recolour) + uExag (relief) are TERRAIN-only; buildings/anatomy keep uGeo=0 (baked + // aColor) + uExag=1. uRuler/uTime are retained-but-0: the Kurvenlineal golden-spiral residue is + // now baked into the mesh (geometry + normals), so the shader applies no ruler — the terrain is a + // static surface, not a shader-animated one. + const uniforms = { uAlpha: { value: 1 }, uGeo: { value: isTerrainScene ? 1 : 0 }, uYMin: { value: yMin }, uYMax: { value: yMax }, uExag: { value: isTerrainScene ? 15 : 1 }, uTime: { value: 0 }, uRuler: { value: 0 } }; const mat = new THREE.ShaderMaterial({ uniforms, vertexShader: VERT, fragmentShader: FRAG, side: THREE.FrontSide }); const mesh = new THREE.Mesh(geom, mat); scene.add(mesh); @@ -480,9 +614,8 @@ function mount(container: HTMLDivElement, d: Decoded, enabled: Float32Array, dir window.addEventListener('resize', onResize); const tick = () => { raf = requestAnimationFrame(tick); - // Kurvenlineal breathing: geo scenes advance uTime and render every frame so the terrain - // breathes; anatomy stays fully on-demand (uTime 0, dirty gating untouched — no idle cost). - if (isGeoScene) { uniforms.uTime.value = performance.now() / 1000; dirty.current = true; } + // Every scene renders on-demand (dirty gating) — no idle redraw. The terrain's Kurvenlineal + // residue is baked into the mesh (static surface), so there is no per-frame shader animation. // server-LOD lifecycle runs even on idle frames (the cascade tracks the static view too); // turning LOD off restores the full geometry. Both are cheap and bounded by the 220 ms poll. const tnow = performance.now(); @@ -553,6 +686,18 @@ async function fetchSoa(): Promise { // separate /osm slippy-map page are both untouched. const scene = new URLSearchParams(window.location.search).get('scene') ?? pathScene(); + // /garmin/ — mod-rewrite style: the SERVER resolves the slug through the + // manifest's garmin_scenes registry (/api/garmin/:location) and serves the bake. + // A 404 body lists the available slugs, so surface it verbatim. + if (scene?.startsWith('garmin:')) { + const loc = scene.slice('garmin:'.length); + const r = await fetch(`/api/garmin/${loc}`); + if (!r.ok) { + const detail = await r.json().then((j) => `${j.error}${j.available ? ` — available: ${j.available.join(', ')}` : ''}`).catch(() => `HTTP ${r.status}`); + throw new Error(`garmin scene "${loc}": ${detail}`); + } + return inflate(r); + } const key = scene ? `${scene}_latest` : 'helix_latest'; const stamped: string | undefined = man?.[key]; if (!stamped) { @@ -609,10 +754,16 @@ export default function BodyHelix() { // Geo scenes: the terrain bake stamps a single layer — present it as "terrain" and drop the // empty anatomy layers, so a map never reads as skin/muscle/skeleton (the layer *id* is kept, so // the show/hide toggle still filters the real geometry). Anatomy keeps the full LAYERS taxonomy. - const geoUI = Boolean(new URLSearchParams(window.location.search).get('scene') ?? pathScene()); + const geoScene = new URLSearchParams(window.location.search).get('scene') ?? pathScene(); + const geoUI = Boolean(geoScene); + // Geo bakes stamp a single layer; present it with a domain-true name (a MAP must never read as + // skin/muscle/skeleton) and drop the empty anatomy layers. Terrain vs buildings gets its own name. + const geoTerrain = geoScene === 'iceland' || Boolean(geoScene?.startsWith('garmin:')); + const geoLayerName = geoTerrain ? 'terrain' : 'buildings'; + const geoLayerColor = geoTerrain ? '#7c8f5c' : '#9aa7b4'; const activeLayers = geoUI && d - ? LAYERS.filter((l) => d.conceptList.some((c) => c.layer === l.id)).map((l) => ({ ...l, name: 'terrain', color: '#7c8f5c' })) + ? LAYERS.filter((l) => d.conceptList.some((c) => c.layer === l.id)).map((l) => ({ ...l, name: geoLayerName, color: geoLayerColor })) : LAYERS; const groups = activeLayers.map((l) => ({ l, items: d ? d.conceptList.filter((c) => c.layer === l.id && (!q || c.name.toLowerCase().includes(q))) : [], @@ -623,14 +774,30 @@ export default function BodyHelix() { const scene = new URLSearchParams(window.location.search).get('scene') ?? pathScene(); const title = scene === 'iceland' ? '/ice — Iceland height-profile terrain' - : scene === 'osm' ? '/geo — OSM terrain' + : scene === 'osm' ? '/geo — OSM buildings' + : scene?.startsWith('garmin:') ? `/garmin/${scene.slice(7)} — Garmin terrain` : scene ? `/helix?scene=${scene}` : '/helix — living anatomy browser'; const subtitle = d ? scene - ? `${d.nVerts.toLocaleString()} verts · ${d.concepts.toLocaleString()} structures · height-profile palette + sky` + ? `${d.nVerts.toLocaleString()} verts · ${d.concepts.toLocaleString()} structures · ${scene === 'iceland' || scene.startsWith('garmin:') ? 'height-profile palette + sky' : 'baked colours + sky'}` : `${d.nVerts.toLocaleString()} verts · ${d.concepts.toLocaleString()} structures · helix::Signed360 normals (Fisher-Z rim)` : 'loading canonical helix bake…'; + // Scene menu — the "menu item" leg of the /garmin arc: built-ins + every slug the + // manifest's garmin_scenes registry names. Plain navigation (each scene is a page). + const [garminScenes, setGarminScenes] = useState([]); + useEffect(() => { + fetch('/body.manifest.json').then((r) => (r.ok ? r.json() : null)) + .then((m) => { if (m?.garmin_scenes) setGarminScenes(Object.keys(m.garmin_scenes)); }) + .catch(() => {}); + }, []); + const sceneOptions = [ + { label: 'body (/helix)', path: '/helix' }, + { label: 'berlin (/geo)', path: '/geo' }, + { label: 'iceland (/ice)', path: '/ice' }, + ...garminScenes.map((s) => ({ label: `garmin: ${s}`, path: `/garmin/${s}` })), + ]; + const herePath = window.location.pathname; return (
@@ -643,6 +810,13 @@ export default function BodyHelix() {
{d && (
+ {activeLayers.map((l) => ( diff --git a/cockpit/src/main.tsx b/cockpit/src/main.tsx index 747d4e98c..5cddbb6ee 100644 --- a/cockpit/src/main.tsx +++ b/cockpit/src/main.tsx @@ -121,6 +121,11 @@ createRoot(document.getElementById('root')!).render( height-profile beautification (ocean/moss/rock/ice terrain palette + procedural sky dome). Client-side render only — no re-bake. */} } /> + {/* /garmin/:location — mod-rewrite style Garmin terrain scenes. The slug + resolves SERVER-side (/api/garmin/:location → manifest garmin_scenes → + embedded bake), so a new scene = one manifest entry + a Dockerfile bake, + no client change. BodyHelix parses the slug from the path (pathScene). */} + } /> {/* /cpic — CPIC pharmacogenomics cockpit (gene-first): {gene, diplotype, drug} → phenotype → recommendation, 2-hop NARS deduction over the real CPIC tables via POST /api/cpic/reason (the standalone cpic crate). Additive, gene-first diff --git a/crates/cockpit-server/src/main.rs b/crates/cockpit-server/src/main.rs index 075eaa5c6..614821fe8 100644 --- a/crates/cockpit-server/src/main.rs +++ b/crates/cockpit-server/src/main.rs @@ -244,6 +244,12 @@ async fn main() { .route("/api/osm/tile/:z/:x/:y", get(osm_tiles::osm_tile_meta_handler)) // The /OSM cockpit page — a slippy map over the tile material, HHTL live. .route("/osm", get(osm::osm_page_handler)) + // Garmin terrain scenes, mod-rewrite style: /api/garmin/:location resolves the + // location slug through the manifest's `garmin_scenes` map to an embedded dist + // bake and serves it. The PAGE /garmin/ is extension-less, so the SPA + // fallback serves index.html and BodyHelix parses the path — one param route, + // N scenes; adding a scene = Dockerfile bake + one manifest entry, no code. + .route("/api/garmin/:location", get(garmin_scene_handler)) // Health .route("/health", get(health_handler)); @@ -315,6 +321,91 @@ async fn shutdown_signal() { tracing::info!("shutting down"); } +// ── Garmin terrain scenes (/api/garmin/:location — mod-rewrite style) ──────── + +/// Resolve a `/garmin/:location` slug through the manifest's `garmin_scenes` +/// map. Pure so it unit-tests without the embedded dist: returns the dist +/// filename for a known slug, `None` for unknown/malformed ones. Slugs are +/// lowercase `[a-z0-9-]` — anything else is rejected before the map lookup +/// (no path traversal into the dist). +fn resolve_garmin_scene(manifest_json: &str, location: &str) -> Option { + if location.is_empty() + || !location + .chars() + .all(|c| c.is_ascii_lowercase() || c.is_ascii_digit() || c == '-') + { + return None; + } + let v: serde_json::Value = serde_json::from_str(manifest_json).ok()?; + v.get("garmin_scenes")? + .get(location)? + .as_str() + .map(String::from) +} + +/// List the available scene slugs (for the self-documenting 404). +fn garmin_scene_slugs(manifest_json: &str) -> Vec { + serde_json::from_str::(manifest_json) + .ok() + .and_then(|v| { + v.get("garmin_scenes").and_then(|m| { + m.as_object() + .map(|o| o.keys().cloned().collect::>()) + }) + }) + .unwrap_or_default() +} + +/// GET /api/garmin/:location — serve the bake registered for `location`. +/// The gz bytes are served raw (application/gzip, no Content-Encoding): +/// BodyHelix inflates client-side via DecompressionStream, same contract as +/// the same-origin body/berlin/iceland wires. +async fn garmin_scene_handler(axum::extract::Path(location): axum::extract::Path) -> Response { + #[cfg(feature = "embed-cockpit")] + { + let Some(man) = COCKPIT_DIST + .get_file("body.manifest.json") + .and_then(|f| std::str::from_utf8(f.contents()).ok()) + else { + return (StatusCode::INTERNAL_SERVER_ERROR, "manifest missing from embedded dist").into_response(); + }; + if let Some(fname) = resolve_garmin_scene(man, &location) { + if let Some(file) = COCKPIT_DIST.get_file(&fname) { + return ( + StatusCode::OK, + [(header::CONTENT_TYPE, "application/gzip")], + file.contents(), + ) + .into_response(); + } + return ( + StatusCode::NOT_FOUND, + Json(serde_json::json!({ + "error": format!("scene '{location}' is registered as '{fname}' but the file is not in the embedded dist — check the Dockerfile bake step"), + })), + ) + .into_response(); + } + return ( + StatusCode::NOT_FOUND, + Json(serde_json::json!({ + "error": format!("unknown garmin scene '{location}'"), + "available": garmin_scene_slugs(man), + })), + ) + .into_response(); + } + #[cfg(not(feature = "embed-cockpit"))] + { + let _ = location; + ( + StatusCode::SERVICE_UNAVAILABLE, + "garmin scenes need the embed-cockpit build", + ) + .into_response() + } +} + // ── Static file handler with SPA fallback ──────────────────────────────────── /// Serves embedded Vite build files. Falls back to index.html for SPA routing. @@ -1137,3 +1228,48 @@ async fn health_handler() -> Json { "cockpit": if cfg!(feature = "embed-cockpit") { "embedded (Vite React)" } else { "fallback shell" }, })) } + +#[cfg(test)] +mod garmin_scene_tests { + use super::{garmin_scene_slugs, resolve_garmin_scene}; + + const MAN: &str = r#"{ + "helix_latest": "body.soa.gz", + "garmin_scenes": { "iceland": "iceland_dem.helix.soa.gz", "canyon": "garmin_canyon.helix.soa.gz" } + }"#; + + #[test] + fn known_slug_resolves_to_dist_filename() { + assert_eq!( + resolve_garmin_scene(MAN, "iceland").as_deref(), + Some("iceland_dem.helix.soa.gz") + ); + assert_eq!( + resolve_garmin_scene(MAN, "canyon").as_deref(), + Some("garmin_canyon.helix.soa.gz") + ); + } + + #[test] + fn unknown_slug_is_none_and_slugs_list_them() { + assert_eq!(resolve_garmin_scene(MAN, "atlantis"), None); + let mut slugs = garmin_scene_slugs(MAN); + slugs.sort(); + assert_eq!(slugs, ["canyon", "iceland"]); + } + + #[test] + fn malformed_slugs_rejected_before_lookup() { + // path traversal / case / separators never reach the map + for bad in ["../index.html", "Iceland", "ice land", "a/b", "", "x_y"] { + assert_eq!(resolve_garmin_scene(MAN, bad), None, "{bad}"); + } + } + + #[test] + fn manifest_without_registry_is_empty() { + assert_eq!(resolve_garmin_scene("{}", "iceland"), None); + assert!(garmin_scene_slugs("{}").is_empty()); + assert_eq!(resolve_garmin_scene("not json", "iceland"), None); + } +} diff --git a/geo/Cargo.lock b/geo/Cargo.lock index 351e8d3f7..232b79d85 100644 --- a/geo/Cargo.lock +++ b/geo/Cargo.lock @@ -180,8 +180,9 @@ dependencies = [ name = "geo-hhtl" version = "0.1.0" dependencies = [ + "helix", "lance-graph-contract", - "ndarray", + "ndarray 0.17.2", "osmpbf", ] @@ -208,6 +209,13 @@ version = "0.17.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "ed5909b6e89a2db4456e54cd5f673791d7eca6732202bbf2a9cc504fe2f9b84a" +[[package]] +name = "helix" +version = "0.1.0" +dependencies = [ + "ndarray 0.17.2 (git+https://github.com/AdaWorldAPI/ndarray.git?branch=master)", +] + [[package]] name = "home" version = "0.5.12" @@ -316,6 +324,22 @@ dependencies = [ "rawpointer", ] +[[package]] +name = "ndarray" +version = "0.17.2" +source = "git+https://github.com/AdaWorldAPI/ndarray.git?branch=master#576028458796724b55ae216a9781d13c44f1781d" +dependencies = [ + "blake3", + "matrixmultiply", + "num-complex", + "num-integer", + "num-traits", + "paste", + "portable-atomic", + "portable-atomic-util", + "rawpointer", +] + [[package]] name = "num-complex" version = "0.4.6" diff --git a/geo/Cargo.toml b/geo/Cargo.toml index 018480ed6..4d0bd83a4 100644 --- a/geo/Cargo.toml +++ b/geo/Cargo.toml @@ -37,12 +37,17 @@ osm = ["dep:osmpbf"] # `helix` a strict superset of `osm`; declare the implication so `--features # helix` (and `cargo test --features helix`) compiles on its own instead of only # when `osm` happens to be co-enabled by the `osm_helix` bin's `required-features`. -helix = ["osm", "dep:ndarray", "dep:lance-graph-contract"] +helix = ["osm", "dep:ndarray", "dep:lance-graph-contract", "dep:helix"] [dependencies] osmpbf = { version = "0.3", optional = true } ndarray = { path = "../../ndarray", default-features = false, features = ["std", "splat3d"], optional = true } lance-graph-contract = { path = "../../lance-graph/crates/lance-graph-contract", features = ["guid-v3-tail"], optional = true } +# THE KURVENLINEAL — the real `helix::CurveRuler` golden-spiral residue (stride-4-over-17), +# the SAME crate the `sculpt` tool and the anatomy `/helix` body use. The DEM baker regenerates +# deterministic fine surface relief from each vertex's lattice address so the terrain carries the +# residue the coarse HHTL grid lacks. Path dep to the sibling clone (like lance-graph-contract). +helix = { path = "../../lance-graph/crates/helix", optional = true } # osm_bake → SPM1 (dep-free lib + osm reader). osm_helix → BSO2 /helix wire. [[bin]] diff --git a/geo/src/bin/iceland_dem.rs b/geo/src/bin/iceland_dem.rs index 483d73e0b..a5661b128 100644 --- a/geo/src/bin/iceland_dem.rs +++ b/geo/src/bin/iceland_dem.rs @@ -16,12 +16,16 @@ //! ``` //! Gzip `.soa` afterwards (BodyHelix fetches a gzip'd artifact). -use geo_hhtl::bso2::{encode_mesh_bso2, MeshConcept, CLASSID_GEO_V3}; +use geo_hhtl::bso2::{encode_grid_bso2, encode_mesh_bso2, MeshConcept, CLASSID_GEO_V3}; use geo_hhtl::hhtl::point_to_hhtl4; use geo_hhtl::osm_read::M_PER_DEG; use lance_graph_contract::canonical_node::{NodeGuid, TailVariant}; +// THE KURVENLINEAL — the real golden-spiral curve-ruler (stride-4-over-17). Same crate the +// `sculpt` tool + the anatomy body use; ported here so the terrain carries the residue. +use helix::CurveRuler; + /// Deep key zoom for the 4-tier HHTL cascade (matches the OSM bake's KEY_ZOOM). const KEY_ZOOM: u32 = 32; /// Terrain basin tier (family) — a natural-area basin, distinct from the OSM @@ -33,19 +37,21 @@ struct Dem { h: usize, west: f64, east: f64, - lats: Vec, // len h, row 0 = north - elev: Vec, // len w*h, row-major, row 0 = north, col 0 = west (metres) + lats: Vec, // len h, row 0 = north + elev: Vec, // len w*h, row-major, row 0 = north, col 0 = west (metres) + rgb: Vec<[u8; 3]>, // len w*h (DEMG v2) or empty (v1) — ESRI imagery drape, same order as elev } fn read_demgrid(path: &str) -> Result> { let b = std::fs::read(path)?; - if b.len() < 24 || &b[0..4] != b"DEMG" { + // Need bytes 0..32 for the fixed header (magic/ver/w/h/west/east); `f64_at(24)` reads 24..32. + if b.len() < 32 || &b[0..4] != b"DEMG" { return Err("not a DEMG file".into()); } let u32_at = |o: usize| u32::from_le_bytes(b[o..o + 4].try_into().unwrap()); let f64_at = |o: usize| f64::from_le_bytes(b[o..o + 8].try_into().unwrap()); let ver = u32_at(4); - if ver != 1 { + if ver != 1 && ver != 2 { return Err(format!("unsupported DEMG version {ver}").into()); } let w = u32_at(8) as usize; @@ -64,15 +70,98 @@ fn read_demgrid(path: &str) -> Result> { elev.push(f32::from_le_bytes(b[o..o + 4].try_into().unwrap())); o += 4; } - Ok(Dem { w, h, west, east, lats, elev }) + // DEMG v2: a per-cell RGB imagery drape (ESRI World Imagery) after the elevation block, + // same row-major order. The baker classifies each cell's colour into a surface KIND. + let mut rgb = Vec::new(); + if ver == 2 { + if b.len() < o + n * 3 { + return Err("DEMG v2 truncated: rgb block short".into()); + } + rgb.reserve(n); + for _ in 0..n { + rgb.push([b[o], b[o + 1], b[o + 2]]); + o += 3; + } + } + Ok(Dem { w, h, west, east, lats, elev, rgb }) +} + +/// Surface KIND — the terrain node's CLASS axis (the "kind" that sits beside the HHTL position +/// address and the Helix residue). Classified from the ESRI imagery colour + elevation; each kind +/// carries a canonical palette so the map reads as clean ocean / green / rock / scree / ice / lava +/// rather than a muddy raw photo. +#[derive(Clone, Copy)] +#[repr(u8)] +enum Kind { + Ocean = 0, + Green = 1, // vegetation / moss / tundra + Rock = 2, // brown highland rock + Scree = 3, // grey bare + Ice = 4, // snow / glacier + Lava = 5, // dark volcanic basalt / black sand +} + +/// Discriminant order for the ver-8 palette block (index = `Kind as u8`). +const KIND_ORDER: [Kind; 6] = [ + Kind::Ocean, + Kind::Green, + Kind::Rock, + Kind::Scree, + Kind::Ice, + Kind::Lava, +]; + +impl Kind { + /// Canonical base colour (sRGB, 0..255) for the kind — a SWEET, clean palette (the /helix + /// conceptColor look), not the muddy raw satellite tone. The imagery classifies the kind; the + /// colour is the kind's own pleasant hue, brightness-varied by the helix residue below. + fn color(self) -> [f32; 3] { + match self { + Kind::Ocean => [26.0, 62.0, 104.0], // clean deep blue + Kind::Green => [104.0, 150.0, 70.0], // fresh grass/moss green + Kind::Rock => [150.0, 120.0, 88.0], // warm tan-brown highland + Kind::Scree => [166.0, 160.0, 150.0], // light warm grey + Kind::Ice => [238.0, 244.0, 252.0], // bright glacier white + Kind::Lava => [66.0, 46.0, 44.0], // dark basalt (not black) + } + } +} + +/// Classify a vertex from its imagery colour + elevation. Elevation gates ocean (y == 0); the +/// imagery luminance/saturation separates ice (bright, low-sat), lava (very dark), green +/// (green-dominant), scree (bright grey) and rock (the mid-tone default). +fn classify_kind(rgb: [u8; 3], elev: f32) -> Kind { + if elev <= 0.0 { + return Kind::Ocean; + } + let (r, g, b) = (rgb[0] as f32, rgb[1] as f32, rgb[2] as f32); + let lum = 0.299 * r + 0.587 * g + 0.114 * b; + let mx = r.max(g).max(b); + let mn = r.min(g).min(b); + let sat = if mx > 0.0 { (mx - mn) / mx } else { 0.0 }; + if lum > 172.0 && sat < 0.18 { + Kind::Ice // bright + desaturated → snow / glacier + } else if lum < 60.0 { + Kind::Lava // very dark → volcanic basalt / black sand + } else if g > r * 1.04 && g > b * 1.10 && g > 58.0 { + Kind::Green // green channel clearly dominant → vegetation / moss + } else if lum > 128.0 && sat < 0.22 { + Kind::Scree // bright-ish grey → scree / bare + } else { + Kind::Rock // mid-tone default → highland rock + } } fn main() { let mut args = std::env::args().skip(1); let (Some(input), Some(output)) = (args.next(), args.next()) else { - eprintln!("usage: iceland_dem "); + eprintln!("usage: iceland_dem [--grid]"); std::process::exit(2); }; + // --grid → the ver-8 radix-grid wire (height F16 + kind u8 only; location / + // connectedness / tilt / residue all client-derived from the address). + // Default stays the ver-7 mesh wire the deployed artifact uses. + let grid_mode = args.next().as_deref() == Some("--grid"); let dem = match read_demgrid(&input) { Ok(d) => d, Err(e) => { @@ -102,7 +191,11 @@ fn main() { let i = r * w + c; mx[i] = x; mz[i] = z; - my[i] = dem.elev[i].max(0.0); // sea level & bathymetry → 0 + // Sea level & bathymetry → 0; clamp the top to Iceland's physical max (Hvannadalshnúkur + // ≈ 2110 m) so a Terrarium nodata sentinel or tile-decode outlier can't become a needle + // spike in the heightfield. A real heightfield is continuous, so the only way to needle + // is an out-of-range vertex — this removes that failure mode at the source. + my[i] = dem.elev[i].clamp(0.0, 2200.0); } } @@ -118,42 +211,47 @@ fn main() { let cz = (loz + hiz) * 0.5; let half = ((hix - lox).max(hiz - loz) * 0.5).max(1.0); let inv = 1.0 / half; - // DISPLAY frame (Dx, Dy=up=elevation, Dz). No vertical exaggeration — the - // shader raises the true-scale island by uExag=10 + the Kurvenlineal. + // DISPLAY frame (Dx, Dy=up=elevation, Dz). No macro exaggeration here — the shader raises the + // true-scale island by uExag; the Kurvenlineal golden-spiral residue is baked in below. let mut pos = vec![[0.0f32; 3]; w * h]; for i in 0..w * h { pos[i] = [(mx[i] - cx) * inv, my[i] * inv, (mz[i] - cz) * inv]; } - // ── Per-vertex normals from the terrain gradient (display frame). ── - // du ≈ +east (col+), dv ≈ +north (row-, since row 0 is north). cross(dv,du) - // points +Dy on flat ground and tilts with the slope. One-sided at edges. let idx = |r: usize, c: usize| r * w + c; - let sub = |a: [f32; 3], b: [f32; 3]| [a[0] - b[0], a[1] - b[1], a[2] - b[2]]; - let cross = |a: [f32; 3], b: [f32; 3]| { - [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]] - }; - let mut nrm = vec![[0.0f32, 1.0, 0.0]; w * h]; - for r in 0..h { - let rn = r.saturating_sub(1); // north neighbour (smaller row) - let rs = (r + 1).min(h - 1); // south neighbour - for c in 0..w { - let cw = c.saturating_sub(1); - let ce = (c + 1).min(w - 1); - let du = sub(pos[idx(r, ce)], pos[idx(r, cw)]); // east tangent - let dv = sub(pos[idx(rn, c)], pos[idx(rs, c)]); // north tangent - let mut n = cross(dv, du); - let m = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt(); - if m > 1e-12 { - n = [n[0] / m, n[1] / m, n[2] / m]; - if n[1] < 0.0 { - n = [-n[0], -n[1], -n[2]]; // keep the up-facing hemisphere - } - } else { - n = [0.0, 1.0, 0.0]; - } - nrm[idx(r, c)] = n; - } + + // Helix CurveRuler detail frequency — used for the within-KIND colour texture below (NOT the + // geometry). At the full 16.5M-vert native resolution the heightfield already carries the fine + // relief, so the golden-spiral residue is NOT displaced into the geometry: `ruler_phase` is a + // per-cell STEP function, and displacing a dense mesh by it strips the surface into stepped + // ridges. The residue instead rides the COLOUR (within-kind variation), where a step reads as + // texture, not spikes. (On the sparse/anatomy meshes the geometry displacement is fine — cells + // are large relative to the mesh; on this dense grid it is not.) + const RES_DETAIL: f32 = 48.0; + let nrm = terrain_normals(&pos, w, h); + + // ── Per-vertex COLOUR = the KIND axis. Classify each vertex's surface class from the ESRI + // imagery colour + elevation, take the kind's canonical palette colour, then texture it by + // (a) the real imagery luminance (light/shadow within the kind) and (b) the Helix residue + // phase (deterministic within-kind variation from the address). Together with the HHTL + // address (per-concept key) and the Helix residue (displacement above), this completes the + // node's three axes: WHERE (HHTL) · SURFACE (Helix) · WHAT (kind). No imagery → a neutral + // moss default so a v1 demgrid still bakes. ── + let has_img = !dem.rgb.is_empty(); + let mut colors: Vec<[u8; 3]> = Vec::with_capacity(w * h); + for i in 0..w * h { + let px = if has_img { dem.rgb[i] } else { [96, 110, 92] }; + let base = classify_kind(px, my[i]).color(); + // SWEET brightness: the kind's clean base, varied WITHIN the kind by the helix residue + // (the "surfel" texture, ±18%) — this is what makes it read like the /helix body, not a + // flat class. Plus a GENTLE, CENTRED imagery light/shadow (±14%, never a muddy darkening): + // bright imagery lifts, dark imagery dips, mid stays put. No harsh luminance multiply. + let sweet = 0.90 + 0.18 * ruler_phase(pos[i], RES_DETAIL); // helix within-kind variation + let lum = (0.299 * px[0] as f32 + 0.587 * px[1] as f32 + 0.114 * px[2] as f32) / 255.0; + let light = 1.0 + 0.28 * (lum - 0.5); // centred real light/shadow, keeps the sweet base + let m = sweet * light; + let q = |c: f32| (c * m).round().clamp(0.0, 255.0) as u8; + colors.push([q(base[0]), q(base[1]), q(base[2])]); } // ── Triangles: two per grid cell, shared vertices (watertight surface). ── @@ -208,13 +306,52 @@ fn main() { } let labels = br#"{"names":["iceland-terrain"]}"#; - let (soa, blocks) = encode_mesh_bso2(&pos, &nrm, &rows, &tris, &concepts, labels); + let (soa, blocks) = if grid_mode { + // ── ver-8 radix-grid: store ONLY height + kind; the client derives + // position (radix x + zrow table), index (grid loop), normals + // (gradient pass) and colour (palette × CurveRuler residue). ── + let palette: Vec<[u8; 3]> = KIND_ORDER + .iter() + .map(|k| { + let c = k.color(); + [c[0] as u8, c[1] as u8, c[2] as u8] + }) + .collect(); + let mut kinds = Vec::with_capacity(w * h); + for i in 0..w * h { + let px = if has_img { dem.rgb[i] } else { [96, 110, 92] }; + kinds.push(classify_kind(px, my[i]) as u8); + } + // Normalize heights to ≤1 so the F16 mantissa carries full precision + // (raw display y spans ~0..0.007 — storing that directly wastes range). + let ymax = pos.iter().map(|p| p[1]).fold(0.0f32, f32::max).max(1e-9); + 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(); + encode_grid_bso2( + w as u32, h as u32, x0, dx, ymax, &zrow, &heights, &kinds, &palette, &concepts, + labels, + ) + } else { + encode_mesh_bso2(&pos, &nrm, &rows, &tris, &concepts, labels, &colors) + }; let nc = u32::from_le_bytes(soa[6..10].try_into().unwrap()); - let nv = u32::from_le_bytes(soa[10..14].try_into().unwrap()); - let nt = u32::from_le_bytes(soa[14..18].try_into().unwrap()); + // ver-8 header carries W,H at the ver-7 nV,nT offsets — derive the real counts. + let (nv, nt) = if grid_mode { + let gw = u32::from_le_bytes(soa[10..14].try_into().unwrap()) as u64; + let gh = u32::from_le_bytes(soa[14..18].try_into().unwrap()) as u64; + (gw * gh, (gw - 1) * (gh - 1) * 2) + } else { + ( + u64::from(u32::from_le_bytes(soa[10..14].try_into().unwrap())), + u64::from(u32::from_le_bytes(soa[14..18].try_into().unwrap())), + ) + }; eprintln!( - "BSO2: {nc} concepts · {nv} verts · {nt} tris · {} B soa · {} B blocks", + "BSO2 ver{}: {nc} concepts · {nv} verts · {nt} tris · {} B soa · {} B blocks", + if grid_mode { 8 } else { 7 }, soa.len(), blocks.len() ); @@ -228,3 +365,73 @@ fn main() { } eprintln!("wrote {output} + {blocks_path}"); } + +/// Per-vertex up-facing gradient normals for a `w×h` heightfield in the display frame. +/// `du ≈ +east` (col+), `dv ≈ +north` (row−, since row 0 is north); `cross(dv, du)` points +/// +Dy on flat ground and tilts with the slope. One-sided at the edges. Called twice: once for +/// the residue displacement direction, once to re-derive shading from the displaced surface. +fn terrain_normals(pos: &[[f32; 3]], w: usize, h: usize) -> Vec<[f32; 3]> { + let idx = |r: usize, c: usize| r * w + c; + let sub = |a: [f32; 3], b: [f32; 3]| [a[0] - b[0], a[1] - b[1], a[2] - b[2]]; + let cross = |a: [f32; 3], b: [f32; 3]| { + [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]] + }; + let mut nrm = vec![[0.0f32, 1.0, 0.0]; w * h]; + for r in 0..h { + let rn = r.saturating_sub(1); // north neighbour (smaller row) + let rs = (r + 1).min(h - 1); // south neighbour + for c in 0..w { + let cw = c.saturating_sub(1); + let ce = (c + 1).min(w - 1); + let du = sub(pos[idx(r, ce)], pos[idx(r, cw)]); // east tangent + let dv = sub(pos[idx(rn, c)], pos[idx(rs, c)]); // north tangent + let mut n = cross(dv, du); + let m = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt(); + if m > 1e-12 { + n = [n[0] / m, n[1] / m, n[2] / m]; + if n[1] < 0.0 { + n = [-n[0], -n[1], -n[2]]; // keep the up-facing hemisphere + } + } else { + n = [0.0, 1.0, 0.0]; + } + nrm[idx(r, c)] = n; + } + } + nrm +} + +// ── THE KURVENLINEAL (ported verbatim from `sculpt::sculpt` — the same helix::CurveRuler the +// anatomy body uses). Deterministic golden-spiral relief; the phase is regenerated from the +// address, NEVER stored. ── +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 a lattice cell into a u64 place anchor (wrapping_mul + xor fold). +fn mix(c: [i64; 3]) -> u64 { + (c[0] as u64).wrapping_mul(MIX_X) + ^ (c[1] as u64).wrapping_mul(MIX_Y) + ^ (c[2] as u64).wrapping_mul(MIX_Z) +} + +/// Deterministic bipolar relief in [-1, 1] from a vertex's lattice address. `cell = +/// floor(pos·detail)` → `place = mix(cell)`; a finer 3× sub-lattice picks `k = mix(sub) mod 17`; +/// value = `(CurveRuler::from_place(place).index(k) / 16)·2 − 1`. Same `pos` + `detail` → same +/// value on every bake, every session (phase is convention, not data — OGAR D-QUANTGATE). +fn ruler_phase(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); + // index ∈ [0,17) → idx/16 ∈ [0,1] → bipolar [-1,1] with 8 → 0. + (idx as f32 / 16.0) * 2.0 - 1.0 +} diff --git a/geo/src/bso2.rs b/geo/src/bso2.rs index fa09b7864..63a6b8f3d 100644 --- a/geo/src/bso2.rs +++ b/geo/src/bso2.rs @@ -343,11 +343,13 @@ pub struct MeshConcept { pub v_count: u32, } -/// Encode an **arbitrary display-frame triangle mesh** into the exact BSO2 ver-6 -/// `/helix` wire (F16 pos + `Signed360` normals + HXFL trailer) that -/// `BodyHelix.tsx` decodes — the generalization of [`encode_bso2`] for meshes -/// that are NOT OSM footprints (e.g. a DEM heightfield grid). The byte layout is -/// identical to [`encode_bso2`]'s so the same reader renders it unchanged. +/// Encode an **arbitrary display-frame triangle mesh** into the BSO2 `/helix` wire +/// (F16 pos + `Signed360` normals + HXFL trailer) that `BodyHelix.tsx` decodes — the +/// generalization of [`encode_bso2`] for meshes that are NOT OSM footprints (e.g. a DEM +/// heightfield grid). With an empty `colors` slice the byte layout is identical to +/// [`encode_bso2`]'s (ver-6); with a per-vertex `colors` drape it is ver-7 — a `rgb[nv]` +/// block after the row column carrying real satellite texture (the reader paints `aColor` +/// from it instead of guessing colour from height). /// /// `pos`/`nrm` are already the normalized DISPLAY frame `(Dx, Dy=up, Dz)`; /// positions are written as `srcPos = (-Dx, Dz, Dy)` (F16) and normals via @@ -363,16 +365,22 @@ pub fn encode_mesh_bso2( tris: &[[u32; 3]], concepts: &[MeshConcept], labels_json: &[u8], + colors: &[[u8; 3]], ) -> (Vec, Vec) { let nv = pos.len(); let nt = tris.len(); let nc = concepts.len(); assert_eq!(nrm.len(), nv, "nrm len must match pos"); assert_eq!(rows.len(), nv, "rows len must match pos"); + // Optional per-vertex true-colour drape. Empty → ver-6 (no colour block, the OSM-mesh + // layout). Non-empty → ver-7 = a `rgb[nv]` block after the row column (before the index), + // carrying the real satellite texture the DEM baker samples per vertex. + assert!(colors.is_empty() || colors.len() == nv, "colors len must be 0 or match pos"); + let ver: u16 = if colors.is_empty() { 6 } else { 7 }; - let mut o = Vec::with_capacity(nc * 40 + nv * 22 + nt * 12 + labels_json.len() + 64); + let mut o = Vec::with_capacity(nc * 40 + nv * 25 + nt * 12 + labels_json.len() + 64); o.extend_from_slice(b"BSO2"); - o.extend_from_slice(&6u16.to_le_bytes()); + o.extend_from_slice(&ver.to_le_bytes()); o.extend_from_slice(&(nc as u32).to_le_bytes()); o.extend_from_slice(&(nv as u32).to_le_bytes()); o.extend_from_slice(&(nt as u32).to_le_bytes()); @@ -408,6 +416,11 @@ pub fn encode_mesh_bso2( for r in rows { o.extend_from_slice(&r.to_le_bytes()); } + // per-vertex true-colour drape (ver-7 only): 3 bytes RGB, same order as pos. The client + // reads this straight into `aColor` for the geo terrain instead of a height-palette guess. + for c in colors { + o.extend_from_slice(c); + } for t in tris { for idx in t { o.extend_from_slice(&idx.to_le_bytes()); @@ -518,3 +531,123 @@ mod tests { assert_eq!(blocks.len(), 16); // one concept × 16 B } } + +/// Encode a **radix-grid heightfield** into the BSO2 **ver-8** wire — the collapse the +/// operator derived (2026-07-07): location, connectedness, tilt, and the kurvenlineal +/// residue are all DETERMINISTIC from the address, so only the two non-deterministic +/// per-vertex facts are stored: +/// +/// height (F16, display units) — the one coordinate the address can't know +/// kind (u8 → header palette) — the surface class +/// +/// The client reconstructs everything else at decode: positions by radix +/// (`i → (row, col) → (x0 + col·dx, z0 + row·dz)`), the triangle index by the grid +/// loop, normals by one gradient pass over the height grid, colour by +/// `palette[kind] × CurveRuler residue` (bit-exact stride-4-over-17). Iceland at +/// 16.5 M verts: 710 MB (ver-7) → ~50 MB raw. Layout: +/// +/// ```text +/// BSO2 | ver=8 u16 | nC u32 | W u32 | H u32 (nV = W·H, nT derived) +/// concepts block (guid | material | layer | label | centroid | v_range) ← ver-6/7 shape +/// grid: x0 f32 | dx f32 | yscale f32 | zrow f32[H] +/// palette: nK u8 | nK × rgb(3×u8) +/// height F16[W·H] | kind u8[W·H] +/// labels_json | materials_json | HXFL trailer +/// ``` +/// +/// x is radix-linear (`x = x0 + col·dx`); z is a per-ROW table (`zrow[row]`) because +/// Mercator-tile rows are not uniform in latitude (~13% spacing drift across Iceland) +/// — still deterministic-from-address, just tabulated once per row (H·4 B ≈ 14 KB). +/// +/// `heights[i]` is display-frame Dy BEFORE `yscale` (the client multiplies), row-major, +/// row 0 = the row at `z0`. The DISPLAY x is stored positive; the client applies its +/// own `(-x, z, y)` source remap convention internally — ver-8 positions are +/// synthesized directly in display space, so no remap round-trip is needed. +/// Returns `(bso2_bytes, blocks_sidecar_bytes)`. +#[must_use] +pub fn encode_grid_bso2( + w: u32, + h: u32, + x0: f32, // display x of col 0 + dx: f32, // display x step per col + yscale: f32, // multiplier applied to decoded F16 heights + zrow: &[f32], // len h — display z per row (Mercator rows are non-uniform) + heights: &[f32], // len w·h, display Dy before yscale (normalized ≤ ~1 for F16) + kinds: &[u8], // len w·h, indices into `palette` + palette: &[[u8; 3]], + concepts: &[MeshConcept], + labels_json: &[u8], +) -> (Vec, Vec) { + let nv = (w as usize) * (h as usize); + let nc = concepts.len(); + assert_eq!(zrow.len(), h as usize, "zrow len must be h"); + assert_eq!(heights.len(), nv, "heights len must be w*h"); + assert_eq!(kinds.len(), nv, "kinds len must be w*h"); + assert!(palette.len() <= 255, "palette must fit u8 count"); + assert!(kinds.iter().all(|&k| (k as usize) < palette.len()), "kind out of palette"); + + let mut o = Vec::with_capacity(nc * 40 + nv * 3 + labels_json.len() + 96); + o.extend_from_slice(b"BSO2"); + o.extend_from_slice(&8u16.to_le_bytes()); + o.extend_from_slice(&(nc as u32).to_le_bytes()); + o.extend_from_slice(&w.to_le_bytes()); + o.extend_from_slice(&h.to_le_bytes()); + // per-concept: guid | material | layer | label | centroid | (v_start,v_count) — + // byte-identical shape to ver-6/7 so the sidebar/LOD reader is shared. + for c in concepts { + o.extend_from_slice(c.key.as_bytes()); + } + o.extend(concepts.iter().map(|_| 4u8)); // material (doppler default) + o.extend(concepts.iter().map(|c| c.layer)); + for c in concepts { + o.extend_from_slice(&c.label.to_le_bytes()); + } + for c in concepts { + // Same (-Dx, Dz, Dy) pre-remap as ver-6/7 focus targets. + for v in [-c.centroid[0], c.centroid[2], c.centroid[1]] { + o.extend_from_slice(&v.to_le_bytes()); + } + } + for c in concepts { + o.extend_from_slice(&c.v_start.to_le_bytes()); + o.extend_from_slice(&c.v_count.to_le_bytes()); + } + // grid header: x radix-linear, z tabulated per row + for v in [x0, dx, yscale] { + o.extend_from_slice(&v.to_le_bytes()); + } + for &z in zrow { + o.extend_from_slice(&z.to_le_bytes()); + } + // palette + o.push(palette.len() as u8); + for p in palette { + o.extend_from_slice(p); + } + // the two stored axes + for &y in heights { + o.extend_from_slice(&F16::from_f32(y).0.to_le_bytes()); + } + o.extend_from_slice(kinds); + // labels_json + materials_json + HXFL (same tail contract as ver-6/7) + o.extend_from_slice(&(labels_json.len() as u32).to_le_bytes()); + o.extend_from_slice(labels_json); + let materials = b"{}"; + o.extend_from_slice(&(materials.len() as u32).to_le_bytes()); + o.extend_from_slice(materials); + o.extend_from_slice(b"HXFL"); + o.extend_from_slice(&(-2.256_794_5f32).to_le_bytes()); + o.extend_from_slice(&11.535_854f32.to_le_bytes()); + + // .blocks sidecar — same per-concept (center3, radius) shape. Radius from the + // concept's grid row span (v_range rows of a w-wide strip). + let mut blocks = Vec::with_capacity(nc * 16); + let half_row = 0.5 * (w as f32) * dx.abs(); + for c in concepts { + for v in c.centroid { + blocks.extend_from_slice(&v.to_le_bytes()); + } + blocks.extend_from_slice(&half_row.max(1e-4).to_le_bytes()); + } + (o, blocks) +} diff --git a/scripts/fetch_iceland_dem.py b/scripts/fetch_iceland_dem.py index 7a3052c22..224153cf7 100755 --- a/scripts/fetch_iceland_dem.py +++ b/scripts/fetch_iceland_dem.py @@ -1,24 +1,30 @@ #!/usr/bin/env python3 -"""Fetch a keyless Iceland DEM heightfield from AWS Terrarium terrain-RGB tiles. +"""Fetch a keyless Iceland DEM heightfield + true-colour imagery drape. -Terrarium tiles (`elevation-tiles-prod`, public, no key) encode elevation in the -PNG's RGB: ``elev_m = R*256 + G + B/256 - 32768``. We fetch the slippy-tile grid -covering Iceland at a chosen zoom, stitch into one elevation raster, block-mean -downsample to a tractable heightfield, and emit a compact ``.demgrid`` binary the -`iceland_dem` Rust baker reads (no HTTP dep in the geo crate). +Two keyless slippy-tile sources on the SAME grid: + * elevation — AWS Terrarium terrain-RGB (`elevation-tiles-prod`, no key): + ``elev_m = R*256 + G + B/256 - 32768``. + * imagery — ESRI World Imagery (`server.arcgisonline.com`, no key, z/y/x + order): real satellite true-colour (green vegetation, white glaciers, black + lava, brown rock) — the TEXTURE the height palette can only approximate. + +We fetch both grids at a chosen zoom, stitch, block-mean downsample, and emit a +compact ``.demgrid`` binary the `iceland_dem` Rust baker reads (no HTTP dep in +the geo crate). The imagery is draped per-vertex into the mesh's colour channel. Grid file layout (little-endian): "DEMG" 4 bytes magic - version u32 = 1 + version u32 = 2 (v1 = elev only; v2 adds the rgb drape below) W u32 columns (west->east, col 0 = west edge) H u32 rows (north->south, row 0 = north edge) west_lon f64 east_lon f64 (lon is LINEAR across columns — WebMercator x is linear in lon) lat[H] f64 latitude of each row centre (row 0 = north); non-linear in row elev[H*W] f32 metres, row-major, row 0 = north, col 0 = west + rgb[H*W*3] u8 v2 ONLY: true-colour drape, same row-major order as elev usage: - python3 scripts/fetch_iceland_dem.py OUT.demgrid [--zoom 9] [--downsample 4] + python3 scripts/fetch_iceland_dem.py OUT.demgrid [--zoom 9] [--downsample 4] [--no-imagery] """ import io import math @@ -35,8 +41,10 @@ LAT_S, LAT_N = 63.0, 67.0 TILE = 256 -BASE = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png" -UA = {"User-Agent": "q2-iceland-dem/1.0 (keyless terrarium bake)"} +DEM_BASE = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png" +# ESRI World Imagery is z/y/x (row before col), unlike the Terrarium z/x/y. +IMG_BASE = "https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}" +UA = {"User-Agent": "q2-iceland-dem/2.0 (keyless terrarium + esri imagery bake)"} def lon_to_tilex(lon, z): @@ -54,19 +62,63 @@ def tiley_to_lat(ty, z): return math.degrees(math.atan(math.sinh(n))) -def fetch_tile(z, x, y, tries=4): - url = BASE.format(z=z, x=x, y=y) +def fetch_rgb_tile(url, tries=4): + """Fetch one tile as an (H, W, 3) uint8 array, or None after `tries` failures.""" last = None for _ in range(tries): try: req = urllib.request.Request(url, headers=UA) data = urllib.request.urlopen(req, timeout=60).read() im = Image.open(io.BytesIO(data)).convert("RGB") - return x, y, np.asarray(im, dtype=np.float32) + return np.asarray(im, dtype=np.uint8) except Exception as e: # noqa: BLE001 last = e - print(f" tile {z}/{x}/{y} FAILED: {last!r}", file=sys.stderr) - return x, y, None + print(f" tile {url} FAILED: {last!r}", file=sys.stderr) + return None + + +def stitch_grid(base_url, zoom, xs, ys, x0, y0, decode, dtype, channels, label): + """Fetch `xs×ys` tiles from `base_url`, decode each, stitch into one raster. + `decode(rgb_uint8) -> array` maps a tile's raw RGB to the stored value. + Aborts (sys.exit) on ANY missing tile after retries — a hole would bake a + synthetic flat patch that looks like real data but isn't.""" + Wtiles, Htiles = len(xs), len(ys) + shape = (Htiles * TILE, Wtiles * TILE) + ((channels,) if channels > 1 else ()) + stitched = np.zeros(shape, dtype=dtype) + jobs = [(x, y) for y in ys for x in xs] + + def one(t): + x, y = t + return x, y, fetch_rgb_tile(base_url.format(z=zoom, x=x, y=y)) + + got, missing = 0, 0 + with ThreadPoolExecutor(max_workers=16) as ex: + for x, y, rgb in ex.map(one, jobs): + gx, gy = x - x0, y - y0 + if rgb is None: + missing += 1 + continue + stitched[gy * TILE:(gy + 1) * TILE, gx * TILE:(gx + 1) * TILE] = decode(rgb) + got += 1 + print(f"{label}: fetched {got} tiles, {missing} missing", file=sys.stderr) + if missing > 0: + print(f"{label}: {missing} tile(s) missing after retries — aborting " + f"(a hole would bake a synthetic patch)", file=sys.stderr) + sys.exit(1) + return stitched + + +def block_mean(a, down): + """Block-mean downsample the first two axes by `down` (crop to a multiple first).""" + if down == 1: + return a + Hpx, Wpx = a.shape[0], a.shape[1] + Hc, Wc = (Hpx // down) * down, (Wpx // down) * down + if a.ndim == 2: + return a[:Hc, :Wc].reshape(Hc // down, down, Wc // down, down).mean(axis=(1, 3)) + # (H, W, C) — mean each channel over the block. + c = a.shape[2] + return a[:Hc, :Wc].reshape(Hc // down, down, Wc // down, down, c).mean(axis=(1, 3)) def main(): @@ -76,11 +128,14 @@ def main(): out = sys.argv[1] zoom = 9 down = 4 + want_imagery = True for i, a in enumerate(sys.argv): if a == "--zoom": zoom = int(sys.argv[i + 1]) if a == "--downsample": down = int(sys.argv[i + 1]) + if a == "--no-imagery": + want_imagery = False x0 = int(math.floor(lon_to_tilex(LON_W, zoom))) x1 = int(math.floor(lon_to_tilex(LON_E, zoom))) @@ -90,41 +145,37 @@ def main(): xs = list(range(x0, x1 + 1)) ys = list(range(y0, y1 + 1)) print(f"zoom {zoom}: x {x0}..{x1} ({len(xs)} tiles), y {y0}..{y1} ({len(ys)} tiles) " - f"= {len(xs) * len(ys)} tiles", file=sys.stderr) + f"= {len(xs) * len(ys)} tiles/layer", file=sys.stderr) - Wtiles, Htiles = len(xs), len(ys) - stitched = np.zeros((Htiles * TILE, Wtiles * TILE), dtype=np.float32) + # ── Elevation (Terrarium terrain-RGB → metres). ── + def dem_decode(rgb): + rgb = rgb.astype(np.float32) + return rgb[:, :, 0] * 256.0 + rgb[:, :, 1] + rgb[:, :, 2] / 256.0 - 32768.0 - jobs = [(zoom, x, y) for y in ys for x in xs] - got, missing = 0, 0 - with ThreadPoolExecutor(max_workers=16) as ex: - for x, y, rgb in ex.map(lambda t: fetch_tile(*t), jobs): - gx = x - x0 - gy = y - y0 - if rgb is None: - missing += 1 - continue - elev = rgb[:, :, 0] * 256.0 + rgb[:, :, 1] + rgb[:, :, 2] / 256.0 - 32768.0 - stitched[gy * TILE:(gy + 1) * TILE, gx * TILE:(gx + 1) * TILE] = elev - got += 1 - print(f"fetched {got} tiles, {missing} missing", file=sys.stderr) - if missing > len(jobs) // 10: - print("too many missing tiles — aborting", file=sys.stderr) - sys.exit(1) + stitched = stitch_grid(DEM_BASE, zoom, xs, ys, x0, y0, dem_decode, np.float32, 1, "dem") + + # ── Imagery (ESRI World Imagery → true-colour RGB). ── + img_stitched = None + if want_imagery: + img_stitched = stitch_grid( + IMG_BASE, zoom, xs, ys, x0, y0, lambda rgb: rgb, np.uint8, 3, "imagery") # Geographic extent of the stitched tile grid (tile edges, exact). west = xs[0] / (1 << zoom) * 360.0 - 180.0 east = (xs[-1] + 1) / (1 << zoom) * 360.0 - 180.0 Hpx, Wpx = stitched.shape - # Block-mean downsample by `down` (crop to a multiple first). - Hc = (Hpx // down) * down - Wc = (Wpx // down) * down - st = stitched[:Hc, :Wc].reshape(Hc // down, down, Wc // down, down).mean(axis=(1, 3)) + st = block_mean(stitched, down) Hout, Wout = st.shape - print(f"stitched {Wpx}x{Hpx} px -> downsampled {Wout}x{Hout} verts " + print(f"stitched {Wpx}x{Hpx} px -> {Wout}x{Hout} verts " f"({Wout * Hout:,} verts)", file=sys.stderr) + img = None + if img_stitched is not None: + img = np.clip(block_mean(img_stitched.astype(np.float32), down), 0, 255).astype(np.uint8) + # Crop imagery to the exact elev grid (block_mean crops to a multiple; match shapes). + img = img[:Hout, :Wout, :] + # Per-row latitude: row r centre sits at tile-y = y0 + (r+0.5)*down/TILE. lats = np.empty(Hout, dtype=np.float64) for r in range(Hout): @@ -134,18 +185,22 @@ def main(): emin, emax = float(st.min()), float(st.max()) land = int((st > 0.5).sum()) print(f"elev range {emin:.1f}..{emax:.1f} m; land verts {land:,} " - f"({100 * land / st.size:.1f}%)", file=sys.stderr) + f"({100 * land / st.size:.1f}%){'; imagery drape' if img is not None else ''}", + file=sys.stderr) + version = 2 if img is not None else 1 with open(out, "wb") as f: f.write(b"DEMG") - f.write(struct.pack(" subfile slices; TRE -> map levels + subdivision tree; RGN -> +points / polylines / polygons with the delta bitstream. Validated visually: +the decoded Grand Canyon tile must LOOK like the canyon. +""" +import struct, sys + +def read_img(path): + b = open(path, 'rb').read() + xor = b[0] + if xor: + b = bytes(x ^ xor for x in b) + assert b[0x10:0x16] == b'DSKIMG', "not a Garmin IMG" + e1, e2 = b[0x61], b[0x62] + blocksize = 1 << (e1 + e2) + # FAT: 512-byte entries from 0x600 (skip the header pseudo-entry ' '*8). + subs = {} + o = 0x600 + while o + 512 <= len(b): + flag = b[o] + name = b[o+1:o+9].decode('ascii', 'replace') + typ = b[o+9:o+12].decode('ascii', 'replace') + if flag != 0x01: + if flag == 0x00: + break + o += 512 + continue + if name.strip() == '' or name[0] in ' \x00': + o += 512 + continue + size = struct.unpack('= 1 << 23: + v -= 1 << 24 + return v * 360.0 / (1 << 24) + +def i24(b, o, signed=True): + v = b[o] | b[o+1] << 8 | b[o+2] << 16 + if signed and v >= 1 << 23: + v -= 1 << 24 + return v + +def u16(b, o): return struct.unpack(' sub_off + sub_size: + break + S = Subdiv() + S.n = n + S.rgn_off = i24(tre, o, signed=False) & 0x0FFFFFFF + S.objtypes = tre[o+3] + S.clon = i24(tre, o+4) + S.clat = i24(tre, o+7) + wraw = u16(tre, o+10) + S.terminate = bool(wraw & 0x8000) + S.width = wraw & 0x7FFF + S.height = u16(tre, o+12) & 0x7FFF + S.next = u16(tre, o+14) if not last else 0 + S.level = li + S.children = [] + subdivs.append(S) + o += rec + n += 1 + return {'bbox': (north, east, south, west), 'levels': levels, 'subdivs': subdivs} + +if __name__ == '__main__': + path = sys.argv[1] if len(sys.argv) > 1 else '/home/user/q2/.claude/maps/garmin-grand-canyon/47505316.img' + subs = read_img(path) + print("subfiles:", {k: len(v) for k, v in subs.items()}) + tre_key = next(k for k in subs if k.endswith('.TRE')) + tre = subs[tre_key] + assert tre[2:12] == b'GARMIN TRE', tre[:16] + T = parse_tre(tre) + n, e, s, w = T['bbox'] + print(f"bbox N{mu2deg(n):.4f} E{mu2deg(e):.4f} S{mu2deg(s):.4f} W{mu2deg(w):.4f}") + for i, L in enumerate(T['levels']): + print(f"level {i}: zoom={L.zoom} bits={L.bits} subdivs={L.nsubdiv} inherited={L.inherited}") + print(f"total subdivs parsed: {len(T['subdivs'])}") + for S in T['subdivs'][:6]: + print(f" sd{S.n} L{S.level} rgn@{S.rgn_off} obj=0x{S.objtypes:02x} " + f"c=({mu2deg(S.clon):.4f},{mu2deg(S.clat):.4f}) wh=({S.width},{S.height}) " + f"term={S.terminate} next={S.next}") + +# ── RGN decode ────────────────────────────────────────────────────────────── + +class BitReader: + def __init__(self, data): + self.d = data + self.pos = 0 # bit position + def take(self, n): + v = 0 + for i in range(n): + byte = self.d[self.pos >> 3] + bit = (byte >> (self.pos & 7)) & 1 + v |= bit << i + self.pos += 1 + return v + def remaining(self): + return len(self.d) * 8 - self.pos + +def base_bits(base): + # imgformat: base 0-9 -> base+2 bits; base >9 -> 2*base-9 bits + return base + 2 if base <= 9 else 2 * base - 9 + +def decode_bitstream(data, lon_base, lat_base, extra_bit): + """Delta bitstream -> [(dlon, dlat)] in level units. + + Exact mirror of QMapShack CGarminPolygon/CShiftReg (the authoritative + open-source decoder): per-axis sign info up front (same-sign flag; if set, + one more bit = constant sign), per-delta two's complement otherwise, and — + the crucial part — the CONTINUATION marker: in signed mode a raw value of + 1<<(n-1) (sign bit only) accumulates (2^(n-1) - 1) and reads on; the final + value closes the sum (negative: tmp - 2^n - acc). + """ + class _Out(Exception): + pass + br = BitReader(data) + def take(n): + if br.remaining() < n: + raise _Out() + return br.take(n) + out = [] + try: + def axis_info(): + if take(1): # same-sign mode + return (False, -1 if take(1) else 1) + return (True, 0) # per-delta signed + x_has_sign, x_const = axis_info() + y_has_sign, y_const = axis_info() + nx = base_bits(lon_base) + (1 if x_has_sign else 0) + ny = base_bits(lat_base) + (1 if y_has_sign else 0) + xsign, ysign = 1 << (nx - 1), 1 << (ny - 1) + while True: + if extra_bit: + take(1) # routing node flag + if x_has_sign: + acc = 0 + while True: + tmp = take(nx) + if tmp != xsign: + break + acc += tmp - 1 # continuation: += 2^(n-1) - 1 + dlon = acc + tmp if tmp < xsign else (tmp - (xsign << 1)) - acc + else: + dlon = take(nx) * x_const + if y_has_sign: + acc = 0 + while True: + tmp = take(ny) + if tmp != ysign: + break + acc += tmp - 1 + dlat = acc + tmp if tmp < ysign else (tmp - (ysign << 1)) - acc + else: + dlat = take(ny) * y_const + out.append((dlon, dlat)) + except _Out: + pass + return out + +def parse_rgn(rgn, tre_info, want_levels=None): + """Yield dicts: {kind, type, level, coords[(lon,lat) mapunits], lbl}.""" + data_off = u32(rgn, 0x15) + data_len = u32(rgn, 0x19) + levels = tre_info['levels'] + # subdivisions that own RGN data, in file order (their rgn_off is ascending) + live = [S for S in tre_info['subdivs'] if S.objtypes] + live.sort(key=lambda S: S.rgn_off) + for i, S in enumerate(live): + S.rgn_end = live[i+1].rgn_off if i + 1 < len(live) else data_len + feats = [] + for S in live: + if want_levels is not None and S.level not in want_levels: + continue + bits = levels[S.level].bits + shift = 24 - bits + block = rgn[data_off + S.rgn_off: data_off + S.rgn_end] + if not block: + continue + kinds = [] # (kind, mask) present, canonical order + for kind, mask in (('point', 0x10), ('ipoint', 0x20), ('line', 0x40), ('poly', 0x80)): + if S.objtypes & mask: + kinds.append(kind) + nptr = len(kinds) - 1 + ptrs = [u16(block, 2*i) for i in range(nptr)] + starts = [2 * nptr] + ptrs + ends = ptrs + [len(block)] + for kind, s, e in zip(kinds, starts, ends): + seg = block[s:e] + o = 0 + if kind in ('point', 'ipoint'): + while o + 9 <= len(seg): + t = seg[o] + lbl = i24(seg, o+1, signed=False) + has_sub = bool(lbl & 0x800000) + dlon = struct.unpack('> 4 + bs = seg[o2+1: o2+1+blen] # blen = bitstream bytes (info byte separate) + o = o2 + 1 + blen + lon = S.clon + (dlon << shift) + lat = S.clat + (dlat << shift) + coords = [(lon, lat)] + for dlo, dla in decode_bitstream(bs, lon_base, lat_base, extra): + lon += dlo << shift + lat += dla << shift + coords.append((lon, lat)) + feats.append({'kind': kind, 'type': t, 'level': S.level, + 'coords': coords, 'lbl': lbl}) + return feats + +# ── LBL decode (6-bit) — enough for numeric contour labels ───────────────── +T6 = " ABCDEFGHIJKLMNOPQRSTUVWXYZ~~~~~0123456789~~~~~~" +def parse_lbl(lbl): + off = u32(lbl, 0x15) + ln = u32(lbl, 0x19) + mult = 1 << lbl[0x1D] + enc = lbl[0x1E] + return {'off': off, 'len': ln, 'mult': mult, 'enc': enc, 'data': lbl} + +def label_text(L, lbloff): + if lbloff == 0: + return '' + o = L['off'] + lbloff * L['mult'] + d = L['data'] + if L['enc'] == 9: # 8-bit + end = d.index(0, o) + return d[o:end].decode('latin1', 'replace') + out = [] + while o + 3 <= len(d): + b0, b1, b2 = d[o], d[o+1], d[o+2] + for c in (b0 >> 2, ((b0 & 3) << 4) | (b1 >> 4), ((b1 & 0xF) << 2) | (b2 >> 6), b2 & 0x3F): + if c > 0x2F: + return ''.join(out) + out.append(T6[c] if c < len(T6) else '?') + o += 3 + return ''.join(out)