diff --git a/cube_bathymetry/include/cube_bathymetry/store_import.h b/cube_bathymetry/include/cube_bathymetry/store_import.h index d797bc9..5070e3c 100644 --- a/cube_bathymetry/include/cube_bathymetry/store_import.h +++ b/cube_bathymetry/include/cube_bathymetry/store_import.h @@ -146,9 +146,23 @@ namespace cube /// uncertainty squared, floored at a small positive epsilon (Node::setPredictedDepth /// requires a finite positive variance; a single-sample CUBE cell can persist a /// zero or non-finite uncertainty, so every finite-depth cell is still seeded). -/// Does NOT insert a CUBE hypothesis and does NOT mark the sheet dirty. +/// Does NOT mark the sheet dirty. +/// +/// @param seed_settled When true (default), ALSO reseed each cell's SETTLED depth +/// as a CUBE hypothesis (`setSettledDepthAt`) so the cell round-trips through +/// `values()` and survives the next whole-tile save -- the lossless warm-start +/// reload of a persisted draft tile (ADR-0001, #21/#70). When false, seed ONLY +/// the predicted surface: the cell carries a **blunder-rejection prior** (it +/// turns on `Node::insert`'s predicted-surface gate so a false-deep sounding +/// below `target - blunder_scalar*sqrt(var)` is rejected) but does NOT +/// fill/settle the cell -- no hypothesis, no `values()` output, the sheet stays +/// clean. This is the `Chart` (contour) prior path (cube#89): settling coarse +/// contour depths would contaminate the survey layer (and its co-estimated +/// backscatter) with non-measured fill, so the prior only gates, it does not +/// fill; survey-falls-through-to-chart gap-filling stays a query-time concern. void primeFromTile( - const marine_bathymetry_store::BathymetryTile & tile, GeoMapSheet & map_sheet); + const marine_bathymetry_store::BathymetryTile & tile, GeoMapSheet & map_sheet, + bool seed_settled = true); /// @brief Load every tile of @p layer from @p store into @p map_sheet, /// priming predicted depths cell-by-cell. @@ -158,10 +172,15 @@ namespace cube /// are seeded. Warm-starts slope correction from persisted draft tiles on /// restart; does not reconstruct CUBE hypothesis state (#21). A no-op if @p layer /// holds no tiles. +/// +/// @param seed_settled Forwarded to @ref primeFromTile (default true = +/// settled+predicted warm-start reload; false = predicted-only blunder-rejection +/// prior, e.g. priming the predicted surface from a `Chart` layer, cube#89). void loadIntoSheet( const marine_bathymetry_store::BathymetryStore & store, marine_bathymetry_store::SourceLayer layer, - GeoMapSheet & map_sheet); + GeoMapSheet & map_sheet, + bool seed_settled = true); } // namespace cube diff --git a/cube_bathymetry/src/import_bag_main.cpp b/cube_bathymetry/src/import_bag_main.cpp index e07b29f..aaf6c40 100644 --- a/cube_bathymetry/src/import_bag_main.cpp +++ b/cube_bathymetry/src/import_bag_main.cpp @@ -73,6 +73,11 @@ "node writes draft. #85)\n"; std::cout << " -o : Output bathymetry-store directory (created if " "needed)\n"; + std::cout << " --prior : seed the CUBE predicted surface from this " + "store's Chart (contour) layer so blunder rejection drops false-deep " + "detections (#89). The chart must be at the survey GGGS level. NOTE: a " + "coarse/shallow-biased chart can also reject LEGITIMATE deeper-than-charted " + "returns; the rejection margin is tunable via the blunder_* params\n"; std::cout << " --bs-store : Also write an MBES backscatter store layer " "(Processed) from the same CUBE pass (optional; surfaces the co-estimated " "intensity -- uncorrected by default, angle-corrected when " @@ -278,6 +283,7 @@ int main(int argc, char * argv[]) std::vector bagfile_names; std::string store_dir; + std::string prior_dir; // optional: Chart-prior store to seed the predicted surface (#89) std::string bs_store_dir; // optional: MBES backscatter store output (#80) std::string bathy_layer_str = "processed"; // bathy target layer (#85): draft|processed std::string detections_topic; // required @@ -318,6 +324,8 @@ int main(int argc, char * argv[]) usage(); } else if (*arg == "-o") { store_dir = next_value("-o"); + } else if (*arg == "--prior") { + prior_dir = next_value("--prior"); } else if (*arg == "--bathy-layer") { bathy_layer_str = next_value("--bathy-layer"); } else if (*arg == "--bs-store") { @@ -489,6 +497,81 @@ int main(int argc, char * argv[]) backscatter_mode, std::move(backscatter_curve.points), backscatter_curve.tl_removed, backscatter_curve.absorption_db_per_m); + // Chart-prior prime (#89): seed the CUBE predicted surface from the prior store's + // Chart (contour) layer BEFORE any soundings are added. primeFromTile lazy-creates + // a node per chart cell carrying a predicted depth, which turns ON CUBE's + // predicted-surface blunder gate (Node::insert is a pass-through when + // predicted_depth_ is NaN); the survey soundings then hit nodes that already carry + // a predicted depth, so a false-deep detection below + // `target - blunder_scalar*sqrt(predicted_var)` is rejected. seed_settled=false: + // seed ONLY the predicted surface, never settle the coarse contour into the survey + // layer (that would contaminate both the bathy and co-estimated backscatter). + // + // ALIGNMENT: load() restores each Chart tile at the GGGS level encoded in its + // filename (the store is multi-level; the fromCellSize() arg below only sets the + // store's default cellIndex() level, it does NOT pin the level tiles load at). A + // chart cell only coincides with -- and thus gates -- a survey node when its level + // matches the survey level, i.e. the Chart layer was imported at the survey's + // resolution. A chart tile at a different level primes a node the survey soundings + // never land on, so it cannot gate. Rather than silently build an un-gated store + // (the operator passed --prior to reject outliers), we validate below and refuse + // when nothing would gate. A future refinement could resample the chart. + if (!prior_dir.empty()) { + marine_bathymetry_store::BathymetryStore prior = + marine_bathymetry_store::BathymetryStore::fromCellSize( + static_cast(geo_map_sheet.nominalCellSizeMeters())); + marine_bathymetry_store::SourceRegistry prior_reg; + try { + marine_bathymetry_store::load(prior, prior_dir, &prior_reg); + } catch (const std::exception & e) { + std::cerr << "ERROR: --prior store '" << prior_dir << "' failed to load: " + << e.what() << std::endl; + return 1; + } + const auto & chart_tiles = + prior.tiles(marine_bathymetry_store::SourceLayer::Chart); + const std::size_t prior_tiles = chart_tiles.size(); + if (prior_tiles == 0) { + std::cerr << "ERROR: --prior store '" << prior_dir << "' has no Chart-layer " + << "tiles; nothing would seed the predicted surface and blunder " + << "rejection would be INACTIVE. Check the path and that the store " + << "has a chart/ layer (refusing to build an un-gated store)." + << std::endl; + return 1; + } + // Only tiles at the survey level can coincide with a survey node and gate + // (see ALIGNMENT above). The prior store's default level is fromCellSize() of + // the survey cell size, so it equals the survey sheet's level. + const uint8_t survey_level = prior.level().level(); + std::size_t mismatched = 0; + for (const auto & grid_tile : chart_tiles) { + if (grid_tile.second.index().level() != survey_level) { + ++mismatched; + } + } + if (mismatched == prior_tiles) { + std::cerr << "ERROR: all " << prior_tiles << " Chart tile(s) in '" << prior_dir + << "' are at a GGGS level other than the survey level " + << static_cast(survey_level) << "; none would coincide with a " + << "survey node, so blunder rejection would be INACTIVE. Re-import " + << "the chart at the survey resolution (refusing to build an " + << "un-gated store)." << std::endl; + return 1; + } + if (mismatched > 0) { + std::cerr << "WARNING: " << mismatched << " of " << prior_tiles << " Chart " + << "tile(s) are not at the survey level " + << static_cast(survey_level) << " and will not gate; only the " + << "level-matched tiles seed the predicted surface." << std::endl; + } + cube::loadIntoSheet( + prior, marine_bathymetry_store::SourceLayer::Chart, geo_map_sheet, + /*seed_settled=*/false); + std::cout << "Primed CUBE predicted surface from " << (prior_tiles - mismatched) + << " of " << prior_tiles << " Chart tile(s) in " << prior_dir + << " (blunder rejection active, #89)." << std::endl; + } + std::cout << "reading messages..." << std::endl; int ping_count = 0; diff --git a/cube_bathymetry/src/store_import.cpp b/cube_bathymetry/src/store_import.cpp index 1614183..5f9e392 100644 --- a/cube_bathymetry/src/store_import.cpp +++ b/cube_bathymetry/src/store_import.cpp @@ -149,7 +149,8 @@ mapSheetToBackscatterCells( } void primeFromTile( - const marine_bathymetry_store::BathymetryTile & tile, GeoMapSheet & map_sheet) + const marine_bathymetry_store::BathymetryTile & tile, GeoMapSheet & map_sheet, + bool seed_settled) { const gggs::GridIndex & grid = tile.index(); const std::vector & depth = tile.depthBand(); @@ -178,6 +179,16 @@ void primeFromTile( map_sheet.setPredictedDepthAt( *it, static_cast(d), static_cast(variance)); + if (!seed_settled) { + // Predicted-only prime (cube#89): seed the slope/blunder-rejection prior + // but do NOT settle the cell. A `Chart` (contour) prior must turn the + // blunder gate on WITHOUT filling the survey layer with coarse contour + // depths (which would also contaminate the co-estimated backscatter with + // non-measured cells). The survey accumulates real data on top; gap-filling + // survey->chart is a query-time overlay, not a settled fill. + continue; + } + // Lossless reload (ADR-0001): also reseed the SETTLED depth as a CUBE // hypothesis so the primed cell round-trips through values() and survives the // next whole-tile save. Without this, a tile primed at startup, partially @@ -191,12 +202,13 @@ void primeFromTile( void loadIntoSheet( const marine_bathymetry_store::BathymetryStore & store, marine_bathymetry_store::SourceLayer layer, - GeoMapSheet & map_sheet) + GeoMapSheet & map_sheet, + bool seed_settled) { // Single fused grid per layer (unh_marine_autonomy#221): iterate the layer's // tiles directly. Empty map -> no-op. for (const auto & grid_tile : store.tiles(layer)) { - primeFromTile(grid_tile.second, map_sheet); + primeFromTile(grid_tile.second, map_sheet, seed_settled); } } diff --git a/cube_bathymetry/test/test_store_import.cpp b/cube_bathymetry/test/test_store_import.cpp index 1e0c160..248e770 100644 --- a/cube_bathymetry/test/test_store_import.cpp +++ b/cube_bathymetry/test/test_store_import.cpp @@ -76,6 +76,25 @@ std::vector makeSoundingsWithIntensity(float intensity) return soundings; } +// The SAME lat/lon footprint as makeSoundings() (so the touchdown cells coincide +// cell-for-cell) but every sounding sits at a single, much deeper depth. Used to +// stage a false-deep "blunder" against a shallow predicted surface (#89). +std::vector makeDeepSoundings(float depth) +{ + std::vector soundings; + const double base_lat = 43.07; + const double base_lon = -70.76; + for (int i = 0; i < 5; ++i) { + gz4d::GeoPointLatLongDegrees point( + base_lat + i * 1e-5, base_lon + i * 1e-5, depth); + GeoSounding s(point); + s.sounding.vertical_error = 0.5f; + s.sounding.horizontal_error = 0.1f; + soundings.push_back(s); + } + return soundings; +} + constexpr int64_t kStamp = 1234567890123456789LL; constexpr uint16_t kSource = 7; constexpr float kIntensity = 42.5f; @@ -370,4 +389,129 @@ TEST(StoreImport, BackscatterNaNPropagation) << "soundings without intensity must surface no backscatter cells"; } +// Predicted-only prime (seed_settled=false, #89): seeds the predicted surface +// (turns the blunder gate on) but must NOT settle the cell -- no CUBE hypothesis, +// no values() output, sheet stays clean. Contrasted in the same test against the +// default seed_settled=true, which DOES settle the cell. +TEST(StoreImport, PredictedOnlyPrimeSeedsNoSettledHypothesis) +{ + // Build a tile with finite depths from the synthetic soundings. + GeoMapSheet source(1.0f); + source.addSoundings(makeSoundings()); + auto grids = source.grids(); + ASSERT_FALSE(grids.empty()); + + const std::vector values = grids.front()->values(); + const marine_bathymetry_store::BathymetryTile tile = + geoGridToTile(*grids.front(), kStamp, kSource); + + // Predicted-only prime: seeds the predicted surface but creates no hypothesis. + GeoMapSheet predicted_only(1.0f); + primeFromTile(tile, predicted_only, /*seed_settled=*/false); + + // Contrast: the default prime (seed_settled=true) DOES settle each cell. + GeoMapSheet settled(1.0f); + primeFromTile(tile, settled, /*seed_settled=*/true); + + auto po_grid = predicted_only.gridAt(tile.index()); + auto st_grid = settled.gridAt(tile.index()); + ASSERT_NE(po_grid, nullptr); + ASSERT_NE(st_grid, nullptr); + + // values() is a const method (it flushes the median pre-filter), so it is safe + // to call through the const grid handle. + const std::vector po_values = po_grid->values(); + const std::vector st_values = st_grid->values(); + + gggs::CellAreaIterator it(tile.index()); + std::size_t k = 0; + std::size_t checked = 0; + for (; it.valid() && k < values.size(); it.next(), ++k) { + if (std::isnan(values[k].depth)) { + continue; + } + ++checked; + // The predicted surface IS seeded for the predicted-only prime. + EXPECT_FLOAT_EQ(po_grid->predictedDepthAt(*it), values[k].depth); + // Predicted-only: NO settled hypothesis -> values() carries no estimate. + EXPECT_TRUE(std::isnan(po_values[k].depth)) + << "predicted-only prime must not settle a hypothesis"; + // Default prime: the SAME cell carries a settled estimate (the distinguishing + // contrast -- predicted-only vs predicted+settled). + EXPECT_FALSE(std::isnan(st_values[k].depth)) + << "default prime (seed_settled=true) must settle the cell"; + } + EXPECT_GT(checked, 0u) << "tile should carry some finite cells to prime"; + + // Priming reproduces persisted data -- it must never mark the sheet dirty. + EXPECT_TRUE(predicted_only.dirtyGrids().empty()); +} + +// The load-bearing behavior (#89): a seeded (shallow) predicted surface turns on +// CUBE's blunder gate so a false-deep sounding is REJECTED, whereas with no prior +// the same deep sounding is accepted. Demonstrates the bathy store is protected +// from false-deep detections by the Chart-prior prime. +TEST(StoreImport, SeededPredictedSurfaceRejectsDeepBlunder) +{ + // A shallow predicted surface (~ -10 m) from synthetic soundings -> the Chart- + // like prior tile we seed from. + GeoMapSheet shallow_src(1.0f); + shallow_src.addSoundings(makeSoundings()); + auto shallow_tiles = mapSheetToTiles(shallow_src, kStamp, kSource); + ASSERT_FALSE(shallow_tiles.empty()); + + // A clearly-too-deep blunder at the SAME locations (identical touchdown cells), + // far below any shallow-surface blunder limit. With target ~ -10 m and the + // defaults (blunder_minimum 10, blunder_percent 0.25, blunder_scalar 3), the + // limit is ~ -20 m, so -150 m is unambiguously a blunder. + constexpr float kDeep = -150.0f; + const std::vector deep = makeDeepSoundings(kDeep); + + // Baseline: WITHOUT a prior, predicted_depth_ stays INVALID_DATA, so Node::insert + // skips the blunder gate and the deep sounding is accepted -> finite deep cells. + GeoMapSheet no_prior(1.0f); + no_prior.addSoundings(deep); + std::size_t accepted_deep = 0; + for (const auto & grid : no_prior.grids()) { + for (const auto & v : grid->values()) { + if (!std::isnan(v.depth) && v.depth < -100.0f) { + ++accepted_deep; + } + } + } + EXPECT_GT(accepted_deep, 0u) + << "without a predicted surface the deep blunder must be accepted"; + + // WITH the shallow predicted surface primed (predicted-only): the blunder gate is + // active at every primed cell, so the deep sounding that touches it is rejected. + GeoMapSheet primed(1.0f); + for (const auto & gt : shallow_tiles) { + primeFromTile(gt.second, primed, /*seed_settled=*/false); + } + primed.addSoundings(deep); + + // Every cell that carries a (shallow) predicted depth must end with NO settled + // estimate: the deep sounding it received was blunder-rejected, and the + // predicted-only prime left no hypothesis of its own. + std::size_t gated_cells = 0; + for (const auto & gt : shallow_tiles) { + const auto & tile = gt.second; + auto grid = primed.gridAt(tile.index()); + ASSERT_NE(grid, nullptr); + const std::vector vals = grid->values(); + const std::vector & depth = tile.depthBand(); + gggs::CellAreaIterator it(tile.index()); + std::size_t k = 0; + for (; it.valid() && k < depth.size() && k < vals.size(); it.next(), ++k) { + if (std::isnan(depth[k])) { + continue; // not a primed cell + } + ++gated_cells; + EXPECT_TRUE(std::isnan(vals[k].depth)) + << "a deep blunder must be rejected where a shallow predicted surface gates"; + } + } + EXPECT_GT(gated_cells, 0u) << "the shallow tile should prime some cells to gate"; +} + } // namespace cube