Tier-2 backscatter correction: remove per-beam 2-way TL (depth-transferable angular response) (#87)

.agent/work-plans/issue-87/progress.md

@@ -0,0 +1,56 @@
+---
+issue: 87
+---
+
+# Issue #87 — Tier-2 TL-removed (depth-transferable) backscatter angular-response curve
+
+## Local Review (Pre-Push)
+**Status**: complete
+**When**: 2026-06-28 20:15 +0000
+**By**: Claude Code Agent (Claude Opus)
+**Verdict**: changes-requested
+
+**Branch**: feature/issue-87 at `24a0187`
+**Mode**: pre-push
+**Depth**: Deep (reason: substantive ADR-0007 addendum + cross-cutting estimator change threading a new field through a pack(1) binary-raster struct)
+**Must-fix**: 2 | **Suggestions**: 3
+**Round**: 1 | **Ship**: continue — one genuine, cross-pass-confirmed correctness concern (live-path R ≠ calibration R) warrants a fix or explicit bounding before the transferability claim holds.
+
+### Findings
+- [ ] (must-fix) Stale inverted-sign formula in the `range` field doc comment (`corrected = raw - (40*log10(R)+2*alpha*R) - residual`); validated code/ADR use `raw + TL(R) - residual`. Commit 24a0187 fixed the sign everywhere except this header — `include/cube_bathymetry/hypothesis.h:55`
+- [ ] (must-fix) Live-path `slant_range = sqrt(x²+y²+z²)` = `R·√(1+sin²tx·sin²rx)` ≠ offline/Python-calibration `R = twtt·c/2` when tx_angle≠0 (cross-pass confirmed, Lens A+B); silently biases the live correction and breaks live/offline equivalence. Comment "norm IS the slant range" is inaccurate. Reconcile R, or bound+document the limitation — `src/cube_bathymetry_node.cpp:973-975` (cf. `sounding.h:52-54`, `import_bag_main.cpp:604`)
+- [ ] (suggestion) `freqs` accumulated unconditionally even in tier-1 (unused); guard on `remove_tl` — `scripts/derive_angular_response.py:168-170`
+- [ ] (suggestion) Python drops non-positive/NaN-range beams from bins entirely while C++ retains them (TL-skipped); note the intentional asymmetry — `scripts/derive_angular_response.py:189-196` / `src/node.cpp:397`
+- [ ] (suggestion) `parse_args` missing a docstring (flake8-docstrings D103) — `scripts/derive_angular_response.py:58`
+
+### Verified correct (both adversarial passes, no action)
+- TL sign (added-back/TVG-style) consistent across Python, CSV, and `node.cpp:401`.
+- Francois-Garrison freshwater A3/P3 coefficients + T>20/≤20 split are standard; 500 kHz/24 °C → ~0.049 dB/m verified by hand. Units (Hz→kHz, dB/m) consistent.
+- Absorption (linear, per-bin via mean R) vs spreading (non-linear log, per-beam) aggregation is identical Python↔C++ by construction; α read verbatim in C++.
+- `DepthAndUncertainty` 16→20 byte growth is layout-safe (only `sizeof()` uses are bag_to_geotiff RasterIO strides that auto-track; all else by-name).
+- `range` threaded through every update/queueEstimate/recordBeam/queueFlush path; no silent drop.
+- Tier-1 backward compatibility preserved (`apply_tl = apply_ara && tl_removed`; defaulted params/args; old loader API retained). Header parser tolerant; tests cover edge cases.
+
+## Local Review (Pre-Push)
+**Status**: complete
+**When**: 2026-06-29 03:56 +0000
+**By**: Claude Code Agent (Claude Opus)
+**Verdict**: approved
+
+**Branch**: feature/issue-87 at `03cf523`
+**Mode**: pre-push
+**Depth**: Deep (reason: substantive ADR-0007 addendum + cross-cutting estimator change threading a new field through a pack(1) binary-raster struct)
+**Must-fix**: 0 | **Suggestions**: 2
+**Round**: 2 | **Ship**: recommended — both round-1 must-fixes resolved (sign comment fixed; live-R bounded+documented for the M3 tx=0 case, the resolution round-1 explicitly permitted); two independent adversarial lenses + static analysis surfaced 0 new must-fix.
+
+### Findings
+- [ ] (suggestion) Optional producer-frame assert on live-path `R = sqrt(x²+y²+z²)` — verified correct for the M3 (sensor-frame norm, tx=0) and honestly documented; an assert would harden against a future source emitting the cloud in a vehicle/base frame — `src/cube_bathymetry_node.cpp:977`
+- [ ] (suggestion) Optional one-shot warn when tier-2 is active but the range-present beam fraction is low; the documented retain-vs-drop asymmetry (C++ keeps NaN/non-positive-range beams as residual-only; Python drops them) is then load-bearing — `src/node.cpp:393-403`
+
+### Verified resolved (round-1 must-fixes)
+- Stale inverted-sign formula in `hypothesis.h` range-field doc now matches the validated `raw + (40*log10(R) + 2*alpha*R) - residual` model — `include/cube_bathymetry/hypothesis.h:54`.
+- Live-path R vs offline-R discrepancy: the inaccurate "norm IS the slant range" comment (present at `24a0187`) replaced with a bounded, M3-tx=0-honest derivation; offline import uses `twtt*c/2` so live/offline match for the M3 — `src/cube_bathymetry_node.cpp:969-977`.
+
+### Verified correct (no action, round 2)
+- Round-1 suggestions also addressed: `freqs` accumulation now guarded on `--remove-tl`; `parse_args` docstring added; retain-vs-drop asymmetry documented.
+- flake8 (ament profile) clean; `py_compile` OK. `pack(1)` 16→20 B growth layout-safe — `bag_to_geotiff.cpp:707` strides by `sizeof(DepthAndUncertainty)`; no whole-struct persistence / on-disk format break. No governance concerns (ADR-0007 addendum extended coherently; consequences handled; tests cover new tier-2 paths).

cube_bathymetry/docs/decisions/0007-mbes-backscatter-store-addendum-phase-b-transition.md

@@ -88,15 +88,71 @@ sonar. The `rx_angles` sign/zero convention was verified against
   symmetry) plus curve-loader/mode-parse tests; the 3 pre-existing default-path
   intensity assertions are unchanged (default `None` ⇒ `corrected == raw`).
 
