surface_align added
This commit is contained in:
@@ -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",
|
||||
]
|
||||
@@ -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
|
||||
@@ -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]
|
||||
@@ -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,
|
||||
|
||||
@@ -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,
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user