From 39b2ff99d58a34576ae6208d02ab1591cd2bc198 Mon Sep 17 00:00:00 2001 From: Guriy Date: Wed, 5 Aug 2026 19:01:26 +0300 Subject: [PATCH] surface_align added --- vna_system/core/processing/__init__.py | 25 + vna_system/core/processing/surface_align.py | 316 ++++++++++++ .../core/processing/surface_alignment.py | 479 ++++++++++++++++++ .../core/processors/configs/bscan_config.json | 13 + .../implementations/bscan_processor.py | 448 +++++++++++++++- 5 files changed, 1269 insertions(+), 12 deletions(-) create mode 100644 vna_system/core/processing/__init__.py create mode 100644 vna_system/core/processing/surface_align.py create mode 100644 vna_system/core/processing/surface_alignment.py diff --git a/vna_system/core/processing/__init__.py b/vna_system/core/processing/__init__.py new file mode 100644 index 0000000..97faaf0 --- /dev/null +++ b/vna_system/core/processing/__init__.py @@ -0,0 +1,25 @@ +"""Signal-processing algorithms shared by the VNA processors.""" + +from vna_system.core.processing.surface_align import ( + AlignConfig, + AlignResult, + apply_alignment, + check_geometry, + estimate_surface, +) +from vna_system.core.processing.surface_alignment import ( + AlignmentState, + AlignOptions, + SurfaceAlignmentEngine, +) + +__all__ = [ + "AlignConfig", + "AlignResult", + "AlignOptions", + "AlignmentState", + "SurfaceAlignmentEngine", + "apply_alignment", + "check_geometry", + "estimate_surface", +] diff --git a/vna_system/core/processing/surface_align.py b/vna_system/core/processing/surface_align.py new file mode 100644 index 0000000..a5de4d7 --- /dev/null +++ b/vna_system/core/processing/surface_align.py @@ -0,0 +1,316 @@ +#!/usr/bin/env python3 +"""Ground-surface alignment for a multi-channel air-coupled SFCW GPR. + +Vendored verbatim from radar_system/python_app/processing/surface_align.py. +Keep it dependency-free (numpy only) and keep the public contract stable -- +`AlignConfig`, `estimate_surface`, `apply_alignment`, `check_geometry` -- so the +two copies stay interchangeable. The application-side glue lives in +`surface_alignment.py`; put integration logic there, not here. + +The estimator does not track each channel separately. It solves for the single +physical quantity that explains every channel at once -- the antenna height h -- +using the known bistatic geometry: + + path_c(h) = sqrt(L_c^2 + (2h)^2) L_c = transmitter-receiver baseline + +A wrong pick in one channel does not lie on the other channels' hyperbolae, so +the joint score suppresses it instead of following it. + +Two smoothing modes: + * "offline" zero-phase low-pass (replay: whole profile available) + * "live" constant-velocity Kalman (only the past is available) + +Background suppression is the operator's choice: a recorded reference, the +global median trace, or a sliding-window median. +""" +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Literal + +import numpy as np + +C_LIGHT = 299_792_458.0 + + +# --------------------------------------------------------------------------- # +# configuration +# --------------------------------------------------------------------------- # +@dataclass +class AlignConfig: + band_hz: tuple[float, float] = (0.3e9, 5.0e9) + pad_factor: int = 16 # zero-padding -> finer peak grid + height_range_m: tuple[float, float] = (0.15, 0.90) + height_step_m: float = 0.001 + track_halfwidth_m: float = 0.09 # search window around the prediction + # first-break: take the first score peak reaching this fraction of the best + first_break_frac: float = 0.45 + min_confidence: float = 2.0 # peak / median of the score curve + background: Literal["reference", "global_median", "sliding_median"] = "global_median" + sliding_window: int = 101 # traces, for background="sliding_median" + smoothing: Literal["offline", "live"] = "offline" + lowpass_hz: float = 1.5 # offline mode + process_noise_m: float = 0.02 # live mode: how fast h may change + measurement_noise_m: float = 0.005 + mad_floor_m: float = 0.002 + outlier_sigmas: float = 4.0 + + +@dataclass +class AlignResult: + height_m: np.ndarray # per cycle, smoothed + height_raw_m: np.ndarray # per cycle, before smoothing + confidence: np.ndarray + accepted: np.ndarray # bool, survived outlier rejection + baselines_m: dict # combo -> L + quality: dict = field(default_factory=dict) + + +# --------------------------------------------------------------------------- # +# helpers +# --------------------------------------------------------------------------- # +def _to_time_domain(spectra: np.ndarray, freq: np.ndarray, pad: int): + """Windowed IFFT -> analytic traces on a zero-padded time grid.""" + n = freq.size + N = n * pad + df = float(freq[1] - freq[0]) + z = np.fft.ifft(spectra * np.hanning(n), n=N, axis=-1) * N + t = np.arange(N) / (N * df) + return z, t + + +def _suppress_background(z: np.ndarray, cfg: AlignConfig, + reference: np.ndarray | None) -> np.ndarray: + """Remove whatever is stationary, by the operator's chosen means.""" + if cfg.background == "reference": + if reference is None: + raise ValueError("background='reference' requires a reference trace") + return z - reference[None, :] + if cfg.background == "global_median": + bg = np.median(z.real, axis=0) + 1j * np.median(z.imag, axis=0) + return z - bg[None, :] + if cfg.background == "sliding_median": + w = max(3, int(cfg.sliding_window) | 1) + half = w // 2 + out = np.empty_like(z) + for i in range(z.shape[0]): + lo, hi = max(0, i - half), min(z.shape[0], i + half + 1) + bg = (np.median(z[lo:hi].real, axis=0) + + 1j * np.median(z[lo:hi].imag, axis=0)) + out[i] = z[i] - bg + return out + raise ValueError(f"unknown background mode {cfg.background!r}") + + +def estimate_baselines(z_by_combo: dict, t: np.ndarray, + search_m: tuple[float, float] = (0.10, 1.20)) -> dict: + """Baseline L from the stationary echo: it sits at one-way delay L/c. + + Returns {combo: L}. Verify with `check_geometry` before trusting it -- if + the stationary echo is a frame reflection rather than the direct path, the + per-channel heights will disagree. + """ + out = {} + for combo, z in z_by_combo.items(): + fixed = np.median(z.real, axis=0) + 1j * np.median(z.imag, axis=0) + d = C_LIGHT * t # one-way path length + m = (d >= search_m[0]) & (d <= search_m[1]) + idx = np.where(m)[0] + k = idx[int(np.argmax(np.abs(fixed[idx])))] + out[combo] = float(d[k]) + return out + + +def _parabolic(y0: float, y1: float, y2: float) -> float: + den = y0 - 2.0 * y1 + y2 + if den == 0.0: + return 0.0 + return float(np.clip(0.5 * (y0 - y2) / den, -1.0, 1.0)) + + +# --------------------------------------------------------------------------- # +# the joint estimator +# --------------------------------------------------------------------------- # +def _score_curve(env_by_combo: dict, baselines: dict, t: np.ndarray, + heights: np.ndarray, cycle: int) -> np.ndarray: + """Sum the channels' envelopes along the bistatic hyperbola for each h.""" + dt = t[1] - t[0] + score = np.zeros_like(heights) + for combo, env in env_by_combo.items(): + L = baselines[combo] + tau = np.sqrt(L * L + 4.0 * heights * heights) / C_LIGHT + idx = tau / dt + i0 = np.floor(idx).astype(int) + frac = idx - i0 + ok = (i0 >= 0) & (i0 + 1 < env.shape[1]) + row = env[cycle] + vals = np.zeros_like(heights) + vals[ok] = row[i0[ok]] * (1 - frac[ok]) + row[i0[ok] + 1] * frac[ok] + # normalise per channel so a loud channel cannot dominate the vote + scale = np.median(row) + 1e-30 + score += vals / scale + return score + + +def _pick_from_score(score: np.ndarray, heights: np.ndarray, + lo: int, hi: int, cfg: AlignConfig) -> tuple[float, float]: + """First significant peak inside [lo, hi) -- nothing subsurface can be + shallower than the surface, so the *first* arrival is the right target.""" + seg = score[lo:hi] + if seg.size < 3: + return float(heights[(lo + hi) // 2]), 0.0 + thr = cfg.first_break_frac * seg.max() + k = None + for i in range(1, seg.size - 1): + if seg[i] >= thr and seg[i] >= seg[i - 1] and seg[i] >= seg[i + 1]: + k = i + break + if k is None: + k = int(np.argmax(seg)) + delta = _parabolic(seg[k - 1], seg[k], seg[k + 1]) if 0 < k < seg.size - 1 else 0.0 + step = heights[1] - heights[0] + h = float(heights[lo + k] + delta * step) + conf = float(seg[k] / (np.median(score) + 1e-30)) + return h, conf + + +def _kalman(meas: np.ndarray, conf: np.ndarray, dt: float, + cfg: AlignConfig) -> np.ndarray: + """Constant-velocity Kalman. Causal: the estimate for cycle i uses only + cycles <= i, so it is usable live, and unlike a causal low-pass it does not + lag (the velocity state extrapolates).""" + x = np.array([meas[0], 0.0]) + P = np.diag([1e-2, 1e-2]) + F = np.array([[1.0, dt], [0.0, 1.0]]) + q = cfg.process_noise_m ** 2 + Q = q * np.array([[dt ** 3 / 3.0, dt ** 2 / 2.0], [dt ** 2 / 2.0, dt]]) + H = np.array([[1.0, 0.0]]) + out = np.empty_like(meas) + for i, (zm, cf) in enumerate(zip(meas, conf)): + x = F @ x + P = F @ P @ F.T + Q + if cf >= cfg.min_confidence: + r = cfg.measurement_noise_m ** 2 * max(1.0, cfg.min_confidence / cf) + S = (H @ P @ H.T).item() + r + K = (P @ H.T / S).ravel() + x = x + K * (zm - (H @ x).item()) + P = (np.eye(2) - np.outer(K, H)) @ P + out[i] = x[0] + return out + + +def _zero_phase_lowpass(y: np.ndarray, fs: float, cut: float) -> np.ndarray: + sp = np.fft.rfft(y - y.mean()) + fr = np.fft.rfftfreq(y.size, 1.0 / fs) + sp[fr > cut] = 0.0 + return np.fft.irfft(sp, n=y.size) + y.mean() + + +def estimate_surface(spectra_by_combo: dict, freq: np.ndarray, + times_s: np.ndarray, cfg: AlignConfig, + baselines: dict | None = None, + references: dict | None = None) -> AlignResult: + """Estimate the antenna height for every cycle, jointly over all channels.""" + z_by_combo, env_by_combo = {}, {} + t = None + for combo, S in spectra_by_combo.items(): + z, t = _to_time_domain(S, freq, cfg.pad_factor) + z_by_combo[combo] = z + if baselines is None: + baselines = estimate_baselines(z_by_combo, t) + for combo, z in z_by_combo.items(): + ref = None if references is None else references.get(combo) + env_by_combo[combo] = np.abs(_suppress_background(z, cfg, ref)) + + heights = np.arange(cfg.height_range_m[0], cfg.height_range_m[1], + cfg.height_step_m) + n_cycles = next(iter(env_by_combo.values())).shape[0] + raw = np.zeros(n_cycles) + conf = np.zeros(n_cycles) + half = int(round(cfg.track_halfwidth_m / cfg.height_step_m)) + + prev = None + for i in range(n_cycles): + score = _score_curve(env_by_combo, baselines, t, heights, i) + if prev is None: + lo, hi = 0, heights.size + else: + c = int(round((prev - heights[0]) / cfg.height_step_m)) + lo, hi = max(0, c - half), min(heights.size, c + half + 1) + h, cf = _pick_from_score(score, heights, lo, hi, cfg) + raw[i], conf[i] = h, cf + if cf >= cfg.min_confidence: + prev = h + elif prev is None: + prev = h + + # robust outlier rejection against a running median + med = np.array([np.median(raw[max(0, i - 4):i + 5]) for i in range(n_cycles)]) + resid = raw - med + mad = max(1.4826 * float(np.median(np.abs(resid - np.median(resid)))), + cfg.mad_floor_m) + accepted = (np.abs(resid) < cfg.outlier_sigmas * mad) & (conf >= cfg.min_confidence) + if accepted.sum() < 3: + accepted = np.ones(n_cycles, dtype=bool) + + filled = np.interp(times_s, times_s[accepted], raw[accepted]) + dt = float(np.median(np.diff(times_s))) + if cfg.smoothing == "offline": + height = _zero_phase_lowpass(filled, 1.0 / dt, cfg.lowpass_hz) + else: + height = _kalman(filled, conf, dt, cfg) + + quality = { + "rejected_fraction": float(1.0 - accepted.mean()), + "confidence_median": float(np.median(conf)), + "mad_m": float(mad), + "height_p2p_m": float(height.max() - height.min()), + "trace_rate_hz": 1.0 / dt, + } + return AlignResult(height, raw, conf, accepted, baselines, quality) + + +def check_geometry(spectra_by_combo: dict, freq: np.ndarray, cfg: AlignConfig, + baselines: dict, height_m: np.ndarray) -> dict: + """Per-channel height implied by each channel's own picked delay. + + All channels see the same platform, so these must agree. Disagreement means + a channel is tracking something that is not the surface -- exactly the + failure a single-channel tracker cannot detect. + """ + out = {} + for combo, S in spectra_by_combo.items(): + z, t = _to_time_domain(S, freq, cfg.pad_factor) + env = np.abs(_suppress_background(z, cfg, None)) + L = baselines[combo] + dt = t[1] - t[0] + hs = [] + for i in range(env.shape[0]): + tau_pred = np.sqrt(L * L + 4 * height_m[i] ** 2) / C_LIGHT + k0 = int(round(tau_pred / dt)) + w = int(round(0.15 / C_LIGHT / dt)) # +-15 cm of path + lo, hi = max(0, k0 - w), min(env.shape[1], k0 + w + 1) + k = lo + int(np.argmax(env[i, lo:hi])) + p = (C_LIGHT * t[k]) ** 2 - L * L + hs.append(np.sqrt(max(p, 0.0)) / 2.0) + out[combo] = np.array(hs) + return out + + +def apply_alignment(spectra_by_combo: dict, freq: np.ndarray, + baselines: dict, height_m: np.ndarray, + height_ref_m: float | None = None) -> dict: + """Shift every channel so the surface lands at a common reference height. + + The shift is a linear phase ramp, which is an *exact* fractional delay: the + magnitude of exp(j2*pi*f*tau) is one, so nothing is filtered or interpolated. + """ + if height_ref_m is None: + height_ref_m = float(np.median(height_m)) + out = {} + for combo, S in spectra_by_combo.items(): + L = baselines[combo] + tau = np.sqrt(L * L + 4.0 * height_m ** 2) / C_LIGHT + tau_ref = np.sqrt(L * L + 4.0 * height_ref_m ** 2) / C_LIGHT + out[combo] = S * np.exp(1j * 2 * np.pi * np.outer(tau - tau_ref, freq)) + return out diff --git a/vna_system/core/processing/surface_alignment.py b/vna_system/core/processing/surface_alignment.py new file mode 100644 index 0000000..1010363 --- /dev/null +++ b/vna_system/core/processing/surface_alignment.py @@ -0,0 +1,479 @@ +"""Glue between the sweep-history model of the processors and `surface_align`. + +`surface_align` works on a whole profile at once: it wants every cycle's +spectrum, a uniform frequency grid and the cycle timestamps. A live VNA hands us +one sweep at a time, and the B-scan can be rebuilt from scratch at any moment +when the operator changes a setting. This module reconciles the two: + +* `plan()` decides whether the cached estimate is still usable, whether it can be + extended with a trailing window (the live case), or whether everything has to + be recomputed (config change, history reload). +* `update()` runs the estimator for the chosen mode and keeps one height per + history index, so a processor can ask "what was the height for cycle i?". +* `shift()` applies the alignment to a single spectrum over the *full* measured + band, even though the height was estimated on a narrower band. + +The estimator itself is untouched -- all physics lives in `surface_align`. +""" +from __future__ import annotations + +from dataclasses import dataclass, field, replace +from typing import Any, Literal + +import numpy as np +from numpy.typing import NDArray + +from vna_system.core.logging.logger import get_component_logger +from vna_system.core.processing.surface_align import ( + AlignConfig, + apply_alignment, + check_geometry, + estimate_surface, +) +# The estimator's own windowed IFFT. Reused rather than reimplemented so a +# recorded reference lands on exactly the time grid the estimator works on. +from vna_system.core.processing.surface_align import _to_time_domain + +logger = get_component_logger(__file__) + +PlanMode = Literal["reuse", "extend", "full"] + +# A single VNA measures one TX/RX pair per sweep, so there is exactly one +# channel. The estimator is dict-based, so adding channels later is a matter of +# passing more keys. +DEFAULT_CHANNEL = "vna" + + +@dataclass(frozen=True) +class AlignOptions: + """Everything the caller can tune, in SI units. + + Frozen so it can be compared as a whole: any change invalidates the cached + estimate, which is exactly the desired behaviour. + """ + band_hz: tuple[float, float] = (0.3e9, 5.0e9) + height_range_m: tuple[float, float] = (0.25, 0.90) + height_step_m: float = 0.001 + background: str = "reference" + smoothing: str = "live" + pad_factor: int = 8 + baseline_m: float | None = None # None -> from the stationary echo + height_ref_m: float | None = None # None or 0 -> median of the first estimate + geometry_check: bool = True + + +@dataclass +class AlignmentState: + """One height per history index, plus the numbers worth watching.""" + heights_m: NDArray[np.floating] + raw_heights_m: NDArray[np.floating] + confidence: NDArray[np.floating] + baseline_m: float + height_ref_m: float + quality: dict[str, float] = field(default_factory=dict) + geometry_residual_m: float | None = None + geometry_spread_m: float | None = None + # True once heights were appended from a trailing window rather than + # estimated over the whole profile in one pass. + incremental: bool = False + + @property + def covered(self) -> int: + """Number of cycles the state has a height for.""" + return int(self.heights_m.size) + + +class SurfaceAlignmentEngine: + """Cached, incrementally extendable surface-height estimate.""" + + # Median-based background needs several traces before it means anything. + MIN_CYCLES = 8 + # With a recorded reference the background is known up front, but the + # smoothers still need a trace interval, i.e. more than one cycle. + MIN_REFERENCE_CYCLES = 3 + # Trailing window used to extend the estimate live. + WINDOW_CYCLES = 128 + # More new cycles than this and a full recompute is cheaper/safer. + MAX_EXTEND = 32 + # Below this many frequency points the band is too narrow to pick a peak. + MIN_BAND_POINTS = 16 + + def __init__(self, channel: str = DEFAULT_CHANNEL) -> None: + self.channel = channel + self._state: AlignmentState | None = None + self._key: tuple[Any, ...] | None = None + # The reference fallback is decided per estimate, i.e. once per sweep + # while streaming; warn about it once instead of on every sweep. + self._warned_reference_missing = False + + # ------------------------------------------------------------------ # + # Cache control + # ------------------------------------------------------------------ # + @property + def state(self) -> AlignmentState | None: + return self._state + + def invalidate(self) -> None: + """Forget the estimate; the next `plan()` asks for a full recompute.""" + self._state = None + self._key = None + self._warned_reference_missing = False + + def refresh_if_incremental(self, smoothing: str) -> None: + """Drop an incrementally built estimate before a full rebuild. + + Live sweeps can only be smoothed causally, so an estimate grown one + window at a time is not what "offline" smoothing promises: that mode + low-passes the whole profile at once. A rebuild is the moment to deliver + it. Cheap modes are left alone -- the Kalman is causal by design, so + re-running it over everything would change little. + """ + if smoothing == "offline" and self._state is not None and self._state.incremental: + logger.debug("Dropping incremental surface estimate for an offline rebuild") + self.invalidate() + + def drop_leading(self, count: int) -> None: + """Drop heights for the `count` oldest cycles (history was trimmed).""" + if count <= 0 or self._state is None: + return + state = self._state + if count >= state.covered: + self.invalidate() + return + self._state = replace( + state, + heights_m=state.heights_m[count:], + raw_heights_m=state.raw_heights_m[count:], + confidence=state.confidence[count:], + ) + + # ------------------------------------------------------------------ # + # Planning + # ------------------------------------------------------------------ # + @staticmethod + def freq_key(freq_hz: NDArray[np.floating]) -> tuple[float, float, int]: + """Identity of a frequency grid: a different grid invalidates the cache.""" + return float(freq_hz[0]), float(freq_hz[-1]), int(freq_hz.size) + + def plan( + self, + n_cycles: int, + options: AlignOptions, + freq_key: tuple[float, float, int], + ) -> tuple[PlanMode, int]: + """Decide what has to be computed for cycles `[0, n_cycles)`. + + Returns + ------- + (mode, start_index) + `start_index` is the first cycle whose spectrum the caller must + supply to `update()`. + """ + key = (options, freq_key, self.channel) + if key != self._key or self._state is None: + return "full", 0 + + covered = self._state.covered + if covered == n_cycles: + return "reuse", 0 + if covered < n_cycles <= covered + self.MAX_EXTEND: + return "extend", max(0, n_cycles - self.WINDOW_CYCLES) + return "full", 0 + + # ------------------------------------------------------------------ # + # Estimation + # ------------------------------------------------------------------ # + def update( + self, + mode: PlanMode, + start_index: int, + spectra: NDArray[np.complexfloating], + freq_hz: NDArray[np.floating], + times_s: NDArray[np.floating], + options: AlignOptions, + reference: NDArray[np.complexfloating] | None = None, + ) -> AlignmentState | None: + """Run the estimator for `mode` and return the updated state. + + `spectra` covers cycles `[start_index, n_cycles)` over the full measured + band; `freq_hz` is that band's grid (uniform, Hz). Returns None when the + data cannot support an estimate yet, in which case the caller should + leave the sweeps unshifted. + """ + if mode == "reuse" and self._state is not None: + return self._state + + band = (freq_hz >= options.band_hz[0]) & (freq_hz <= options.band_hz[1]) + if int(band.sum()) < self.MIN_BAND_POINTS: + logger.warning( + "Surface alignment band too narrow; skipping", + points=int(band.sum()), + band_hz=options.band_hz, + ) + return None + + fb = np.asarray(freq_hz, dtype=float)[band] + window = np.ascontiguousarray(np.asarray(spectra, dtype=np.complex128)[:, band]) + times = self._sanitize_times(np.asarray(times_s, dtype=float)) + n_window = window.shape[0] + n_cycles = start_index + n_window + + cfg, refs = self._build_config(options, fb, times, reference, band, n_cycles) + if cfg is None: + return None + + if mode == "extend" and self._state is not None: + return self._extend(cfg, refs, fb, window, times, n_cycles) + + state = self._full(cfg, refs, fb, window, times, options) + if state is not None: + self._key = (options, self.freq_key(freq_hz), self.channel) + self._state = state + return state + + def _build_config( + self, + options: AlignOptions, + fb: NDArray[np.floating], + times: NDArray[np.floating], + reference: NDArray[np.complexfloating] | None, + band: NDArray[np.bool_], + n_cycles: int, + ) -> tuple[AlignConfig | None, dict[str, NDArray[np.complexfloating]] | None]: + """Translate `AlignOptions` into an `AlignConfig` plus reference traces.""" + background = options.background + refs: dict[str, NDArray[np.complexfloating]] | None = None + + if background == "reference": + if reference is None: + if not self._warned_reference_missing: + logger.warning( + "Alignment background='reference' but no matching reference sweep; " + "falling back to the global median" + ) + self._warned_reference_missing = True + background = "global_median" + else: + ref_band = np.asarray(reference, dtype=np.complex128)[band] + z_ref, _ = _to_time_domain(ref_band[None, :], fb, options.pad_factor) + refs = {self.channel: z_ref[0]} + + # A median background is built across cycles, so it needs a few of them. + min_cycles = self.MIN_REFERENCE_CYCLES if refs is not None else self.MIN_CYCLES + if n_cycles < min_cycles: + logger.debug("Not enough cycles for surface alignment yet", cycles=n_cycles, needed=min_cycles) + return None, None + + dt = self._median_dt(times) + trace_rate_hz = 1.0 / dt + # The estimator's 1.5 Hz default assumes a fast towed profile. A VNA + # sweep takes about a second, so tie the cutoff to the actual trace rate + # instead -- otherwise it sits above Nyquist and smooths nothing. + lowpass_hz = float(min(1.5, max(0.05, 0.25 * trace_rate_hz))) + + cfg = AlignConfig( + band_hz=options.band_hz, + pad_factor=int(options.pad_factor), + height_range_m=options.height_range_m, + height_step_m=float(options.height_step_m), + background=background, # type: ignore[arg-type] + smoothing=options.smoothing, # type: ignore[arg-type] + lowpass_hz=lowpass_hz, + ) + return cfg, refs + + def _full( + self, + cfg: AlignConfig, + refs: dict[str, NDArray[np.complexfloating]] | None, + fb: NDArray[np.floating], + spectra: NDArray[np.complexfloating], + times: NDArray[np.floating], + options: AlignOptions, + ) -> AlignmentState | None: + """Estimate every cycle from scratch.""" + baselines = None if options.baseline_m is None else {self.channel: float(options.baseline_m)} + + try: + res = estimate_surface( + {self.channel: spectra}, fb, times, cfg, + baselines=baselines, references=refs, + ) + except Exception as exc: # noqa: BLE001 + logger.error("Surface estimation failed", error=repr(exc)) + return None + + baseline_m = float(res.baselines_m[self.channel]) + height_ref = ( + float(options.height_ref_m) + if options.height_ref_m + else float(np.median(res.height_m)) + ) + + residual: float | None = None + spread: float | None = None + if options.geometry_check: + residual, spread = self._geometry_check( + cfg, fb, spectra, res.baselines_m, res.height_m, + ) + + logger.info( + "Surface alignment estimated", + cycles=int(res.height_m.size), + baseline_cm=round(baseline_m * 100.0, 2), + height_cm=round(float(np.median(res.height_m)) * 100.0, 2), + height_ref_cm=round(height_ref * 100.0, 2), + rejected_fraction=round(res.quality.get("rejected_fraction", 0.0), 3), + confidence_median=round(res.quality.get("confidence_median", 0.0), 2), + geometry_residual_mm=None if residual is None else round(residual * 1000.0, 1), + ) + + return AlignmentState( + heights_m=np.asarray(res.height_m, dtype=float), + raw_heights_m=np.asarray(res.height_raw_m, dtype=float), + confidence=np.asarray(res.confidence, dtype=float), + baseline_m=baseline_m, + height_ref_m=height_ref, + quality=dict(res.quality), + geometry_residual_m=residual, + geometry_spread_m=spread, + ) + + def _extend( + self, + cfg: AlignConfig, + refs: dict[str, NDArray[np.complexfloating]] | None, + fb: NDArray[np.floating], + window: NDArray[np.complexfloating], + times: NDArray[np.floating], + n_cycles: int, + ) -> AlignmentState | None: + """Estimate only the new cycles from a trailing window. + + The baseline and the reference height stay frozen -- they define the + common target every already-drawn column was shifted to. Smoothing is + forced to the causal Kalman: inside a trailing window only the past is + available, and a zero-phase low-pass would distort exactly the newest + samples we are after. + """ + state = self._state + assert state is not None + new_count = n_cycles - state.covered + + try: + res = estimate_surface( + {self.channel: window}, fb, times, replace(cfg, smoothing="live"), + baselines={self.channel: state.baseline_m}, references=refs, + ) + except Exception as exc: # noqa: BLE001 + logger.error("Surface estimation failed for trailing window", error=repr(exc)) + return None + + updated = replace( + state, + heights_m=np.concatenate([state.heights_m, res.height_m[-new_count:]]), + raw_heights_m=np.concatenate([state.raw_heights_m, res.height_raw_m[-new_count:]]), + confidence=np.concatenate([state.confidence, res.confidence[-new_count:]]), + quality=dict(res.quality), + incremental=True, + ) + self._state = updated + logger.debug( + "Surface alignment extended", + new_cycles=new_count, + covered=updated.covered, + height_cm=round(float(updated.heights_m[-1]) * 100.0, 2), + ) + return updated + + def _geometry_check( + self, + cfg: AlignConfig, + fb: NDArray[np.floating], + spectra: NDArray[np.complexfloating], + baselines: dict[str, float], + height_m: NDArray[np.floating], + ) -> tuple[float | None, float | None]: + """Re-derive the height from each channel's own pick and compare. + + With several channels the spread of the per-channel means is the metric + that catches a confident mistracking. With a single channel there is no + spread to look at, so the useful number is the residual between the + joint height and what that channel alone implies: it grows when the + tracked peak is not where the geometry says the surface should be. + """ + # check_geometry subtracts a background itself and is not handed the + # reference trace, so a "reference" config would make it raise. + probe_cfg = cfg if cfg.background != "reference" else replace(cfg, background="global_median") + try: + per_channel = check_geometry({self.channel: spectra}, fb, probe_cfg, baselines, height_m) + except Exception as exc: # noqa: BLE001 + logger.warning("Geometry check failed", error=repr(exc)) + return None, None + + means = [float(np.mean(h)) for h in per_channel.values()] + spread = float(np.max(means) - np.min(means)) if len(means) > 1 else 0.0 + implied = per_channel[self.channel] + residual = float(np.sqrt(np.mean((implied - height_m) ** 2))) + return residual, spread + + # ------------------------------------------------------------------ # + # Application + # ------------------------------------------------------------------ # + def shift( + self, + spectrum: NDArray[np.complexfloating], + freq_hz: NDArray[np.floating], + state: AlignmentState, + index: int, + ) -> NDArray[np.complexfloating]: + """Move cycle `index` so its surface lands at the reference height. + + Applied over the full measured band: the shift is a linear phase ramp + derived from geometry, so it is an exact fractional delay and does not + care which sub-band the height came from. + """ + if index < 0 or index >= state.covered: + return spectrum + + aligned = apply_alignment( + {self.channel: np.asarray(spectrum, dtype=np.complex128)[None, :]}, + np.asarray(freq_hz, dtype=float), + {self.channel: state.baseline_m}, + np.asarray([state.heights_m[index]], dtype=float), + height_ref_m=state.height_ref_m, + ) + return aligned[self.channel][0] + + def surface_path_m(self, state: AlignmentState) -> float: + """Two-way path length of the aligned surface echo, in metres. + + This is where the surface sits after alignment: `sqrt(L^2 + (2h)^2)`. + """ + return float(np.hypot(state.baseline_m, 2.0 * state.height_ref_m)) + + # ------------------------------------------------------------------ # + # Helpers + # ------------------------------------------------------------------ # + @staticmethod + def _median_dt(times_s: NDArray[np.floating]) -> float: + if times_s.size < 2: + return 1.0 + dt = float(np.median(np.diff(times_s))) + return dt if dt > 0.0 else 1.0 + + @staticmethod + def _sanitize_times(times_s: NDArray[np.floating]) -> NDArray[np.floating]: + """Return a strictly increasing time axis. + + The estimator interpolates over rejected cycles and derives the trace + rate from these values, both of which need monotonic times. Imported + histories can carry missing or duplicated timestamps, so fall back to a + unit-spaced index in that case. + """ + t = np.asarray(times_s, dtype=float) + if t.size == 0: + return t + if not np.all(np.isfinite(t)) or (t.size > 1 and np.any(np.diff(t) <= 0.0)): + return np.arange(t.size, dtype=float) + return t - t[0] diff --git a/vna_system/core/processors/configs/bscan_config.json b/vna_system/core/processors/configs/bscan_config.json index 6db4813..2081d20 100644 --- a/vna_system/core/processors/configs/bscan_config.json +++ b/vna_system/core/processors/configs/bscan_config.json @@ -19,6 +19,19 @@ "apply_eps_correction": true, "eps_r": 3.7, "eps_boundary_m": 0.0, + "surface_align": true, + "align_start_freq": 300.0, + "align_stop_freq": 6010.0, + "align_height_min": 0.25, + "align_height_max": 0.9, + "align_background": "reference", + "align_smoothing": "offline", + "align_baseline_mode": "manual", + "align_baseline_m": 0.36, + "align_ref_m": 0.0, + "align_pad_factor": 8, + "align_geometry_check": false, + "align_show_line": false, "data_limit": 500, "y_min": -50, "y_max": 40, diff --git a/vna_system/core/processors/implementations/bscan_processor.py b/vna_system/core/processors/implementations/bscan_processor.py index 7ae50ec..323be5e 100644 --- a/vna_system/core/processors/implementations/bscan_processor.py +++ b/vna_system/core/processors/implementations/bscan_processor.py @@ -10,6 +10,11 @@ from numpy.typing import NDArray from scipy.ndimage import gaussian_filter1d from vna_system.core.logging.logger import get_component_logger +from vna_system.core.processing.surface_alignment import ( + AlignmentState, + AlignOptions, + SurfaceAlignmentEngine, +) from vna_system.core.processors.base_processor import BaseProcessor, UIParameter, ProcessedResult from vna_system.core.acquisition.sweep_buffer import SweepData from vna_system.core.config import SPEED_OF_LIGHT_M_S @@ -53,6 +58,8 @@ class BScanProcessor(BaseProcessor): - IFFT with frequency range filtering - Depth windowing with gain shaping - Plot history accumulation for multi-sweep heatmaps + - Optional surface alignment: every sweep is delayed so the ground echo + lands at a common height, which flattens the surface across the B-scan """ @@ -64,6 +71,9 @@ class BScanProcessor(BaseProcessor): # Local plot history (separate from sweep history maintained by BaseProcessor) self._plot_history: list[dict[str, Any]] = [] + + # Surface-height tracker; holds one height per sweep-history index + self._aligner = SurfaceAlignmentEngine() self._ach_norm_curve: NDArray[np.complex128] | None = None self._ach_norm_mtime: float | None = None self._s11_norm_curve_1: NDArray[np.complex128] | None = None @@ -98,9 +108,30 @@ class BScanProcessor(BaseProcessor): "if_normalize" : False, "if_draw_level" : False, "detection_level" : 5, - "apply_eps_correction": False, - "eps_r": 4.0, - "eps_boundary_m": 0.0, + "apply_eps_correction": False, + "eps_r": 4.0, + "eps_boundary_m": 0.0, + # Surface alignment (see vna_system.core.processing.surface_align) + "surface_align": False, # Master toggle + "align_start_freq": 300.0, # Estimation band start (MHz) + "align_stop_freq": 5000.0, # Estimation band stop (MHz) + "align_height_min": 0.25, # Height search range (m) + "align_height_max": 0.90, + # A recorded reference is the honest background for an air-coupled + # setup: a median background also removes a surface that is not + # moving, which biases the height. Falls back to the global median + # when no matching reference is loaded. + "align_background": "reference", + # Causal Kalman by default: it costs nothing live and, unlike the + # brick-wall low-pass, it does not ring at the ends of a profile + # whose height drifts. + "align_smoothing": "live", + "align_baseline_mode": "auto", # TX-RX baseline: measured or manual + "align_baseline_m": 0.30, # Used when mode == "manual" + "align_ref_m": 0.0, # Target height; 0 -> median + "align_pad_factor": 8, # Zero-padding of the estimator IFFT + "align_geometry_check": True, # Cross-check the tracked height + "align_show_line": True, # Draw the aligned-surface marker } def get_ui_parameters(self) -> list[UIParameter]: @@ -252,6 +283,96 @@ class BScanProcessor(BaseProcessor): value=cfg["eps_boundary_m"], options={"min": 0.0, "max": 2.5, "step": 0.01, "dtype": "float"}, ), + + # --- Surface alignment --- + UIParameter( + name="surface_align", + label="Выравнивание по поверхности", + type="toggle", + value=cfg["surface_align"], + ), + UIParameter( + name="align_start_freq", + label="Выравнивание: нач. частота (МГц)", + type="slider", + value=cfg["align_start_freq"], + options={"min": 100.0, "max": 8800.0, "step": 10.0, "dtype": "float"}, + ), + UIParameter( + name="align_stop_freq", + label="Выравнивание: кон. частота (МГц)", + type="slider", + value=cfg["align_stop_freq"], + options={"min": 100.0, "max": 8800.0, "step": 10.0, "dtype": "float"}, + ), + UIParameter( + name="align_height_min", + label="Выравнивание: мин. высота (м)", + type="slider", + value=cfg["align_height_min"], + options={"min": 0.05, "max": 2.0, "step": 0.01, "dtype": "float"}, + ), + UIParameter( + name="align_height_max", + label="Выравнивание: макс. высота (м)", + type="slider", + value=cfg["align_height_max"], + options={"min": 0.10, "max": 3.0, "step": 0.01, "dtype": "float"}, + ), + UIParameter( + name="align_background", + label="Выравнивание: фон", + type="select", + value=cfg["align_background"], + options={"choices": ["global_median", "sliding_median", "reference"]}, + ), + UIParameter( + name="align_smoothing", + label="Выравнивание: сглаживание", + type="select", + value=cfg["align_smoothing"], + options={"choices": ["offline", "live"]}, + ), + UIParameter( + name="align_baseline_mode", + label="Выравнивание: база TX-RX", + type="select", + value=cfg["align_baseline_mode"], + options={"choices": ["auto", "manual"]}, + ), + UIParameter( + name="align_baseline_m", + label="Выравнивание: база вручную (м)", + type="slider", + value=cfg["align_baseline_m"], + options={"min": 0.0, "max": 2.0, "step": 0.01, "dtype": "float"}, + ), + UIParameter( + name="align_ref_m", + label="Выравнивание: опорная высота (м, 0 = медиана)", + type="slider", + value=cfg["align_ref_m"], + options={"min": 0.0, "max": 3.0, "step": 0.01, "dtype": "float"}, + ), + UIParameter( + name="align_pad_factor", + label="Выравнивание: паддинг", + type="select", + value=cfg["align_pad_factor"], + options={"choices": [4, 8, 16]}, + ), + UIParameter( + name="align_geometry_check", + label="Выравнивание: контроль геометрии", + type="toggle", + value=cfg["align_geometry_check"], + ), + UIParameter( + name="align_show_line", + label="Выравнивание: линия поверхности", + type="toggle", + value=cfg["align_show_line"], + ), ] def update_config(self, updates: dict[str, Any]) -> None: @@ -269,8 +390,37 @@ class BScanProcessor(BaseProcessor): with self._lock: self._sweep_history.clear() self._plot_history.clear() + self._aligner.invalidate() logger.info("Plot and sweep history cleared completely", processor_id=self.processor_id) + # ------------------------------------------------------------------------- # + # History bookkeeping (surface alignment is indexed by sweep history) + # ------------------------------------------------------------------------- # + + def _trim_history(self) -> None: + """Trim sweep history, keeping the surface-height estimate index-aligned.""" + before = len(self._sweep_history) + super()._trim_history() + dropped = before - len(self._sweep_history) + if dropped > 0: + self._aligner.drop_leading(dropped) + + def import_history_data(self, history_data: list[dict[str, Any]]) -> None: + """Replace sweep history; the surface estimate is rebuilt from scratch.""" + with self._lock: + self._aligner.invalidate() + super().import_history_data(history_data) + + def append_history(self, history_data: list[dict[str, Any]]) -> ProcessedResult | None: + """Append sweep history; force a full re-estimate of the surface. + + A bulk append brings in timestamps from another recording, so the + trailing-window shortcut used for live sweeps does not apply. + """ + with self._lock: + self._aligner.invalidate() + return super().append_history(history_data) + def delete_column(self, column_index: int) -> bool: """ Delete a specific column (sweep) from the plot history. @@ -305,6 +455,9 @@ class BScanProcessor(BaseProcessor): if array_index < len(self._sweep_history): del self._sweep_history[array_index] + # Heights are indexed by sweep history, so the mapping is gone + self._aligner.invalidate() + logger.info( "Column deleted successfully", column_index=column_index, @@ -364,10 +517,19 @@ class BScanProcessor(BaseProcessor): sweep_data: SweepData, calibrated_data: SweepData | None, vna_config: dict[str, Any], + history_index: int | None = None, ) -> dict[str, Any]: """ Process a single sweep and prepare B-scan data. + Parameters + ---------- + history_index: + Position of this sweep in `_sweep_history`, used to look up its + surface height. None means "the newest entry", which is what the + live path needs; `recalculate` passes the explicit index while + replaying the history. + Returns ------- dict @@ -408,6 +570,12 @@ class BScanProcessor(BaseProcessor): complex_data = self._apply_ach_normalization(complex_data) complex_data = self._apply_s11_normalization(complex_data) + # Optional surface alignment: delay this sweep so its ground echo + # sits at the same height as every other sweep in the B-scan + complex_data, align_info = self._apply_surface_alignment( + complex_data, vna_config, history_index + ) + # Keep frequency controls in sync with the current VNA config self._update_frequency_ranges(vna_config) @@ -453,6 +621,8 @@ class BScanProcessor(BaseProcessor): "all_distance_data": all_distance, "all_sweep_numbers": all_sweep_numbers, "all_timestamps": all_timestamps, + # Surface alignment (None when disabled or not yet available) + "surface_alignment": align_info, } except Exception as exc: # noqa: BLE001 @@ -588,6 +758,18 @@ class BScanProcessor(BaseProcessor): if processed_data.get("reference_used", False): config_info += " | Открытый воздух: ВКЛ" + alignment = processed_data.get("surface_alignment") + if alignment: + quality = alignment.get("quality") or {} + config_info += ( + f" | Выравнивание: h={alignment['height_m'] * 100:.1f} см" + f" (опора {alignment['height_ref_m'] * 100:.1f} см," + f" база {alignment['baseline_m'] * 100:.1f} см," + f" брак {quality.get('rejected_fraction', 0.0) * 100:.0f}%)" + ) + elif self._config.get("surface_align", False): + config_info += " | Выравнивание: нет оценки" + # if self._config["data_limitation"]: # config_info += f" | Limit: {self._config['data_limitation']}" @@ -642,6 +824,25 @@ class BScanProcessor(BaseProcessor): } ] + if alignment and self._config.get("align_show_line", True): + # Where the surface now sits: every column was shifted to this depth + layout["shapes"] = layout.get("shapes", []) + [ + { + "type": "line", + "xref": "paper", + "yref": "y", + "x0": 0, + "x1": 1, + "y0": alignment["surface_depth_m"], + "y1": alignment["surface_depth_m"], + "line": { + "width": 2, + "dash": "dashdot", + "color": "#FFB300", + }, + } + ] + if detected_trace is not None: return {"data": [heatmap_trace,detected_trace], "layout": layout} return {"data": [heatmap_trace], "layout": layout} @@ -673,6 +874,7 @@ class BScanProcessor(BaseProcessor): "all_distance_data": [], "all_sweep_numbers": [], "all_timestamps": [], + "surface_alignment": None, } plotly_conf = self.generate_plotly_config(empty_data, {}) ui_params = self.get_ui_parameters() @@ -689,16 +891,21 @@ class BScanProcessor(BaseProcessor): # Clear existing plot history to rebuild from scratch self._plot_history.clear() + # A rebuild has the whole profile in hand, so let the surface + # estimate be smoothed the way the operator asked for + if self._config.get("surface_align", False): + self._aligner.refresh_if_incremental(str(self._config["align_smoothing"])) + # Process all sweeps in history with current config last_processed = None last_vna_config = {} - for entry in self._sweep_history: + for index, entry in enumerate(self._sweep_history): sweep_data = entry["sweep_data"] calibrated_data = entry["calibrated_data"] vna_config = entry["vna_config"] # Use process_sweep to handle the processing logic - processed = self.process_sweep(sweep_data, calibrated_data, vna_config) + processed = self.process_sweep(sweep_data, calibrated_data, vna_config, index) # Skip if processing failed if "error" not in processed: @@ -985,6 +1192,209 @@ class BScanProcessor(BaseProcessor): logger.error("Failed to load S11 normalization file", file=str(file_path), error=repr(exc)) return None + # ------------------------------------------------------------------------- + # Surface alignment + # ------------------------------------------------------------------------- + + def _frequency_bounds(self, vna_config: dict[str, Any]) -> tuple[float, float]: + """Measured band of a sweep in Hz; the device config wins over the sliders.""" + if vna_config: + return ( + float(vna_config.get("start_freq", 100e6)), + float(vna_config.get("stop_freq", 8.8e9)), + ) + return ( + float(self._config["start_freq"]) * 1e6, + float(self._config["stop_freq"]) * 1e6, + ) + + def _frequency_axis(self, vna_config: dict[str, Any], n_points: int) -> NDArray[np.floating]: + """Uniform frequency grid (Hz) for a sweep of `n_points` samples.""" + start_hz, stop_hz = self._frequency_bounds(vna_config) + return np.linspace(start_hz, stop_hz, n_points, dtype=float) + + def _align_options(self) -> AlignOptions: + """Translate the UI config into estimator options (SI units).""" + cfg = self._config + band_lo, band_hi = sorted( + (float(cfg["align_start_freq"]) * 1e6, float(cfg["align_stop_freq"]) * 1e6) + ) + height_lo, height_hi = sorted( + (float(cfg["align_height_min"]), float(cfg["align_height_max"])) + ) + step_m = 0.001 + if height_hi - height_lo < 4.0 * step_m: + # A degenerate range leaves the search grid with no room to pick a peak + height_hi = height_lo + 4.0 * step_m + + manual_baseline = str(cfg["align_baseline_mode"]) == "manual" + height_ref = float(cfg["align_ref_m"]) + + return AlignOptions( + band_hz=(band_lo, band_hi), + height_range_m=(height_lo, height_hi), + height_step_m=step_m, + background=str(cfg["align_background"]), + smoothing=str(cfg["align_smoothing"]), + pad_factor=int(cfg["align_pad_factor"]), + baseline_m=float(cfg["align_baseline_m"]) if manual_baseline else None, + height_ref_m=height_ref if height_ref > 0.0 else None, + geometry_check=bool(cfg["align_geometry_check"]), + ) + + def _apply_surface_alignment( + self, + complex_data: NDArray[np.complex128], + vna_config: dict[str, Any], + history_index: int | None, + ) -> tuple[NDArray[np.complex128], dict[str, Any] | None]: + """ + Delay one sweep so its ground echo lands at the common reference height. + + Returns the (possibly unchanged) spectrum and, when a height was + available, the numbers describing the decision for this sweep. + + Runs under the processor lock: the live path calls this outside it, while + a websocket-triggered rebuild holds it for the whole replay, and both + read and update the same height cache. + """ + if not self._config.get("surface_align", False): + return complex_data, None + + with self._lock: + state = self._ensure_alignment_state(vna_config) + if state is None: + return complex_data, None + + index = ( + len(self._sweep_history) - 1 if history_index is None else int(history_index) + ) + if index < 0 or index >= state.covered: + logger.debug( + "No surface height for this sweep; leaving it unshifted", + index=index, + covered=state.covered, + ) + return complex_data, None + + freq_axis = self._frequency_axis(vna_config, complex_data.size) + shifted = self._aligner.shift(complex_data, freq_axis, state, index) + + info = { + "height_m": float(state.heights_m[index]), + "height_raw_m": float(state.raw_heights_m[index]), + "confidence": float(state.confidence[index]), + "height_ref_m": float(state.height_ref_m), + "baseline_m": float(state.baseline_m), + "surface_depth_m": self._surface_depth_m(state), + "cycles": state.covered, + "quality": dict(state.quality), + "geometry_residual_m": state.geometry_residual_m, + "geometry_spread_m": state.geometry_spread_m, + } + return shifted, info + + def _surface_depth_m(self, state: AlignmentState) -> float: + """Where the aligned surface lands on the plotted depth axis (metres).""" + two_way_path_m = self._aligner.surface_path_m(state) + depth_m = (two_way_path_m - 2.0 * float(self._config["cut"])) / 2.0 + corrected = self._apply_eps_depth_correction(np.asarray([depth_m], dtype=float)) + return float(corrected[0]) + + def _ensure_alignment_state(self, vna_config: dict[str, Any]) -> AlignmentState | None: + """ + Return an up-to-date surface-height estimate for the whole sweep history. + + Cheap when nothing changed: the engine reuses its cached heights, extends + them from a trailing window when new sweeps arrived, and only recomputes + everything when the settings or the frequency grid changed. + """ + with self._lock: + n_cycles = len(self._sweep_history) + latest = self._sweep_history[-1] if n_cycles else None + + if latest is None: + return None + + source = latest.get("calibrated_data") or latest.get("sweep_data") + n_points = len(getattr(source, "points", None) or ()) + if n_points < 2: + return None + + freq_axis = self._frequency_axis(vna_config, n_points) + options = self._align_options() + mode, start_index = self._aligner.plan( + n_cycles, options, self._aligner.freq_key(freq_axis) + ) + if mode == "reuse": + return self._aligner.state + + collected = self._collect_alignment_spectra(start_index, n_cycles, n_points) + if collected is None: + return None + + spectra, times = collected + reference = ( + self._alignment_reference(n_points) + if options.background == "reference" + else None + ) + return self._aligner.update( + mode, start_index, spectra, freq_axis, times, options, reference + ) + + def _collect_alignment_spectra( + self, + start: int, + stop: int, + n_points: int, + ) -> tuple[NDArray[np.complex128], NDArray[np.floating]] | None: + """ + Stack calibrated spectra of history entries `[start, stop)` for the estimator. + + Deliberately the calibrated data rather than the display pipeline output: + open-air subtraction and the amplitude normalizations are display choices, + while the estimator brings its own background model. What it does need is + a uniform grid, so a history with mixed sweep lengths disables alignment + instead of silently mixing bands. + """ + with self._lock: + entries = self._sweep_history[start:stop] + + rows: list[NDArray[np.complex128]] = [] + times: list[float] = [] + for entry in entries: + source = entry.get("calibrated_data") or entry.get("sweep_data") + spectrum = self._get_complex_s11(source) if source is not None else None + if spectrum is None or spectrum.size != n_points: + logger.warning( + "Surface alignment skipped: sweep history is not uniform", + expected_points=n_points, + got=None if spectrum is None else int(spectrum.size), + ) + return None + rows.append(spectrum) + # Acquisition time when we have it: the smoothing derives the trace + # rate from these, and the entry timestamp is when the sweep was + # processed, which is not the same thing during a replay. + timestamp = getattr(source, "timestamp", None) or entry.get("timestamp") + times.append(float(timestamp or 0.0)) + + if not rows: + return None + return np.asarray(rows, dtype=np.complex128), np.asarray(times, dtype=float) + + def _alignment_reference(self, n_points: int) -> NDArray[np.complex128] | None: + """Calibrated open-air reference used by background='reference'.""" + with self._lock: + latest = self._sweep_history[-1] if self._sweep_history else None + + reference_data = None if latest is None else latest.get("reference_data") + spectrum = self._get_complex_s11(reference_data) if reference_data is not None else None + if spectrum is None or spectrum.size != n_points: + return None + return spectrum + def _update_frequency_ranges(self, vna_config: dict[str, Any]) -> None: """Clamp configured frequency sliders to VNA limits.""" if not vna_config: @@ -1067,12 +1477,7 @@ class BScanProcessor(BaseProcessor): """Full analysis pipeline: limit -> IFFT -> depth shaping.""" try: # Determine effective frequency range (Hz) - if vna_config: - freq_start = float(vna_config.get("start_freq", 100e6)) - freq_stop = float(vna_config.get("stop_freq", 8.8e9)) - else: - freq_start = self._config["start_freq"] * 1e6 - freq_stop = self._config["stop_freq"] * 1e6 + freq_start, freq_stop = self._frequency_bounds(vna_config) # Determine sigma for smoothing if vna_config: @@ -1083,7 +1488,7 @@ class BScanProcessor(BaseProcessor): sigma = self._config["sigma"] # Frequency vector over current data length - freq_axis = np.linspace(freq_start, freq_stop, complex_data.size, dtype=float) + freq_axis = self._frequency_axis(vna_config, complex_data.size) # Hardcoded frequency-domain notch filter (manual on/off via comments) complex_data, _ = self._apply_hardcoded_notch_filter(complex_data, freq_axis) @@ -1369,5 +1774,24 @@ class BScanProcessor(BaseProcessor): "all_sweep_numbers": all_sweep_numbers, "all_timestamps": all_timestamps, } + state["surface_alignment"] = self._alignment_summary() return state + + def _alignment_summary(self) -> dict[str, Any] | None: + """Full surface-height track, for export and offline inspection.""" + align_state = self._aligner.state + if align_state is None: + return None + + return { + "enabled": bool(self._config.get("surface_align", False)), + "heights_m": align_state.heights_m.tolist(), + "raw_heights_m": align_state.raw_heights_m.tolist(), + "confidence": align_state.confidence.tolist(), + "baseline_m": align_state.baseline_m, + "height_ref_m": align_state.height_ref_m, + "quality": dict(align_state.quality), + "geometry_residual_m": align_state.geometry_residual_m, + "geometry_spread_m": align_state.geometry_spread_m, + }