+## Tier-2 addendum — TL-removed, depth-transferable curve (cube_bathymetry#87)
+
+The tier-1 curve above bins backscatter by `|rx_angle|` only, so per-beam range
+(spreading + absorption) is baked into the curve and it is **depth-regime
+specific** (must be re-derived per survey). Tier-2 removes the per-beam **2-way
+transmission loss** by **compensating (adding it back)** before the angular
+response is characterized and applied — a distant return lost more energy, so it
+is boosted to recover range-independent backscatter (TVG-style) — so the residual
+curve becomes **range/depth transferable**:
+
+`corrected = raw + TL(R) − residualCurve(|beam_angle|)`,
+`TL(R) = 40·log₁₀(R) + 2·α·R` (R = per-beam slant range `twtt·c/2`, m).
+
+> Sign note: TL is **added back**, not subtracted. An earlier draft subtracted it
+> (making far beams dimmer → steeper curve + larger cross-bag spread); the
+> validated fix compensates the loss, which flattens the curve and shrinks the
+> spread.
+
+Design decisions:
+
+1. **The curve file is self-describing — one correction mode, not two.** The
+   `Empirical` mode is unchanged; the estimator applies TL **iff the loaded curve
+   says so**. The curve CSV header gains `# tl_removed: true` and
+   `# absorption_db_per_m: <α>` (plus `# water_temp_c` / `# tl_model` provenance).
+   Tier-1 curves carry `tl_removed: false` (or omit the lines) and keep working
+   unchanged — fully backward compatible.
+2. **α lives only in Python.** The Francois-Garrison freshwater absorption is
+   computed once by `derive_angular_response.py` (`--remove-tl --water-temp-c`)
+   and written into the header as a scalar. The C++ estimator reads that scalar
+   verbatim and never recomputes α, so Python and C++ apply an **identical** TL by
+   construction (a divergence would silently corrupt the correction). At 500 kHz /
+   24 °C / fresh water α ≈ 0.049 dB/m.
+3. **Per-beam slant range R is a new sufficient statistic**, threaded exactly like
+   `beam_angle`: `Sounding::slant_range` → `DepthAndUncertainty::range` (the
+   pack(1) raster struct grows to 20 bytes; layout-safe because every raster read
+   strides by `sizeof(DepthAndUncertainty)`) → `BeamIntensitySample::range`,
+   recorded on the winning depth hypothesis. The offline importer threads R from
+   the `Sounding` detections ctor; the live node recovers R as the norm of the
+   sensor-frame `/soundings` point (no new cloud field). A NaN / non-positive R
+   skips the TL term (no log of a non-positive range).
+4. **Fresh water only.** The salinity (boric-acid + MgSO₄) seawater absorption
+   terms are out of scope for this issue; `--salinity > 0` is rejected by the tool.
+5. **Multi-bag derivation.** `derive_angular_response.py` now accepts multiple
+   bags (`nargs='+'`) and merges the per-bin sums — a real survey calibration
+   spans many bags.
+
+This still does **not** implement the full GeoCoder (insonified-area +
+beam-pattern), which is the deferred tier-3. The TVG/absorption/frequency state
+ultimately belongs in `SonarInfo` (unh_marine_autonomy#240); α is computed
+locally meanwhile.
+
 ## Deferred (explicit follow-ups)
 
-- **Full radiometric GeoCoder** — insonified-area, beam-pattern, TVG residual,
-  and the depth/slope incidence term (ADR-0007 D3 / cube_bathymetry#15 / #59).
-  The empirical curve absorbs the aggregate angular falloff for a flat bottom;
-  the slope-aware incidence correction remains future work, gated on the
+- **Full radiometric GeoCoder (tier-3)** — insonified-area, beam-pattern, TVG
+  residual, and the depth/slope incidence term (ADR-0007 D3 / cube_bathymetry#15
+  / #59). The empirical curve (tier-1) absorbs the aggregate angular falloff for a
+  flat bottom; tier-2 makes it range/depth transferable; the slope-aware incidence
+  correction and the area/beam-pattern terms remain future work, gated on the
   predicted-surface producer (#59). **Follow-up to file:** author the full,
   canonical `docs/decisions/0007-mbes-backscatter-store.md` document (this
   addendum folds into it).
+- **TVG/absorption/frequency state in SonarInfo** (unh_marine_autonomy#240): the
+  M3's already-applied TVG is the one genuine remaining unknown; the tier-1-vs-2
+  comparison resolves it empirically. Until then α is computed locally from the
+  per-ping frequency + water temperature.
 - **Intensity domain (dB vs linear) as a first-class sonar property.** The
   correction assumes dB (a subtraction). Sonars reporting linear intensity would
   need a divide, or a domain conversion at ingest. **Follow-up to file:** model

cube_bathymetry/include/cube_bathymetry/angular_response_curve.h

@@ -37,15 +37,41 @@ namespace cube
   bool parseBackscatterAngleCorrection(
     const std::string & text, BackscatterAngleCorrection & out);
 
+/// A loaded angular-response curve plus its self-describing TL provenance
+/// (cube_bathymetry#87). `points` is the ascending {abs_angle_deg,
+/// db_relative_to_nadir} curve. When `tl_removed` is true the curve is a TL-
+/// REMOVED residual (tier-2) and the estimator must remove the per-beam 2-way
+/// transmission loss `40*log10(R) + 2*alpha*R` (alpha == `absorption_db_per_m`)
+/// before subtracting the residual. A tier-1 curve has `tl_removed == false`
+/// and `absorption_db_per_m == 0`, fully backward compatible.
+  struct AngularResponseCurve
+  {
+    std::vector < std::pair < float, float >> points;
+    bool tl_removed = false;
+    float absorption_db_per_m = 0.0f;
+  };
+
 /// Load an empirical angular-response curve from a CSV file matching the seed at
 /// ~/data/logs/analysis/m3_angular_response_curve.csv: a header line then rows
 /// `abs_angle_deg_center,mean_bs_db,n,db_relative_to_nadir`. Reads the first
 /// (abs_angle_deg_center) and fourth (db_relative_to_nadir) columns into ascending
 /// {abs_angle_deg, db_relative_to_nadir} pairs (sorted by angle).
 ///
+/// Also parses the optional self-describing TL header comments written by the
+/// derive tool (cube_bathymetry#87):
+///   `# tl_removed: true|false`
+///   `# absorption_db_per_m: <float>`
+/// Their absence yields tier-1 defaults (tl_removed=false, absorption=0), so a
+/// tier-1 curve keeps loading unchanged. The scalar `absorption_db_per_m` is read
+/// verbatim (the C++ estimator never recomputes alpha).
+///
 /// Comment lines (`#`...), blank lines, the header row, and rows that fail to
-/// parse are skipped. A missing/unreadable file yields an empty vector (the
+/// parse are skipped. A missing/unreadable file yields an empty curve (the
 /// Empirical correction then degrades to a no-op -- callers should warn).
+  AngularResponseCurve loadAngularResponseCurveWithHeader(const std::string & path);
+
+/// Backward-compatible convenience: the curve `points` only (drops the TL
+/// provenance). Equivalent to `loadAngularResponseCurveWithHeader(path).points`.
   std::vector < std::pair < float, float >> loadAngularResponseCurve(const std::string & path);
 
 }  // namespace cube

cube_bathymetry/include/cube_bathymetry/common.h

@@ -85,17 +85,28 @@ namespace cube
 
   /// Per-beam receive/steering angle (radians) accompanying the intensity, the
   /// {raw intensity, angle} sufficient-statistics pair (ADR-0007 D3). NaN when
-  /// not reported. NOTE: extending this #pragma pack(push,1) struct from 8 to
-  /// 16 bytes is layout-safe -- no caller depends on sizeof(DepthAndUncertainty)
-  /// (the only sizeof uses are raster-band strides, 2*sizeof(float)).
+  /// not reported.
     float beam_angle;
 
+  /// Per-beam slant range R (meters) accompanying the intensity, carried through
+  /// the median queue bound to its depth so the tier-2 backscatter TL correction
+  /// (cube_bathymetry#87) reads the same beam's range. NaN when not reported.
+  ///
+  /// LAYOUT NOTE: this #pragma pack(push,1) struct is now 20 bytes (5 floats).
+  /// The only sizeof uses are raster-band strides that read the first two floats
+  /// (depth, uncertainty) as bands using sizeof(DepthAndUncertainty) as the
+  /// element stride (e.g. bag_to_geotiff.cpp RasterIO). Appending floats at the
+  /// END is layout-safe ONLY because every raster read uses
+  /// sizeof(DepthAndUncertainty) (NOT a hardcoded 8/16); verified for #87.
+    float range;
+
     DepthAndUncertainty(float depth = std::numeric_limits < float > ::quiet_NaN(),
       float uncertainty = std::numeric_limits < float > ::quiet_NaN(),
       float intensity = std::numeric_limits < float > ::quiet_NaN(),
-      float beam_angle = std::numeric_limits < float > ::quiet_NaN())
+      float beam_angle = std::numeric_limits < float > ::quiet_NaN(),
+      float range = std::numeric_limits < float > ::quiet_NaN())
       : depth(depth), uncertainty(uncertainty), intensity(intensity),
-      beam_angle(beam_angle) {
+      beam_angle(beam_angle), range(range) {
     }
   };
 #pragma pack(pop)

cube_bathymetry/include/cube_bathymetry/geo_map_sheet.h

@@ -57,9 +57,18 @@ namespace cube
   /// every owned GeoGrid holds by const reference -- so the setting reaches all
   /// existing and future grids. Call AFTER construction and BEFORE/at processing.
   /// An Empirical mode with an empty curve is a no-op (the caller should warn).
+  ///
+  /// @p tl_removed / @p absorption_db_per_m carry the tier-2 TL provenance
+  /// (cube_bathymetry#87): when @p tl_removed is true the @p curve is a TL-
+  /// removed residual and the estimator removes `40*log10(R) + 2*alpha*R` per
+  /// beam (alpha == @p absorption_db_per_m) before subtracting the residual.
+  /// Both default to the tier-1 values (false / 0), so existing callers are
+  /// unchanged.
     void setBackscatterCorrection(
       BackscatterAngleCorrection mode,
-      std::vector < std::pair < float, float >> curve);
+      std::vector < std::pair < float, float >> curve,
+      bool tl_removed = false,
+      float absorption_db_per_m = 0.0f);
 
   /// Return the grids within the bounds, creating new ones if necessary
     std::vector < std::shared_ptr <

cube_bathymetry/include/cube_bathymetry/hypothesis.h

@@ -24,6 +24,7 @@
 #define CUBE_BATHYMETRY__HYPOTHESIS_H_
 
 #include <cstdint>
+#include <limits>
 #include <memory>
 #include <vector>
 #include "cube_bathymetry/parameters.h"
@@ -48,6 +49,15 @@ namespace cube
   /// treats a NaN angle as "no angle correction available" and emits that beam
   /// uncorrected.
     float beam_angle;
+
+  /// Per-beam slant range R from the sonar head to the touchdown, in meters,
+  /// for the tier-2 backscatter 2-way transmission-loss correction
+  /// (cube_bathymetry#87): `corrected = raw + (40*log10(R) + 2*alpha*R) -
+  /// residualCurve(|angle|)` (the TL is ADDED BACK to compensate the loss). NaN
+  /// (or non-positive) when not reported; the TL
+  /// term is then skipped (identity) so the beam is corrected by the residual
+  /// angular-response curve alone (tier-1 behavior).
+    float range = std::numeric_limits < float > ::quiet_NaN();
   };
 
 /// Depth hypothesis structure used to maintain a current track on the depth
@@ -115,8 +125,12 @@ namespace cube
   ///
   /// A NaN raw_intensity is skipped (a source that omits intensities must never
   /// inject a phantom sample); a NaN beam_angle is retained (the beam is
-  /// still a valid intensity sample, merely uncorrectable for angle).
-    void recordBeam(float raw_intensity, float beam_angle);
+  /// still a valid intensity sample, merely uncorrectable for angle). A NaN /
+  /// non-positive range is also retained (the beam is a valid intensity sample,
+  /// merely uncorrectable for the tier-2 TL term, cube_bathymetry#87).
+    void recordBeam(
+      float raw_intensity, float beam_angle,
+      float range = std::numeric_limits < float > ::quiet_NaN());
 
   /// Current depth mean estimate
     double current_estimate;

cube_bathymetry/include/cube_bathymetry/node.h

@@ -113,9 +113,12 @@ namespace cube
   ///     backscatter association tracks the depth association (ADR-0007 D2/D3).
   ///   beam_angle: per-beam receive/steering angle (radians, NaN when absent),
   ///     the angle half of the {raw intensity, angle} sufficient-stats pair.
+  ///   range: per-beam slant range (m, NaN when absent), for the tier-2 TL
+  ///     correction (cube_bathymetry#87).
     bool update(
       float depth, float variance, const Parameters & parameters,
-      float intensity = std::nan(""), float beam_angle = std::nan(""));
+      float intensity = std::nan(""), float beam_angle = std::nan(""),
+      float range = std::nan(""));
 
   /// Find the closest matching hypothesis in the current linked list.
   /// This computes the normalised absolute error between one-step
