facett-map3d-flight/tests/ground_loss_probe.rs
//! **MEASUREMENT: where the GROUND is lost between korp's `map_ways` and the 3D mesh.**
//!
//! An instrument, not a guard. It reads a REAL region-scale GeoParquet ways extract off
//! disk (no server, no S3, no GPU), runs it through the EXACT functions korp's
//! `rebuild_map3d` runs it through —
//!
//! ```text
//! facett_osm::read_ways_with(.., MapPalette::cartographic(day)) // == korp's classify
//! -> facett_map3d::extrude::densest_cluster_subset // == from_geometry step 1
//! -> facett_map3d::extrude::build_with_water_terrain_styled // == from_geometry step 2
//! ```
//!
//! — and prints the way/triangle census **per `WayKind`** at each hop. The question it
//! exists to answer is the one four local fixes did not: the 2D pane paints 16.89 %
//! green out of these same ways, and the 3D pane paints none, so which hop deletes the
//! cover?
//!
//! `#[ignore]` — it needs a local extract. Run:
//!
//! ```text
//! GROUND_PROBE_WAYS=/path/to/ways.parquet \
//! cargo test -p facett-map3d-flight --test ground_loss_probe -- --ignored --nocapture
//! ```
use facett_map::{MapPalette, Way};
use facett_map3d::extrude;
use facett_map3d::extrude::M_PER_DEG_LAT;
/// The census key: the classified `WayKind` string every downstream gate reads.
const KINDS: [&str; 7] = ["road", "path", "coastline", "water", "park", "building", "other"];
fn census(ways: &[Way]) -> Vec<(&'static str, usize, usize)> {
KINDS
.iter()
.map(|k| {
let sel: Vec<&Way> = ways.iter().filter(|w| w.kind == *k).collect();
let rings = sel.iter().filter(|w| w.fill.is_some()).count();
(*k, sel.len(), rings)
})
.collect()
}
fn print_census(label: &str, ways: &[Way]) {
let n = ways.len();
println!(" {label}: {n} ways");
for (k, c, rings) in census(ways) {
if c == 0 {
continue;
}
println!(
" {k:<10} {c:>7} ({:>5.2} % of set, {rings} filled rings)",
100.0 * c as f64 / n.max(1) as f64
);
}
}
/// Triangles in the opaque mesh whose source colour matches a way class' fill, counted
/// straight off `Mesh3::colors` — the only honest "did this class reach the mesh".
fn tri_by_color(mesh: &extrude::Mesh3, pal: &MapPalette) -> Vec<(&'static str, usize)> {
let want: [(&'static str, egui::Color32); 4] = [
("green (park/landuse)", pal.green),
("building", pal.building),
("land/terrain-quad", egui::Color32::from_rgb(38, 44, 40)),
("water fill", pal.water),
];
want.iter()
.map(|(n, c)| (*n, mesh.colors.iter().filter(|k| *k == c).count()))
.collect()
}
#[test]
#[ignore = "needs a local ways.parquet — set GROUND_PROBE_WAYS"]
fn where_the_ground_is_lost() {
let path = std::env::var("GROUND_PROBE_WAYS")
.expect("set GROUND_PROBE_WAYS to a GeoParquet ways file");
let cap: usize = std::env::var("GROUND_PROBE_MAX")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(320_000);
// The palette korp classifies AND paints with (cartographic DAY) — the same object
// `Korp::map_theme().palette` hands `map_paint::classify_ways`.
let pal = MapPalette::cartographic(true);
let ways = facett_osm::read_ways_with(&path, cap, &pal).expect("read ways");
println!("\n== HOP 1: korp `map_ways` (classified, what BOTH 2D and 3D are handed) ==");
print_census("map_ways", &ways);
let span = extrude::world_span_metres(&ways);
println!(" world span {:?} m", span);
println!(" metres_per_design_unit {:.1}", facett_map3d::metres_per_design_unit(&ways));
println!("\n== HOP 2: `densest_cluster_subset` (from_geometry step 1) ==");
let town = extrude::densest_cluster_subset(&ways);
print_census("town", &town);
let town_span = extrude::world_span_metres(&town);
println!(" world span {:?} m", town_span);
println!(" metres_per_design_unit {:.1}", facett_map3d::metres_per_design_unit(&town));
// The per-kind SURVIVAL rate — the number the whole document turns on.
println!(" survival by kind:");
let before = census(&ways);
let after = census(&town);
for (i, (k, c0, r0)) in before.iter().enumerate() {
let (_, c1, r1) = after[i];
if *c0 == 0 {
continue;
}
println!(
" {k:<10} {c0:>7} -> {c1:>7} ({:>6.2} % survive; filled rings {r0} -> {r1})",
100.0 * c1 as f64 / *c0 as f64
);
}
// WHY a ring dies: the three exits in `densest_cluster_subset`'s keep loop, counted
// on the same window that function computes. Re-derived here (read-only) because the
// function returns only the survivors.
println!("\n== HOP 2b: why the rings die (window classification) ==");
classify_ring_loss(&ways, &town);
println!("\n== HOP 3: `build_with_water_terrain_styled` (from_geometry step 2) ==");
let (mesh, tiles, water, _roof) = extrude::build_with_water_terrain_styled(
&town,
0.22,
extrude::GroundTerrain::Flat,
extrude::RoofStyle::Flat,
);
println!(" tiles(footprints) {tiles}");
println!(" opaque triangles {}", mesh.triangle_count());
println!(" water triangles {}", water.triangle_count());
for (n, c) in tri_by_color(&mesh, &pal) {
println!(
" tris coloured {n:<22} {c:>7} ({:>5.2} % of mesh)",
100.0 * c as f64 / mesh.triangle_count().max(1) as f64
);
}
let du_per_m = town_span.map(|s| extrude::DESIGN_FILL / s);
println!(" du_per_m {:?}", du_per_m);
println!(" ground_band {:e}", extrude::ground_band(du_per_m));
println!(" ground_top {:e}", extrude::ground_top(extrude::ground_band(du_per_m)));
println!(" building_rise_eps {:e}", extrude::building_rise_eps(du_per_m));
// The AREA the cover actually covers, in design units² — a triangle count says
// nothing about how much paper a class paints. This is the 2D "16.89 % green"
// question asked of the 3D mesh.
println!("\n== HOP 4: design-plane AREA by class (what a pixel measure would see) ==");
let areas = area_by_color(&mesh, &pal);
// ── THE VERDICT. Everything above is a table; this is the sentence. ───────────
//
// This file had ZERO assertions — `grep -c assert` = 0 — in the instrument written
// for "the 3D pane was reported at 0.00 % green". It printed its census and exited
// 0, so it would have printed `green 0.00 %` and STILL passed. A heading carries a
// verdict and a heading cannot go red; only an assertion can. These restate, as
// assertions, exactly the losses the probe exists to detect — each one on the
// measured quantity, not on infrastructure.
let share = |name: &str| areas.iter().find(|(n, _)| *n == name).map(|(_, s)| *s).unwrap_or(0.0);
let green = share("green (park/landuse)");
let building = share("building");
println!(
"\n VERDICT: cover {green:.2} % of the design box, buildings {building:.2} % \
(floors: {MIN_COVER_AREA_PCT} / {MIN_BUILDING_AREA_PCT})"
);
assert!(
town.iter().any(|w| w.kind == "park" || w.kind == "water"),
"the town window kept no cover ways at all — this is the 0.00 % report, one hop earlier"
);
assert!(
green > MIN_COVER_AREA_PCT,
"the extruded ground cover paints {green:.2} % of the design box — this IS the \
`0.00 % green` report, restated as a number the build can fail on"
);
assert!(
building > MIN_BUILDING_AREA_PCT,
"the buildings paint {building:.2} % of the design box — a 3D city with no city"
);
}
/// The floor under the cover class' share of the design box (percent). The 2D pane
/// paints **16.89 %** green out of these same ways and the mesh census reads 16.6 %, so
/// a real extract that clears anything at all clears this by an order of magnitude; it
/// is set to catch the 0.00 %, not to pin the exact figure of one extract.
const MIN_COVER_AREA_PCT: f64 = 1.0;
/// The same, for building footprints — a town window whose buildings vanish is the
/// other half of the same defect.
const MIN_BUILDING_AREA_PCT: f64 = 0.5;
/// Re-derive `densest_cluster_subset`'s town window from `ways` (same constants) and
/// count, for every filled ring, which of the three exits it takes.
fn classify_ring_loss(ways: &[Way], town: &[Way]) {
const WINDOW_M: f64 = 1500.0;
if town.len() == ways.len() {
println!(" (subset returned the input unchanged — no window was applied)");
return;
}
let mut pts: Vec<(f64, f64)> = Vec::new();
for w in ways {
if w.kind != "building" || w.lonlat.len() < 3 {
continue;
}
let (mut sx, mut sy, mut n) = (0.0, 0.0, 0.0);
for &(lon, lat) in &w.lonlat {
sx += lon;
sy += lat;
n += 1.0;
}
pts.push((sx / n, sy / n));
}
let (mut lon0, mut lat0, mut lon1, mut lat1) = (f64::MAX, f64::MAX, f64::MIN, f64::MIN);
for &(lon, lat) in &pts {
lon0 = lon0.min(lon);
lat0 = lat0.min(lat);
lon1 = lon1.max(lon);
lat1 = lat1.max(lat);
}
let lat_mid = (lat0 + lat1) * 0.5;
let m_per_deg_lon = M_PER_DEG_LAT * lat_mid.to_radians().cos().abs().max(1e-6);
let cell_deg_lat = WINDOW_M / M_PER_DEG_LAT;
let cell_deg_lon = WINDOW_M / m_per_deg_lon;
let nx = (((lon1 - lon0) / cell_deg_lon).ceil() as usize + 1).max(1);
let ny = (((lat1 - lat0) / cell_deg_lat).ceil() as usize + 1).max(1);
let cx = |lon: f64| (((lon - lon0) / cell_deg_lon) as usize).min(nx - 1);
let cy = |lat: f64| (((lat - lat0) / cell_deg_lat) as usize).min(ny - 1);
let mut bins = vec![0u32; nx * ny];
for &(lon, lat) in &pts {
bins[cy(lat) * nx + cx(lon)] += 1;
}
let best = bins.iter().enumerate().max_by_key(|(_, c)| **c).map(|(i, _)| i).unwrap_or(0);
let (bcy, bcx) = (best / nx, best % nx);
let town_lon = lon0 + (bcx as f64 + 0.5) * cell_deg_lon;
let town_lat = lat0 + (bcy as f64 + 0.5) * cell_deg_lat;
let (wlon, wlat) = (WINDOW_M / m_per_deg_lon, WINDOW_M / M_PER_DEG_LAT);
let (klon0, klon1) = (town_lon - wlon, town_lon + wlon);
let (klat0, klat1) = (town_lat - wlat, town_lat + wlat);
println!(" window lon [{klon0:.5}, {klon1:.5}] lat [{klat0:.5}, {klat1:.5}]");
let inside = |&(lon, lat): &(f64, f64)| {
lon >= klon0 && lon <= klon1 && lat >= klat0 && lat <= klat1
};
// **"No vertex inside" is TWO different things and they must not share a bucket.**
// A ring far from the town and a ring that ENCLOSES the town both have zero
// vertices in the window. The first is genuinely absent; the second is the bay the
// city sits on, and it is the one a vertex test can never see. Split by BBOX
// overlap, which can.
let bbox_of = |pts: &[(f64, f64)]| {
let (mut a0, mut b0, mut a1, mut b1) = (f64::MAX, f64::MAX, f64::MIN, f64::MIN);
for &(lon, lat) in pts {
a0 = a0.min(lon);
b0 = b0.min(lat);
a1 = a1.max(lon);
b1 = b1.max(lat);
}
(a0, b0, a1, b1)
};
println!(
" {:<9}{:>9}{:>10}{:>11}{:>11}{:>11}",
"kind", "rings", "elsewhere", "encloses", "straddles", "inside"
);
for k in ["park", "water", "building", "coastline"] {
let (mut elsewhere, mut encloses, mut straddle_ring, mut all_in) = (0usize, 0usize, 0usize, 0usize);
for w in ways.iter().filter(|w| w.kind == k && w.fill.is_some()) {
if w.lonlat.is_empty() {
elsewhere += 1;
continue;
}
let (a0, b0, a1, b1) = bbox_of(&w.lonlat);
let overlaps = a0 <= klon1 && a1 >= klon0 && b0 <= klat1 && b1 >= klat0;
if !w.lonlat.iter().any(inside) {
// No vertex in the window: either nowhere near it, or wrapped around it.
if overlaps {
encloses += 1;
} else {
elsewhere += 1;
}
} else if w.lonlat.iter().all(inside) {
all_in += 1;
} else {
straddle_ring += 1;
}
}
let tot = elsewhere + encloses + straddle_ring + all_in;
println!(" {k:<9}{tot:>9}{elsewhere:>10}{encloses:>11}{straddle_ring:>11}{all_in:>11}");
}
println!(
" `straddles` + `encloses` are what `clip_ring_to_window` recovers; `elsewhere`\n \
is genuinely not in this town and no clip can bring it back."
);
}
/// Total `(x, z)` area of the mesh triangles carrying each class colour, as a share of
/// the `[-1,1]` design box — the closest CPU-side answer to "what fraction of the pane
/// would this class paint if the camera looked straight down".
fn area_by_color(mesh: &extrude::Mesh3, pal: &MapPalette) -> Vec<(&'static str, f64)> {
let want: [(&'static str, egui::Color32); 4] = [
("green (park/landuse)", pal.green),
("building", pal.building),
("land/terrain-quad", egui::Color32::from_rgb(38, 44, 40)),
("water fill", pal.water),
];
let mut acc = vec![0.0f64; want.len()];
for (t, col) in mesh.colors.iter().enumerate() {
let Some(i) = want.iter().position(|(_, c)| c == col) else { continue };
let (a, b, c) = (
mesh.positions[mesh.indices[t * 3] as usize],
mesh.positions[mesh.indices[t * 3 + 1] as usize],
mesh.positions[mesh.indices[t * 3 + 2] as usize],
);
// Planar (x, z) area — the footprint the triangle casts on the ground.
let ar = 0.5
* ((b[0] - a[0]) as f64 * (c[2] - a[2]) as f64
- (c[0] - a[0]) as f64 * (b[2] - a[2]) as f64)
.abs();
acc[i] += ar;
}
let box_area = 4.0f64; // [-1,1]²
let mut shares = Vec::with_capacity(want.len());
for (i, (n, _)) in want.iter().enumerate() {
println!(" {n:<24} {:>8.4} du² ({:>6.2} % of the design box)", acc[i], 100.0 * acc[i] / box_area);
shares.push((*n, 100.0 * acc[i] / box_area));
}
shares
}