@@ -164,9 +167,12 @@ namespace cube
   ///   intensity/beam_angle: the {raw intensity, angle} pair for this beam,
   ///     carried on the queue entry bound to its depth so the median sort never
   ///     mismatches a depth with a foreign intensity (ADR-0007 D3).
+  ///   range: per-beam slant range (m, NaN when absent), carried alongside so the
+  ///     tier-2 TL correction reads the same beam's range (cube_bathymetry#87).
     bool queueEstimate(
       float depth, float variance, const Parameters & parameters,
-      float intensity = std::nan(""), float beam_angle = std::nan(""));
+      float intensity = std::nan(""), float beam_angle = std::nan(""),
+      float range = std::nan(""));
 
   /* Routine: cube_node_extract_depth_unct
   * Purpose:  Extract depth and uncertainty of current best estimate

cube_bathymetry/include/cube_bathymetry/parameters.h

@@ -190,6 +190,24 @@ namespace cube
   /// CSV. Empty -> the Empirical correction is a no-op (logged as a warning at
   /// configure time). Used only when backscatter_angle_correction == Empirical.
     std::vector < std::pair < float, float >> angular_response_curve;
+
+  /// Tier-2 backscatter correction (cube_bathymetry#87): when true, the loaded
+  /// angular_response_curve is a TL-REMOVED residual, so the estimator removes
+  /// the per-beam 2-way transmission loss `40*log10(R) + 2*alpha*R` (R = per-beam
+  /// slant range, m) BEFORE subtracting the residual curve, making the curve
+  /// depth/range transferable. Set by the curve loader from the CSV header
+  /// (`# tl_removed: <bool>`). Default false = tier-1 (no TL term), fully
+  /// backward compatible. Used only when backscatter_angle_correction ==
+  /// Empirical.
+    bool backscatter_tl_removed = false;
+
+  /// Tier-2 absorption coefficient alpha in dB/m (cube_bathymetry#87). The scalar
+  /// the estimator multiplies into the TL term `2*alpha*R`; read verbatim from
+  /// the curve CSV header (`# absorption_db_per_m: <float>`) so the C++ estimator
+  /// NEVER recomputes the (Francois-Garrison) absorption -- the Python derive tool
+  /// is the single source of truth, guaranteeing Python/C++ consistency. Default
+  /// 0 = no absorption term (tier-1, or fresh water at negligible alpha).
+    float backscatter_absorption_db_per_m = 0.0f;
   };
 
 }  // namespace cube