From f9ea8817c0123bf0f5849812fd91af839b79cebe Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Mon, 24 Nov 2025 17:57:28 +0100 Subject: [PATCH 1/6] Refine single target detection implementation --- echopype/mask/__init__.py | 16 +- echopype/mask/api.py | 41 ++ .../single_target_detection/detect_esp3.py | 625 ++++++++++++++++++ .../single_target_detection/detect_matecho.py | 584 ++++++++++++++++ 4 files changed, 1264 insertions(+), 2 deletions(-) create mode 100644 echopype/mask/single_target_detection/detect_esp3.py create mode 100644 echopype/mask/single_target_detection/detect_matecho.py diff --git a/echopype/mask/__init__.py b/echopype/mask/__init__.py index b45d84b33..a9a3a207b 100644 --- a/echopype/mask/__init__.py +++ b/echopype/mask/__init__.py @@ -1,3 +1,15 @@ -from .api import apply_mask, detect_seafloor, detect_shoal, frequency_differencing +from .api import ( + apply_mask, + detect_seafloor, + detect_shoal, + detect_single_targets, + frequency_differencing, +) -__all__ = ["frequency_differencing", "apply_mask", "detect_seafloor", "detect_shoal"] +__all__ = [ + "frequency_differencing", + "apply_mask", + "detect_seafloor", + "detect_shoal", + "detect_single_targets", +] diff --git a/echopype/mask/api.py b/echopype/mask/api.py index 48c3d88c0..71a6d241c 100644 --- a/echopype/mask/api.py +++ b/echopype/mask/api.py @@ -13,6 +13,10 @@ from echopype.mask.seafloor_detection.bottom_basic import bottom_basic from echopype.mask.seafloor_detection.bottom_blackwell import bottom_blackwell +# for single_target_detection +from echopype.mask.single_target_detection.detect_esp3 import detect_esp3 +from echopype.mask.single_target_detection.detect_matecho import detect_matecho + from ..utils.io import validate_source from ..utils.prov import add_processing_level, echopype_prov_attrs, insert_input_processing_level from .freq_diff import _check_freq_diff_source_Sv, _parse_freq_diff_eq @@ -799,3 +803,40 @@ def detect_shoal( raise ValueError(f"Unsupported shoal detection method: {method}") return METHODS_SHOAL[method](ds, **params) + + +# Registry of supported methods for single_target_detection +METHODS_SINGLE_TARGET = { + "esp3": detect_esp3, + "matecho": detect_matecho, +} + + +def detect_single_targets( + ds: xr.Dataset, + method: str, + params: dict | None = None, +): + """ + Run single-target detection using the selected method. + + Parameters + ---------- + ds : xr.Dataset + Sv dataset (must contain what's needed by the method). + method : str + Name of the detection method to use (e.g., "esp3"). + params : dict, optional + Parameters for the detection function (method-specific). + + Returns + ------- + Whatever the method returns (e.g., dict or xr.Dataset). + """ + if method not in METHODS_SINGLE_TARGET: + raise ValueError(f"Unsupported single-target method: {method}") + + if params is None: + raise ValueError("No parameters given.") + + return METHODS_SINGLE_TARGET[method](ds, params) diff --git a/echopype/mask/single_target_detection/detect_esp3.py b/echopype/mask/single_target_detection/detect_esp3.py new file mode 100644 index 000000000..af5a82092 --- /dev/null +++ b/echopype/mask/single_target_detection/detect_esp3.py @@ -0,0 +1,625 @@ +import math + +import numpy as np +import xarray as xr + +DEFAULTS = { + "SoundSpeed": 1500.0, + "TS_threshold": -50.0, # dB + "PLDL": 6.0, # dB below peak + "MinNormPL": 0.7, + "MaxNormPL": 1.5, + "DataType": "CW", # CW only for now + "block_len": 1e7 / 3, # ~3.33M cells per block + "MaxBeamComp": 4.0, # dB (like default in MATLAB) + "MaxStdMinAxisAngle": 0.6, # rad (≈ 34°) + "MaxStdMajAxisAngle": 0.6, # rad +} + + +def init_st_struct(): + """Return a dict mirroring the MATLAB single_targets_tot struct.""" + return { + "TS_comp": [], + "TS_uncomp": [], + "Target_range": [], + "Target_range_disp": [], + "Target_range_min": [], + "Target_range_max": [], + "idx_r": [], + "StandDev_Angles_Minor_Axis": [], + "StandDev_Angles_Major_Axis": [], + "Angle_minor_axis": [], + "Angle_major_axis": [], + "Ping_number": [], + "Time": [], + "nb_valid_targets": 0, + "idx_target_lin": [], + "pulse_env_before_lin": [], + "pulse_env_after_lin": [], + "PulseLength_Normalized_PLDL": [], + "Transmitted_pulse_length": [], + "Heave": [], + "Roll": [], + "Pitch": [], + "Heading": [], + "Dist": [], + } + + +def _validate_params(p: dict): + if p.get("DataType", "CW") not in ("CW", "FM"): + raise ValueError("DataType must be 'CW' or 'FM'.") + if not (-120.0 <= float(p["TS_threshold"]) <= -20.0): + raise ValueError("TS_threshold must be in [-120, -20] dB.") + if not (1.0 <= float(p["PLDL"]) <= 30.0): + raise ValueError("PLDL must be in [1, 30] dB.") + if not (0.0 <= float(p["MinNormPL"]) <= 10.0): + raise ValueError("MinNormPL must be in [0, 10].") + if not (0.0 <= float(p["MaxNormPL"]) <= 10.0): + raise ValueError("MaxNormPL must be in [0, 10].") + if float(p["block_len"]) <= 0: + raise ValueError("block_len must be > 0.") + if not (0.0 <= float(p["MaxBeamComp"]) <= 18.0): + raise ValueError("MaxBeamComp must be in [0, 18] dB.") + if not (0.0 <= float(p["MaxStdMinAxisAngle"]) <= 45.0): + # keep same numeric limits as MATLAB; here we assume radians + raise ValueError("MaxStdMinAxisAngle must be in [0, 45].") + if not (0.0 <= float(p["MaxStdMajAxisAngle"]) <= 45.0): + raise ValueError("MaxStdMajAxisAngle must be in [0, 45].") + + +def _pull_nav(ds, name_candidates, ping_time_ref): + """Return a 1D float array aligned to ping_time_ref (or None).""" + for nm in name_candidates: + if nm in ds: + da = ds[nm] + for tdim in ("ping_time", "time", "time1", "time2"): + if tdim in da.dims: + if tdim != "ping_time": + da = da.rename({tdim: "ping_time"}) + da = da.reindex( + ping_time=ping_time_ref, + method="nearest", + tolerance=np.timedelta64(500, "ms"), + ) + return np.asarray(da.values, dtype=float) + if da.ndim == 1 and da.size == ping_time_ref.size: + return np.asarray(da.values, dtype=float) + return None + + +def simrad_beam_compensation( + along_angles_rad: np.ndarray, + athwart_angles_rad: np.ndarray, + bw_along_rad: float, + bw_athwart_rad: float, +) -> np.ndarray: + """ + Simrad beam pattern compensation [dB], same polynomial as in ESP3/Matecho. + + Parameters + ---------- + along_angles_rad : array + Alongship split-beam angles in radians. + athwart_angles_rad : array + Athwartship split-beam angles in radians. + bw_along_rad : float + 3 dB beamwidth alongship in radians. + bw_athwart_rad : float + 3 dB beamwidth athwartship in radians. + + Returns + ------- + comp : array + Beam compensation in dB (one-way). + """ + # same form as used in Matecho & ESP3: + # x = 2 * (phi_along / BW_along) + # y = 2 * (phi_athwart / BW_athwart) + # comp = 6.0206 * (x^2 + y^2 - 0.18 * x^2 * y^2) + x = 2.0 * (along_angles_rad / bw_along_rad) + y = 2.0 * (athwart_angles_rad / bw_athwart_rad) + comp = 6.0206 * (x**2 + y**2 - 0.18 * (x**2) * (y**2)) + return comp + + +def detect_esp3(ds: xr.Dataset, params: dict): + """ + ESP3-style single-target detection (CW branch) translated to Python. + + Parameters + ---------- + ds : xarray.Dataset + Calibrated Sv dataset (channel, ping_time, range_sample). + Must contain Sv (TS-like), depth, and optionally angles. + + params : dict + Must include: + - "channel": channel name + May include: + - "bottom_da" : bottom depth vs ping_time (DataArray) + - TS_threshold, PLDL, MinNormPL, MaxNormPL, SoundSpeed, + MaxBeamComp, MaxStdMinAxisAngle, MaxStdMajAxisAngle, + block_len, DataType + + Returns + ------- + out : dict + ESP3-like structure with TS_comp, TS_uncomp, Target_range, etc. + """ + if params is None: + params = {} + channel = params.get("channel") + if channel is None: + raise ValueError("params['channel'] is required.") + bottom_da = params.get("bottom_da", None) + + # Merge defaults + p = { + **DEFAULTS, + **{k: v for k, v in params.items() if k not in ("channel", "bottom_da")}, + } + p["DataType"] = p.get("DataType", "CW") + _validate_params(p) + + # ------------------------------------------------------------------ + # 1) Select data + # ------------------------------------------------------------------ + Sv = ds["Sv"].sel(channel=channel).transpose("ping_time", "range_sample") + Depth = ds["depth"].sel(channel=channel).transpose("ping_time", "range_sample") + + along = ds.get("angle_alongship") + if along is not None: + if "channel" in along.dims: + along = along.sel(channel=channel) + along = along.transpose("ping_time", "range_sample") + + athwt = ds.get("angle_athwartship") + if athwt is not None: + if "channel" in athwt.dims: + athwt = athwt.sel(channel=channel) + athwt = athwt.transpose("ping_time", "range_sample") + + nb_pings_tot = Sv.sizes["ping_time"] + nb_samples_tot = Sv.sizes["range_sample"] + idx_pings_tot = np.arange(nb_pings_tot, dtype=int) + idx_r_tot = np.arange(nb_samples_tot, dtype=int) + + # Optional bottom cropping (like idx_r_tot>max(idx_bot) removal) + if bottom_da is not None: + bb_global = bottom_da + if "channel" in bb_global.dims: + bb_global = bb_global.sel(channel=channel) + bb_global = bb_global.sel(ping_time=Sv["ping_time"]) + + D_all = Depth.values # (ping, range_sample) + b_all = bb_global.values # (ping,) + idx_bot = np.empty(nb_pings_tot, dtype=int) + for ip in range(nb_pings_tot): + under = np.nonzero(np.isfinite(D_all[ip]) & (D_all[ip] >= b_all[ip]))[0] + idx_bot[ip] = under[0] if under.size else (nb_samples_tot - 1) + max_idx = int(np.nanmax(idx_bot)) + idx_r_tot = np.arange(max_idx + 1, dtype=int) + + # Region/bad-data mask (placeholder: none => all False) + mask_inter_tot = np.zeros((idx_r_tot.size, idx_pings_tot.size), dtype=bool) + + # Block size + cells_per_ping = max(1, idx_r_tot.size) + block_size = int(min(math.ceil(float(p["block_len"]) / cells_per_ping), idx_pings_tot.size)) + num_ite = int(math.ceil(idx_pings_tot.size / max(block_size, 1))) + + # ------------------------------------------------------------------ + # 2) Nav / attitude + # ------------------------------------------------------------------ + pt_ref = Sv["ping_time"] + heading_arr = _pull_nav(ds, ["heading", "Heading"], pt_ref) + pitch_arr = _pull_nav(ds, ["pitch", "Pitch"], pt_ref) + roll_arr = _pull_nav(ds, ["roll", "Roll"], pt_ref) + heave_arr = _pull_nav(ds, ["heave", "Heave", "vertical_offset"], pt_ref) + dist_arr = _pull_nav(ds, ["dist", "Dist", "distance"], pt_ref) + + def _fallback(a, n): + return a if a is not None else np.zeros(n, dtype=float) + + heading = _fallback(heading_arr, nb_pings_tot) + pitch = _fallback(pitch_arr, nb_pings_tot) + roll = _fallback(roll_arr, nb_pings_tot) + heave = _fallback(heave_arr, nb_pings_tot) + dist = _fallback(dist_arr, nb_pings_tot) + + times = Sv["ping_time"].values + + # ------------------------------------------------------------------ + # 3) Prepare constants (alpha, c, T, Np) + # ------------------------------------------------------------------ + alpha = 0.0 + if "sound_absorption" in ds: + sa = ds["sound_absorption"] + try: + alpha = ( + float(sa.sel(channel=channel).values.item()) + if "channel" in sa.dims + else float(sa.values.item()) + ) + except Exception: + alpha = 0.0 + + c = float(p["SoundSpeed"]) + + # Pulse length & number of samples + if "Np" in p and p["Np"] is not None and int(p["Np"]) > 2: + Np = int(p["Np"]) + T = float(p.get("pulse_length", 1e-3)) + else: + dstep_da = Depth.diff("range_sample").median(skipna=True) + dstep = ( + float(dstep_da.values.item()) + if getattr(dstep_da, "size", 1) == 1 + else float(dstep_da.values) + ) + dt = (2.0 * dstep) / c if dstep > 0 else 1e-4 + T = float(p.get("pulse_length", 1e-3)) + Np = max(3, int(round(T / dt))) + + TS_thr = float(p["TS_threshold"]) + PLDL = float(p["PLDL"]) + min_len = max(1, int(Np * float(p["MinNormPL"]))) + max_len = max(1, int(math.ceil(Np * float(p["MaxNormPL"])))) + MaxBeamComp = float(p["MaxBeamComp"]) + MaxStdMin = float(p["MaxStdMinAxisAngle"]) + MaxStdMaj = float(p["MaxStdMajAxisAngle"]) + + # Beamwidths (taken from metadata if present; else fallback) + bw_along_rad = None + bw_athwart_rad = None + if "beamwidth_alongship" in ds: + bw_al = ds["beamwidth_alongship"].sel(channel=channel).values + bw_al = float(np.asarray(bw_al).item()) + if abs(bw_al) > np.pi: # assume deg + bw_al = np.deg2rad(bw_al) + bw_along_rad = bw_al + if "beamwidth_athwartship" in ds: + bw_at = ds["beamwidth_athwartship"].sel(channel=channel).values + bw_at = float(np.asarray(bw_at).item()) + if abs(bw_at) > np.pi: + bw_at = np.deg2rad(bw_at) + bw_athwart_rad = bw_at + + # Fallback beamwidths if missing + if bw_along_rad is None: + bw_along_rad = np.deg2rad(7.0) + if bw_athwart_rad is None: + bw_athwart_rad = np.deg2rad(7.0) + + # ------------------------------------------------------------------ + # 4) Output struct (like single_targets_tot) + # ------------------------------------------------------------------ + out = init_st_struct() + + # Bottom DA per-ping, sliced once + bb = None + if bottom_da is not None: + bb = bottom_da + if "channel" in bb.dims: + bb = bb.sel(channel=channel) + bb = bb.sel(ping_time=Sv["ping_time"]) + + # ------------------------------------------------------------------ + # 5) Block loop + # ------------------------------------------------------------------ + for ui in range(num_ite): + start = ui * block_size + stop = min((ui + 1) * block_size, idx_pings_tot.size) + idx_pings = idx_pings_tot[start:stop] + + # Per-block rows (respect bottom) + idx_r = idx_r_tot.copy() + if bb is not None and idx_pings.size > 0: + D_block = Depth.isel(ping_time=idx_pings).values + b_block = bb.isel(ping_time=idx_pings).values + idx_bot = np.empty(idx_pings.size, dtype=int) + for k in range(idx_pings.size): + under = np.nonzero(np.isfinite(D_block[k]) & (D_block[k] >= b_block[k]))[0] + idx_bot[k] = under[0] if under.size else (nb_samples_tot - 1) + max_idx = int(np.nanmax(idx_bot)) + idx_r = idx_r[idx_r <= max_idx] + + # Mask (regions) + mask = np.zeros((idx_r.size, idx_pings.size), dtype=bool) + if mask_inter_tot.size: + r_pos = np.searchsorted(idx_r_tot, idx_r) + p_pos = np.searchsorted(idx_pings_tot, idx_pings) + mask |= mask_inter_tot[np.ix_(r_pos, p_pos)] + + # Extract submatrices (TS ~ Sv here) + TS = ( + Sv.isel(ping_time=idx_pings, range_sample=idx_r) + .transpose("range_sample", "ping_time") + .values.copy() + ) + DEP = ( + Depth.isel(ping_time=idx_pings, range_sample=idx_r) + .transpose("range_sample", "ping_time") + .values + ) + + nb_samples, nb_pings = TS.shape + + # Angles sub-block (in radians) + along_block = None + athwt_block = None + if along is not None: + along_block = ( + along.isel(ping_time=idx_pings, range_sample=idx_r) + .transpose("range_sample", "ping_time") + .data + ) + along_block = np.deg2rad(along_block) + if athwt is not None: + athwt_block = ( + athwt.isel(ping_time=idx_pings, range_sample=idx_r) + .transpose("range_sample", "ping_time") + .data + ) + athwt_block = np.deg2rad(athwt_block) + + # Under-bottom mask + if bb is not None and idx_pings.size > 0: + bcol = bb.isel(ping_time=idx_pings).values.reshape(1, -1) + mask |= DEP >= bcol + + TS[mask] = -999.0 + if not np.any(TS > -999.0): + continue + + # Remove trailing all-masked rows + valid_rows = np.where((TS > -999.0).any(axis=1))[0] + last_row = valid_rows.max() + row_sel = slice(0, last_row + 1) + TS = TS[row_sel, :] + DEP = DEP[row_sel, :] + if along_block is not None: + along_block = along_block[row_sel, :] + if athwt_block is not None: + athwt_block = athwt_block[row_sel, :] + idx_r = idx_r[row_sel] + nb_samples, nb_pings = TS.shape + + # TVG matrix and Power (as in MATLAB) + r_eff = DEP - c * T / 4.0 + TVG_mat = np.where( + r_eff > 0.0, + 40.0 * np.log10(r_eff) + 2.0 * alpha * r_eff, + np.nan, + ) + Power = TS - TVG_mat + + # Peak detection on smoothed TS (or Power), like ESP3 + peak_calc = "TS" # ESP3 uses 'TS' by default in your snippet + if peak_calc.lower() == "power": + base_mat = Power + else: + base_mat = TS + + # moving average over Np/2 samples in linear domain + win = max(int(np.floor(Np / 2)), 1) + kernel = np.ones(win, dtype=float) / float(win) + + base_lin = 10.0 ** (base_mat / 10.0) + smooth_lin = np.apply_along_axis( + lambda col: np.convolve(col, kernel, mode="same"), + axis=0, + arr=base_lin, + ) + peak_mat = 10.0 * np.log10(smooth_lin) + peak_mat[TS <= -999.0] = -999.0 + + # --- Find peaks with MinSeparation ≈ Np (along range axis) --- + idx_peaks_bool = np.zeros_like(peak_mat, dtype=bool) + for jp in range(nb_pings): + z = peak_mat[:, jp] + valid = z > -999.0 + if np.count_nonzero(valid) < 3: + continue + cand = np.where(valid[1:-1] & (z[1:-1] > z[:-2]) & (z[1:-1] >= z[2:]))[0] + 1 + + # MinSeparation = Np + last = -(10**9) + keep_idx = [] + for k in cand: + if k - last >= Np: + keep_idx.append(k) + last = k + if keep_idx: + idx_peaks_bool[np.array(keep_idx, dtype=int), jp] = True + + idx_peaks_bool[TS <= -999.0] = False + + # linear indices of peaks (like idx_peaks_lin) + i_peaks_lin, j_peaks_lin = np.where(idx_peaks_bool) + nb_peaks = len(i_peaks_lin) + if nb_peaks == 0: + continue + + # Level of local maxima minus PLDL + peak_vals = peak_mat[idx_peaks_bool] # 1D (same order as i_peaks_lin) + pulse_level = peak_vals - PLDL + + # To mimic MATLAB's pulse-env scan, we will grow left/right until + # we fall below pulse_level or reach max_len. + pulse_env_before = np.zeros(nb_peaks, dtype=int) + pulse_env_after = np.zeros(nb_peaks, dtype=int) + + for idx_peak in range(nb_peaks): + ip = i_peaks_lin[idx_peak] + jp = j_peaks_lin[idx_peak] + thr = pulse_level[idx_peak] + zcol = peak_mat[:, jp] + + # grow up (before) + left = 0 + i = ip - 1 + while i >= 0 and left < max_len and np.isfinite(zcol[i]) and zcol[i] >= thr: + left += 1 + i -= 1 + + # grow down (after) + right = 0 + i = ip + 1 + nP = zcol.size + while i < nP and right < max_len and np.isfinite(zcol[i]) and zcol[i] >= thr: + right += 1 + i += 1 + + pulse_env_before[idx_peak] = left + pulse_env_after[idx_peak] = right + + pulse_length_lin = pulse_env_before + pulse_env_after + 1 + # good pulses according to [MinNormPL, MaxNormPL] + good_pulses = (pulse_length_lin >= min_len) & (pulse_length_lin <= max_len) + + # Filter peaks to targets + i_targets = i_peaks_lin[good_pulses] + j_targets = j_peaks_lin[good_pulses] + pulse_before = pulse_env_before[good_pulses] + pulse_after = pulse_env_after[good_pulses] + pulse_len = pulse_length_lin[good_pulses] + nb_targets = len(i_targets) + if nb_targets == 0: + continue + + # ------------------------------------------------------------------ + # For each target, compute: + # - samples_targets_power / range / angles + # - target_range (weighted by Power) + # - beam comp at peak (target_comp) + # - TS_uncomp = target_peak_power + TVG(target_range) + # - TS_comp = TS_uncomp + target_comp + # - apply angle std filter and MaxBeamComp & TS_threshold + # ------------------------------------------------------------------ + for it in range(nb_targets): + ip = int(i_targets[it]) + jp = int(j_targets[it]) + left = int(pulse_before[it]) + right = int(pulse_after[it]) + + i0 = max(ip - left, 0) + i1 = min(ip + right, nb_samples - 1) + seg_slice = slice(i0, i1 + 1) + + pow_seg = Power[seg_slice, jp] + dep_seg = DEP[seg_slice, jp] + + # If all NaNs or -inf => skip + if not np.isfinite(pow_seg).any(): + continue + + # angles segment + if along_block is not None and athwt_block is not None: + al_seg = along_block[seg_slice, jp] + at_seg = athwt_block[seg_slice, jp] + else: + al_seg = np.full_like(pow_seg, np.nan, dtype=float) + at_seg = np.full_like(pow_seg, np.nan, dtype=float) + + # std + mean of angles (for filtering & output) + std_al = float(np.nanstd(al_seg)) + std_at = float(np.nanstd(at_seg)) + phi_al = float(np.nanmean(al_seg)) + phi_at = float(np.nanmean(at_seg)) + + # Angle stability filter (CW only) + if p["DataType"] == "CW": + if std_al > MaxStdMin or std_at > MaxStdMaj: + continue + + # Power-weighted range (like target_range) + w = 10.0 ** (pow_seg / 10.0) + if not np.isfinite(w).any() or np.nansum(w) <= 0: + continue + target_range = float(np.nansum(w * dep_seg) / np.nansum(w) - c * T / 4.0) + if target_range < 0: + target_range = 0.0 + + target_range_min = float(np.nanmin(dep_seg)) + target_range_max = float(np.nanmax(dep_seg)) + + # Peak power within pulse window + # (we take max of pow_seg and its index) + idx_loc_peak = int(np.nanargmax(pow_seg)) + target_peak_power = float(pow_seg[idx_loc_peak]) + if not np.isfinite(target_peak_power): + continue + + # Beam comp at peak sample + if along_block is not None and athwt_block is not None: + al_peak = al_seg[idx_loc_peak] + at_peak = at_seg[idx_loc_peak] + if np.isfinite(al_peak) and np.isfinite(at_peak): + target_comp = float( + simrad_beam_compensation( + np.array([al_peak]), + np.array([at_peak]), + bw_along_rad, + bw_athwart_rad, + )[0] + ) + else: + target_comp = 0.0 + else: + target_comp = 0.0 + + # TVG at target_range (like ESP3: TVG = 40 log10(R) + 2 alpha R) + if target_range > 0: + TVG = 40.0 * np.log10(target_range) + 2.0 * alpha * target_range + else: + TVG = np.nan + + if not np.isfinite(TVG): + continue + + target_TS_uncomp = target_peak_power + TVG + target_TS_comp = target_TS_uncomp + target_comp + + # Filters on TS_comp and beam comp + if target_TS_comp <= TS_thr: + continue + if abs(target_comp) > MaxBeamComp: + continue + + # Store detection + ping_glob = int(idx_pings[jp]) + range_idx_glob = int(idx_r[ip]) + + out["TS_comp"].append(target_TS_comp) + out["TS_uncomp"].append(target_TS_uncomp) + out["Target_range"].append(target_range) + out["Target_range_disp"].append(target_range + c * T / 4.0) + out["Target_range_min"].append(target_range_min) + out["Target_range_max"].append(target_range_max) + + out["idx_r"].append(range_idx_glob) + out["Ping_number"].append(ping_glob) + out["Time"].append(times[ping_glob]) + out["idx_target_lin"].append(int(ping_glob * nb_samples_tot + range_idx_glob)) + + out["pulse_env_before_lin"].append(int(left)) + out["pulse_env_after_lin"].append(int(right)) + out["PulseLength_Normalized_PLDL"].append(pulse_len[it] / float(Np)) + out["Transmitted_pulse_length"].append(int(pulse_len[it])) + + out["StandDev_Angles_Minor_Axis"].append(std_al) + out["StandDev_Angles_Major_Axis"].append(std_at) + out["Angle_minor_axis"].append(phi_al) + out["Angle_major_axis"].append(phi_at) + + out["Heave"].append(float(heave[ping_glob])) + out["Roll"].append(float(roll[ping_glob])) + out["Pitch"].append(float(pitch[ping_glob])) + out["Heading"].append(float(heading[ping_glob])) + out["Dist"].append(float(dist[ping_glob])) + + out["nb_valid_targets"] = len(out["TS_comp"]) + return out diff --git a/echopype/mask/single_target_detection/detect_matecho.py b/echopype/mask/single_target_detection/detect_matecho.py new file mode 100644 index 000000000..aed7b277e --- /dev/null +++ b/echopype/mask/single_target_detection/detect_matecho.py @@ -0,0 +1,584 @@ +# matecho.py +import math + +import numpy as np +import xarray as xr + +MATECHO_DEFAULTS = { + "SoundSpeed": 1500.0, + "TS_threshold": -50.0, # dB + "MaxAngleOneWayCompression": 6.0, # dB (one-way; test uses 2x) + "MaxPhaseDeviation": 8.0, # phase steps (128/pi units), Matecho GUI + "MinEchoLength": 0.8, # in pulse lengths + "MaxEchoLength": 1.8, + "MinEchoSpace": 1.0, # in pulse lengths + "MinEchoDepthM": 3.0, + "MaxEchoDepthM": 38.0, + "tvg_start_sample": 3, # EK60=3, EK80=1 + "block_len": 1e7 / 3, + # Beam pattern (fallbacks if no metadata): + "beamwidth_along_3dB_rad": np.deg2rad(7.0), + "beamwidth_athwart_3dB_rad": np.deg2rad(7.0), + "steer_along_rad": 0.0, + "steer_athwart_rad": 0.0, + # Angle sensitivities (phase steps), overwritten from ds: + "angle_sens_al": 1.0, + "angle_sens_at": 1.0, + # Sv->TS constant terms; keep 0 if unknown: + "psi_two_way": 0.0, + "Sa_correction": 0.0, + "Sa_EK80_nominal": 0.0, +} + + +def _fallback_1d_aligned(ds, names, ping_time): + """ + Try to get a 1D variable aligned to ping_time from ds, using + several possible names and time dims, with nearest-neighbour + reindexing and 500 ms tolerance. + """ + for nm in names: + if nm in ds: + da = ds[nm] + for tdim in ("ping_time", "time", "time1", "time2"): + if tdim in da.dims: + if tdim != "ping_time": + da = da.rename({tdim: "ping_time"}) + da = da.reindex( + ping_time=ping_time, + method="nearest", + tolerance=np.timedelta64(500, "ms"), + ) + return np.asarray(da.values, dtype=float) + if da.ndim == 1 and da.size == ping_time.size: + return np.asarray(da.values, dtype=float) + return None + + +def detect_matecho(ds: xr.Dataset, params: dict): + """ + Matecho-inspired single-target detector (CW path only). + + Parameters + ---------- + ds : xarray.Dataset + Calibrated Sv dataset (e.g., from echopype.calibrate.compute_Sv), + with dimensions (channel, ping_time, range_sample) and variables: + Sv, depth, angle_alongship, angle_athwartship, and instrument + metadata (beamwidths, offsets, angle sensitivities). + params : dict + Parameters including at least: + - "channel": channel name (string) + Optional: + - "bottom_da": xarray.DataArray of bottom depth vs ping_time + - Matecho parameters: TS_threshold, MaxAngleOneWayCompression, + MaxPhaseDeviation, MinEchoLength, MaxEchoLength, MinEchoSpace, + MinEchoDepthM, MaxEchoDepthM, pulse_length, tvg_start_sample, etc. + + Returns + ------- + out : dict + Matecho-like single-target "struct" (Python dict), with keys such as + TS_comp, TS_uncomp, Target_range, Ping_number, Time, etc. + """ + if params is None: + raise ValueError("params is required.") + channel = params.get("channel") + if channel is None: + raise ValueError("params['channel'] is required.") + bottom_da = params.get("bottom_da", None) + + # --- Compact metadata warnings (key variables) + _missing = [] + + for v in [ + "angle_alongship", + "angle_athwartship", + "beamwidth_alongship", + "beamwidth_athwartship", + "angle_offset_alongship", + "angle_offset_athwartship", + "angle_sensitivity_alongship", + "angle_sensitivity_athwartship", + "sound_speed", + ]: + if v not in ds: + _missing.append(v) + + if _missing: + print( + f"Warning: the following variables are missing for channel '{channel}': " + + ", ".join(_missing) + + ". Defaults will be used where applicable." + ) + + # --- Select main variables + Sv = ds["Sv"].sel(channel=channel).transpose("ping_time", "range_sample") + Depth = ds["depth"].sel(channel=channel).transpose("ping_time", "range_sample") + pt = Sv["ping_time"] + + along = ds.get("angle_alongship") + if along is not None: + if "channel" in along.dims: + along = along.sel(channel=channel) + along = along.transpose("ping_time", "range_sample") + athwt = ds.get("angle_athwartship") + if athwt is not None: + if "channel" in athwt.dims: + athwt = athwt.sel(channel=channel) + athwt = athwt.transpose("ping_time", "range_sample") + + # --- Build local defaults and override with ds-based geometry + defaults = MATECHO_DEFAULTS.copy() + + # Beamwidths (3 dB) from ds_Sv (per channel) + if "beamwidth_alongship" in ds: + bw_al_val = ds["beamwidth_alongship"].sel(channel=channel).values + bw_al_val = float(np.asarray(bw_al_val).item()) + if abs(bw_al_val) > np.pi: # assume degrees + bw_al_val = np.deg2rad(bw_al_val) + defaults["beamwidth_along_3dB_rad"] = bw_al_val + + if "beamwidth_athwartship" in ds: + bw_at_val = ds["beamwidth_athwartship"].sel(channel=channel).values + bw_at_val = float(np.asarray(bw_at_val).item()) + if abs(bw_at_val) > np.pi: # assume degrees + bw_at_val = np.deg2rad(bw_at_val) + defaults["beamwidth_athwart_3dB_rad"] = bw_at_val + + # Beam steering offsets + if "angle_offset_alongship" in ds: + off_al_val = ds["angle_offset_alongship"].sel(channel=channel).values + off_al_val = float(np.asarray(off_al_val).item()) + if abs(off_al_val) > np.pi: # assume degrees + off_al_val = np.deg2rad(off_al_val) + defaults["steer_along_rad"] = off_al_val + + if "angle_offset_athwartship" in ds: + off_at_val = ds["angle_offset_athwartship"].sel(channel=channel).values + off_at_val = float(np.asarray(off_at_val).item()) + if abs(off_at_val) > np.pi: # assume degrees + off_at_val = np.deg2rad(off_at_val) + defaults["steer_athwart_rad"] = off_at_val + + # Angle sensitivities (phase steps) + if "angle_sensitivity_alongship" in ds: + sens_al_val = ds["angle_sensitivity_alongship"].sel(channel=channel).values + defaults["angle_sens_al"] = float(np.asarray(sens_al_val).item()) + + if "angle_sensitivity_athwartship" in ds: + sens_at_val = ds["angle_sensitivity_athwartship"].sel(channel=channel).values + defaults["angle_sens_at"] = float(np.asarray(sens_at_val).item()) + + # Sound speed from data (median), if present + if "sound_speed" in ds: + try: + c_med = float(ds["sound_speed"].median().values) + defaults["SoundSpeed"] = c_med + except Exception: + pass + + # --- Merge defaults and user params + p = { + **defaults, + **{k: v for k, v in params.items() if k not in ("channel", "bottom_da")}, + } + + # --- Compact warnings for optional-but-important metadata + _missing2 = [] + + # navigation / platform variables (checked via fallback) + _nav_vars = { + "heading": ["heading", "Heading"], + "pitch": ["pitch", "Pitch"], + "roll": ["roll", "Roll"], + "heave": ["heave", "Heave", "vertical_offset"], + "dist": ["dist", "Dist", "distance"], + } + for logical_name, candidates in _nav_vars.items(): + if _fallback_1d_aligned(ds, candidates, ds["Sv"].sel(channel=channel)["ping_time"]) is None: + _missing2.append(logical_name) + + # direct dataset variables + _direct_vars = [ + "sound_absorption", + "transducer_depth", + ] + for v in _direct_vars: + if v not in ds: + _missing2.append(v) + + # bottom line if provided + if bottom_da is None: + _missing2.append("bottom_da") + + # pulse parameters from user / defaults + if "pulse_length" not in p: + _missing2.append("pulse_length") + + # emit a single clean warning + if _missing2: + print( + f"[matecho] Warning: optional metadata missing or falling back to defaults " + f"for channel '{channel}': {', '.join(_missing2)}." + ) + + ############## + + c = float(p["SoundSpeed"]) + TS_thr = float(p["TS_threshold"]) + + # --- Nav / platform + heading = _fallback_1d_aligned(ds, ["heading", "Heading"], pt) + pitch = _fallback_1d_aligned(ds, ["pitch", "Pitch"], pt) + roll = _fallback_1d_aligned(ds, ["roll", "Roll"], pt) + heave = _fallback_1d_aligned(ds, ["heave", "Heave", "vertical_offset"], pt) + dist = _fallback_1d_aligned(ds, ["dist", "Dist", "distance"], pt) + n_ping = Sv.sizes["ping_time"] + zeros = np.zeros(n_ping, dtype=float) + heading = heading if heading is not None else zeros + pitch = pitch if pitch is not None else zeros + roll = roll if roll is not None else zeros + heave = heave if heave is not None else zeros + dist = dist if dist is not None else zeros + + # --- Per-channel absorption (optional) + alpha = 0.0 + if "sound_absorption" in ds: + sa = ds["sound_absorption"] + try: + alpha = ( + float(sa.sel(channel=channel).values.item()) + if "channel" in sa.dims + else float(sa.values.item()) + ) + except Exception: + alpha = 0.0 + + # --- Transducer depth (fallback 0) + TD = 0.0 + if "transducer_depth" in ds: + try: + td = ds["transducer_depth"].sel(channel=channel) + TD = float(np.asarray(td).item()) if td.ndim == 0 else float(td.values) + except Exception: + TD = 0.0 + + # --- Range from transducer face (remove TD + heave) + R = Depth - (TD + heave.reshape(-1, 1)) + R = R.clip(min=0.0) + + # --- Vertical step & timing + dstep_da = Depth.isel(ping_time=0).diff("range_sample") # 1D + dstep_arr = np.asarray(dstep_da) + + if np.isfinite(dstep_arr).any(): + dstep = float(np.nanmedian(dstep_arr)) + else: + dstep = np.nan + + if not np.isfinite(dstep) or dstep <= 0: + dstep_da_fb = Depth.diff("range_sample") # 2D + dstep_arr_fb = np.asarray(dstep_da_fb) + if np.isfinite(dstep_arr_fb).any(): + dstep = float(np.nanmedian(dstep_arr_fb)) + else: + dstep = np.nan + + dt = (2.0 * dstep) / c if np.isfinite(dstep) and dstep > 0 else 1e-4 + + T = float(p.get("pulse_length", 1e-3)) + # if "Np" in p and p["Np"] is not None and int(p["Np"]) > 2: + # Np = int(p["Np"]) + # else: + # Np = max(3, int(round(T / max(dt, 1e-6)))) + + # --- TS constants + sv2ts = ( + 10.0 * np.log10(c * T / 2.0) + + float(p["psi_two_way"]) + + 2.0 * float(p["Sa_correction"]) + + 2.0 * float(p["Sa_EK80_nominal"]) + ) + + # --- Index domains & bottom cropping + nb_pings_tot = Sv.sizes["ping_time"] + nb_samples_tot = Sv.sizes["range_sample"] + idx_pings_tot = np.arange(nb_pings_tot, dtype=int) + idx_r_tot = np.arange(nb_samples_tot, dtype=int) + + bb = None + if bottom_da is not None: + bb = bottom_da + if "channel" in bb.dims: + bb = bb.sel(channel=channel) + bb = bb.sel(ping_time=Sv["ping_time"]) + + D = Depth.values # (ping, range_sample) + b = bb.values # (ping,) + idx_bot = np.empty(nb_pings_tot, dtype=int) + for ip in range(nb_pings_tot): + under = np.nonzero(np.isfinite(D[ip]) & (D[ip] >= b[ip]))[0] + if under.size: + idx_bot[ip] = max(under[0] - 1, 0) + else: + idx_bot[ip] = nb_samples_tot - 1 + max_idx = int(np.nanmax(idx_bot)) + idx_r_tot = np.arange(max_idx + 1, dtype=int) + + # --- Blocking in ping dimension + cells_per_ping = max(1, idx_r_tot.size) + block_size = int(min(math.ceil(float(p["block_len"]) / cells_per_ping), idx_pings_tot.size)) + num_ite = int(math.ceil(idx_pings_tot.size / max(block_size, 1))) + + # --- TVG start & helper vectors + ind_range0 = int(p["tvg_start_sample"]) # 3 (EK60) or 1 (EK80) + # Matecho: ts_range = max(delta, ((1:NbSamp)-ind_range0)*delta) + ts_range = np.maximum( + dstep, + (np.arange(nb_samples_tot) + 1 - ind_range0) * dstep, + ) + NechP = max(3, int(round(T / max(dt, 1e-6)))) # samples per pulse + + # --- Output struct (compatible with existing ESP3-style code) + try: + from .esp3 import init_st_struct # if helper exists + except Exception: + + def init_st_struct(): + return { + "TS_comp": [], + "TS_uncomp": [], + "Target_range": [], + "Target_range_disp": [], + "Target_range_min": [], + "Target_range_max": [], + "idx_r": [], + "StandDev_Angles_Minor_Axis": [], + "StandDev_Angles_Major_Axis": [], + "Angle_minor_axis": [], + "Angle_major_axis": [], + "Ping_number": [], + "Time": [], + "nb_valid_targets": 0, + "idx_target_lin": [], + "pulse_env_before_lin": [], + "pulse_env_after_lin": [], + "PulseLength_Normalized_PLDL": [], + "Transmitted_pulse_length": [], + "Heave": [], + "Roll": [], + "Pitch": [], + "Heading": [], + "Dist": [], + } + + out = init_st_struct() + times = Sv["ping_time"].values + + # --- Parameters for selection + MaxAngleOW = float(p["MaxAngleOneWayCompression"]) + MinEchoLen = float(p["MinEchoLength"]) + MaxEchoLen = float(p["MaxEchoLength"]) + MinEchoSpace = float(p["MinEchoSpace"]) + MinEchoDepthM = float(p["MinEchoDepthM"]) + MaxEchoDepthM = float(p["MaxEchoDepthM"]) + max_phase = float(p.get("MaxPhaseDeviation", np.inf)) + dec_tir = 8 # ignore first samples near TX + + # Beam pattern params + bw_al = float(p["beamwidth_along_3dB_rad"]) + bw_at = float(p["beamwidth_athwart_3dB_rad"]) + dec_al = float(p["steer_along_rad"]) + dec_at = float(p["steer_athwart_rad"]) + + # Phase sensitivities + sens_al = float(p.get("angle_sens_al", 1.0)) + sens_at = float(p.get("angle_sens_at", 1.0)) + + for ui in range(num_ite): + start = ui * block_size + stop = min((ui + 1) * block_size, idx_pings_tot.size) + idx_pings = idx_pings_tot[start:stop] + + idx_r = idx_r_tot.copy() + if bb is not None and idx_pings.size > 0: + D_block = Depth.isel(ping_time=idx_pings).values + b_block = bb.isel(ping_time=idx_pings).values + idx_bot = np.empty(idx_pings.size, dtype=int) + for k in range(idx_pings.size): + under = np.nonzero(np.isfinite(D_block[k]) & (D_block[k] >= b_block[k]))[0] + if under.size: + idx_bot[k] = max(under[0] - 1, 0) + else: + idx_bot[k] = nb_samples_tot - 1 + max_idx = int(np.nanmax(idx_bot)) + idx_r = idx_r[idx_r <= max_idx] + + # Subset arrays -> (samples, pings) + Sv_blk = ( + Sv.isel(ping_time=idx_pings, range_sample=idx_r) + .transpose("range_sample", "ping_time") + .values + ) + R_blk = ( + R.isel(ping_time=idx_pings, range_sample=idx_r) + .transpose("range_sample", "ping_time") + .values + ) + # DEPblk = ( + # Depth.isel(ping_time=idx_pings, range_sample=idx_r) + # .transpose("range_sample", "ping_time") + # .values + # ) + tsr_blk = ts_range[idx_r].reshape(-1, 1) + + # Angles: convert to radians here (echopype angles usually in degrees) + along_blk = None + athwt_blk = None + if along is not None: + along_blk = ( + along.isel(ping_time=idx_pings, range_sample=idx_r) + .transpose("range_sample", "ping_time") + .data + ) + along_blk = np.deg2rad(along_blk) + if athwt is not None: + athwt_blk = ( + athwt.isel(ping_time=idx_pings, range_sample=idx_r) + .transpose("range_sample", "ping_time") + .data + ) + athwt_blk = np.deg2rad(athwt_blk) + + # Sv->TSu and Plike + r_eff = np.maximum(tsr_blk, 1e-6) + tsu = Sv_blk + 20.0 * np.log10(r_eff) + sv2ts + Plike = tsu - 40.0 * np.log10(r_eff) - 2.0 * alpha * tsr_blk + + # Beam compensation (one-way approx); if no angles, comp=0 + if along_blk is not None and athwt_blk is not None: + x = 2.0 * ((along_blk - dec_al) / bw_al) + y = 2.0 * ((athwt_blk - dec_at) / bw_at) + one_way_comp = 6.0206 * (x**2 + y**2 - 0.18 * (x**2) * (y**2)) + else: + one_way_comp = np.zeros_like(tsu) + + ts = tsu + one_way_comp + + # Per-ping detection + for j in range(Plike.shape[1]): + zP = Plike[:, j] + zTS = ts[:, j] + zTSU = tsu[:, j] + rcol = R_blk[:, j] + + if not np.isfinite(zP).any(): + continue + + # local maxima of Plike, skipping first samples + valid = np.isfinite(zP) + valid[:dec_tir] = False + loc = np.where(valid[1:-1] & (zP[1:-1] > zP[:-2]) & (zP[1:-1] >= zP[2:]))[0] + 1 + + if loc.size == 0: + continue + + # Matecho-like: process stronger echoes first (Condition 8) + order = sorted(loc.tolist(), key=lambda kk: zTS[kk], reverse=True) + keep = [] + + for k in order: + # TS threshold + if not np.isfinite(zTS[k]) or zTS[k] <= TS_thr: + continue + + # Angle-comp guard: Matecho uses 2 * MaxAngleOW (one-way) + if (zTS[k] - zTSU[k]) > 2.0 * MaxAngleOW: + continue + + # Phase deviation test (Condition 6) + if max_phase < np.inf and along_blk is not None and athwt_blk is not None: + i0_phase = max(k - 2, 0) + i1_phase = min(k + 2, Plike.shape[0] - 1) + + al_win = along_blk[i0_phase : i1_phase + 1, j] + at_win = athwt_blk[i0_phase : i1_phase + 1, j] + + al_steps = al_win * (128.0 / np.pi) * sens_al + at_steps = at_win * (128.0 / np.pi) * sens_at + + if np.nanstd(al_steps) > max_phase or np.nanstd(at_steps) > max_phase: + continue + + # 6-dB width around Plike peak (Condition 7) + base = zP[k] + i0 = k + while i0 > 0 and np.isfinite(zP[i0 - 1]) and (zP[i0 - 1] >= base - 6.0): + i0 -= 1 + i1 = k + nP = zP.size + while i1 + 1 < nP and np.isfinite(zP[i1 + 1]) and (zP[i1 + 1] >= base - 6.0): + i1 += 1 + plen = i1 - i0 + 1 + + # echo length (in samples) bounds + if plen < round(NechP * MinEchoLen) or plen > round(NechP * MaxEchoLen): + continue + + # minimal spacing between echoes in samples (Matecho Condition 8) + if keep and min(abs(k - kk) for kk in keep) < int(round(MinEchoSpace * NechP)): + continue + + # depth band test (surface-referenced) + depth_here = float(rcol[k] + TD + heave[idx_pings[j]]) + if not (MinEchoDepthM <= depth_here <= MaxEchoDepthM): + continue + + # extra guard: do not keep targets below bottom + if bb is not None: + b_here = float(bb.values[idx_pings[j]]) + if depth_here > b_here: + continue + + keep.append(k) + + # Save detection + r_seg = R_blk[i0 : i1 + 1, j] + r_seg = r_seg[np.isfinite(r_seg)] + if r_seg.size == 0: + continue + r_min = float(r_seg.min()) + r_max = float(r_seg.max()) + r_peak = float(rcol[k]) + + out["TS_comp"].append(float(zTS[k])) + out["TS_uncomp"].append(float(zTSU[k])) + out["Target_range"].append(r_peak) + out["Target_range_disp"].append(r_peak + c * T / 4.0) + out["Target_range_min"].append(r_min) + out["Target_range_max"].append(r_max) + + out["idx_r"].append(int(idx_r[k])) + out["Ping_number"].append(int(idx_pings[j])) + out["Time"].append(times[idx_pings[j]]) + out["idx_target_lin"].append(int(idx_pings[j] * nb_samples_tot + idx_r[k])) + + out["pulse_env_before_lin"].append(int(k - i0)) + out["pulse_env_after_lin"].append(int(i1 - k)) + out["PulseLength_Normalized_PLDL"].append(plen / float(NechP)) + out["Transmitted_pulse_length"].append(int(plen)) + + out["StandDev_Angles_Minor_Axis"].append(np.nan) + out["StandDev_Angles_Major_Axis"].append(np.nan) + out["Angle_minor_axis"].append(np.nan) + out["Angle_major_axis"].append(np.nan) + + out["Heave"].append(float(heave[idx_pings[j]])) + out["Roll"].append(float(roll[idx_pings[j]])) + out["Pitch"].append(float(pitch[idx_pings[j]])) + out["Heading"].append(float(heading[idx_pings[j]])) + out["Dist"].append(float(dist[idx_pings[j]])) + + out["nb_valid_targets"] = len(out["TS_comp"]) + return out From d4c74c3b3a7add8601b0c9ae049fae2b6ecb2cc5 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Fri, 19 Dec 2025 22:01:01 +0100 Subject: [PATCH 2/6] Add matecho single-target detection stub and API scaffold - Remove preliminary Matecho algorithm logic - Introduce a stub implementation for single-target detection - Define and validate the public API (inputs, params, outputs) - Return a zero-target xr.Dataset with a stable schema - Add parameter and dataset validation checks - Prepare ground for future algorithmic implementation --- echopype/mask/api.py | 31 +- .../single_target_detection/detect_esp3.py | 625 ----------------- .../single_target_detection/detect_matecho.py | 646 ++++-------------- echopype/tests/mask/test_mask.py | 159 +++++ 4 files changed, 307 insertions(+), 1154 deletions(-) delete mode 100644 echopype/mask/single_target_detection/detect_esp3.py diff --git a/echopype/mask/api.py b/echopype/mask/api.py index 71a6d241c..1f71f49ec 100644 --- a/echopype/mask/api.py +++ b/echopype/mask/api.py @@ -14,7 +14,6 @@ from echopype.mask.seafloor_detection.bottom_blackwell import bottom_blackwell # for single_target_detection -from echopype.mask.single_target_detection.detect_esp3 import detect_esp3 from echopype.mask.single_target_detection.detect_matecho import detect_matecho from ..utils.io import validate_source @@ -807,7 +806,6 @@ def detect_shoal( # Registry of supported methods for single_target_detection METHODS_SINGLE_TARGET = { - "esp3": detect_esp3, "matecho": detect_matecho, } @@ -815,23 +813,31 @@ def detect_shoal( def detect_single_targets( ds: xr.Dataset, method: str, - params: dict | None = None, -): + params: dict, +) -> xr.Dataset: """ Run single-target detection using the selected method. Parameters ---------- ds : xr.Dataset - Sv dataset (must contain what's needed by the method). + Acoustic dataset containing the fields required by the selected method. + Typical dimensions include ``ping_time`` and ``range_sample`` (or ``depth``). + Depending on the method, this may require TS, Sv, and optionally split-beam + angles (e.g. alongship / athwartship angles). method : str - Name of the detection method to use (e.g., "esp3"). - params : dict, optional - Parameters for the detection function (method-specific). + Name of the detection method to use (e.g., ``"matecho"``, ``"esp3"``, ...). + params : dict + Method-specific parameters. This argument is required and no defaults + are assumed. Returns ------- - Whatever the method returns (e.g., dict or xr.Dataset). + xr.Dataset + Per-target results with dimension ``target`` and coordinates + ``ping_time`` and ``ping_number``. Variables include compensated and + uncompensated TS, range metrics, and optional + navigation/attitude fields when available. """ if method not in METHODS_SINGLE_TARGET: raise ValueError(f"Unsupported single-target method: {method}") @@ -839,4 +845,9 @@ def detect_single_targets( if params is None: raise ValueError("No parameters given.") - return METHODS_SINGLE_TARGET[method](ds, params) + out = METHODS_SINGLE_TARGET[method](ds, params) + + if not isinstance(out, xr.Dataset) or "target" not in out.dims: + raise TypeError(f"{method} must return an xr.Dataset with a 'target' dimension.") + + return out diff --git a/echopype/mask/single_target_detection/detect_esp3.py b/echopype/mask/single_target_detection/detect_esp3.py deleted file mode 100644 index af5a82092..000000000 --- a/echopype/mask/single_target_detection/detect_esp3.py +++ /dev/null @@ -1,625 +0,0 @@ -import math - -import numpy as np -import xarray as xr - -DEFAULTS = { - "SoundSpeed": 1500.0, - "TS_threshold": -50.0, # dB - "PLDL": 6.0, # dB below peak - "MinNormPL": 0.7, - "MaxNormPL": 1.5, - "DataType": "CW", # CW only for now - "block_len": 1e7 / 3, # ~3.33M cells per block - "MaxBeamComp": 4.0, # dB (like default in MATLAB) - "MaxStdMinAxisAngle": 0.6, # rad (≈ 34°) - "MaxStdMajAxisAngle": 0.6, # rad -} - - -def init_st_struct(): - """Return a dict mirroring the MATLAB single_targets_tot struct.""" - return { - "TS_comp": [], - "TS_uncomp": [], - "Target_range": [], - "Target_range_disp": [], - "Target_range_min": [], - "Target_range_max": [], - "idx_r": [], - "StandDev_Angles_Minor_Axis": [], - "StandDev_Angles_Major_Axis": [], - "Angle_minor_axis": [], - "Angle_major_axis": [], - "Ping_number": [], - "Time": [], - "nb_valid_targets": 0, - "idx_target_lin": [], - "pulse_env_before_lin": [], - "pulse_env_after_lin": [], - "PulseLength_Normalized_PLDL": [], - "Transmitted_pulse_length": [], - "Heave": [], - "Roll": [], - "Pitch": [], - "Heading": [], - "Dist": [], - } - - -def _validate_params(p: dict): - if p.get("DataType", "CW") not in ("CW", "FM"): - raise ValueError("DataType must be 'CW' or 'FM'.") - if not (-120.0 <= float(p["TS_threshold"]) <= -20.0): - raise ValueError("TS_threshold must be in [-120, -20] dB.") - if not (1.0 <= float(p["PLDL"]) <= 30.0): - raise ValueError("PLDL must be in [1, 30] dB.") - if not (0.0 <= float(p["MinNormPL"]) <= 10.0): - raise ValueError("MinNormPL must be in [0, 10].") - if not (0.0 <= float(p["MaxNormPL"]) <= 10.0): - raise ValueError("MaxNormPL must be in [0, 10].") - if float(p["block_len"]) <= 0: - raise ValueError("block_len must be > 0.") - if not (0.0 <= float(p["MaxBeamComp"]) <= 18.0): - raise ValueError("MaxBeamComp must be in [0, 18] dB.") - if not (0.0 <= float(p["MaxStdMinAxisAngle"]) <= 45.0): - # keep same numeric limits as MATLAB; here we assume radians - raise ValueError("MaxStdMinAxisAngle must be in [0, 45].") - if not (0.0 <= float(p["MaxStdMajAxisAngle"]) <= 45.0): - raise ValueError("MaxStdMajAxisAngle must be in [0, 45].") - - -def _pull_nav(ds, name_candidates, ping_time_ref): - """Return a 1D float array aligned to ping_time_ref (or None).""" - for nm in name_candidates: - if nm in ds: - da = ds[nm] - for tdim in ("ping_time", "time", "time1", "time2"): - if tdim in da.dims: - if tdim != "ping_time": - da = da.rename({tdim: "ping_time"}) - da = da.reindex( - ping_time=ping_time_ref, - method="nearest", - tolerance=np.timedelta64(500, "ms"), - ) - return np.asarray(da.values, dtype=float) - if da.ndim == 1 and da.size == ping_time_ref.size: - return np.asarray(da.values, dtype=float) - return None - - -def simrad_beam_compensation( - along_angles_rad: np.ndarray, - athwart_angles_rad: np.ndarray, - bw_along_rad: float, - bw_athwart_rad: float, -) -> np.ndarray: - """ - Simrad beam pattern compensation [dB], same polynomial as in ESP3/Matecho. - - Parameters - ---------- - along_angles_rad : array - Alongship split-beam angles in radians. - athwart_angles_rad : array - Athwartship split-beam angles in radians. - bw_along_rad : float - 3 dB beamwidth alongship in radians. - bw_athwart_rad : float - 3 dB beamwidth athwartship in radians. - - Returns - ------- - comp : array - Beam compensation in dB (one-way). - """ - # same form as used in Matecho & ESP3: - # x = 2 * (phi_along / BW_along) - # y = 2 * (phi_athwart / BW_athwart) - # comp = 6.0206 * (x^2 + y^2 - 0.18 * x^2 * y^2) - x = 2.0 * (along_angles_rad / bw_along_rad) - y = 2.0 * (athwart_angles_rad / bw_athwart_rad) - comp = 6.0206 * (x**2 + y**2 - 0.18 * (x**2) * (y**2)) - return comp - - -def detect_esp3(ds: xr.Dataset, params: dict): - """ - ESP3-style single-target detection (CW branch) translated to Python. - - Parameters - ---------- - ds : xarray.Dataset - Calibrated Sv dataset (channel, ping_time, range_sample). - Must contain Sv (TS-like), depth, and optionally angles. - - params : dict - Must include: - - "channel": channel name - May include: - - "bottom_da" : bottom depth vs ping_time (DataArray) - - TS_threshold, PLDL, MinNormPL, MaxNormPL, SoundSpeed, - MaxBeamComp, MaxStdMinAxisAngle, MaxStdMajAxisAngle, - block_len, DataType - - Returns - ------- - out : dict - ESP3-like structure with TS_comp, TS_uncomp, Target_range, etc. - """ - if params is None: - params = {} - channel = params.get("channel") - if channel is None: - raise ValueError("params['channel'] is required.") - bottom_da = params.get("bottom_da", None) - - # Merge defaults - p = { - **DEFAULTS, - **{k: v for k, v in params.items() if k not in ("channel", "bottom_da")}, - } - p["DataType"] = p.get("DataType", "CW") - _validate_params(p) - - # ------------------------------------------------------------------ - # 1) Select data - # ------------------------------------------------------------------ - Sv = ds["Sv"].sel(channel=channel).transpose("ping_time", "range_sample") - Depth = ds["depth"].sel(channel=channel).transpose("ping_time", "range_sample") - - along = ds.get("angle_alongship") - if along is not None: - if "channel" in along.dims: - along = along.sel(channel=channel) - along = along.transpose("ping_time", "range_sample") - - athwt = ds.get("angle_athwartship") - if athwt is not None: - if "channel" in athwt.dims: - athwt = athwt.sel(channel=channel) - athwt = athwt.transpose("ping_time", "range_sample") - - nb_pings_tot = Sv.sizes["ping_time"] - nb_samples_tot = Sv.sizes["range_sample"] - idx_pings_tot = np.arange(nb_pings_tot, dtype=int) - idx_r_tot = np.arange(nb_samples_tot, dtype=int) - - # Optional bottom cropping (like idx_r_tot>max(idx_bot) removal) - if bottom_da is not None: - bb_global = bottom_da - if "channel" in bb_global.dims: - bb_global = bb_global.sel(channel=channel) - bb_global = bb_global.sel(ping_time=Sv["ping_time"]) - - D_all = Depth.values # (ping, range_sample) - b_all = bb_global.values # (ping,) - idx_bot = np.empty(nb_pings_tot, dtype=int) - for ip in range(nb_pings_tot): - under = np.nonzero(np.isfinite(D_all[ip]) & (D_all[ip] >= b_all[ip]))[0] - idx_bot[ip] = under[0] if under.size else (nb_samples_tot - 1) - max_idx = int(np.nanmax(idx_bot)) - idx_r_tot = np.arange(max_idx + 1, dtype=int) - - # Region/bad-data mask (placeholder: none => all False) - mask_inter_tot = np.zeros((idx_r_tot.size, idx_pings_tot.size), dtype=bool) - - # Block size - cells_per_ping = max(1, idx_r_tot.size) - block_size = int(min(math.ceil(float(p["block_len"]) / cells_per_ping), idx_pings_tot.size)) - num_ite = int(math.ceil(idx_pings_tot.size / max(block_size, 1))) - - # ------------------------------------------------------------------ - # 2) Nav / attitude - # ------------------------------------------------------------------ - pt_ref = Sv["ping_time"] - heading_arr = _pull_nav(ds, ["heading", "Heading"], pt_ref) - pitch_arr = _pull_nav(ds, ["pitch", "Pitch"], pt_ref) - roll_arr = _pull_nav(ds, ["roll", "Roll"], pt_ref) - heave_arr = _pull_nav(ds, ["heave", "Heave", "vertical_offset"], pt_ref) - dist_arr = _pull_nav(ds, ["dist", "Dist", "distance"], pt_ref) - - def _fallback(a, n): - return a if a is not None else np.zeros(n, dtype=float) - - heading = _fallback(heading_arr, nb_pings_tot) - pitch = _fallback(pitch_arr, nb_pings_tot) - roll = _fallback(roll_arr, nb_pings_tot) - heave = _fallback(heave_arr, nb_pings_tot) - dist = _fallback(dist_arr, nb_pings_tot) - - times = Sv["ping_time"].values - - # ------------------------------------------------------------------ - # 3) Prepare constants (alpha, c, T, Np) - # ------------------------------------------------------------------ - alpha = 0.0 - if "sound_absorption" in ds: - sa = ds["sound_absorption"] - try: - alpha = ( - float(sa.sel(channel=channel).values.item()) - if "channel" in sa.dims - else float(sa.values.item()) - ) - except Exception: - alpha = 0.0 - - c = float(p["SoundSpeed"]) - - # Pulse length & number of samples - if "Np" in p and p["Np"] is not None and int(p["Np"]) > 2: - Np = int(p["Np"]) - T = float(p.get("pulse_length", 1e-3)) - else: - dstep_da = Depth.diff("range_sample").median(skipna=True) - dstep = ( - float(dstep_da.values.item()) - if getattr(dstep_da, "size", 1) == 1 - else float(dstep_da.values) - ) - dt = (2.0 * dstep) / c if dstep > 0 else 1e-4 - T = float(p.get("pulse_length", 1e-3)) - Np = max(3, int(round(T / dt))) - - TS_thr = float(p["TS_threshold"]) - PLDL = float(p["PLDL"]) - min_len = max(1, int(Np * float(p["MinNormPL"]))) - max_len = max(1, int(math.ceil(Np * float(p["MaxNormPL"])))) - MaxBeamComp = float(p["MaxBeamComp"]) - MaxStdMin = float(p["MaxStdMinAxisAngle"]) - MaxStdMaj = float(p["MaxStdMajAxisAngle"]) - - # Beamwidths (taken from metadata if present; else fallback) - bw_along_rad = None - bw_athwart_rad = None - if "beamwidth_alongship" in ds: - bw_al = ds["beamwidth_alongship"].sel(channel=channel).values - bw_al = float(np.asarray(bw_al).item()) - if abs(bw_al) > np.pi: # assume deg - bw_al = np.deg2rad(bw_al) - bw_along_rad = bw_al - if "beamwidth_athwartship" in ds: - bw_at = ds["beamwidth_athwartship"].sel(channel=channel).values - bw_at = float(np.asarray(bw_at).item()) - if abs(bw_at) > np.pi: - bw_at = np.deg2rad(bw_at) - bw_athwart_rad = bw_at - - # Fallback beamwidths if missing - if bw_along_rad is None: - bw_along_rad = np.deg2rad(7.0) - if bw_athwart_rad is None: - bw_athwart_rad = np.deg2rad(7.0) - - # ------------------------------------------------------------------ - # 4) Output struct (like single_targets_tot) - # ------------------------------------------------------------------ - out = init_st_struct() - - # Bottom DA per-ping, sliced once - bb = None - if bottom_da is not None: - bb = bottom_da - if "channel" in bb.dims: - bb = bb.sel(channel=channel) - bb = bb.sel(ping_time=Sv["ping_time"]) - - # ------------------------------------------------------------------ - # 5) Block loop - # ------------------------------------------------------------------ - for ui in range(num_ite): - start = ui * block_size - stop = min((ui + 1) * block_size, idx_pings_tot.size) - idx_pings = idx_pings_tot[start:stop] - - # Per-block rows (respect bottom) - idx_r = idx_r_tot.copy() - if bb is not None and idx_pings.size > 0: - D_block = Depth.isel(ping_time=idx_pings).values - b_block = bb.isel(ping_time=idx_pings).values - idx_bot = np.empty(idx_pings.size, dtype=int) - for k in range(idx_pings.size): - under = np.nonzero(np.isfinite(D_block[k]) & (D_block[k] >= b_block[k]))[0] - idx_bot[k] = under[0] if under.size else (nb_samples_tot - 1) - max_idx = int(np.nanmax(idx_bot)) - idx_r = idx_r[idx_r <= max_idx] - - # Mask (regions) - mask = np.zeros((idx_r.size, idx_pings.size), dtype=bool) - if mask_inter_tot.size: - r_pos = np.searchsorted(idx_r_tot, idx_r) - p_pos = np.searchsorted(idx_pings_tot, idx_pings) - mask |= mask_inter_tot[np.ix_(r_pos, p_pos)] - - # Extract submatrices (TS ~ Sv here) - TS = ( - Sv.isel(ping_time=idx_pings, range_sample=idx_r) - .transpose("range_sample", "ping_time") - .values.copy() - ) - DEP = ( - Depth.isel(ping_time=idx_pings, range_sample=idx_r) - .transpose("range_sample", "ping_time") - .values - ) - - nb_samples, nb_pings = TS.shape - - # Angles sub-block (in radians) - along_block = None - athwt_block = None - if along is not None: - along_block = ( - along.isel(ping_time=idx_pings, range_sample=idx_r) - .transpose("range_sample", "ping_time") - .data - ) - along_block = np.deg2rad(along_block) - if athwt is not None: - athwt_block = ( - athwt.isel(ping_time=idx_pings, range_sample=idx_r) - .transpose("range_sample", "ping_time") - .data - ) - athwt_block = np.deg2rad(athwt_block) - - # Under-bottom mask - if bb is not None and idx_pings.size > 0: - bcol = bb.isel(ping_time=idx_pings).values.reshape(1, -1) - mask |= DEP >= bcol - - TS[mask] = -999.0 - if not np.any(TS > -999.0): - continue - - # Remove trailing all-masked rows - valid_rows = np.where((TS > -999.0).any(axis=1))[0] - last_row = valid_rows.max() - row_sel = slice(0, last_row + 1) - TS = TS[row_sel, :] - DEP = DEP[row_sel, :] - if along_block is not None: - along_block = along_block[row_sel, :] - if athwt_block is not None: - athwt_block = athwt_block[row_sel, :] - idx_r = idx_r[row_sel] - nb_samples, nb_pings = TS.shape - - # TVG matrix and Power (as in MATLAB) - r_eff = DEP - c * T / 4.0 - TVG_mat = np.where( - r_eff > 0.0, - 40.0 * np.log10(r_eff) + 2.0 * alpha * r_eff, - np.nan, - ) - Power = TS - TVG_mat - - # Peak detection on smoothed TS (or Power), like ESP3 - peak_calc = "TS" # ESP3 uses 'TS' by default in your snippet - if peak_calc.lower() == "power": - base_mat = Power - else: - base_mat = TS - - # moving average over Np/2 samples in linear domain - win = max(int(np.floor(Np / 2)), 1) - kernel = np.ones(win, dtype=float) / float(win) - - base_lin = 10.0 ** (base_mat / 10.0) - smooth_lin = np.apply_along_axis( - lambda col: np.convolve(col, kernel, mode="same"), - axis=0, - arr=base_lin, - ) - peak_mat = 10.0 * np.log10(smooth_lin) - peak_mat[TS <= -999.0] = -999.0 - - # --- Find peaks with MinSeparation ≈ Np (along range axis) --- - idx_peaks_bool = np.zeros_like(peak_mat, dtype=bool) - for jp in range(nb_pings): - z = peak_mat[:, jp] - valid = z > -999.0 - if np.count_nonzero(valid) < 3: - continue - cand = np.where(valid[1:-1] & (z[1:-1] > z[:-2]) & (z[1:-1] >= z[2:]))[0] + 1 - - # MinSeparation = Np - last = -(10**9) - keep_idx = [] - for k in cand: - if k - last >= Np: - keep_idx.append(k) - last = k - if keep_idx: - idx_peaks_bool[np.array(keep_idx, dtype=int), jp] = True - - idx_peaks_bool[TS <= -999.0] = False - - # linear indices of peaks (like idx_peaks_lin) - i_peaks_lin, j_peaks_lin = np.where(idx_peaks_bool) - nb_peaks = len(i_peaks_lin) - if nb_peaks == 0: - continue - - # Level of local maxima minus PLDL - peak_vals = peak_mat[idx_peaks_bool] # 1D (same order as i_peaks_lin) - pulse_level = peak_vals - PLDL - - # To mimic MATLAB's pulse-env scan, we will grow left/right until - # we fall below pulse_level or reach max_len. - pulse_env_before = np.zeros(nb_peaks, dtype=int) - pulse_env_after = np.zeros(nb_peaks, dtype=int) - - for idx_peak in range(nb_peaks): - ip = i_peaks_lin[idx_peak] - jp = j_peaks_lin[idx_peak] - thr = pulse_level[idx_peak] - zcol = peak_mat[:, jp] - - # grow up (before) - left = 0 - i = ip - 1 - while i >= 0 and left < max_len and np.isfinite(zcol[i]) and zcol[i] >= thr: - left += 1 - i -= 1 - - # grow down (after) - right = 0 - i = ip + 1 - nP = zcol.size - while i < nP and right < max_len and np.isfinite(zcol[i]) and zcol[i] >= thr: - right += 1 - i += 1 - - pulse_env_before[idx_peak] = left - pulse_env_after[idx_peak] = right - - pulse_length_lin = pulse_env_before + pulse_env_after + 1 - # good pulses according to [MinNormPL, MaxNormPL] - good_pulses = (pulse_length_lin >= min_len) & (pulse_length_lin <= max_len) - - # Filter peaks to targets - i_targets = i_peaks_lin[good_pulses] - j_targets = j_peaks_lin[good_pulses] - pulse_before = pulse_env_before[good_pulses] - pulse_after = pulse_env_after[good_pulses] - pulse_len = pulse_length_lin[good_pulses] - nb_targets = len(i_targets) - if nb_targets == 0: - continue - - # ------------------------------------------------------------------ - # For each target, compute: - # - samples_targets_power / range / angles - # - target_range (weighted by Power) - # - beam comp at peak (target_comp) - # - TS_uncomp = target_peak_power + TVG(target_range) - # - TS_comp = TS_uncomp + target_comp - # - apply angle std filter and MaxBeamComp & TS_threshold - # ------------------------------------------------------------------ - for it in range(nb_targets): - ip = int(i_targets[it]) - jp = int(j_targets[it]) - left = int(pulse_before[it]) - right = int(pulse_after[it]) - - i0 = max(ip - left, 0) - i1 = min(ip + right, nb_samples - 1) - seg_slice = slice(i0, i1 + 1) - - pow_seg = Power[seg_slice, jp] - dep_seg = DEP[seg_slice, jp] - - # If all NaNs or -inf => skip - if not np.isfinite(pow_seg).any(): - continue - - # angles segment - if along_block is not None and athwt_block is not None: - al_seg = along_block[seg_slice, jp] - at_seg = athwt_block[seg_slice, jp] - else: - al_seg = np.full_like(pow_seg, np.nan, dtype=float) - at_seg = np.full_like(pow_seg, np.nan, dtype=float) - - # std + mean of angles (for filtering & output) - std_al = float(np.nanstd(al_seg)) - std_at = float(np.nanstd(at_seg)) - phi_al = float(np.nanmean(al_seg)) - phi_at = float(np.nanmean(at_seg)) - - # Angle stability filter (CW only) - if p["DataType"] == "CW": - if std_al > MaxStdMin or std_at > MaxStdMaj: - continue - - # Power-weighted range (like target_range) - w = 10.0 ** (pow_seg / 10.0) - if not np.isfinite(w).any() or np.nansum(w) <= 0: - continue - target_range = float(np.nansum(w * dep_seg) / np.nansum(w) - c * T / 4.0) - if target_range < 0: - target_range = 0.0 - - target_range_min = float(np.nanmin(dep_seg)) - target_range_max = float(np.nanmax(dep_seg)) - - # Peak power within pulse window - # (we take max of pow_seg and its index) - idx_loc_peak = int(np.nanargmax(pow_seg)) - target_peak_power = float(pow_seg[idx_loc_peak]) - if not np.isfinite(target_peak_power): - continue - - # Beam comp at peak sample - if along_block is not None and athwt_block is not None: - al_peak = al_seg[idx_loc_peak] - at_peak = at_seg[idx_loc_peak] - if np.isfinite(al_peak) and np.isfinite(at_peak): - target_comp = float( - simrad_beam_compensation( - np.array([al_peak]), - np.array([at_peak]), - bw_along_rad, - bw_athwart_rad, - )[0] - ) - else: - target_comp = 0.0 - else: - target_comp = 0.0 - - # TVG at target_range (like ESP3: TVG = 40 log10(R) + 2 alpha R) - if target_range > 0: - TVG = 40.0 * np.log10(target_range) + 2.0 * alpha * target_range - else: - TVG = np.nan - - if not np.isfinite(TVG): - continue - - target_TS_uncomp = target_peak_power + TVG - target_TS_comp = target_TS_uncomp + target_comp - - # Filters on TS_comp and beam comp - if target_TS_comp <= TS_thr: - continue - if abs(target_comp) > MaxBeamComp: - continue - - # Store detection - ping_glob = int(idx_pings[jp]) - range_idx_glob = int(idx_r[ip]) - - out["TS_comp"].append(target_TS_comp) - out["TS_uncomp"].append(target_TS_uncomp) - out["Target_range"].append(target_range) - out["Target_range_disp"].append(target_range + c * T / 4.0) - out["Target_range_min"].append(target_range_min) - out["Target_range_max"].append(target_range_max) - - out["idx_r"].append(range_idx_glob) - out["Ping_number"].append(ping_glob) - out["Time"].append(times[ping_glob]) - out["idx_target_lin"].append(int(ping_glob * nb_samples_tot + range_idx_glob)) - - out["pulse_env_before_lin"].append(int(left)) - out["pulse_env_after_lin"].append(int(right)) - out["PulseLength_Normalized_PLDL"].append(pulse_len[it] / float(Np)) - out["Transmitted_pulse_length"].append(int(pulse_len[it])) - - out["StandDev_Angles_Minor_Axis"].append(std_al) - out["StandDev_Angles_Major_Axis"].append(std_at) - out["Angle_minor_axis"].append(phi_al) - out["Angle_major_axis"].append(phi_at) - - out["Heave"].append(float(heave[ping_glob])) - out["Roll"].append(float(roll[ping_glob])) - out["Pitch"].append(float(pitch[ping_glob])) - out["Heading"].append(float(heading[ping_glob])) - out["Dist"].append(float(dist[ping_glob])) - - out["nb_valid_targets"] = len(out["TS_comp"]) - return out diff --git a/echopype/mask/single_target_detection/detect_matecho.py b/echopype/mask/single_target_detection/detect_matecho.py index aed7b277e..026285d04 100644 --- a/echopype/mask/single_target_detection/detect_matecho.py +++ b/echopype/mask/single_target_detection/detect_matecho.py @@ -1,10 +1,9 @@ -# matecho.py -import math +import warnings import numpy as np import xarray as xr -MATECHO_DEFAULTS = { +PARAM_DEFAULTS = { "SoundSpeed": 1500.0, "TS_threshold": -50.0, # dB "MaxAngleOneWayCompression": 6.0, # dB (one-way; test uses 2x) @@ -31,66 +30,110 @@ } -def _fallback_1d_aligned(ds, names, ping_time): +def _matecho_struct_to_dataset(out: dict, *, channel: str | None = None) -> xr.Dataset: """ - Try to get a 1D variable aligned to ping_time from ds, using - several possible names and time dims, with nearest-neighbour - reindexing and 500 ms tolerance. + Convert struct-of-lists into an xr.Dataset with dim 'target'. """ - for nm in names: - if nm in ds: - da = ds[nm] - for tdim in ("ping_time", "time", "time1", "time2"): - if tdim in da.dims: - if tdim != "ping_time": - da = da.rename({tdim: "ping_time"}) - da = da.reindex( - ping_time=ping_time, - method="nearest", - tolerance=np.timedelta64(500, "ms"), - ) - return np.asarray(da.values, dtype=float) - if da.ndim == 1 and da.size == ping_time.size: - return np.asarray(da.values, dtype=float) - return None + n = int(out.get("nb_valid_targets", len(out.get("TS_comp", [])))) + target = np.arange(n, dtype=int) + + def _arr(key, dtype=float): + vals = out.get(key, []) + if n == 0: + return np.asarray([], dtype=dtype) + a = np.asarray(vals, dtype=dtype) + if a.size != n: + raise ValueError(f"Output field '{key}' has length {a.size}, expected {n}.") + return a + + ds_out = xr.Dataset( + data_vars=dict( + TS_comp=("target", _arr("TS_comp", float), {"units": "dB"}), + TS_uncomp=("target", _arr("TS_uncomp", float), {"units": "dB"}), + target_range=("target", _arr("Target_range", float), {"units": "m"}), + target_range_disp=("target", _arr("Target_range_disp", float), {"units": "m"}), + target_range_min=("target", _arr("Target_range_min", float), {"units": "m"}), + target_range_max=("target", _arr("Target_range_max", float), {"units": "m"}), + idx_r=("target", _arr("idx_r", int)), + idx_target_lin=("target", _arr("idx_target_lin", int)), + pulse_env_before_lin=("target", _arr("pulse_env_before_lin", int)), + pulse_env_after_lin=("target", _arr("pulse_env_after_lin", int)), + pulse_length_normalized_pldl=("target", _arr("PulseLength_Normalized_PLDL", float)), + transmitted_pulse_length=("target", _arr("Transmitted_pulse_length", int)), + angle_minor_axis=("target", _arr("Angle_minor_axis", float), {"units": "rad"}), + angle_major_axis=("target", _arr("Angle_major_axis", float), {"units": "rad"}), + std_angle_minor_axis=( + "target", + _arr("StandDev_Angles_Minor_Axis", float), + {"units": "phase_steps"}, + ), + std_angle_major_axis=( + "target", + _arr("StandDev_Angles_Major_Axis", float), + {"units": "phase_steps"}, + ), + heave=("target", _arr("Heave", float), {"units": "m"}), + roll=("target", _arr("Roll", float), {"units": "rad"}), + pitch=("target", _arr("Pitch", float), {"units": "rad"}), + heading=("target", _arr("Heading", float), {"units": "rad"}), + dist=("target", _arr("Dist", float), {"units": "m"}), + ), + coords=dict( + target=target, + ping_number=("target", _arr("Ping_number", int)), + ping_time=("target", _arr("Time", "datetime64[ns]")), + ), + attrs=dict( + method="matecho", + channel=str(channel) if channel is not None else "", + nb_valid_targets=n, + ), + ) + + if channel is not None: + # scalar coordinate + ds_out = ds_out.assign_coords(channel=np.asarray(channel)) + + return ds_out -def detect_matecho(ds: xr.Dataset, params: dict): +def detect_matecho(ds: xr.Dataset, params: dict) -> xr.Dataset: """ - Matecho-inspired single-target detector (CW path only). + Matecho-inspired single-target detector (CW only). - Parameters - ---------- - ds : xarray.Dataset - Calibrated Sv dataset (e.g., from echopype.calibrate.compute_Sv), - with dimensions (channel, ping_time, range_sample) and variables: - Sv, depth, angle_alongship, angle_athwartship, and instrument - metadata (beamwidths, offsets, angle sensitivities). - params : dict - Parameters including at least: - - "channel": channel name (string) - Optional: - - "bottom_da": xarray.DataArray of bottom depth vs ping_time - - Matecho parameters: TS_threshold, MaxAngleOneWayCompression, - MaxPhaseDeviation, MinEchoLength, MaxEchoLength, MinEchoSpace, - MinEchoDepthM, MaxEchoDepthM, pulse_length, tvg_start_sample, etc. + This is a placeholder implementation that demonstrates: + - required signature + - required parameter checks + - required output structure (xr.Dataset with dim 'target') - Returns - ------- - out : dict - Matecho-like single-target "struct" (Python dict), with keys such as - TS_comp, TS_uncomp, Target_range, Ping_number, Time, etc. + No detection is performed. """ if params is None: raise ValueError("params is required.") + channel = params.get("channel") if channel is None: raise ValueError("params['channel'] is required.") - bottom_da = params.get("bottom_da", None) - # --- Compact metadata warnings (key variables) - _missing = [] + var_name = params.get("var_name", "Sv") + + if var_name not in ds: + raise ValueError(f"var_name '{var_name}' not found in input dataset.") + da = ds[var_name] + + required_dims = {"channel", "ping_time", "range_sample"} + if not required_dims.issubset(set(da.dims)): + raise ValueError(f"{var_name} must have dims {sorted(required_dims)}. Got: {da.dims}.") + + # channel must exist in Sv(channel, ping_time, range_sample) + if "channel" not in da.coords: + raise ValueError(f"{var_name} has no 'channel' coordinate.") + if channel not in da["channel"].values: + raise ValueError(f"Channel '{channel}' not found in {var_name}.") + + # --- optional metadata warnings (kept for API clarity) + _missing = [] for v in [ "angle_alongship", "angle_athwartship", @@ -106,479 +149,44 @@ def detect_matecho(ds: xr.Dataset, params: dict): _missing.append(v) if _missing: - print( - f"Warning: the following variables are missing for channel '{channel}': " + warnings.warn( + f"The following variables are missing for channel '{channel}': " + ", ".join(_missing) - + ". Defaults will be used where applicable." + + ". Defaults will be used where applicable.", + UserWarning, + stacklevel=2, ) - # --- Select main variables - Sv = ds["Sv"].sel(channel=channel).transpose("ping_time", "range_sample") - Depth = ds["depth"].sel(channel=channel).transpose("ping_time", "range_sample") - pt = Sv["ping_time"] - - along = ds.get("angle_alongship") - if along is not None: - if "channel" in along.dims: - along = along.sel(channel=channel) - along = along.transpose("ping_time", "range_sample") - athwt = ds.get("angle_athwartship") - if athwt is not None: - if "channel" in athwt.dims: - athwt = athwt.sel(channel=channel) - athwt = athwt.transpose("ping_time", "range_sample") - - # --- Build local defaults and override with ds-based geometry - defaults = MATECHO_DEFAULTS.copy() - - # Beamwidths (3 dB) from ds_Sv (per channel) - if "beamwidth_alongship" in ds: - bw_al_val = ds["beamwidth_alongship"].sel(channel=channel).values - bw_al_val = float(np.asarray(bw_al_val).item()) - if abs(bw_al_val) > np.pi: # assume degrees - bw_al_val = np.deg2rad(bw_al_val) - defaults["beamwidth_along_3dB_rad"] = bw_al_val - - if "beamwidth_athwartship" in ds: - bw_at_val = ds["beamwidth_athwartship"].sel(channel=channel).values - bw_at_val = float(np.asarray(bw_at_val).item()) - if abs(bw_at_val) > np.pi: # assume degrees - bw_at_val = np.deg2rad(bw_at_val) - defaults["beamwidth_athwart_3dB_rad"] = bw_at_val - - # Beam steering offsets - if "angle_offset_alongship" in ds: - off_al_val = ds["angle_offset_alongship"].sel(channel=channel).values - off_al_val = float(np.asarray(off_al_val).item()) - if abs(off_al_val) > np.pi: # assume degrees - off_al_val = np.deg2rad(off_al_val) - defaults["steer_along_rad"] = off_al_val - - if "angle_offset_athwartship" in ds: - off_at_val = ds["angle_offset_athwartship"].sel(channel=channel).values - off_at_val = float(np.asarray(off_at_val).item()) - if abs(off_at_val) > np.pi: # assume degrees - off_at_val = np.deg2rad(off_at_val) - defaults["steer_athwart_rad"] = off_at_val - - # Angle sensitivities (phase steps) - if "angle_sensitivity_alongship" in ds: - sens_al_val = ds["angle_sensitivity_alongship"].sel(channel=channel).values - defaults["angle_sens_al"] = float(np.asarray(sens_al_val).item()) - - if "angle_sensitivity_athwartship" in ds: - sens_at_val = ds["angle_sensitivity_athwartship"].sel(channel=channel).values - defaults["angle_sens_at"] = float(np.asarray(sens_at_val).item()) - - # Sound speed from data (median), if present - if "sound_speed" in ds: - try: - c_med = float(ds["sound_speed"].median().values) - defaults["SoundSpeed"] = c_med - except Exception: - pass - - # --- Merge defaults and user params - p = { - **defaults, - **{k: v for k, v in params.items() if k not in ("channel", "bottom_da")}, + # ------------------------------------------------------------------ + # ## algo ## + # Intentionally empty: no detection performed + # ------------------------------------------------------------------ + + out = { + "TS_comp": [], + "TS_uncomp": [], + "Target_range": [], + "Target_range_disp": [], + "Target_range_min": [], + "Target_range_max": [], + "idx_r": [], + "StandDev_Angles_Minor_Axis": [], + "StandDev_Angles_Major_Axis": [], + "Angle_minor_axis": [], + "Angle_major_axis": [], + "Ping_number": [], + "Time": [], + "idx_target_lin": [], + "pulse_env_before_lin": [], + "pulse_env_after_lin": [], + "PulseLength_Normalized_PLDL": [], + "Transmitted_pulse_length": [], + "Heave": [], + "Roll": [], + "Pitch": [], + "Heading": [], + "Dist": [], + "nb_valid_targets": 0, } - # --- Compact warnings for optional-but-important metadata - _missing2 = [] - - # navigation / platform variables (checked via fallback) - _nav_vars = { - "heading": ["heading", "Heading"], - "pitch": ["pitch", "Pitch"], - "roll": ["roll", "Roll"], - "heave": ["heave", "Heave", "vertical_offset"], - "dist": ["dist", "Dist", "distance"], - } - for logical_name, candidates in _nav_vars.items(): - if _fallback_1d_aligned(ds, candidates, ds["Sv"].sel(channel=channel)["ping_time"]) is None: - _missing2.append(logical_name) - - # direct dataset variables - _direct_vars = [ - "sound_absorption", - "transducer_depth", - ] - for v in _direct_vars: - if v not in ds: - _missing2.append(v) - - # bottom line if provided - if bottom_da is None: - _missing2.append("bottom_da") - - # pulse parameters from user / defaults - if "pulse_length" not in p: - _missing2.append("pulse_length") - - # emit a single clean warning - if _missing2: - print( - f"[matecho] Warning: optional metadata missing or falling back to defaults " - f"for channel '{channel}': {', '.join(_missing2)}." - ) - - ############## - - c = float(p["SoundSpeed"]) - TS_thr = float(p["TS_threshold"]) - - # --- Nav / platform - heading = _fallback_1d_aligned(ds, ["heading", "Heading"], pt) - pitch = _fallback_1d_aligned(ds, ["pitch", "Pitch"], pt) - roll = _fallback_1d_aligned(ds, ["roll", "Roll"], pt) - heave = _fallback_1d_aligned(ds, ["heave", "Heave", "vertical_offset"], pt) - dist = _fallback_1d_aligned(ds, ["dist", "Dist", "distance"], pt) - n_ping = Sv.sizes["ping_time"] - zeros = np.zeros(n_ping, dtype=float) - heading = heading if heading is not None else zeros - pitch = pitch if pitch is not None else zeros - roll = roll if roll is not None else zeros - heave = heave if heave is not None else zeros - dist = dist if dist is not None else zeros - - # --- Per-channel absorption (optional) - alpha = 0.0 - if "sound_absorption" in ds: - sa = ds["sound_absorption"] - try: - alpha = ( - float(sa.sel(channel=channel).values.item()) - if "channel" in sa.dims - else float(sa.values.item()) - ) - except Exception: - alpha = 0.0 - - # --- Transducer depth (fallback 0) - TD = 0.0 - if "transducer_depth" in ds: - try: - td = ds["transducer_depth"].sel(channel=channel) - TD = float(np.asarray(td).item()) if td.ndim == 0 else float(td.values) - except Exception: - TD = 0.0 - - # --- Range from transducer face (remove TD + heave) - R = Depth - (TD + heave.reshape(-1, 1)) - R = R.clip(min=0.0) - - # --- Vertical step & timing - dstep_da = Depth.isel(ping_time=0).diff("range_sample") # 1D - dstep_arr = np.asarray(dstep_da) - - if np.isfinite(dstep_arr).any(): - dstep = float(np.nanmedian(dstep_arr)) - else: - dstep = np.nan - - if not np.isfinite(dstep) or dstep <= 0: - dstep_da_fb = Depth.diff("range_sample") # 2D - dstep_arr_fb = np.asarray(dstep_da_fb) - if np.isfinite(dstep_arr_fb).any(): - dstep = float(np.nanmedian(dstep_arr_fb)) - else: - dstep = np.nan - - dt = (2.0 * dstep) / c if np.isfinite(dstep) and dstep > 0 else 1e-4 - - T = float(p.get("pulse_length", 1e-3)) - # if "Np" in p and p["Np"] is not None and int(p["Np"]) > 2: - # Np = int(p["Np"]) - # else: - # Np = max(3, int(round(T / max(dt, 1e-6)))) - - # --- TS constants - sv2ts = ( - 10.0 * np.log10(c * T / 2.0) - + float(p["psi_two_way"]) - + 2.0 * float(p["Sa_correction"]) - + 2.0 * float(p["Sa_EK80_nominal"]) - ) - - # --- Index domains & bottom cropping - nb_pings_tot = Sv.sizes["ping_time"] - nb_samples_tot = Sv.sizes["range_sample"] - idx_pings_tot = np.arange(nb_pings_tot, dtype=int) - idx_r_tot = np.arange(nb_samples_tot, dtype=int) - - bb = None - if bottom_da is not None: - bb = bottom_da - if "channel" in bb.dims: - bb = bb.sel(channel=channel) - bb = bb.sel(ping_time=Sv["ping_time"]) - - D = Depth.values # (ping, range_sample) - b = bb.values # (ping,) - idx_bot = np.empty(nb_pings_tot, dtype=int) - for ip in range(nb_pings_tot): - under = np.nonzero(np.isfinite(D[ip]) & (D[ip] >= b[ip]))[0] - if under.size: - idx_bot[ip] = max(under[0] - 1, 0) - else: - idx_bot[ip] = nb_samples_tot - 1 - max_idx = int(np.nanmax(idx_bot)) - idx_r_tot = np.arange(max_idx + 1, dtype=int) - - # --- Blocking in ping dimension - cells_per_ping = max(1, idx_r_tot.size) - block_size = int(min(math.ceil(float(p["block_len"]) / cells_per_ping), idx_pings_tot.size)) - num_ite = int(math.ceil(idx_pings_tot.size / max(block_size, 1))) - - # --- TVG start & helper vectors - ind_range0 = int(p["tvg_start_sample"]) # 3 (EK60) or 1 (EK80) - # Matecho: ts_range = max(delta, ((1:NbSamp)-ind_range0)*delta) - ts_range = np.maximum( - dstep, - (np.arange(nb_samples_tot) + 1 - ind_range0) * dstep, - ) - NechP = max(3, int(round(T / max(dt, 1e-6)))) # samples per pulse - - # --- Output struct (compatible with existing ESP3-style code) - try: - from .esp3 import init_st_struct # if helper exists - except Exception: - - def init_st_struct(): - return { - "TS_comp": [], - "TS_uncomp": [], - "Target_range": [], - "Target_range_disp": [], - "Target_range_min": [], - "Target_range_max": [], - "idx_r": [], - "StandDev_Angles_Minor_Axis": [], - "StandDev_Angles_Major_Axis": [], - "Angle_minor_axis": [], - "Angle_major_axis": [], - "Ping_number": [], - "Time": [], - "nb_valid_targets": 0, - "idx_target_lin": [], - "pulse_env_before_lin": [], - "pulse_env_after_lin": [], - "PulseLength_Normalized_PLDL": [], - "Transmitted_pulse_length": [], - "Heave": [], - "Roll": [], - "Pitch": [], - "Heading": [], - "Dist": [], - } - - out = init_st_struct() - times = Sv["ping_time"].values - - # --- Parameters for selection - MaxAngleOW = float(p["MaxAngleOneWayCompression"]) - MinEchoLen = float(p["MinEchoLength"]) - MaxEchoLen = float(p["MaxEchoLength"]) - MinEchoSpace = float(p["MinEchoSpace"]) - MinEchoDepthM = float(p["MinEchoDepthM"]) - MaxEchoDepthM = float(p["MaxEchoDepthM"]) - max_phase = float(p.get("MaxPhaseDeviation", np.inf)) - dec_tir = 8 # ignore first samples near TX - - # Beam pattern params - bw_al = float(p["beamwidth_along_3dB_rad"]) - bw_at = float(p["beamwidth_athwart_3dB_rad"]) - dec_al = float(p["steer_along_rad"]) - dec_at = float(p["steer_athwart_rad"]) - - # Phase sensitivities - sens_al = float(p.get("angle_sens_al", 1.0)) - sens_at = float(p.get("angle_sens_at", 1.0)) - - for ui in range(num_ite): - start = ui * block_size - stop = min((ui + 1) * block_size, idx_pings_tot.size) - idx_pings = idx_pings_tot[start:stop] - - idx_r = idx_r_tot.copy() - if bb is not None and idx_pings.size > 0: - D_block = Depth.isel(ping_time=idx_pings).values - b_block = bb.isel(ping_time=idx_pings).values - idx_bot = np.empty(idx_pings.size, dtype=int) - for k in range(idx_pings.size): - under = np.nonzero(np.isfinite(D_block[k]) & (D_block[k] >= b_block[k]))[0] - if under.size: - idx_bot[k] = max(under[0] - 1, 0) - else: - idx_bot[k] = nb_samples_tot - 1 - max_idx = int(np.nanmax(idx_bot)) - idx_r = idx_r[idx_r <= max_idx] - - # Subset arrays -> (samples, pings) - Sv_blk = ( - Sv.isel(ping_time=idx_pings, range_sample=idx_r) - .transpose("range_sample", "ping_time") - .values - ) - R_blk = ( - R.isel(ping_time=idx_pings, range_sample=idx_r) - .transpose("range_sample", "ping_time") - .values - ) - # DEPblk = ( - # Depth.isel(ping_time=idx_pings, range_sample=idx_r) - # .transpose("range_sample", "ping_time") - # .values - # ) - tsr_blk = ts_range[idx_r].reshape(-1, 1) - - # Angles: convert to radians here (echopype angles usually in degrees) - along_blk = None - athwt_blk = None - if along is not None: - along_blk = ( - along.isel(ping_time=idx_pings, range_sample=idx_r) - .transpose("range_sample", "ping_time") - .data - ) - along_blk = np.deg2rad(along_blk) - if athwt is not None: - athwt_blk = ( - athwt.isel(ping_time=idx_pings, range_sample=idx_r) - .transpose("range_sample", "ping_time") - .data - ) - athwt_blk = np.deg2rad(athwt_blk) - - # Sv->TSu and Plike - r_eff = np.maximum(tsr_blk, 1e-6) - tsu = Sv_blk + 20.0 * np.log10(r_eff) + sv2ts - Plike = tsu - 40.0 * np.log10(r_eff) - 2.0 * alpha * tsr_blk - - # Beam compensation (one-way approx); if no angles, comp=0 - if along_blk is not None and athwt_blk is not None: - x = 2.0 * ((along_blk - dec_al) / bw_al) - y = 2.0 * ((athwt_blk - dec_at) / bw_at) - one_way_comp = 6.0206 * (x**2 + y**2 - 0.18 * (x**2) * (y**2)) - else: - one_way_comp = np.zeros_like(tsu) - - ts = tsu + one_way_comp - - # Per-ping detection - for j in range(Plike.shape[1]): - zP = Plike[:, j] - zTS = ts[:, j] - zTSU = tsu[:, j] - rcol = R_blk[:, j] - - if not np.isfinite(zP).any(): - continue - - # local maxima of Plike, skipping first samples - valid = np.isfinite(zP) - valid[:dec_tir] = False - loc = np.where(valid[1:-1] & (zP[1:-1] > zP[:-2]) & (zP[1:-1] >= zP[2:]))[0] + 1 - - if loc.size == 0: - continue - - # Matecho-like: process stronger echoes first (Condition 8) - order = sorted(loc.tolist(), key=lambda kk: zTS[kk], reverse=True) - keep = [] - - for k in order: - # TS threshold - if not np.isfinite(zTS[k]) or zTS[k] <= TS_thr: - continue - - # Angle-comp guard: Matecho uses 2 * MaxAngleOW (one-way) - if (zTS[k] - zTSU[k]) > 2.0 * MaxAngleOW: - continue - - # Phase deviation test (Condition 6) - if max_phase < np.inf and along_blk is not None and athwt_blk is not None: - i0_phase = max(k - 2, 0) - i1_phase = min(k + 2, Plike.shape[0] - 1) - - al_win = along_blk[i0_phase : i1_phase + 1, j] - at_win = athwt_blk[i0_phase : i1_phase + 1, j] - - al_steps = al_win * (128.0 / np.pi) * sens_al - at_steps = at_win * (128.0 / np.pi) * sens_at - - if np.nanstd(al_steps) > max_phase or np.nanstd(at_steps) > max_phase: - continue - - # 6-dB width around Plike peak (Condition 7) - base = zP[k] - i0 = k - while i0 > 0 and np.isfinite(zP[i0 - 1]) and (zP[i0 - 1] >= base - 6.0): - i0 -= 1 - i1 = k - nP = zP.size - while i1 + 1 < nP and np.isfinite(zP[i1 + 1]) and (zP[i1 + 1] >= base - 6.0): - i1 += 1 - plen = i1 - i0 + 1 - - # echo length (in samples) bounds - if plen < round(NechP * MinEchoLen) or plen > round(NechP * MaxEchoLen): - continue - - # minimal spacing between echoes in samples (Matecho Condition 8) - if keep and min(abs(k - kk) for kk in keep) < int(round(MinEchoSpace * NechP)): - continue - - # depth band test (surface-referenced) - depth_here = float(rcol[k] + TD + heave[idx_pings[j]]) - if not (MinEchoDepthM <= depth_here <= MaxEchoDepthM): - continue - - # extra guard: do not keep targets below bottom - if bb is not None: - b_here = float(bb.values[idx_pings[j]]) - if depth_here > b_here: - continue - - keep.append(k) - - # Save detection - r_seg = R_blk[i0 : i1 + 1, j] - r_seg = r_seg[np.isfinite(r_seg)] - if r_seg.size == 0: - continue - r_min = float(r_seg.min()) - r_max = float(r_seg.max()) - r_peak = float(rcol[k]) - - out["TS_comp"].append(float(zTS[k])) - out["TS_uncomp"].append(float(zTSU[k])) - out["Target_range"].append(r_peak) - out["Target_range_disp"].append(r_peak + c * T / 4.0) - out["Target_range_min"].append(r_min) - out["Target_range_max"].append(r_max) - - out["idx_r"].append(int(idx_r[k])) - out["Ping_number"].append(int(idx_pings[j])) - out["Time"].append(times[idx_pings[j]]) - out["idx_target_lin"].append(int(idx_pings[j] * nb_samples_tot + idx_r[k])) - - out["pulse_env_before_lin"].append(int(k - i0)) - out["pulse_env_after_lin"].append(int(i1 - k)) - out["PulseLength_Normalized_PLDL"].append(plen / float(NechP)) - out["Transmitted_pulse_length"].append(int(plen)) - - out["StandDev_Angles_Minor_Axis"].append(np.nan) - out["StandDev_Angles_Major_Axis"].append(np.nan) - out["Angle_minor_axis"].append(np.nan) - out["Angle_major_axis"].append(np.nan) - - out["Heave"].append(float(heave[idx_pings[j]])) - out["Roll"].append(float(roll[idx_pings[j]])) - out["Pitch"].append(float(pitch[idx_pings[j]])) - out["Heading"].append(float(heading[idx_pings[j]])) - out["Dist"].append(float(dist[idx_pings[j]])) - - out["nb_valid_targets"] = len(out["TS_comp"]) - return out + return _matecho_struct_to_dataset(out, channel=channel) diff --git a/echopype/tests/mask/test_mask.py b/echopype/tests/mask/test_mask.py index 5b212c78f..cd45e620c 100644 --- a/echopype/tests/mask/test_mask.py +++ b/echopype/tests/mask/test_mask.py @@ -21,6 +21,10 @@ # for schoals from echopype.mask import detect_shoal + +# for single targets +from echopype.mask import detect_single_targets + from scipy import ndimage as ndi from typing import List, Union, Optional @@ -2006,3 +2010,158 @@ def test_echoview_mincan_no_linking(): # Label the mask and confirm they are separate components (not connected) _, nlab = ndi.label(mask.values, structure=np.ones((3, 3), dtype=bool)) assert nlab == 2 + +# test for single target detections + +def _make_ds_single_target_stub( + n_ping=5, n_range=10, channels=("chan1",), var_name="Sv" +) -> xr.Dataset: + coords = {"ping_time": np.arange(n_ping), "range_sample": np.arange(n_range)} + da = xr.DataArray( + np.full((len(channels), n_ping, n_range), -90.0, dtype=float), + dims=("channel", "ping_time", "range_sample"), + coords={**coords, "channel": list(channels)}, + name=var_name, + ) + return da.to_dataset() + +def _make_ds_single_target_wrong_dims( + n_ping=5, n_range=10, channels=("chan1",), var_name="Sv" +) -> xr.Dataset: + # Missing 'range_sample' on purpose + coords = {"ping_time": np.arange(n_ping)} + da = xr.DataArray( + np.full((len(channels), n_ping), -90.0, dtype=float), + dims=("channel", "ping_time"), + coords={**coords, "channel": list(channels)}, + name=var_name, + ) + return da.to_dataset() + +@pytest.mark.unit +def test_detect_single_targets_var_name_missing_in_ds_raises(): + ds = _make_ds_single_target_stub(var_name="Sv") # dataset has Sv + with pytest.raises(ValueError, match=r"var_name 'Sv_corrected' not found"): + detect_single_targets(ds, method="matecho", params={"channel": "chan1", "var_name": "Sv_corrected"}) + + +@pytest.mark.unit +def test_detect_single_targets_default_var_name_Sv_ok(): + ds = _make_ds_single_target_stub(var_name="Sv") + out = detect_single_targets(ds, method="matecho", params={"channel": "chan1"}) # no var_name -> default Sv + assert isinstance(out, xr.Dataset) + assert "target" in out.dims + assert out.sizes["target"] == 0 + + +@pytest.mark.unit +def test_detect_single_targets_var_name_wrong_dims_raises(): + ds = _make_ds_single_target_wrong_dims(var_name="Sv") + with pytest.raises(ValueError, match=r"Sv must have dims"): + detect_single_targets(ds, method="matecho", params={"channel": "chan1", "var_name": "Sv"}) + + +@pytest.mark.unit +def test_detect_single_targets_channel_not_in_ds_raises(): + ds = _make_ds_single_target_stub(channels=("chan1",)) + with pytest.raises(ValueError, match=r"Channel 'chanX' not found"): + detect_single_targets(ds, method="matecho", params={"channel": "chanX", "var_name": "Sv"}) + + +@pytest.mark.unit +def test_detect_single_targets_unknown_method_raises(): + ds = _make_ds_single_target_stub() + with pytest.raises(ValueError, match="Unsupported single-target method"): + detect_single_targets(ds, method="__bad__", params={"channel": "chan1"}) + + +@pytest.mark.unit +def test_detect_single_targets_params_required_raises(): + ds = _make_ds_single_target_stub() + with pytest.raises(ValueError, match="No parameters given"): + detect_single_targets(ds, method="matecho", params=None) + + +@pytest.mark.unit +def test_detect_single_targets_missing_channel_raises(): + ds = _make_ds_single_target_stub() + with pytest.raises(ValueError, match=r"params\['channel'\] is required"): + detect_single_targets(ds, method="matecho", params={}) + + +@pytest.mark.unit +def test_detect_single_targets_returns_dataset_with_target_dim_zero(): + ds = _make_ds_single_target_stub(channels=("chan1",)) + out = detect_single_targets(ds, method="matecho", params={"channel": "chan1"}) + + assert isinstance(out, xr.Dataset) + assert "target" in out.dims + assert out.sizes["target"] == 0 + + # required coords exist and have the right dims + assert "ping_time" in out.coords + assert "ping_number" in out.coords + assert out["ping_time"].dims == ("target",) + assert out["ping_number"].dims == ("target",) + + # scalar channel coord is OK (what your stub does) + assert "channel" in out.coords + assert out["channel"].ndim == 0 # scalar coordinate + + +@pytest.mark.unit +def test_detect_single_targets_schema_variables_present_even_if_empty(): + ds = _make_ds_single_target_stub() + out = detect_single_targets(ds, method="matecho", params={"channel": "chan1"}) + + # These are the variables your _matecho_struct_to_dataset currently creates + expected_vars = [ + "TS_comp", + "TS_uncomp", + "target_range", + "target_range_disp", + "target_range_min", + "target_range_max", + "idx_r", + "idx_target_lin", + "pulse_env_before_lin", + "pulse_env_after_lin", + "pulse_length_normalized_pldl", + "transmitted_pulse_length", + "angle_minor_axis", + "angle_major_axis", + "std_angle_minor_axis", + "std_angle_major_axis", + "heave", + "roll", + "pitch", + "heading", + "dist", + ] + for v in expected_vars: + assert v in out.data_vars + assert out[v].dims == ("target",) + assert out[v].sizes["target"] == 0 + + +@pytest.mark.unit +def test_detect_single_targets_concat_empty_is_stable(): + ds = _make_ds_single_target_stub() + out1 = detect_single_targets(ds, method="matecho", params={"channel": "chan1"}) + out2 = detect_single_targets(ds, method="matecho", params={"channel": "chan1"}) + + outc = xr.concat([out1, out2], dim="target") + assert "target" in outc.dims + assert outc.sizes["target"] == 0 + + +@pytest.mark.unit +def test_detect_single_targets_missing_optional_metadata_warns(): + ds = _make_ds_single_target_stub() + + with pytest.warns(UserWarning, match="Defaults will be used"): + detect_single_targets( + ds, + method="matecho", + params={"channel": "chan1"}, + ) From 63038e4f731728ae220d1656c2bc8f6d893503ae Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Fri, 19 Dec 2025 22:11:11 +0100 Subject: [PATCH 3/6] edit comments --- echopype/mask/single_target_detection/detect_matecho.py | 2 +- echopype/tests/mask/test_mask.py | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/echopype/mask/single_target_detection/detect_matecho.py b/echopype/mask/single_target_detection/detect_matecho.py index 026285d04..f10547d7c 100644 --- a/echopype/mask/single_target_detection/detect_matecho.py +++ b/echopype/mask/single_target_detection/detect_matecho.py @@ -132,7 +132,7 @@ def detect_matecho(ds: xr.Dataset, params: dict) -> xr.Dataset: if channel not in da["channel"].values: raise ValueError(f"Channel '{channel}' not found in {var_name}.") - # --- optional metadata warnings (kept for API clarity) + # --- metadata warnings _missing = [] for v in [ "angle_alongship", diff --git a/echopype/tests/mask/test_mask.py b/echopype/tests/mask/test_mask.py index cd45e620c..32fabcd8f 100644 --- a/echopype/tests/mask/test_mask.py +++ b/echopype/tests/mask/test_mask.py @@ -2098,13 +2098,11 @@ def test_detect_single_targets_returns_dataset_with_target_dim_zero(): assert "target" in out.dims assert out.sizes["target"] == 0 - # required coords exist and have the right dims assert "ping_time" in out.coords assert "ping_number" in out.coords assert out["ping_time"].dims == ("target",) assert out["ping_number"].dims == ("target",) - # scalar channel coord is OK (what your stub does) assert "channel" in out.coords assert out["channel"].ndim == 0 # scalar coordinate @@ -2114,7 +2112,7 @@ def test_detect_single_targets_schema_variables_present_even_if_empty(): ds = _make_ds_single_target_stub() out = detect_single_targets(ds, method="matecho", params={"channel": "chan1"}) - # These are the variables your _matecho_struct_to_dataset currently creates + # Compare with the variables _matecho_struct_to_dataset currently creates expected_vars = [ "TS_comp", "TS_uncomp", From 04c5ca2ac777bcce7138e24c054286cc59648df1 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Mon, 22 Dec 2025 20:20:03 +0100 Subject: [PATCH 4/6] Update detect_matecho.py --- .../single_target_detection/detect_matecho.py | 1035 +++++++++++++++-- 1 file changed, 924 insertions(+), 111 deletions(-) diff --git a/echopype/mask/single_target_detection/detect_matecho.py b/echopype/mask/single_target_detection/detect_matecho.py index f10547d7c..8fece863f 100644 --- a/echopype/mask/single_target_detection/detect_matecho.py +++ b/echopype/mask/single_target_detection/detect_matecho.py @@ -1,39 +1,42 @@ import warnings +from typing import Any, Dict, List, Tuple import numpy as np import xarray as xr +# --- +# Defaults PARAM_DEFAULTS = { "SoundSpeed": 1500.0, - "TS_threshold": -50.0, # dB - "MaxAngleOneWayCompression": 6.0, # dB (one-way; test uses 2x) - "MaxPhaseDeviation": 8.0, # phase steps (128/pi units), Matecho GUI - "MinEchoLength": 0.8, # in pulse lengths + "MinThreshold": -60.0, + "MaxAngleOneWayCompression": 6.0, + "MaxPhaseDeviation": 8.0, + "MinEchoLength": 0.8, "MaxEchoLength": 1.8, - "MinEchoSpace": 1.0, # in pulse lengths + "MinEchoSpace": 1.0, "MinEchoDepthM": 3.0, "MaxEchoDepthM": 38.0, - "tvg_start_sample": 3, # EK60=3, EK80=1 - "block_len": 1e7 / 3, - # Beam pattern (fallbacks if no metadata): - "beamwidth_along_3dB_rad": np.deg2rad(7.0), - "beamwidth_athwart_3dB_rad": np.deg2rad(7.0), - "steer_along_rad": 0.0, - "steer_athwart_rad": 0.0, - # Angle sensitivities (phase steps), overwritten from ds: - "angle_sens_al": 1.0, - "angle_sens_at": 1.0, - # Sv->TS constant terms; keep 0 if unknown: + "FloorFilterTolerance": 0.0, + "tvg_start_sample": 3, # (EK60=3, EK80=1) + # Sv->TS "psi_two_way": 0.0, "Sa_correction": 0.0, "Sa_EK80_nominal": 0.0, + # Required for CW range/TS conversion + "pulse_length": None, # seconds (pd in MATLAB) + "sound_speed": None, # m/s (c in MATLAB) + "alpha": 0.0, # dB/m (alpha(kf)) + # TD/CH + "transducer_depth": None, # TD (m), from surface + "heave_compensation": 0.0, # CH (m) per ping or scalar + # (optional) + "bottom_da": None, } +# --- +# Output conversion def _matecho_struct_to_dataset(out: dict, *, channel: str | None = None) -> xr.Dataset: - """ - Convert struct-of-lists into an xr.Dataset with dim 'target'. - """ n = int(out.get("nb_valid_targets", len(out.get("TS_comp", [])))) target = np.arange(n, dtype=int) @@ -41,47 +44,62 @@ def _arr(key, dtype=float): vals = out.get(key, []) if n == 0: return np.asarray([], dtype=dtype) - a = np.asarray(vals, dtype=dtype) + a = np.asarray(vals) + if dtype is not None: + a = a.astype(dtype) if a.size != n: raise ValueError(f"Output field '{key}' has length {a.size}, expected {n}.") return a + tsfreq = np.asarray(out.get("TSfreq_matrix", np.empty((0, 23))), dtype=float) + if tsfreq.ndim == 1: + tsfreq = tsfreq.reshape(1, -1) + if tsfreq.size == 0: + tsfreq = np.empty((0, 23), dtype=float) + if tsfreq.shape[0] != n: + raise ValueError( + f"TSfreq_matrix has {tsfreq.shape[0]} rows but expected {n} (nb_valid_targets)." + ) + ds_out = xr.Dataset( data_vars=dict( + nb_valid_targets=((), n), TS_comp=("target", _arr("TS_comp", float), {"units": "dB"}), TS_uncomp=("target", _arr("TS_uncomp", float), {"units": "dB"}), - target_range=("target", _arr("Target_range", float), {"units": "m"}), - target_range_disp=("target", _arr("Target_range_disp", float), {"units": "m"}), - target_range_min=("target", _arr("Target_range_min", float), {"units": "m"}), - target_range_max=("target", _arr("Target_range_max", float), {"units": "m"}), + Target_range=("target", _arr("Target_range", float), {"units": "m"}), + Target_range_disp=("target", _arr("Target_range_disp", float), {"units": "m"}), + Target_range_min=("target", _arr("Target_range_min", float), {"units": "m"}), + Target_range_max=("target", _arr("Target_range_max", float), {"units": "m"}), idx_r=("target", _arr("idx_r", int)), idx_target_lin=("target", _arr("idx_target_lin", int)), - pulse_env_before_lin=("target", _arr("pulse_env_before_lin", int)), - pulse_env_after_lin=("target", _arr("pulse_env_after_lin", int)), - pulse_length_normalized_pldl=("target", _arr("PulseLength_Normalized_PLDL", float)), - transmitted_pulse_length=("target", _arr("Transmitted_pulse_length", int)), - angle_minor_axis=("target", _arr("Angle_minor_axis", float), {"units": "rad"}), - angle_major_axis=("target", _arr("Angle_major_axis", float), {"units": "rad"}), - std_angle_minor_axis=( + pulse_env_before_lin=("target", _arr("pulse_env_before_lin", float)), + pulse_env_after_lin=("target", _arr("pulse_env_after_lin", float)), + PulseLength_Normalized_PLDL=("target", _arr("PulseLength_Normalized_PLDL", float)), + Transmitted_pulse_length=("target", _arr("Transmitted_pulse_length", float)), + Angle_minor_axis=("target", _arr("Angle_minor_axis", float), {"units": "rad"}), + Angle_major_axis=("target", _arr("Angle_major_axis", float), {"units": "rad"}), + StandDev_Angles_Minor_Axis=( "target", _arr("StandDev_Angles_Minor_Axis", float), {"units": "phase_steps"}, ), - std_angle_major_axis=( + StandDev_Angles_Major_Axis=( "target", _arr("StandDev_Angles_Major_Axis", float), {"units": "phase_steps"}, ), - heave=("target", _arr("Heave", float), {"units": "m"}), - roll=("target", _arr("Roll", float), {"units": "rad"}), - pitch=("target", _arr("Pitch", float), {"units": "rad"}), - heading=("target", _arr("Heading", float), {"units": "rad"}), - dist=("target", _arr("Dist", float), {"units": "m"}), + Heave=("target", _arr("Heave", float), {"units": "m"}), + Roll=("target", _arr("Roll", float), {"units": "rad"}), + Pitch=("target", _arr("Pitch", float), {"units": "rad"}), + Heading=("target", _arr("Heading", float), {"units": "rad"}), + Dist=("target", _arr("Dist", float), {"units": "m"}), + TSfreq_matrix=(("target", "tsfreq_col"), tsfreq, {}), ), coords=dict( target=target, - ping_number=("target", _arr("Ping_number", int)), - ping_time=("target", _arr("Time", "datetime64[ns]")), + tsfreq_col=np.arange(tsfreq.shape[1], dtype=int), + Ping_number=("target", _arr("Ping_number", int)), + Time=("target", _arr("Time", np.dtype("datetime64[ns]"))), ), attrs=dict( method="matecho", @@ -89,104 +107,899 @@ def _arr(key, dtype=float): nb_valid_targets=n, ), ) - if channel is not None: - # scalar coordinate - ds_out = ds_out.assign_coords(channel=np.asarray(channel)) - + ds_out = ds_out.assign_coords(channel=channel) return ds_out -def detect_matecho(ds: xr.Dataset, params: dict) -> xr.Dataset: - """ - Matecho-inspired single-target detector (CW only). +def _empty_out(channel: str) -> xr.Dataset: + out = { + "nb_valid_targets": 0, + "TS_comp": [], + "TS_uncomp": [], + "Target_range": [], + "Target_range_disp": [], + "Target_range_min": [], + "Target_range_max": [], + "idx_r": [], + "idx_target_lin": [], + "pulse_env_before_lin": [], + "pulse_env_after_lin": [], + "PulseLength_Normalized_PLDL": [], + "Transmitted_pulse_length": [], + "Angle_minor_axis": [], + "Angle_major_axis": [], + "StandDev_Angles_Minor_Axis": [], + "StandDev_Angles_Major_Axis": [], + "Heave": [], + "Roll": [], + "Pitch": [], + "Heading": [], + "Dist": [], + "Ping_number": [], + "Time": [], + "TSfreq_matrix": np.empty((0, 23)), + } + return _matecho_struct_to_dataset(out, channel=channel) - This is a placeholder implementation that demonstrates: - - required signature - - required parameter checks - - required output structure (xr.Dataset with dim 'target') - No detection is performed. - """ +# --- +# Helpers +def _make_dbg() -> dict: # to print debug + return { + "enable": True, + "print_header": True, + "print_block": True, + "print_cond_1to5": True, + "print_cond6": True, + "print_cond7": True, + "print_per_ping": False, + "print_examples": True, + "n_examples": 10, + } + + +def _nanmin(a): + return float(np.nanmin(a)) if np.any(np.isfinite(a)) else np.nan + + +def _nanmax(a): + return float(np.nanmax(a)) if np.any(np.isfinite(a)) else np.nan + + +def _nanmed(a): + return float(np.nanmedian(a)) if np.any(np.isfinite(a)) else np.nan + + +def _merge_params(params: dict) -> dict: if params is None: raise ValueError("params is required.") + return {**PARAM_DEFAULTS, **params} - channel = params.get("channel") + +def _validate_inputs(ds: xr.Dataset, Param: dict) -> Tuple[str, xr.DataArray]: + channel = Param.get("channel") if channel is None: raise ValueError("params['channel'] is required.") - var_name = params.get("var_name", "Sv") - + var_name = Param.get("var_name", "Sv") if var_name not in ds: raise ValueError(f"var_name '{var_name}' not found in input dataset.") - da = ds[var_name] - + Sv = ds[var_name] required_dims = {"channel", "ping_time", "range_sample"} - if not required_dims.issubset(set(da.dims)): - raise ValueError(f"{var_name} must have dims {sorted(required_dims)}. Got: {da.dims}.") + if not required_dims.issubset(set(Sv.dims)): + raise ValueError(f"{var_name} must have dims {sorted(required_dims)}. Got: {Sv.dims}.") - # channel must exist in Sv(channel, ping_time, range_sample) - if "channel" not in da.coords: - raise ValueError(f"{var_name} has no 'channel' coordinate.") - if channel not in da["channel"].values: + if channel not in Sv["channel"].values: raise ValueError(f"Channel '{channel}' not found in {var_name}.") - # --- metadata warnings - _missing = [] - for v in [ - "angle_alongship", - "angle_athwartship", - "beamwidth_alongship", - "beamwidth_athwartship", - "angle_offset_alongship", - "angle_offset_athwartship", - "angle_sensitivity_alongship", - "angle_sensitivity_athwartship", - "sound_speed", - ]: - if v not in ds: - _missing.append(v) - - if _missing: + return channel, Sv + + +def _extract_matrices(ds: xr.Dataset, Sv: xr.DataArray, channel: str) -> Dict[str, Any]: + sv_mat = Sv.sel(channel=channel).transpose("ping_time", "range_sample").values + npings, NbSamp = sv_mat.shape + + if "angle_alongship" in ds: + al_mat = ( + ds["angle_alongship"].sel(channel=channel).transpose("ping_time", "range_sample").values + ) + else: + al_mat = np.full_like(sv_mat, np.nan, dtype=float) + + if "angle_athwartship" in ds: + ath_mat = ( + ds["angle_athwartship"] + .sel(channel=channel) + .transpose("ping_time", "range_sample") + .values + ) + else: + ath_mat = np.full_like(sv_mat, np.nan, dtype=float) + + if "depth" in ds: + depth = ds["depth"].sel(channel=channel).median("ping_time", skipna=True).values + else: + depth = np.arange(NbSamp, dtype=float) + + if len(depth) < 2: + raise ValueError("Need at least 2 samples to compute delta.") + + # trim NaN tail + valid = np.isfinite(depth) + if not np.any(valid): + raise ValueError("Depth vector is all-NaN for this channel.") + last = np.where(valid)[0][-1] + 1 + + depth = depth[:last] + sv_mat = sv_mat[:, :last] + al_mat = al_mat[:, :last] + ath_mat = ath_mat[:, :last] + NbSamp = last + + return { + "sv_mat": sv_mat, + "al_mat": al_mat, + "ath_mat": ath_mat, + "depth": depth, + "npings": npings, + "NbSamp": NbSamp, + } + + +def _compute_delta(depth: np.ndarray) -> float: + dd = np.diff(depth) + dd = dd[np.isfinite(dd)] + if dd.size == 0: + raise ValueError("Cannot compute delta from depth: all diffs are NaN.") + return float(np.nanmedian(dd)) + + +def _extract_metadata( + ds: xr.Dataset, Param: dict, params: dict, channel: str, npings: int, DBG: dict +) -> Dict[str, Any]: + # c, pd + c = Param["sound_speed"] + if c is None: + if "sound_speed" in ds: + c = float(ds["sound_speed"].values) + else: + c = float(PARAM_DEFAULTS["SoundSpeed"]) + + pd = Param["pulse_length"] + if pd is None: + raise ValueError("params['pulse_length'] is required for CW (pd in MATLAB).") + + alpha = float(Param.get("alpha", 0.0)) + psi = float(Param.get("psi_two_way", 0.0)) + sa_cor = float(Param.get("Sa_correction", 0.0)) + SacorEK80Nominal = float(Param.get("Sa_EK80_nominal", 0.0)) + + # beam meta + ouv_al = ( + float(ds.get("beamwidth_alongship", xr.DataArray(np.nan)).sel(channel=channel).values) + if "beamwidth_alongship" in ds + else np.nan + ) + ouv_at = ( + float(ds.get("beamwidth_athwartship", xr.DataArray(np.nan)).sel(channel=channel).values) + if "beamwidth_athwartship" in ds + else np.nan + ) + dec_al = ( + float(ds.get("angle_offset_alongship", xr.DataArray(0.0)).sel(channel=channel).values) + if "angle_offset_alongship" in ds + else 0.0 + ) + dec_at = ( + float(ds.get("angle_offset_athwartship", xr.DataArray(0.0)).sel(channel=channel).values) + if "angle_offset_athwartship" in ds + else 0.0 + ) + al_sens = ( + float(ds.get("angle_sensitivity_alongship", xr.DataArray(1.0)).sel(channel=channel).values) + if "angle_sensitivity_alongship" in ds + else 1.0 + ) + ath_sens = ( + float( + ds.get("angle_sensitivity_athwartship", xr.DataArray(1.0)).sel(channel=channel).values + ) + if "angle_sensitivity_athwartship" in ds + else 1.0 + ) + + missing_meta = [] + if np.isnan(ouv_al): + missing_meta.append("beamwidth_alongship") + if np.isnan(ouv_at): + missing_meta.append("beamwidth_athwartship") + if missing_meta: warnings.warn( - f"The following variables are missing for channel '{channel}': " - + ", ".join(_missing) - + ". Defaults will be used where applicable.", + f"Missing beam metadata for channel '{channel}': {', '.join(missing_meta)}. ", UserWarning, stacklevel=2, ) - # ------------------------------------------------------------------ - # ## algo ## - # Intentionally empty: no detection performed - # ------------------------------------------------------------------ + # TD/CH + TD = Param.get("transducer_depth", None) + if TD is None: + if "transducer_depth" in ds: + TD = float(ds["transducer_depth"].sel(channel=channel).values) + else: + TD = 0.0 + TD = float(TD) + + CH = Param.get("heave_compensation", 0.0) + if np.ndim(CH) == 0: + CH_vec = np.full((npings,), float(CH), dtype=float) + else: + CH_vec = np.asarray(CH, dtype=float) + if CH_vec.size != npings: + raise ValueError("heave_compensation length must match number of pings.") + + # bottom + bottom_da = Param.get("bottom_da", None) + if bottom_da is None: + bot = np.full((npings,), np.inf, dtype=float) + else: + bot = bottom_da.sel(ping_time=ds["ping_time"]).values.astype(float) + + # timing + timeSampleInterval_usec = params.get("timeSampleInterval_usec", None) + if timeSampleInterval_usec is None: + raise ValueError("params['timeSampleInterval_usec'] is required.") + + # angle units detection + beamwidth_units = "" + if "beamwidth_alongship" in ds and hasattr(ds["beamwidth_alongship"], "attrs"): + beamwidth_units = str(ds["beamwidth_alongship"].attrs.get("units", "")).lower() + angles_are_degrees = "deg" in beamwidth_units + + if DBG["enable"]: + warnings.warn( + f"beamwidth_alongship units='{beamwidth_units}' -> \ + angles_are_degrees = {angles_are_degrees}", + category=UserWarning, + stacklevel=2, + ) + + return { + "c": float(c), + "pd": float(pd), + "alpha": float(alpha), + "psi": float(psi), + "sa_cor": float(sa_cor), + "SacorEK80Nominal": float(SacorEK80Nominal), + "ouv_al": float(ouv_al), + "ouv_at": float(ouv_at), + "dec_al": float(dec_al), + "dec_at": float(dec_at), + "al_sens": float(al_sens), + "ath_sens": float(ath_sens), + "TD": float(TD), + "CH_vec": CH_vec, + "bot": bot, + "timeSampleInterval_usec": float(timeSampleInterval_usec), + "angles_are_degrees": bool(angles_are_degrees), + } + + +def _normalise_angles(al_mat, ath_mat, meta: dict, DBG: dict): + if not meta["angles_are_degrees"]: + return al_mat, ath_mat, meta + + DEG2RAD = np.pi / 180.0 + if DBG["enable"]: + warnings.warn( + "Converting angles, beamwidths, and offsets from degrees to radians.", + category=UserWarning, + stacklevel=2, + ) + + al_mat = al_mat * DEG2RAD + ath_mat = ath_mat * DEG2RAD + meta = dict(meta) + meta["dec_al"] *= DEG2RAD + meta["dec_at"] *= DEG2RAD + meta["ouv_al"] *= DEG2RAD + meta["ouv_at"] *= DEG2RAD + return al_mat, ath_mat, meta + + +def _compute_pulse_derived(meta: dict, Param: dict, delta: float) -> Dict[str, float]: + c = meta["c"] + pd = meta["pd"] + pulse_len_m = c * pd / 2.0 + + n_min = int(np.round((pulse_len_m * float(Param["MinEchoLength"])) / delta)) + n_max = int(np.round((pulse_len_m * float(Param["MaxEchoLength"])) / delta)) + n_min = max(n_min, 1) + n_max = max(n_max, n_min) + + min_space_m = float(Param["MinEchoSpace"]) * pulse_len_m + return { + "pulse_len_m": float(pulse_len_m), + "n_min": int(n_min), + "n_max": int(n_max), + "min_space_m": float(min_space_m), + } + + +# --- +# Helpers: core matrices (ts_range, tsu, Plike, ts) +def _compute_core_matrices( + sv_mat: np.ndarray, + al_mat: np.ndarray, + ath_mat: np.ndarray, + depth: np.ndarray, + delta: float, + meta: dict, + Param: dict, + channel: str, +) -> Dict[str, Any]: + npings, NbSamp = sv_mat.shape + + ind_range0 = int(Param["tvg_start_sample"]) + dec_tir = 8 + + c = meta["c"] + pd = meta["pd"] + alpha = meta["alpha"] + psi = meta["psi"] + sa_cor = meta["sa_cor"] + SacorEK80Nominal = meta["SacorEK80Nominal"] + + # CW conversion constant + sv2ts = 10.0 * np.log10(c * pd / 2.0) + psi + 2.0 * sa_cor + 2.0 * SacorEK80Nominal + + # NechP (dt-based) + NechP = pd / (meta["timeSampleInterval_usec"] * 1e-6) + + # range vector in metres (one-way) + jj = np.arange(1, NbSamp + 1, dtype=float) + ts_range_1d = np.maximum(delta, (jj - ind_range0) * delta) + ts_range = np.tile(ts_range_1d[None, :], (npings, 1)) + + # bottom index constraint + bot = meta["bot"] + TD = meta["TD"] + CH_vec = meta["CH_vec"] + ir = (bot - TD - CH_vec - float(Param["FloorFilterTolerance"])) / delta + ir = np.minimum(ir, NbSamp).astype(float) + + # tsu / Plike / ts + tsu_mat = sv_mat + 20.0 * np.log10(ts_range) + sv2ts + Plike_mat = tsu_mat - 40.0 * np.log10(ts_range) - 2.0 * alpha * ts_range + + ouv_al = meta["ouv_al"] + ouv_at = meta["ouv_at"] + dec_al = meta["dec_al"] + dec_at = meta["dec_at"] + + x = 2.0 * (al_mat - dec_al) / ouv_al + y = 2.0 * (ath_mat - dec_at) / ouv_at + ts_mat = tsu_mat + 6.0206 * (x * x + y * y - 0.18 * (x * x) * (y * y)) + + return { + "npings": npings, + "NbSamp": NbSamp, + "ind_range0": ind_range0, + "dec_tir": dec_tir, + "NechP": float(NechP), + "sv2ts": float(sv2ts), + "ts_range": ts_range, + "ir": ir, + "tsu_mat": tsu_mat, + "Plike_mat": Plike_mat, + "ts_mat": ts_mat, + "depth": depth, + "channel": channel, + } + + +# --- +# Helpers: Conditions 1–5 +def _cond_1to5(core: dict, Param: dict, DBG: dict) -> Tuple[np.ndarray, np.ndarray]: + NbSamp = core["NbSamp"] + ts_mat = core["ts_mat"] + tsu_mat = core["tsu_mat"] + Plike_mat = core["Plike_mat"] + ir = core["ir"] + dec_tir = core["dec_tir"] + + cols = np.arange(1, NbSamp - 1, dtype=int) # never evaluate edge samples + + cond1 = ts_mat[:, 1:-1] > float(Param["MinThreshold"]) + cond2 = (ts_mat[:, 1:-1] - tsu_mat[:, 1:-1]) <= 2.0 * float(Param["MaxAngleOneWayCompression"]) + + left = Plike_mat[:, 0:-2] + mid = Plike_mat[:, 1:-1] + right = Plike_mat[:, 2:] + cond3 = mid >= np.nanmax(np.stack([left, right]), axis=0) + + col_1based = cols + 1 + cond4 = col_1based[None, :] < ir[:, None] + cond5 = col_1based[None, :] > dec_tir + + mask_1to5 = cond1 & cond2 & cond3 & cond4 & cond5 + i1, ind0 = np.where(mask_1to5) + ind = ind0 + 1 # back to full-matrix python col index + + if DBG["enable"] and DBG["print_cond_1to5"]: + warnings.warn( + "[Cond 1–5] Counts (pixel-level):\n" + f" c1(ts > thr): {int(np.count_nonzero(cond1))}\n" + f" c2(comp <=): {int(np.count_nonzero(cond2))}\n" + f" c3(loc max): {int(np.count_nonzero(cond3))}\n" + f" c4(above bottom): {int(np.count_nonzero(cond4))}\n" + f" c5(idx > dec_tir): {int(np.count_nonzero(cond5))}\n" + f" ALL: {int(np.count_nonzero(mask_1to5))}\n" + f" detections after Cond 1–5: {ind.size}", + category=UserWarning, + stacklevel=2, + ) + + return i1, ind + + +# --- +# Helpers: Condition 6 +def _cond_6_phase( + i1: np.ndarray, + ind: np.ndarray, + al_mat: np.ndarray, + ath_mat: np.ndarray, + core: dict, + meta: dict, + Param: dict, + DBG: dict, +) -> np.ndarray: + NbSamp = core["NbSamp"] + al_sens = meta["al_sens"] + ath_sens = meta["ath_sens"] + + std_ph_al = np.full((ind.size,), np.nan, dtype=float) + std_ph_ath = np.full((ind.size,), np.nan, dtype=float) + + for i in range(ind.size): + r0 = max(ind[i] - 2, 0) + r1 = min(ind[i] + 2, NbSamp - 1) + sl = slice(r0, r1 + 1) + std_ph_al[i] = np.nanstd(al_mat[i1[i], sl] * 128.0 / np.pi * al_sens) + std_ph_ath[i] = np.nanstd(ath_mat[i1[i], sl] * 128.0 / np.pi * ath_sens) + + ind2 = np.where( + (std_ph_al <= float(Param["MaxPhaseDeviation"])) + & (std_ph_ath <= float(Param["MaxPhaseDeviation"])) + )[0] + + if DBG["enable"] and DBG["print_cond6"]: + warnings.warn( + "[Cond 6] Phase deviation:\n" + f" std_ph_al: min={_nanmin(std_ph_al):.3f} " + f"med={_nanmed(std_ph_al):.3f} " + f"max={_nanmax(std_ph_al):.3f}\n" + f" std_ph_ath: min={_nanmin(std_ph_ath):.3f} " + f"med={_nanmed(std_ph_ath):.3f} " + f"max={_nanmax(std_ph_ath):.3f}\n" + f" Passing (<= MaxPhaseDeviation={float(Param['MaxPhaseDeviation']):.3f}): " + f"{ind2.size} / {ind.size}", + category=UserWarning, + stacklevel=2, + ) + + return ind2 + + +# --- +# Helpers: Condition 7 +def _cond_7_echo_length( + i1: np.ndarray, + ind: np.ndarray, + ind2: np.ndarray, + core: dict, + derived: dict, + Param: dict, + DBG: dict, +) -> Dict[str, Any]: + NbSamp = core["NbSamp"] + Plike_mat = core["Plike_mat"] + ir = core["ir"] + + iid = ind[ind2] # python col index (0-based) + iip = i1[ind2] # ping indices + necho = iid.size + + isup = np.full((necho,), np.nan, dtype=float) + iinf = np.full((necho,), np.nan, dtype=float) + + for i in range(necho): + ip = iip[i] + ii = iid[i] + + ir_floor = int(np.floor(ir[ip])) + ir_floor = max(ir_floor, 1) + ir_floor = min(ir_floor, NbSamp) + + thr = Plike_mat[ip, ii] - 6.0 + + forward = Plike_mat[ip, ii:ir_floor] < thr + if np.any(forward): + umin = np.where(forward)[0].min() + isup[i] = ii + umin + 1.0 + else: + isup[i] = ii + 1.0 + + backward = Plike_mat[ip, 0 : ii + 1] < thr + if np.any(backward): + umax = np.where(backward)[0].max() + iinf[i] = umax + 2.0 + else: + iinf[i] = ii + 1.0 + + L6dB = isup - iinf + 1.0 + + n_min = derived["n_min"] + n_max = derived["n_max"] + ind4 = np.where((L6dB >= n_min) & (L6dB <= n_max))[0] + + if DBG["enable"] and DBG["print_cond7"]: + warnings.warn( + "[Cond 7] Echo length (L6dB):\n" + f" L6dB: min={_nanmin(L6dB):.0f} " + f"med={_nanmed(L6dB):.0f} " + f"max={_nanmax(L6dB):.0f}\n" + f" Bounds depth-based [{n_min}..{n_max}] samples\n" + f" Passing depth-based: {ind4.size} / {necho}", + category=UserWarning, + stacklevel=2, + ) + + return {"iid": iid, "iip": iip, "iinf": iinf, "isup": isup, "L6dB": L6dB, "ind4": ind4} + + +# --- +# Helpers: Condition 8 (refinement + spacing + output build) +def _refine_ir_fin_1b(Plike_row: np.ndarray, iinf_1b: int, isup_1b: int) -> float: + il6dB = np.arange(iinf_1b - 1, isup_1b, dtype=int) + w = 10.0 ** (Plike_row[il6dB] / 10.0) + return float(np.nansum((il6dB + 1) * w) / np.nansum(w)) + + +def _cond_8_build_outputs( + ds: xr.Dataset, + params: dict, + i1: np.ndarray, + ind2: np.ndarray, + c7: dict, + core: dict, + meta: dict, + derived: dict, + al_mat: np.ndarray, + ath_mat: np.ndarray, + depth: np.ndarray, + Param: dict, + DBG: dict, + delta: float, +) -> Tuple[dict, List[List[float]]]: + + ts_mat = core["ts_mat"] + tsu_mat = core["tsu_mat"] + Plike_mat = core["Plike_mat"] + ind_range0 = core["ind_range0"] + + Ping_number_full = np.arange(1, core["npings"] + 1, dtype=int) + ping_time = ds["ping_time"].values + + # candidates per ping after Cond7 + iid = c7["iid"] + iinf = c7["iinf"] + isup = c7["isup"] + ind4 = c7["ind4"] + + ip0 = i1[ind2[ind4]] + ipu = np.unique(ip0) out = { - "TS_comp": [], - "TS_uncomp": [], - "Target_range": [], - "Target_range_disp": [], - "Target_range_min": [], - "Target_range_max": [], - "idx_r": [], - "StandDev_Angles_Minor_Axis": [], - "StandDev_Angles_Major_Axis": [], - "Angle_minor_axis": [], - "Angle_major_axis": [], - "Ping_number": [], - "Time": [], - "idx_target_lin": [], - "pulse_env_before_lin": [], - "pulse_env_after_lin": [], - "PulseLength_Normalized_PLDL": [], - "Transmitted_pulse_length": [], - "Heave": [], - "Roll": [], - "Pitch": [], - "Heading": [], - "Dist": [], - "nb_valid_targets": 0, + k: [] + for k in [ + "TS_comp", + "TS_uncomp", + "Target_range", + "Target_range_disp", + "Target_range_min", + "Target_range_max", + "idx_r", + "StandDev_Angles_Minor_Axis", + "StandDev_Angles_Major_Axis", + "Angle_minor_axis", + "Angle_major_axis", + "Ping_number", + "Time", + "idx_target_lin", + "pulse_env_before_lin", + "pulse_env_after_lin", + "PulseLength_Normalized_PLDL", + "Transmitted_pulse_length", + "Heave", + "Roll", + "Pitch", + "Heading", + "Dist", + ] } + TSfreq_rows: List[List[float]] = [] + + alpha = meta["alpha"] + TD = meta["TD"] + CH_vec = meta["CH_vec"] + bot = meta["bot"] + min_space_m = derived["min_space_m"] + + for ip in ipu: + ip2 = np.where(ip0 == ip)[0] + sel = ind4[ip2] + + # compute fine range ir_fin + ir_fin = np.full((sel.size,), np.nan, dtype=float) + for j, idx in enumerate(sel): + ir_fin[j] = _refine_ir_fin_1b(Plike_mat[ip, :], int(iinf[idx]), int(isup[idx])) + + i_ini = iid[sel] + i_ini_1b = i_ini + 1.0 + + ts_fin = ts_mat[ip, i_ini] + ( + 40.0 * np.log10((ir_fin - ind_range0) / (i_ini_1b - ind_range0)) + + 2.0 * alpha * delta * (ir_fin - i_ini_1b) + ) + tsu_fin = tsu_mat[ip, i_ini] + ( + 40.0 * np.log10((ir_fin - ind_range0) / (i_ini_1b - ind_range0)) + + 2.0 * alpha * delta * (ir_fin - i_ini_1b) + ) + + # convert to metres once + range_fin_m = (ir_fin - ind_range0) * delta + + # spacing in metres + if sel.size > 1: + inds = np.argsort(-ts_fin) + keep = [inds[0]] + for k in inds[1:]: + if np.nanmin(np.abs(range_fin_m[k] - range_fin_m[np.array(keep)])) >= min_space_m: + keep.append(k) + ind5 = np.array(keep, dtype=int) + else: + ind5 = np.array([0], dtype=int) + + targetRange = range_fin_m[ind5] # metres (one-way range from transducer) + D = targetRange + TD + CH_vec[ip] # metres (depth from surface) + + ok = (D >= float(Param["MinEchoDepthM"])) & (D <= float(Param["MaxEchoDepthM"])) + if not np.any(ok): + continue + + targetRange = targetRange[ok] + D = D[ok] + ind5_ok = ind5[ok] + + compensatedTS = ts_fin[ind5_ok] + unCompensatedTS = tsu_fin[ind5_ok] + + samp_idx = i_ini[ind5_ok] + AlongShipAngleRad = al_mat[ip, samp_idx] + AthwartShipAngleRad = ath_mat[ip, samp_idx] + + for pidx in range(targetRange.size): + positionWorldCoord = [np.nan, np.nan, np.nan] # to add + positionBeamCoord = [ + targetRange[pidx] * np.tan(AlongShipAngleRad[pidx]), + targetRange[pidx] * np.tan(AthwartShipAngleRad[pidx]), + targetRange[pidx], + ] + + DepthIndex = int(np.argmin(np.abs(depth - D[pidx])) + 1) + IndexPingClean = int(Ping_number_full[ip]) + DepthIndexFromTransducer = int(np.round(targetRange[pidx] / delta)) + + kf = int(params.get("Frequency_index", 1)) + + TSfreq_rows.append( + [ + float(Ping_number_full[ip]), + float(D[pidx]), + float(DepthIndex), + float(compensatedTS[pidx]), + float(unCompensatedTS[pidx]), + float(AlongShipAngleRad[pidx]), + float(AthwartShipAngleRad[pidx]), + 0.0, + float(positionWorldCoord[0]), + float(positionWorldCoord[1]), + float(positionWorldCoord[2]), + float(positionBeamCoord[0]), + float(positionBeamCoord[1]), + float(positionBeamCoord[2]), + float(Param["MinEchoDepthM"]), + float(Param["MaxEchoDepthM"]), + float(bot[ip]) if np.isfinite(bot[ip]) else np.nan, + float(IndexPingClean), + float(DepthIndexFromTransducer), + float(targetRange.size), + 1.0, + 0.0, + float(kf), + ] + ) + + out["TS_comp"].append(float(compensatedTS[pidx])) + out["TS_uncomp"].append(float(unCompensatedTS[pidx])) + + out["Target_range"].append(float(D[pidx])) + out["Target_range_disp"].append(float(D[pidx])) + out["Target_range_min"].append(float(Param["MinEchoDepthM"])) + out["Target_range_max"].append(float(Param["MaxEchoDepthM"])) + + out["idx_r"].append(int(samp_idx[pidx] + 1)) + out["idx_target_lin"].append(int(samp_idx[pidx] + 1)) + + out["pulse_env_before_lin"].append(np.nan) + out["pulse_env_after_lin"].append(np.nan) + out["PulseLength_Normalized_PLDL"].append(np.nan) + out["Transmitted_pulse_length"].append(float(meta["pd"])) + + out["Angle_minor_axis"].append(float(AlongShipAngleRad[pidx])) + out["Angle_major_axis"].append(float(AthwartShipAngleRad[pidx])) + + out["StandDev_Angles_Minor_Axis"].append(float(np.nan)) + out["StandDev_Angles_Major_Axis"].append(float(np.nan)) + + out["Heave"].append(float(CH_vec[ip])) + out["Roll"].append( + float(ds.get("roll", xr.DataArray(np.nan)).isel(ping_time=ip).values) + if "roll" in ds + else np.nan + ) + out["Pitch"].append( + float(ds.get("pitch", xr.DataArray(np.nan)).isel(ping_time=ip).values) + if "pitch" in ds + else np.nan + ) + out["Heading"].append( + float(ds.get("heading", xr.DataArray(np.nan)).isel(ping_time=ip).values) + if "heading" in ds + else np.nan + ) + out["Dist"].append(np.nan) + out["Ping_number"].append(int(Ping_number_full[ip])) + out["Time"].append(ping_time[ip]) + + return out, TSfreq_rows + + +# --- +# main +def detect_matecho(ds: xr.Dataset, params: dict) -> xr.Dataset: + """ + Single-target detector (split-beam style) on CW Sv data. + + Overview + -------- + Detect candidate single targets ping-by-ping using a peak-based workflow: + + 1) Build core matrices + - Convert Sv (dB) to TS-like quantities (uncompensated and compensated). + - Compute one-way range (m) and apply bottom/near-surface constraints. + + 2) Phase I: candidate peak selection (Conditions 1–5) + - C1: compensated TS above `MinThreshold` + - C2: beam compensation not exceeding `MaxAngleOneWayCompression` + - C3: local maximum in power-like metric (Plike) + - C4: sample index above bottom (optionally provided) with `FloorFilterTolerance` + - C5: exclude early samples (near transducer) using a fixed guard index + + Note: edge samples are not evaluated (requires left/right neighbours). + + 3) Condition 6: angle stability + - Reject peaks whose along/athwart angle standard deviation within ±2 samples + exceeds `MaxPhaseDeviation` (in phase-step units after sensitivity scaling). + + 4) Condition 7: echo-length + - Define pulse envelope at -6 dB relative to the peak in Plike. + - Keep candidates whose envelope length is within + [`MinEchoLength`, `MaxEchoLength`] expressed in pulse-length units. + + 5) Phase II: refinement and deconfliction (Condition 8) + - Refine target position using a weighted centroid within the -6 dB envelope. + - Convert refined index to range (m) and enforce minimum spacing + `MinEchoSpace` (pulse-length units). + - Apply depth window [`MinEchoDepthM`, `MaxEchoDepthM`] using + `transducer_depth` and `heave_compensation`. (to check further) + + Parameters + ---------- + ds : xr.Dataset + Dataset containing Sv (dB) and optional angle/metadata variables. + params : dict + Detection parameters (see PARAM_DEFAULTS). Must include: + - channel, pulse_length, timeSampleInterval_usec + Optional: sound_speed, alpha, bottom_da, transducer_depth, + heave_compensation, etc. (add all) + + Returns + ------- + xr.Dataset + Per-target outputs (TS_comp, TS_uncomp, angles, indices, time/ping number, etc.). + Empty dataset is returned when no targets are detected. + """ + + Param = _merge_params(params) + DBG = _make_dbg() + + channel, Sv = _validate_inputs(ds, Param) + + # extract variables + ex = _extract_matrices(ds, Sv, channel) + sv_mat = ex["sv_mat"] + al_mat = ex["al_mat"] + ath_mat = ex["ath_mat"] + depth = ex["depth"] + npings = ex["npings"] + + # extract Delta + metadata + angle normalisation + delta = _compute_delta(depth) + meta = _extract_metadata(ds, Param, params, channel, npings, DBG) + al_mat, ath_mat, meta = _normalise_angles(al_mat, ath_mat, meta, DBG) + + # compute derived pulse quantities + derived = _compute_pulse_derived(meta, Param, delta) + + # get core matrices (ts_range, tsu, Plike, ts) + constraints + core = _compute_core_matrices( + sv_mat=sv_mat, + al_mat=al_mat, + ath_mat=ath_mat, + depth=depth, + delta=delta, + meta=meta, + Param=Param, + channel=channel, + ) + + # pass conditions + i1, ind = _cond_1to5(core, Param, DBG) + if ind.size == 0: + return _empty_out(channel) + + ind2 = _cond_6_phase(i1, ind, al_mat, ath_mat, core, meta, Param, DBG) + if ind2.size == 0: + return _empty_out(channel) + + c7 = _cond_7_echo_length(i1, ind, ind2, core, derived, Param, DBG) + if c7["ind4"].size == 0: + return _empty_out(channel) + + out, TSfreq_rows = _cond_8_build_outputs( + ds=ds, + params=params, + i1=i1, + ind2=ind2, + c7=c7, + core=core, + meta=meta, + derived=derived, + al_mat=al_mat, + ath_mat=ath_mat, + depth=depth, + Param=Param, + DBG=DBG, + delta=delta, + ) + + out["nb_valid_targets"] = len(out["TS_comp"]) + out["TSfreq_matrix"] = ( + np.asarray(TSfreq_rows, dtype=float) if TSfreq_rows else np.empty((0, 23)) + ) return _matecho_struct_to_dataset(out, channel=channel) From a6a8282e4c7a5b66512d81d58ee5f466fc2471a9 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Wed, 24 Dec 2025 08:26:27 +0100 Subject: [PATCH 5/6] Update test_mask.py fix tests to match updated single-target detection code --- echopype/tests/mask/test_mask.py | 102 +++++++++++++++++++------------ 1 file changed, 63 insertions(+), 39 deletions(-) diff --git a/echopype/tests/mask/test_mask.py b/echopype/tests/mask/test_mask.py index 32fabcd8f..a57c42b10 100644 --- a/echopype/tests/mask/test_mask.py +++ b/echopype/tests/mask/test_mask.py @@ -2013,6 +2013,12 @@ def test_echoview_mincan_no_linking(): # test for single target detections +REQ = { + "channel": "chan1", + "pulse_length": 1e-3, # 1 ms dummy + "timeSampleInterval_usec": 40.0, # dummy +} + def _make_ds_single_target_stub( n_ping=5, n_range=10, channels=("chan1",), var_name="Sv" ) -> xr.Dataset: @@ -2042,13 +2048,21 @@ def _make_ds_single_target_wrong_dims( def test_detect_single_targets_var_name_missing_in_ds_raises(): ds = _make_ds_single_target_stub(var_name="Sv") # dataset has Sv with pytest.raises(ValueError, match=r"var_name 'Sv_corrected' not found"): - detect_single_targets(ds, method="matecho", params={"channel": "chan1", "var_name": "Sv_corrected"}) + detect_single_targets( + ds, + method="matecho", + params={**REQ, "var_name": "Sv_corrected"}, + ) @pytest.mark.unit def test_detect_single_targets_default_var_name_Sv_ok(): ds = _make_ds_single_target_stub(var_name="Sv") - out = detect_single_targets(ds, method="matecho", params={"channel": "chan1"}) # no var_name -> default Sv + out = detect_single_targets( + ds, + method="matecho", + params=REQ, # no var_name -> default Sv + ) assert isinstance(out, xr.Dataset) assert "target" in out.dims assert out.sizes["target"] == 0 @@ -2058,21 +2072,29 @@ def test_detect_single_targets_default_var_name_Sv_ok(): def test_detect_single_targets_var_name_wrong_dims_raises(): ds = _make_ds_single_target_wrong_dims(var_name="Sv") with pytest.raises(ValueError, match=r"Sv must have dims"): - detect_single_targets(ds, method="matecho", params={"channel": "chan1", "var_name": "Sv"}) + detect_single_targets( + ds, + method="matecho", + params={**REQ, "var_name": "Sv"}, + ) @pytest.mark.unit def test_detect_single_targets_channel_not_in_ds_raises(): ds = _make_ds_single_target_stub(channels=("chan1",)) with pytest.raises(ValueError, match=r"Channel 'chanX' not found"): - detect_single_targets(ds, method="matecho", params={"channel": "chanX", "var_name": "Sv"}) + detect_single_targets( + ds, + method="matecho", + params={**REQ, "channel": "chanX", "var_name": "Sv"}, + ) @pytest.mark.unit def test_detect_single_targets_unknown_method_raises(): ds = _make_ds_single_target_stub() with pytest.raises(ValueError, match="Unsupported single-target method"): - detect_single_targets(ds, method="__bad__", params={"channel": "chan1"}) + detect_single_targets(ds, method="__bad__", params={**REQ}) @pytest.mark.unit @@ -2086,22 +2108,24 @@ def test_detect_single_targets_params_required_raises(): def test_detect_single_targets_missing_channel_raises(): ds = _make_ds_single_target_stub() with pytest.raises(ValueError, match=r"params\['channel'\] is required"): - detect_single_targets(ds, method="matecho", params={}) + p = dict(REQ) + p.pop("channel", None) + detect_single_targets(ds, method="matecho", params=p) @pytest.mark.unit def test_detect_single_targets_returns_dataset_with_target_dim_zero(): ds = _make_ds_single_target_stub(channels=("chan1",)) - out = detect_single_targets(ds, method="matecho", params={"channel": "chan1"}) + out = detect_single_targets(ds, method="matecho", params=REQ) assert isinstance(out, xr.Dataset) assert "target" in out.dims assert out.sizes["target"] == 0 - assert "ping_time" in out.coords - assert "ping_number" in out.coords - assert out["ping_time"].dims == ("target",) - assert out["ping_number"].dims == ("target",) + assert "Time" in out.coords + assert "Ping_number" in out.coords + assert out["Time"].dims == ("target",) + assert out["Ping_number"].dims == ("target",) assert "channel" in out.coords assert out["channel"].ndim == 0 # scalar coordinate @@ -2110,43 +2134,48 @@ def test_detect_single_targets_returns_dataset_with_target_dim_zero(): @pytest.mark.unit def test_detect_single_targets_schema_variables_present_even_if_empty(): ds = _make_ds_single_target_stub() - out = detect_single_targets(ds, method="matecho", params={"channel": "chan1"}) + out = detect_single_targets(ds, method="matecho", params=REQ) - # Compare with the variables _matecho_struct_to_dataset currently creates expected_vars = [ "TS_comp", "TS_uncomp", - "target_range", - "target_range_disp", - "target_range_min", - "target_range_max", + "Target_range", + "Target_range_disp", + "Target_range_min", + "Target_range_max", "idx_r", "idx_target_lin", "pulse_env_before_lin", "pulse_env_after_lin", - "pulse_length_normalized_pldl", - "transmitted_pulse_length", - "angle_minor_axis", - "angle_major_axis", - "std_angle_minor_axis", - "std_angle_major_axis", - "heave", - "roll", - "pitch", - "heading", - "dist", + "PulseLength_Normalized_PLDL", + "Transmitted_pulse_length", + "Angle_minor_axis", + "Angle_major_axis", + "StandDev_Angles_Minor_Axis", + "StandDev_Angles_Major_Axis", + "Heave", + "Roll", + "Pitch", + "Heading", + "Dist", + "TSfreq_matrix", ] for v in expected_vars: assert v in out.data_vars - assert out[v].dims == ("target",) - assert out[v].sizes["target"] == 0 + # TSfreq_matrix is 2D even if empty; the rest are 1D on target + if v == "TSfreq_matrix": + assert out[v].dims[0] == "target" + assert out[v].sizes["target"] == 0 + else: + assert out[v].dims == ("target",) + assert out[v].sizes["target"] == 0 @pytest.mark.unit def test_detect_single_targets_concat_empty_is_stable(): ds = _make_ds_single_target_stub() - out1 = detect_single_targets(ds, method="matecho", params={"channel": "chan1"}) - out2 = detect_single_targets(ds, method="matecho", params={"channel": "chan1"}) + out1 = detect_single_targets(ds, method="matecho", params=REQ) + out2 = detect_single_targets(ds, method="matecho", params=REQ) outc = xr.concat([out1, out2], dim="target") assert "target" in outc.dims @@ -2156,10 +2185,5 @@ def test_detect_single_targets_concat_empty_is_stable(): @pytest.mark.unit def test_detect_single_targets_missing_optional_metadata_warns(): ds = _make_ds_single_target_stub() - - with pytest.warns(UserWarning, match="Defaults will be used"): - detect_single_targets( - ds, - method="matecho", - params={"channel": "chan1"}, - ) + with pytest.warns(UserWarning, match=r"Missing beam metadata"): + detect_single_targets(ds, method="matecho", params=REQ) \ No newline at end of file From 75732ce1e178c7dec7fa2fadb0136d3a11d7f394 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Tue, 24 Feb 2026 00:37:57 -0800 Subject: [PATCH 6/6] Reshape Matecho detector, add EV Method 2, propagate tau_effective - Refactor / reshape Matecho CW detector interface and outputs - Propagate tau_effective from compute_Sv (complex encoding path) (Power encoding support still pending) - Standardise detector return structure (ping_time, range_sample, frequency_nominal - Additional diagnostic variables (TS, etc.) temporarily retained - Add minimal Echoview split-beam Method 2 implementation - Update and refactor associated tests --- echopype/calibrate/calibrate_ek.py | 14 + echopype/mask/api.py | 47 +- .../detect_echoview_split_method2.py | 508 +++++ .../single_target_detection/detect_matecho.py | 1658 +++++++++-------- echopype/tests/mask/test_mask.py | 468 +++-- 5 files changed, 1747 insertions(+), 948 deletions(-) create mode 100644 echopype/mask/single_target_detection/detect_echoview_split_method2.py diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index c75d873b9..6159e77c2 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -581,6 +581,20 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: # Attach calculated range (with units meter) into data set out = out.to_dataset().merge(range_meter) + # Add tau_effective to output (Sv only) + if cal_type == "Sv": + out["tau_effective"] = tau_effective + out["tau_effective"].attrs.update( + { + "long_name": "Effective pulse length", + "units": "s", + "description": ( + "Effective pulse length used in Sv calibration. " + "For GPT channels, transmit_duration_nominal is used." + ), + } + ) + # Add frequency_nominal to data set out["frequency_nominal"] = beam["frequency_nominal"] diff --git a/echopype/mask/api.py b/echopype/mask/api.py index 1f71f49ec..7b1e5d12a 100644 --- a/echopype/mask/api.py +++ b/echopype/mask/api.py @@ -12,6 +12,9 @@ # for seafloor detection from echopype.mask.seafloor_detection.bottom_basic import bottom_basic from echopype.mask.seafloor_detection.bottom_blackwell import bottom_blackwell +from echopype.mask.single_target_detection.detect_echoview_split_method2 import ( + detect_echoview_split_method2, +) # for single_target_detection from echopype.mask.single_target_detection.detect_matecho import detect_matecho @@ -807,6 +810,7 @@ def detect_shoal( # Registry of supported methods for single_target_detection METHODS_SINGLE_TARGET = { "matecho": detect_matecho, + "echoview_split_method2": detect_echoview_split_method2, } @@ -822,11 +826,11 @@ def detect_single_targets( ---------- ds : xr.Dataset Acoustic dataset containing the fields required by the selected method. - Typical dimensions include ``ping_time`` and ``range_sample`` (or ``depth``). - Depending on the method, this may require TS, Sv, and optionally split-beam + Typical dimensions include ``ping_time`` and ``range_sample``. + Depending on the method, this may require TS, Sv and split-beam angles (e.g. alongship / athwartship angles). method : str - Name of the detection method to use (e.g., ``"matecho"``, ``"esp3"``, ...). + Name of the detection method to use (e.g., ``"matecho"``, ...). params : dict Method-specific parameters. This argument is required and no defaults are assumed. @@ -834,20 +838,49 @@ def detect_single_targets( Returns ------- xr.Dataset - Per-target results with dimension ``target`` and coordinates - ``ping_time`` and ``ping_number``. Variables include compensated and - uncompensated TS, range metrics, and optional - navigation/attitude fields when available. + Per-target detection results with dimension ``target``. + + Coordinates (defined on ``target``) are: + - ``ping_time`` + - ``range_sample`` + - ``frequency_nominal`` + + Each row corresponds to a single detection occurring at one + time and one range sample. The frequency may be constant + (CW data) or vary per target (FM data). + """ + if method not in METHODS_SINGLE_TARGET: raise ValueError(f"Unsupported single-target method: {method}") if params is None: raise ValueError("No parameters given.") + if "beam_type" not in ds: + raise ValueError("beam_type variable is missing from dataset.") + + beam_vals = np.unique(ds["beam_type"].values) + + if not np.all(np.isin(beam_vals, [1, 65])): + raise ValueError( + f"Only split-beam data supported (beam_type 1 or 65). " f"Found: {beam_vals}" + ) + out = METHODS_SINGLE_TARGET[method](ds, params) if not isinstance(out, xr.Dataset) or "target" not in out.dims: raise TypeError(f"{method} must return an xr.Dataset with a 'target' dimension.") + required = ("ping_time", "range_sample", "frequency_nominal") + missing = [v for v in required if v not in out] + if missing: + raise ValueError( + f"{method} output missing required field(s): {missing} (expected {list(required)})." + ) + + bad_dims = [v for v in required if out[v].dims != ("target",)] + if bad_dims: + raise ValueError(f"{method} field(s) must have dims ('target',): {bad_dims}.") + return out diff --git a/echopype/mask/single_target_detection/detect_echoview_split_method2.py b/echopype/mask/single_target_detection/detect_echoview_split_method2.py new file mode 100644 index 000000000..c56a7e5c9 --- /dev/null +++ b/echopype/mask/single_target_detection/detect_echoview_split_method2.py @@ -0,0 +1,508 @@ +from __future__ import annotations + +import numpy as np +import xarray as xr + +# ## +# Params validation +# ## + +REQUIRED_PARAMS = { + "ts_threshold_db", + "pldl_db", + "min_norm_pulse", + "max_norm_pulse", + "beam_comp_model", + "max_beam_comp_db", + "max_sd_minor_deg", + "max_sd_major_deg", +} + +OPTIONAL_PARAMS = { + "dec_tir_samples", + "bottom_offset_m", + "exclude_above_m", + "exclude_below_m", + "allow_nans_inside_envelope", +} + + +def _validate_params(params: dict) -> dict: + if params is None: + raise ValueError("No parameters given.") + unknown = set(params.keys()) - (REQUIRED_PARAMS | OPTIONAL_PARAMS) + if unknown: + raise ValueError(f"Unknown parameters: {sorted(unknown)}") + missing = REQUIRED_PARAMS - set(params.keys()) + if missing: + raise ValueError(f"Missing required parameters: {sorted(missing)}") + if params["min_norm_pulse"] > params["max_norm_pulse"]: + raise ValueError("min_norm_pulse must be <= max_norm_pulse") + return params + + +def _validate_ds_m2(ds_m2: xr.Dataset) -> xr.Dataset: + must = [ + "TS", + "echo_range", + "angle_alongship", + "angle_athwartship", + "sound_absorption", + "sample_interval", + "tau_effective", + ] + for v in must: + if v not in ds_m2: + raise ValueError(f"ds_m2 missing required variable: {v}") + + ts_mat = ds_m2["TS"] + if ts_mat.dims != ("ping_time", "range_sample"): + raise ValueError("Expected ds_m2['TS'] dims exactly ('ping_time','range_sample').") + for v in ["echo_range", "angle_alongship", "angle_athwartship"]: + if ds_m2[v].dims != ("ping_time", "range_sample"): + raise ValueError(f"Expected ds_m2['{v}'] dims exactly ('ping_time','range_sample').") + + # Broadcast alpha to 2D if needed (for simple math) + alpha = ds_m2["sound_absorption"] + if alpha.ndim == 1: + ds_m2 = ds_m2.assign(sound_absorption=alpha.broadcast_like(ts_mat)) + elif alpha.ndim == 0: + ds_m2 = ds_m2.assign(sound_absorption=(xr.zeros_like(ts_mat) + alpha)) + elif alpha.ndim == 2: + pass + else: + raise ValueError("sound_absorption must be scalar, 1D(ping_time), or 2D like TS.") + + return ds_m2 + + +# ## +# Core algorithm equations +# ## + + +def _plike_from_ts( + ts_db: xr.DataArray, r_m: xr.DataArray, alpha_db_m: xr.DataArray +) -> xr.DataArray: + # Plike = TS - 40log10(r) - 2 alpha r + r = xr.where(r_m > 0, r_m, np.nan) + return ts_db - 40.0 * xr.apply_ufunc(np.log10, r) - 2.0 * alpha_db_m * r_m + + +def _local_max_first_plateau(plike_mat: xr.DataArray) -> xr.DataArray: + prev = plike_mat.shift(range_sample=1) + nxt = plike_mat.shift(range_sample=-1) + peak = (plike_mat > prev) & (plike_mat >= nxt) & ~(plike_mat == prev) + peak = peak & xr.apply_ufunc(np.isfinite, prev) & xr.apply_ufunc(np.isfinite, nxt) + peak = peak & xr.apply_ufunc(np.isfinite, plike_mat) + return peak + + +def _nech_p_samples(ds_m2: xr.Dataset) -> np.ndarray: + # NechP (samples) = tau_effective / sample_interval + tau = ds_m2["tau_effective"] + dt = ds_m2["sample_interval"] + + # ---- tau -> 1D per ping_time ---- + if tau.ndim == 0: + tau_vec = np.full(ds_m2.sizes["ping_time"], float(tau.values), dtype=float) + elif tau.ndim == 1 and tau.dims == ("ping_time",): + tau_vec = tau.values.astype(float) + else: + tau_vec = tau.isel(range_sample=0).values.astype(float) + + if dt.ndim == 0: + dt_vec = np.full(ds_m2.sizes["ping_time"], float(dt.values), dtype=float) + elif dt.ndim == 1 and dt.dims == ("ping_time",): + dt_vec = dt.values.astype(float) + else: + dt_vec = dt.isel(range_sample=0).values.astype(float) + + sample_interval_sec = dt_vec + + return tau_vec / sample_interval_sec + + +def _envelope_bounds_1d( + plike_row: np.ndarray, p: int, thr: float, allow_nans: bool +) -> tuple[int | None, int | None]: + # expand left + m = p + while m > 0: + v = plike_row[m - 1] + if not np.isfinite(v): + return (None, None) if not allow_nans else (m, p) + if v >= thr: + m -= 1 + else: + break + + # expand right + last = p + while last < plike_row.size - 1: + v = plike_row[last + 1] + if not np.isfinite(v): + return (None, None) if not allow_nans else (m, last) + if v >= thr: + last += 1 + else: + break + + return m, last + + +# ## +# Beam compensation +# ## + + +def _beam_comp_db(ds_m2: xr.Dataset, params: dict) -> xr.DataArray: + + model = params["beam_comp_model"] + + if model == "none": + return xr.zeros_like(ds_m2["TS"]) + + elif model == "simrad_lobe": + + th_al = ds_m2["angle_alongship"] + th_at = ds_m2["angle_athwartship"] + + bw_al = ds_m2["beamwidth_alongship"].broadcast_like(th_al) + bw_at = ds_m2["beamwidth_athwartship"].broadcast_like(th_at) + + off_al = ds_m2["angle_offset_alongship"].broadcast_like(th_al) + off_at = ds_m2["angle_offset_athwartship"].broadcast_like(th_at) + + x = 2 * (th_al - off_al) / bw_al + y = 2 * (th_at - off_at) / bw_at + + beam_comp_db = 6.0206 * (x**2 + y**2 - 0.18 * x**2 * y**2) + + return beam_comp_db.broadcast_like(ds_m2["TS"]) + + elif model == "provided": + if "beam_comp_db" not in ds_m2: + raise ValueError("beam_comp_db must exist if beam_comp_model='provided'") + return ds_m2["beam_comp_db"].broadcast_like(ds_m2["TS"]) + + else: + raise ValueError(f"Unknown beam_comp_model: {model}") + + +# ## +# Phase I +# ## + + +def _phase1_simple(ds_m2: xr.Dataset, params: dict, beam_comp_db: xr.DataArray) -> xr.Dataset: + """ + Returns a compact xr.Dataset with dim 'target' and vars: + ping_index, range_sample, iinf, isup, pulse_len_samples, norm_pulse_len, plike_peak + """ + plike_mat = _plike_from_ts(ds_m2["TS"], ds_m2["echo_range"], ds_m2["sound_absorption"]) + + cand_mask = _local_max_first_plateau(plike_mat) + cand_mask = cand_mask & (beam_comp_db <= float(params["max_beam_comp_db"])) + + # optional gates + if params.get("dec_tir_samples") is not None: + dec_tir = int(params["dec_tir_samples"]) + idx = xr.DataArray( + np.arange(plike_mat.sizes["range_sample"]), + dims=("range_sample",), + coords={"range_sample": plike_mat["range_sample"]}, + ) + cand_mask = cand_mask & (idx >= dec_tir) + + if params.get("exclude_above_m") is not None: + cand_mask = cand_mask & (ds_m2["echo_range"] >= float(params["exclude_above_m"])) + if params.get("exclude_below_m") is not None: + cand_mask = cand_mask & (ds_m2["echo_range"] <= float(params["exclude_below_m"])) + + if "bottom" in ds_m2 and params.get("bottom_offset_m") is not None: + off = float(params["bottom_offset_m"]) + bottom2d = ds_m2["bottom"].broadcast_like(ds_m2["TS"]) + cand_mask = cand_mask & (ds_m2["echo_range"] <= (bottom2d - off)) + + # numpy loop (envelopes) + plike_np = plike_mat.values + cand_np = cand_mask.values + al_np = ds_m2["angle_alongship"].values + ath_np = ds_m2["angle_athwartship"].values + + nech_p = _nech_p_samples(ds_m2) + pldl_db = float(params["pldl_db"]) + min_norm_pulse = float(params["min_norm_pulse"]) + max_norm_pulse = float(params["max_norm_pulse"]) + max_sd_minor_deg = float(params["max_sd_minor_deg"]) + max_sd_major_deg = float(params["max_sd_major_deg"]) + allow_nans = bool(params.get("allow_nans_inside_envelope", False)) + + ( + ping_index_list, + range_sample_list, + iinf_list, + isup_list, + pulse_len_samples_list, + norm_pulse_len_list, + plike_peak_list, + ) = ([] for _ in range(7)) + + for it in range(plike_np.shape[0]): + peaks = np.where(cand_np[it])[0] + if peaks.size == 0: + continue + + nch = nech_p[it] + if not np.isfinite(nch) or nch <= 0: + continue + + plike_row = plike_np[it] + ali = al_np[it] + athi = ath_np[it] + + for p in peaks: + plike_peak = plike_row[p] + if not np.isfinite(plike_peak): + continue + + iinf, isup = _envelope_bounds_1d( + plike_row, int(p), plike_peak - pldl_db, allow_nans=allow_nans + ) + if iinf is None: + continue + + pulse_len_samples = isup - iinf + 1 + norm_pulse_len = pulse_len_samples / nch + if norm_pulse_len < min_norm_pulse or norm_pulse_len > max_norm_pulse: + continue + + seg_al = ali[iinf : isup + 1] * 180.0 / np.pi + seg_ath = athi[iinf : isup + 1] * 180.0 / np.pi + if np.nanstd(seg_ath) > max_sd_minor_deg: + continue + if np.nanstd(seg_al) > max_sd_major_deg: + continue + + ping_index_list.append(it) + range_sample_list.append(int(p)) + iinf_list.append(int(iinf)) + isup_list.append(int(isup)) + pulse_len_samples_list.append(int(pulse_len_samples)) + norm_pulse_len_list.append(float(norm_pulse_len)) + plike_peak_list.append(float(plike_peak)) + + if len(range_sample_list) == 0: + return xr.Dataset(coords={"target": np.arange(0, dtype=np.int64)}) + + return xr.Dataset( + data_vars=dict( + ping_index=("target", np.array(ping_index_list, dtype=np.int64)), + range_sample=("target", np.array(range_sample_list, dtype=np.int64)), + iinf=("target", np.array(iinf_list, dtype=np.int64)), + isup=("target", np.array(isup_list, dtype=np.int64)), + pulse_len_samples=("target", np.array(pulse_len_samples_list, dtype=np.int64)), + norm_pulse_len=("target", np.array(norm_pulse_len_list, dtype=np.float64)), + plike_peak=("target", np.array(plike_peak_list, dtype=np.float64)), + ), + coords=dict(target=np.arange(len(range_sample_list), dtype=np.int64)), + ) + + +# ## +# Phase II (TS computation + threshold + overlap rejection) +# ## + + +def _compute_ts_at_peaks( + feats: xr.Dataset, ds_m2: xr.Dataset, beam_comp_db: xr.DataArray +) -> xr.Dataset: + it = feats["ping_index"].values.astype(np.int64) + p = feats["range_sample"].values.astype(np.int64) + + r = ds_m2["echo_range"].values[it, p] + a = ds_m2["sound_absorption"].values[it, p] + plike_peak = feats["plike_peak"].values + beam_comp_p = beam_comp_db.values[it, p] + + ts_uncomp = plike_peak + 40.0 * np.log10(r) + 2.0 * a * r + ts_comp = ts_uncomp + beam_comp_p + + return feats.assign( + target_range=("target", r.astype(float)), + ts_uncomp=("target", ts_uncomp.astype(float)), + ts_comp=("target", ts_comp.astype(float)), + ) + + +def _reject_overlaps_per_ping(feats: xr.Dataset) -> xr.Dataset: + if feats.dims.get("target", 0) <= 1: + return feats + + ping_idx = feats["ping_index"].values + iinf = feats["iinf"].values + isup = feats["isup"].values + ts_comp = feats["ts_comp"].values + + keep = np.ones(feats.dims["target"], dtype=bool) + + for it in np.unique(ping_idx): + ii = np.where(ping_idx == it)[0] + if ii.size <= 1: + continue + + order = ii[np.argsort(iinf[ii])] + accepted = [] + + for j in order: + if not accepted: + accepted.append(j) + continue + + k = accepted[-1] + if iinf[j] > isup[k]: + accepted.append(j) + continue + + # overlap: reject lower ts_comp + if ts_comp[j] >= ts_comp[k]: + keep[k] = False + accepted[-1] = j + else: + keep[j] = False + + return feats.isel(target=keep) + + +def _phase2( + ds_m2: xr.Dataset, params: dict, beam_comp_db: xr.DataArray, feats: xr.Dataset +) -> xr.Dataset: + if feats.dims.get("target", 0) == 0: + return feats + + feats = _compute_ts_at_peaks(feats, ds_m2, beam_comp_db) + + # EV Method2: threshold applied to compensated TS + keep0 = feats["ts_comp"] >= float(params["ts_threshold_db"]) + feats = feats.isel(target=keep0) + + if feats.dims.get("target", 0) == 0: + return feats + + feats = _reject_overlaps_per_ping(feats) + return feats + + +# # ## +# Output packing +# ## + + +def _pack_targets(feats: xr.Dataset, ds_m2: xr.Dataset) -> xr.Dataset: + if feats.sizes.get("target", 0) == 0: + return xr.Dataset(coords={"target": np.arange(0, dtype=np.int64)}) + + it = feats["ping_index"].values.astype(np.int64) + p = feats["range_sample"].values.astype(np.int64) + + ping_time = ds_m2["ping_time"].values[it] + angle_major_deg = ds_m2["angle_alongship"].values[it, p] * 180.0 / np.pi + angle_minor_deg = ds_m2["angle_athwartship"].values[it, p] * 180.0 / np.pi + + # add constant per target for single-channel ds_m2 + fn = ds_m2["frequency_nominal"] + if fn.ndim == 0: + freq_val = float(fn.values) + else: + freq_val = float(fn.values[0]) + frequency_nominal = np.full(it.shape[0], freq_val, dtype=np.float64) + + return xr.Dataset( + data_vars=dict( + ping_time=("target", ping_time), + range_sample=("target", p), + frequency_nominal=("target", frequency_nominal), + ping_index=("target", it), + iinf=("target", feats["iinf"].values.astype(np.int64)), + isup=("target", feats["isup"].values.astype(np.int64)), + pulse_len_samples=("target", feats["pulse_len_samples"].values.astype(np.int64)), + norm_pulse_len=("target", feats["norm_pulse_len"].values.astype(np.float64)), + target_range=("target", feats["target_range"].values.astype(np.float64)), + ts_uncomp=("target", feats["ts_uncomp"].values.astype(np.float64)), + ts_comp=("target", feats["ts_comp"].values.astype(np.float64)), + angle_major_deg=("target", angle_major_deg.astype(np.float64)), + angle_minor_deg=("target", angle_minor_deg.astype(np.float64)), + ), + coords=dict(target=np.arange(feats.sizes["target"], dtype=np.int64)), + attrs=dict(method="echoview_split_method2"), + ) + + +# ## +# Public API +# ## + + +def detect_echoview_split_method2(ds_m2: xr.Dataset, params: dict) -> xr.Dataset: + """ + Echoview split-beam single-target detection (Method 2) port. + + Implements in python the EV Method 2 workflow on precomputed TS (dB) and split-beam angles: + (i) compute Plike from TS + range + absorption, (ii) find local Plike maxima and build + a -PLDL envelope, (iii) apply Phase I gates (beam-comp limit, echo-length / normalised + pulse length, angle std-dev), then (iv) compute TS at peaks, apply the compensated TS + threshold, and reject overlaps within each ping. + + Notes: + - Input TS must be computed upstream. + - Angles/beam geometry are converted from degrees to radians for the beam compensation model. + + Parameters + ---------- + ds_m2 : xr.Dataset + Single-channel dataset with 2D fields (ping_time, range_sample), at minimum: + TS, echo_range, angle_alongship, angle_athwartship, sound_absorption, + sample_interval, tau_effective, plus beam geometry terms if beam compensation is used. + + params : dict + Required keys: + - ts_threshold_db, pldl_db + - min_norm_pulse, max_norm_pulse + - beam_comp_model ('none'|'simrad_lobe'|'provided'), max_beam_comp_db + - max_sd_minor_deg, max_sd_major_deg + Optional keys: + - dec_tir_samples, bottom_offset_m, exclude_above_m, exclude_below_m, + allow_nans_inside_envelope + + Returns + ------- + xr.Dataset + Dataset with `target` detections including ts_comp/ts_uncomp, target_range, + ping_time, range_sample, angles, and method metadata. + """ + params = _validate_params(params) + ds_m2 = _validate_ds_m2(ds_m2) + + # ADAPTATION: degrees -> radians + deg2rad = np.pi / 180.0 + + ds_m2 = ds_m2.copy() + + ds_m2["angle_alongship"] = ds_m2["angle_alongship"] * deg2rad + ds_m2["angle_athwartship"] = ds_m2["angle_athwartship"] * deg2rad + + ds_m2["beamwidth_alongship"] = ds_m2["beamwidth_alongship"] * deg2rad + ds_m2["beamwidth_athwartship"] = ds_m2["beamwidth_athwartship"] * deg2rad + + ds_m2["angle_offset_alongship"] = ds_m2["angle_offset_alongship"] * deg2rad + ds_m2["angle_offset_athwartship"] = ds_m2["angle_offset_athwartship"] * deg2rad + ### + + beam_comp_db = _beam_comp_db(ds_m2, params) + + feats = _phase1_simple(ds_m2, params, beam_comp_db) + feats = _phase2(ds_m2, params, beam_comp_db, feats) + + out = _pack_targets(feats, ds_m2) + return out diff --git a/echopype/mask/single_target_detection/detect_matecho.py b/echopype/mask/single_target_detection/detect_matecho.py index 8fece863f..0e8086c1a 100644 --- a/echopype/mask/single_target_detection/detect_matecho.py +++ b/echopype/mask/single_target_detection/detect_matecho.py @@ -1,551 +1,441 @@ -import warnings -from typing import Any, Dict, List, Tuple +from __future__ import annotations + +from typing import Any, Dict import numpy as np import xarray as xr -# --- -# Defaults -PARAM_DEFAULTS = { - "SoundSpeed": 1500.0, - "MinThreshold": -60.0, - "MaxAngleOneWayCompression": 6.0, - "MaxPhaseDeviation": 8.0, - "MinEchoLength": 0.8, - "MaxEchoLength": 1.8, - "MinEchoSpace": 1.0, - "MinEchoDepthM": 3.0, - "MaxEchoDepthM": 38.0, - "FloorFilterTolerance": 0.0, - "tvg_start_sample": 3, # (EK60=3, EK80=1) - # Sv->TS - "psi_two_way": 0.0, - "Sa_correction": 0.0, - "Sa_EK80_nominal": 0.0, - # Required for CW range/TS conversion - "pulse_length": None, # seconds (pd in MATLAB) - "sound_speed": None, # m/s (c in MATLAB) - "alpha": 0.0, # dB/m (alpha(kf)) - # TD/CH - "transducer_depth": None, # TD (m), from surface - "heave_compensation": 0.0, # CH (m) per ping or scalar - # (optional) - "bottom_da": None, -} - - -# --- -# Output conversion -def _matecho_struct_to_dataset(out: dict, *, channel: str | None = None) -> xr.Dataset: - n = int(out.get("nb_valid_targets", len(out.get("TS_comp", [])))) - target = np.arange(n, dtype=int) - def _arr(key, dtype=float): - vals = out.get(key, []) - if n == 0: - return np.asarray([], dtype=dtype) - a = np.asarray(vals) - if dtype is not None: - a = a.astype(dtype) - if a.size != n: - raise ValueError(f"Output field '{key}' has length {a.size}, expected {n}.") - return a +def _compute_core_matrices(sv_mat, al_mat, ath_mat, meta, params, derived): + sv = np.asarray(sv_mat, dtype=np.float64) + al = np.asarray(al_mat, dtype=np.float64) + ath = np.asarray(ath_mat, dtype=np.float64) - tsfreq = np.asarray(out.get("TSfreq_matrix", np.empty((0, 23))), dtype=float) - if tsfreq.ndim == 1: - tsfreq = tsfreq.reshape(1, -1) - if tsfreq.size == 0: - tsfreq = np.empty((0, 23), dtype=float) - if tsfreq.shape[0] != n: - raise ValueError( - f"TSfreq_matrix has {tsfreq.shape[0]} rows but expected {n} (nb_valid_targets)." - ) + _, n_samp = sv.shape - ds_out = xr.Dataset( - data_vars=dict( - nb_valid_targets=((), n), - TS_comp=("target", _arr("TS_comp", float), {"units": "dB"}), - TS_uncomp=("target", _arr("TS_uncomp", float), {"units": "dB"}), - Target_range=("target", _arr("Target_range", float), {"units": "m"}), - Target_range_disp=("target", _arr("Target_range_disp", float), {"units": "m"}), - Target_range_min=("target", _arr("Target_range_min", float), {"units": "m"}), - Target_range_max=("target", _arr("Target_range_max", float), {"units": "m"}), - idx_r=("target", _arr("idx_r", int)), - idx_target_lin=("target", _arr("idx_target_lin", int)), - pulse_env_before_lin=("target", _arr("pulse_env_before_lin", float)), - pulse_env_after_lin=("target", _arr("pulse_env_after_lin", float)), - PulseLength_Normalized_PLDL=("target", _arr("PulseLength_Normalized_PLDL", float)), - Transmitted_pulse_length=("target", _arr("Transmitted_pulse_length", float)), - Angle_minor_axis=("target", _arr("Angle_minor_axis", float), {"units": "rad"}), - Angle_major_axis=("target", _arr("Angle_major_axis", float), {"units": "rad"}), - StandDev_Angles_Minor_Axis=( - "target", - _arr("StandDev_Angles_Minor_Axis", float), - {"units": "phase_steps"}, - ), - StandDev_Angles_Major_Axis=( - "target", - _arr("StandDev_Angles_Major_Axis", float), - {"units": "phase_steps"}, - ), - Heave=("target", _arr("Heave", float), {"units": "m"}), - Roll=("target", _arr("Roll", float), {"units": "rad"}), - Pitch=("target", _arr("Pitch", float), {"units": "rad"}), - Heading=("target", _arr("Heading", float), {"units": "rad"}), - Dist=("target", _arr("Dist", float), {"units": "m"}), - TSfreq_matrix=(("target", "tsfreq_col"), tsfreq, {}), - ), - coords=dict( - target=target, - tsfreq_col=np.arange(tsfreq.shape[1], dtype=int), - Ping_number=("target", _arr("Ping_number", int)), - Time=("target", _arr("Time", np.dtype("datetime64[ns]"))), - ), - attrs=dict( - method="matecho", - channel=str(channel) if channel is not None else "", - nb_valid_targets=n, - ), - ) - if channel is not None: - ds_out = ds_out.assign_coords(channel=channel) - return ds_out + ind_range0 = int(derived["ind_range0"]) + delta = float(derived["delta"]) + sv2ts_f32 = np.float32(derived["sv2ts_f32"]) + alpha = float(meta["alpha"]) + # --- range vectors in float64 (Matecho parity) --- + idx = np.arange(1, n_samp + 1, dtype=np.float64) + ts_range64 = np.maximum(np.float64(delta), (idx - ind_range0) * np.float64(delta)) + log20 = 20.0 * np.log10(ts_range64) + log40 = 40.0 * np.log10(ts_range64) -def _empty_out(channel: str) -> xr.Dataset: - out = { - "nb_valid_targets": 0, - "TS_comp": [], - "TS_uncomp": [], - "Target_range": [], - "Target_range_disp": [], - "Target_range_min": [], - "Target_range_max": [], - "idx_r": [], - "idx_target_lin": [], - "pulse_env_before_lin": [], - "pulse_env_after_lin": [], - "PulseLength_Normalized_PLDL": [], - "Transmitted_pulse_length": [], - "Angle_minor_axis": [], - "Angle_major_axis": [], - "StandDev_Angles_Minor_Axis": [], - "StandDev_Angles_Major_Axis": [], - "Heave": [], - "Roll": [], - "Pitch": [], - "Heading": [], - "Dist": [], - "Ping_number": [], - "Time": [], - "TSfreq_matrix": np.empty((0, 23)), - } - return _matecho_struct_to_dataset(out, channel=channel) + # --- TSU: compute in float64 then cast to float32 (Matecho single) --- + tsu = (sv + log20[None, :] + np.float64(sv2ts_f32)).astype(np.float32) + + # --- Plike: do the casts once --- + ts_range32 = ts_range64.astype(np.float32) + log40_32 = log40.astype(np.float32) + two_alpha32 = np.float32(2.0 * alpha) + + plike = tsu - log40_32[None, :] - two_alpha32 * ts_range32[None, :] + # --- beam compensation (float32) --- + x = (2.0 * (al - float(meta["dec_al"])) / float(meta["ouv_al"])).astype(np.float32) + y = (2.0 * (ath - float(meta["dec_at"])) / float(meta["ouv_at"])).astype(np.float32) + ts = tsu + np.float32(6.0206) * (x * x + y * y - np.float32(0.18) * x * x * y * y) + + # --- bottom gate --- + ir = ( + meta["bot"] + - meta["TD_vec"] + - meta["heave_vec"] + - float(params.get("floor_filter_tolerance_m", 0.0)) + ) / delta + ir = np.clip(ir, 1.0, float(n_samp)) -# --- -# Helpers -def _make_dbg() -> dict: # to print debug return { - "enable": True, - "print_header": True, - "print_block": True, - "print_cond_1to5": True, - "print_cond6": True, - "print_cond7": True, - "print_per_ping": False, - "print_examples": True, - "n_examples": 10, + "tsu_mat": tsu, + "ts_mat": ts, + "plike_mat": plike, + "ir": ir, + "nbSamp": int(n_samp), + "dec_tir": int(params.get("dec_tir", derived.get("dec_tir", 8))), } -def _nanmin(a): - return float(np.nanmin(a)) if np.any(np.isfinite(a)) else np.nan +def _matlab_round_pos(x: float) -> int: + # MATLAB round for positive numbers: halves away from zero + return int(np.floor(x + 0.5)) -def _nanmax(a): - return float(np.nanmax(a)) if np.any(np.isfinite(a)) else np.nan +def _require(ds: xr.Dataset, name: str) -> xr.DataArray: + if name not in ds: + raise ValueError(f"Missing required variable: '{name}'") + return ds[name] -def _nanmed(a): - return float(np.nanmedian(a)) if np.any(np.isfinite(a)) else np.nan +def _require_dims(da: xr.DataArray, name: str, dims: tuple[str, ...]) -> None: + got = tuple(da.dims) + for d in dims: + if d not in got: + raise ValueError(f"Variable '{name}' must have dim '{d}'. Got dims={got}") -def _merge_params(params: dict) -> dict: +def validate_matecho_inputs_strict(ds: xr.Dataset, params: dict) -> tuple[str, xr.DataArray]: if params is None: raise ValueError("params is required.") - return {**PARAM_DEFAULTS, **params} + if "channel" not in params or params["channel"] is None: + raise ValueError("Missing required parameter: params['channel']") + channel = str(params["channel"]) + + if "channel" not in ds.coords: + raise ValueError("Missing required coordinate: 'channel'") + if "ping_time" not in ds.coords and "ping_time" not in ds.dims: + raise ValueError("Missing required coordinate/dimension: 'ping_time'") + if "range_sample" not in ds.coords and "range_sample" not in ds.dims: + raise ValueError("Missing required coordinate/dimension: 'range_sample'") + + if channel not in ds["channel"].values: + raise ValueError(f"Channel '{channel}' not found in ds['channel'].") + + required_3d = ["Sv", "angle_alongship", "angle_athwartship"] + for v in required_3d: + da = _require(ds, v) + _require_dims(da, v, ("channel", "ping_time", "range_sample")) + + if "echo_range" in ds: + axis_name = "echo_range" + else: + raise ValueError("Missing required variable: 'echo_range'") + axis = _require(ds, axis_name) + if "range_sample" not in axis.dims: + raise ValueError( + f"Variable '{axis_name}' must include dim 'range_sample'. Got dims={axis.dims}" + ) -def _validate_inputs(ds: xr.Dataset, Param: dict) -> Tuple[str, xr.DataArray]: - channel = Param.get("channel") - if channel is None: - raise ValueError("params['channel'] is required.") + required_meta = [ + "sound_speed", + "transmit_duration_nominal", + "tau_effective", + "sample_interval", + "sound_absorption", + "equivalent_beam_angle", + "sa_correction", + "beamwidth_alongship", + "beamwidth_athwartship", + "angle_offset_alongship", + "angle_offset_athwartship", + "angle_sensitivity_alongship", + "angle_sensitivity_athwartship", + "transducer_depth", + "heave_compensation", + ] + for v in required_meta: + da = _require(ds, v) + if v in ( + "beamwidth_alongship", + "beamwidth_athwartship", + "angle_offset_alongship", + "angle_offset_athwartship", + "angle_sensitivity_alongship", + "angle_sensitivity_athwartship", + ): + _require_dims(da, v, ("channel",)) + if v in ( + "sound_absorption", + "transmit_duration_nominal", + "transducer_depth", + "heave_compensation", + ): + _require_dims(da, v, ("ping_time", "channel")) + + if params.get("seafloor") is None: + raise ValueError("Need bottom information: provide params['seafloor'].") + + return channel, ds["Sv"] + + +def _extract_matrices( + ds: xr.Dataset, + Sv: xr.DataArray, + channel: str, +) -> Dict[str, Any]: + sv_mat = Sv.sel(channel=channel).transpose("ping_time", "range_sample").values + al_mat = ( + ds["angle_alongship"].sel(channel=channel).transpose("ping_time", "range_sample").values + ) + ath_mat = ( + ds["angle_athwartship"].sel(channel=channel).transpose("ping_time", "range_sample").values + ) - var_name = Param.get("var_name", "Sv") - if var_name not in ds: - raise ValueError(f"var_name '{var_name}' not found in input dataset.") + r_da = ds["echo_range"] if "echo_range" in ds else ds["depth"] + if "channel" in r_da.dims: + r_da = r_da.sel(channel=channel) + axis = np.asarray(r_da.median("ping_time", skipna=True).values, dtype=float).squeeze() - Sv = ds[var_name] - required_dims = {"channel", "ping_time", "range_sample"} - if not required_dims.issubset(set(Sv.dims)): - raise ValueError(f"{var_name} must have dims {sorted(required_dims)}. Got: {Sv.dims}.") + last_idx = int(np.where(np.isfinite(axis))[0].max()) + 1 + sv_mat, al_mat, ath_mat = sv_mat[:, :last_idx], al_mat[:, :last_idx], ath_mat[:, :last_idx] + axis = axis[:last_idx] - if channel not in Sv["channel"].values: - raise ValueError(f"Channel '{channel}' not found in {var_name}.") + delta = float(axis[1] - axis[0]) + depth1 = float(axis[0]) + depth2 = float(axis[1]) - return channel, Sv + # SHIFT (Matecho StartDepthSample alignment) + if "start_depth_m" not in ds: + raise ValueError("Matecho alignment requires ds['start_depth_m'], but it is missing.") + sdm_da = ds["start_depth_m"] + if "channel" in sdm_da.dims: + sdm_da = sdm_da.sel(channel=channel) + sdm_vec = np.asarray(sdm_da.values, dtype=float).squeeze() -def _extract_matrices(ds: xr.Dataset, Sv: xr.DataArray, channel: str) -> Dict[str, Any]: - sv_mat = Sv.sel(channel=channel).transpose("ping_time", "range_sample").values - npings, NbSamp = sv_mat.shape + npings, nbSamp0 = sv_mat.shape - if "angle_alongship" in ds: - al_mat = ( - ds["angle_alongship"].sel(channel=channel).transpose("ping_time", "range_sample").values - ) - else: - al_mat = np.full_like(sv_mat, np.nan, dtype=float) - - if "angle_athwartship" in ds: - ath_mat = ( - ds["angle_athwartship"] - .sel(channel=channel) - .transpose("ping_time", "range_sample") - .values - ) - else: - ath_mat = np.full_like(sv_mat, np.nan, dtype=float) + start_samp_1b = np.floor(sdm_vec / delta).astype(int) + start_samp_1b = np.clip(start_samp_1b, 1, nbSamp0) - if "depth" in ds: - depth = ds["depth"].sel(channel=channel).median("ping_time", skipna=True).values - else: - depth = np.arange(NbSamp, dtype=float) + lens = nbSamp0 - (start_samp_1b - 1) + nbSamp_new = int(np.nanmax(lens)) - if len(depth) < 2: - raise ValueError("Need at least 2 samples to compute delta.") + sv_shift = np.full((npings, nbSamp_new), np.nan, dtype=np.float64) + al_shift = np.full((npings, nbSamp_new), np.nan, dtype=np.float64) + ath_shift = np.full((npings, nbSamp_new), np.nan, dtype=np.float64) - # trim NaN tail - valid = np.isfinite(depth) - if not np.any(valid): - raise ValueError("Depth vector is all-NaN for this channel.") - last = np.where(valid)[0][-1] + 1 + for k in range(npings): + s0 = int(start_samp_1b[k] - 1) + L = int(lens[k]) + sv_shift[k, :L] = sv_mat[k, s0 : s0 + L] + al_shift[k, :L] = al_mat[k, s0 : s0 + L] + ath_shift[k, :L] = ath_mat[k, s0 : s0 + L] - depth = depth[:last] - sv_mat = sv_mat[:, :last] - al_mat = al_mat[:, :last] - ath_mat = ath_mat[:, :last] - NbSamp = last + sv_mat, al_mat, ath_mat = sv_shift, al_shift, ath_shift return { "sv_mat": sv_mat, "al_mat": al_mat, "ath_mat": ath_mat, - "depth": depth, - "npings": npings, - "NbSamp": NbSamp, + "delta": delta, + "nbSamp": sv_mat.shape[1], + "depth1": depth1, + "depth2": depth2, + "start_samp_1b": start_samp_1b, + "start_depth_m": sdm_vec, } -def _compute_delta(depth: np.ndarray) -> float: - dd = np.diff(depth) - dd = dd[np.isfinite(dd)] - if dd.size == 0: - raise ValueError("Cannot compute delta from depth: all diffs are NaN.") - return float(np.nanmedian(dd)) +def _extract_metadata(ds: xr.Dataset, params: dict) -> dict: + channel = params["channel"] + + def _need(name: str) -> xr.DataArray: + if name not in ds: + raise ValueError(f"Missing required variable in ds: '{name}'") + return ds[name] + + def _sel_channel(da: xr.DataArray) -> xr.DataArray: + if "channel" in da.dims: + return da.sel(channel=channel) + return da + + def _scalar_from_da(da: xr.DataArray, name: str) -> float: + da = _sel_channel(da) + if "ping_time" in da.dims: + da = da.median("ping_time", skipna=True) + val = float(np.asarray(da.values).squeeze()) + if not np.isfinite(val): + raise ValueError(f"'{name}' is present but not finite for channel='{channel}'.") + return val + + def _ping_vector_from_da(da: xr.DataArray, name: str, npings: int) -> np.ndarray: + da = _sel_channel(da) + if "ping_time" not in da.dims: + val = float(np.asarray(da.values).squeeze()) + return np.full((npings,), val, dtype=float) + v = np.asarray(da.values, dtype=float).squeeze() + if v.shape[0] != npings: + raise ValueError(f"'{name}' ping_time length mismatch: {v.shape[0]} != {npings}") + return v + + if "channel" not in ds.coords: + raise ValueError("ds must have a 'channel' coordinate.") + if channel not in ds["channel"].values: + raise ValueError(f"Channel '{channel}' not found in ds['channel'].") + + if "ping_time" in ds.coords: + npings = int(ds.sizes["ping_time"]) + else: + raise ValueError("ds must have 'ping_time' coordinate/dimension.") + c = _scalar_from_da(ds["sound_speed"], "sound_speed") -def _extract_metadata( - ds: xr.Dataset, Param: dict, params: dict, channel: str, npings: int, DBG: dict -) -> Dict[str, Any]: - # c, pd - c = Param["sound_speed"] - if c is None: - if "sound_speed" in ds: - c = float(ds["sound_speed"].values) - else: - c = float(PARAM_DEFAULTS["SoundSpeed"]) - - pd = Param["pulse_length"] - if pd is None: - raise ValueError("params['pulse_length'] is required for CW (pd in MATLAB).") - - alpha = float(Param.get("alpha", 0.0)) - psi = float(Param.get("psi_two_way", 0.0)) - sa_cor = float(Param.get("Sa_correction", 0.0)) - SacorEK80Nominal = float(Param.get("Sa_EK80_nominal", 0.0)) - - # beam meta - ouv_al = ( - float(ds.get("beamwidth_alongship", xr.DataArray(np.nan)).sel(channel=channel).values) - if "beamwidth_alongship" in ds - else np.nan - ) - ouv_at = ( - float(ds.get("beamwidth_athwartship", xr.DataArray(np.nan)).sel(channel=channel).values) - if "beamwidth_athwartship" in ds - else np.nan - ) - dec_al = ( - float(ds.get("angle_offset_alongship", xr.DataArray(0.0)).sel(channel=channel).values) - if "angle_offset_alongship" in ds - else 0.0 - ) - dec_at = ( - float(ds.get("angle_offset_athwartship", xr.DataArray(0.0)).sel(channel=channel).values) - if "angle_offset_athwartship" in ds - else 0.0 - ) - al_sens = ( - float(ds.get("angle_sensitivity_alongship", xr.DataArray(1.0)).sel(channel=channel).values) - if "angle_sensitivity_alongship" in ds - else 1.0 - ) - ath_sens = ( - float( - ds.get("angle_sensitivity_athwartship", xr.DataArray(1.0)).sel(channel=channel).values - ) - if "angle_sensitivity_athwartship" in ds - else 1.0 - ) + if "transmit_duration_nominal" in ds: + pd = _scalar_from_da(ds["transmit_duration_nominal"], "transmit_duration_nominal") + else: + raise ValueError("Need ds['transmit_duration_nominal'] (seconds).") - missing_meta = [] - if np.isnan(ouv_al): - missing_meta.append("beamwidth_alongship") - if np.isnan(ouv_at): - missing_meta.append("beamwidth_athwartship") - if missing_meta: - warnings.warn( - f"Missing beam metadata for channel '{channel}': {', '.join(missing_meta)}. ", - UserWarning, - stacklevel=2, - ) + if "sample_interval" in ds: + dt_sec = _scalar_from_da(ds["sample_interval"], "sample_interval") + dt_usec = dt_sec * 1e6 + else: + raise ValueError("Need ds['sample_interval'].") - # TD/CH - TD = Param.get("transducer_depth", None) - if TD is None: - if "transducer_depth" in ds: - TD = float(ds["transducer_depth"].sel(channel=channel).values) - else: - TD = 0.0 - TD = float(TD) + alpha = _scalar_from_da(ds["sound_absorption"], "sound_absorption") + psi = _scalar_from_da(ds["equivalent_beam_angle"], "equivalent_beam_angle") + sa_cor = _scalar_from_da(ds["sa_correction"], "sa_correction") + + effpd = float(np.asarray(_sel_channel(ds["tau_effective"]).values).squeeze()) + + # EK80 nominal Sa term: compute-only (like Matecho) + sa_cor_ek80_nominal = 0.0 - CH = Param.get("heave_compensation", 0.0) - if np.ndim(CH) == 0: - CH_vec = np.full((npings,), float(CH), dtype=float) + if np.isfinite(effpd) and (effpd > 0.0) and (pd > 0.0): + sa_cor_ek80_nominal = float(np.float32(5.0 * np.log10(effpd / pd))) else: - CH_vec = np.asarray(CH, dtype=float) - if CH_vec.size != npings: - raise ValueError("heave_compensation length must match number of pings.") - - # bottom - bottom_da = Param.get("bottom_da", None) - if bottom_da is None: - bot = np.full((npings,), np.inf, dtype=float) + sa_cor_ek80_nominal = 0.0 + + ## + + ouv_al = _scalar_from_da(_need("beamwidth_alongship"), "beamwidth_alongship") + ouv_at = _scalar_from_da(_need("beamwidth_athwartship"), "beamwidth_athwartship") + + dec_al = _scalar_from_da(_need("angle_offset_alongship"), "angle_offset_alongship") + dec_at = _scalar_from_da(_need("angle_offset_athwartship"), "angle_offset_athwartship") + al_sens = _scalar_from_da(_need("angle_sensitivity_alongship"), "angle_sensitivity_alongship") + ath_sens = _scalar_from_da( + _need("angle_sensitivity_athwartship"), "angle_sensitivity_athwartship" + ) + + if "transducer_depth" in ds: + TD_vec = _ping_vector_from_da(ds["transducer_depth"], "transducer_depth", npings) else: - bot = bottom_da.sel(ping_time=ds["ping_time"]).values.astype(float) - - # timing - timeSampleInterval_usec = params.get("timeSampleInterval_usec", None) - if timeSampleInterval_usec is None: - raise ValueError("params['timeSampleInterval_usec'] is required.") - - # angle units detection - beamwidth_units = "" - if "beamwidth_alongship" in ds and hasattr(ds["beamwidth_alongship"], "attrs"): - beamwidth_units = str(ds["beamwidth_alongship"].attrs.get("units", "")).lower() - angles_are_degrees = "deg" in beamwidth_units - - if DBG["enable"]: - warnings.warn( - f"beamwidth_alongship units='{beamwidth_units}' -> \ - angles_are_degrees = {angles_are_degrees}", - category=UserWarning, - stacklevel=2, - ) + raise ValueError("Need ds['transducer_depth'].") + + if "heave_compensation" in ds: + CH_vec = _ping_vector_from_da(ds["heave_compensation"], "heave_compensation", npings) + else: + raise ValueError("Need ds['heave_compensation'].") + + # seafloor + sf = params["seafloor"] + + if isinstance(sf, xr.DataArray): + sf = sf.sel(ping_time=ds["ping_time"], method="nearest") + bot = np.asarray(sf.values, dtype=float).squeeze() + else: + bot = np.asarray(sf, dtype=float).squeeze() + + if bot.shape[0] != npings: + raise ValueError(f"seafloor length mismatch: {bot.shape[0]} != {npings}") return { "c": float(c), "pd": float(pd), + "sample_interval_usec": float(dt_usec), "alpha": float(alpha), "psi": float(psi), "sa_cor": float(sa_cor), - "SacorEK80Nominal": float(SacorEK80Nominal), + "sa_cor_ek80_nominal": float(sa_cor_ek80_nominal), + "tau_effective": float(effpd), "ouv_al": float(ouv_al), "ouv_at": float(ouv_at), "dec_al": float(dec_al), "dec_at": float(dec_at), "al_sens": float(al_sens), "ath_sens": float(ath_sens), - "TD": float(TD), - "CH_vec": CH_vec, + "TD_vec": TD_vec, + "heave_vec": CH_vec, "bot": bot, - "timeSampleInterval_usec": float(timeSampleInterval_usec), - "angles_are_degrees": bool(angles_are_degrees), } -def _normalise_angles(al_mat, ath_mat, meta: dict, DBG: dict): - if not meta["angles_are_degrees"]: - return al_mat, ath_mat, meta +def _compute_derived_for_parity(meta: dict, ex: dict, params: dict) -> dict: + c = float(meta["c"]) + pd = float(meta["pd"]) + dt_usec = float(meta["sample_interval_usec"]) + dt_s = dt_usec * 1e-6 - DEG2RAD = np.pi / 180.0 - if DBG["enable"]: - warnings.warn( - "Converting angles, beamwidths, and offsets from degrees to radians.", - category=UserWarning, - stacklevel=2, - ) + effpd = float(meta.get("tau_effective", np.nan)) + sa_cor_nominal = float(meta.get("sa_cor_ek80_nominal", np.nan)) + psi = float(meta["psi"]) + sa_cor = float(meta["sa_cor"]) - al_mat = al_mat * DEG2RAD - ath_mat = ath_mat * DEG2RAD - meta = dict(meta) - meta["dec_al"] *= DEG2RAD - meta["dec_at"] *= DEG2RAD - meta["ouv_al"] *= DEG2RAD - meta["ouv_at"] *= DEG2RAD - return al_mat, ath_mat, meta + ind_range0 = 1 if np.isfinite(effpd) else 3 + nech_p = pd / dt_s + sv2ts_raw64 = 10.0 * np.log10(c * pd / 2.0) + psi + 2.0 * sa_cor + 2.0 * sa_cor_nominal + sv2ts_f32 = float(np.float32(sv2ts_raw64)) -def _compute_pulse_derived(meta: dict, Param: dict, delta: float) -> Dict[str, float]: - c = meta["c"] - pd = meta["pd"] - pulse_len_m = c * pd / 2.0 + if "min_echo_length_pulse" not in params or params["min_echo_length_pulse"] is None: + raise ValueError("Missing required parameter: params['min_echo_length_pulse']") + if "max_echo_length_pulse" not in params or params["max_echo_length_pulse"] is None: + raise ValueError("Missing required parameter: params['max_echo_length_pulse']") + + min_echo_len = float(params["min_echo_length_pulse"]) + max_echo_len = float(params["max_echo_length_pulse"]) + + n_min = _matlab_round_pos(nech_p * min_echo_len) + n_max = _matlab_round_pos(nech_p * max_echo_len) - n_min = int(np.round((pulse_len_m * float(Param["MinEchoLength"])) / delta)) - n_max = int(np.round((pulse_len_m * float(Param["MaxEchoLength"])) / delta)) n_min = max(n_min, 1) - n_max = max(n_max, n_min) + n_max = max(n_max, 1) - min_space_m = float(Param["MinEchoSpace"]) * pulse_len_m return { - "pulse_len_m": float(pulse_len_m), + "timeSampleInterval": float(dt_s), + "tau_effective": float(effpd), + "ind_range0": int(ind_range0), + "nech_p": float(nech_p), + "dec_tir": int(params.get("dec_tir", 8)), + "depth1": float(ex.get("depth1", np.nan)), + "depth2": float(ex.get("depth2", np.nan)), + "delta": float(ex["delta"]), + "n_depth": int(ex["nbSamp"]), + "sv2ts_f32": float(sv2ts_f32), + "sv2ts_raw64": float(sv2ts_raw64), "n_min": int(n_min), "n_max": int(n_max), - "min_space_m": float(min_space_m), } -# --- -# Helpers: core matrices (ts_range, tsu, Plike, ts) -def _compute_core_matrices( - sv_mat: np.ndarray, - al_mat: np.ndarray, - ath_mat: np.ndarray, - depth: np.ndarray, - delta: float, - meta: dict, - Param: dict, - channel: str, -) -> Dict[str, Any]: - npings, NbSamp = sv_mat.shape - - ind_range0 = int(Param["tvg_start_sample"]) - dec_tir = 8 - - c = meta["c"] - pd = meta["pd"] - alpha = meta["alpha"] - psi = meta["psi"] - sa_cor = meta["sa_cor"] - SacorEK80Nominal = meta["SacorEK80Nominal"] - - # CW conversion constant - sv2ts = 10.0 * np.log10(c * pd / 2.0) + psi + 2.0 * sa_cor + 2.0 * SacorEK80Nominal - - # NechP (dt-based) - NechP = pd / (meta["timeSampleInterval_usec"] * 1e-6) - - # range vector in metres (one-way) - jj = np.arange(1, NbSamp + 1, dtype=float) - ts_range_1d = np.maximum(delta, (jj - ind_range0) * delta) - ts_range = np.tile(ts_range_1d[None, :], (npings, 1)) - - # bottom index constraint - bot = meta["bot"] - TD = meta["TD"] - CH_vec = meta["CH_vec"] - ir = (bot - TD - CH_vec - float(Param["FloorFilterTolerance"])) / delta - ir = np.minimum(ir, NbSamp).astype(float) - - # tsu / Plike / ts - tsu_mat = sv_mat + 20.0 * np.log10(ts_range) + sv2ts - Plike_mat = tsu_mat - 40.0 * np.log10(ts_range) - 2.0 * alpha * ts_range - - ouv_al = meta["ouv_al"] - ouv_at = meta["ouv_at"] - dec_al = meta["dec_al"] - dec_at = meta["dec_at"] - - x = 2.0 * (al_mat - dec_al) / ouv_al - y = 2.0 * (ath_mat - dec_at) / ouv_at - ts_mat = tsu_mat + 6.0206 * (x * x + y * y - 0.18 * (x * x) * (y * y)) - - return { - "npings": npings, - "NbSamp": NbSamp, - "ind_range0": ind_range0, - "dec_tir": dec_tir, - "NechP": float(NechP), - "sv2ts": float(sv2ts), - "ts_range": ts_range, - "ir": ir, - "tsu_mat": tsu_mat, - "Plike_mat": Plike_mat, - "ts_mat": ts_mat, - "depth": depth, - "channel": channel, - } +# conditions -# --- -# Helpers: Conditions 1–5 -def _cond_1to5(core: dict, Param: dict, DBG: dict) -> Tuple[np.ndarray, np.ndarray]: - NbSamp = core["NbSamp"] +def _cond_1to5(core: dict, params: dict): ts_mat = core["ts_mat"] tsu_mat = core["tsu_mat"] - Plike_mat = core["Plike_mat"] + plike_mat = core["plike_mat"] ir = core["ir"] + nbSamp = core["nbSamp"] dec_tir = core["dec_tir"] - cols = np.arange(1, NbSamp - 1, dtype=int) # never evaluate edge samples - - cond1 = ts_mat[:, 1:-1] > float(Param["MinThreshold"]) - cond2 = (ts_mat[:, 1:-1] - tsu_mat[:, 1:-1]) <= 2.0 * float(Param["MaxAngleOneWayCompression"]) - - left = Plike_mat[:, 0:-2] - mid = Plike_mat[:, 1:-1] - right = Plike_mat[:, 2:] - cond3 = mid >= np.nanmax(np.stack([left, right]), axis=0) - - col_1based = cols + 1 - cond4 = col_1based[None, :] < ir[:, None] - cond5 = col_1based[None, :] > dec_tir - - mask_1to5 = cond1 & cond2 & cond3 & cond4 & cond5 - i1, ind0 = np.where(mask_1to5) - ind = ind0 + 1 # back to full-matrix python col index - - if DBG["enable"] and DBG["print_cond_1to5"]: - warnings.warn( - "[Cond 1–5] Counts (pixel-level):\n" - f" c1(ts > thr): {int(np.count_nonzero(cond1))}\n" - f" c2(comp <=): {int(np.count_nonzero(cond2))}\n" - f" c3(loc max): {int(np.count_nonzero(cond3))}\n" - f" c4(above bottom): {int(np.count_nonzero(cond4))}\n" - f" c5(idx > dec_tir): {int(np.count_nonzero(cond5))}\n" - f" ALL: {int(np.count_nonzero(mask_1to5))}\n" - f" detections after Cond 1–5: {ind.size}", - category=UserWarning, - stacklevel=2, - ) - - return i1, ind - - -# --- -# Helpers: Condition 6 + # MATLAB: never evaluate 1st and last sample + cols0 = np.arange(1, nbSamp - 1) # Python 0-based + cols1 = cols0 + 1 # MATLAB 1-based + + ts = ts_mat[:, 1:-1] + tsu = tsu_mat[:, 1:-1] + pmid = plike_mat[:, 1:-1] + pleft = plike_mat[:, 0:-2] + pright = plike_mat[:, 2:] + + # Cond-1 : TS threshold + # ts_mat(:,2:end-1) > MinThreshold + cond1 = ts > float(params["min_threshold_db"]) + # Cond-2 : angle compression + # (ts − tsu) <= 2 * MaxAngleOneWayCompression + cond2 = (ts - tsu) <= 2.0 * float(params["max_angle_oneway_compression_db"]) + # Cond-3 : local max of Plike + cond3 = pmid >= np.maximum(pleft, pright) + # Cond-4 : above bottom gate + # col < ir(ping) + cond4 = cols1[None, :] < ir[:, None] + # Cond-5 : index > dec_tir + cond5 = cols1[None, :] > dec_tir + mask = cond1 & cond2 & cond3 & cond4 & cond5 + + # Pixel indices of valid detections + iping, icol0 = np.where(mask) + icol = icol0 + 1 + + return iping, icol + + +# --- Condition 6 def _cond_6_phase( i1: np.ndarray, ind: np.ndarray, @@ -553,453 +443,567 @@ def _cond_6_phase( ath_mat: np.ndarray, core: dict, meta: dict, - Param: dict, - DBG: dict, -) -> np.ndarray: - NbSamp = core["NbSamp"] - al_sens = meta["al_sens"] - ath_sens = meta["ath_sens"] - - std_ph_al = np.full((ind.size,), np.nan, dtype=float) - std_ph_ath = np.full((ind.size,), np.nan, dtype=float) - - for i in range(ind.size): - r0 = max(ind[i] - 2, 0) - r1 = min(ind[i] + 2, NbSamp - 1) - sl = slice(r0, r1 + 1) - std_ph_al[i] = np.nanstd(al_mat[i1[i], sl] * 128.0 / np.pi * al_sens) - std_ph_ath[i] = np.nanstd(ath_mat[i1[i], sl] * 128.0 / np.pi * ath_sens) + params: dict, +): + """ + Returns: + ind2: indices into (i1, ind) that pass Cond-6 + std_ph_al: per-candidate std (same length as i1/ind) + std_ph_ath: per-candidate std (same length as i1/ind) + """ + nbSamp = int(core["nbSamp"]) + max_phase_deviation_deg = float(params["max_phase_deviation_deg"]) + + # MATLAB scaling: * 128/pi * sens + al_sens = float(meta["al_sens"]) + ath_sens = float(meta["ath_sens"]) + scale_al = (128.0 / np.pi) * al_sens + scale_ath = (128.0 / np.pi) * ath_sens + + necho = ind.size + std_ph_al = np.full(necho, np.nan, dtype=float) + std_ph_ath = np.full(necho, np.nan, dtype=float) + + for k in range(necho): + p = int(i1[k]) + j = int(ind[k]) + + lo = j - 2 + hi = j + 2 + if hi >= nbSamp: + hi = nbSamp - 1 + if lo < 0: + lo = 0 + + a = al_mat[p, lo : hi + 1] * scale_al + b = ath_mat[p, lo : hi + 1] * scale_ath + + std_ph_al[k] = np.nanstd(a, ddof=1) + std_ph_ath[k] = np.nanstd(b, ddof=1) ind2 = np.where( - (std_ph_al <= float(Param["MaxPhaseDeviation"])) - & (std_ph_ath <= float(Param["MaxPhaseDeviation"])) + (std_ph_al <= max_phase_deviation_deg) & (std_ph_ath <= max_phase_deviation_deg) )[0] - if DBG["enable"] and DBG["print_cond6"]: - warnings.warn( - "[Cond 6] Phase deviation:\n" - f" std_ph_al: min={_nanmin(std_ph_al):.3f} " - f"med={_nanmed(std_ph_al):.3f} " - f"max={_nanmax(std_ph_al):.3f}\n" - f" std_ph_ath: min={_nanmin(std_ph_ath):.3f} " - f"med={_nanmed(std_ph_ath):.3f} " - f"max={_nanmax(std_ph_ath):.3f}\n" - f" Passing (<= MaxPhaseDeviation={float(Param['MaxPhaseDeviation']):.3f}): " - f"{ind2.size} / {ind.size}", - category=UserWarning, - stacklevel=2, - ) - - return ind2 + return ind2, std_ph_al, std_ph_ath -# --- -# Helpers: Condition 7 def _cond_7_echo_length( i1: np.ndarray, ind: np.ndarray, ind2: np.ndarray, core: dict, derived: dict, - Param: dict, - DBG: dict, -) -> Dict[str, Any]: - NbSamp = core["NbSamp"] - Plike_mat = core["Plike_mat"] + params: dict, +) -> dict: + """ + Matecho Cond-7 (CW): echo length gate based on Plike -6 dB bounds. + + Uses echopype port parameter names: + - params["min_echo_length_pulse"], params["max_echo_length_pulse"] + """ + + nbSamp = int(core["nbSamp"]) + plike_mat = core["plike_mat"] ir = core["ir"] - iid = ind[ind2] # python col index (0-based) - iip = i1[ind2] # ping indices - necho = iid.size + n_min = int(derived["n_min"]) + n_max = int(derived["n_max"]) + + # survivors after Cond-6 + iip = i1[ind2].astype(np.int64) + iid0 = ind[ind2].astype(np.int64) - isup = np.full((necho,), np.nan, dtype=float) - iinf = np.full((necho,), np.nan, dtype=float) + necho = iid0.size + iinf = np.empty(necho, dtype=np.int64) + isup = np.empty(necho, dtype=np.int64) - for i in range(necho): - ip = iip[i] - ii = iid[i] + for k in range(necho): + ip = int(iip[k]) + iid_1b = int(iid0[k] + 1) + # MATLAB: ir_floor = floor(ir(ip)) ir_floor = int(np.floor(ir[ip])) ir_floor = max(ir_floor, 1) - ir_floor = min(ir_floor, NbSamp) - - thr = Plike_mat[ip, ii] - 6.0 - - forward = Plike_mat[ip, ii:ir_floor] < thr - if np.any(forward): - umin = np.where(forward)[0].min() - isup[i] = ii + umin + 1.0 + ir_floor = min(ir_floor, nbSamp) + + peak = float(plike_mat[ip, iid_1b - 1]) + thr = peak - 6.0 + + # forward: iid .. ir_floor + fwd = plike_mat[ip, (iid_1b - 1) : ir_floor] < thr + if np.any(fwd): + umin0 = int(np.where(fwd)[0].min()) + # FIX: MATLAB isup = iid - 2 + min(u) + # Here umin0 = min(u) - 1 => isup = iid_1b - 1 + umin0 + isup[k] = iid_1b - 1 + umin0 else: - isup[i] = ii + 1.0 - - backward = Plike_mat[ip, 0 : ii + 1] < thr - if np.any(backward): - umax = np.where(backward)[0].max() - iinf[i] = umax + 2.0 + isup[k] = iid_1b + + # backward: 1 .. iid + bwd = plike_mat[ip, :iid_1b] < thr + if np.any(bwd): + umax0 = int(np.where(bwd)[0].max()) + # MATLAB: iinf = 1 + max(u) where max(u)=umax0+1 + iinf[k] = umax0 + 2 else: - iinf[i] = ii + 1.0 + iinf[k] = iid_1b - L6dB = isup - iinf + 1.0 - - n_min = derived["n_min"] - n_max = derived["n_max"] + L6dB = (isup - iinf + 1).astype(np.int64) ind4 = np.where((L6dB >= n_min) & (L6dB <= n_max))[0] - if DBG["enable"] and DBG["print_cond7"]: - warnings.warn( - "[Cond 7] Echo length (L6dB):\n" - f" L6dB: min={_nanmin(L6dB):.0f} " - f"med={_nanmed(L6dB):.0f} " - f"max={_nanmax(L6dB):.0f}\n" - f" Bounds depth-based [{n_min}..{n_max}] samples\n" - f" Passing depth-based: {ind4.size} / {necho}", - category=UserWarning, - stacklevel=2, - ) - - return {"iid": iid, "iip": iip, "iinf": iinf, "isup": isup, "L6dB": L6dB, "ind4": ind4} + return { + "iip": iip, + "iid0": iid0, + "iinf_1b": iinf, + "isup_1b": isup, + "L6dB": L6dB, + "ind4": ind4, + } -# --- -# Helpers: Condition 8 (refinement + spacing + output build) -def _refine_ir_fin_1b(Plike_row: np.ndarray, iinf_1b: int, isup_1b: int) -> float: - il6dB = np.arange(iinf_1b - 1, isup_1b, dtype=int) - w = 10.0 ** (Plike_row[il6dB] / 10.0) - return float(np.nansum((il6dB + 1) * w) / np.nansum(w)) +# cond 8 -def _cond_8_build_outputs( - ds: xr.Dataset, - params: dict, - i1: np.ndarray, - ind2: np.ndarray, +def _cond_8_refine_and_spacing_matecho( c7: dict, core: dict, meta: dict, derived: dict, - al_mat: np.ndarray, - ath_mat: np.ndarray, - depth: np.ndarray, - Param: dict, - DBG: dict, - delta: float, -) -> Tuple[dict, List[List[float]]]: - + params: dict, +) -> dict: + """ + Cond-8 (Matecho CW): refine detections, enforce within-ping spacing, and apply a depth gate. + + Steps (per ping) + - Refine range index `ir_fin` as a power-weighted centroid of samples within the -6 dB bounds + (`w = 10**(Plike/10)`). + - Correct TS/TSU to `ir_fin` using Matecho’s range/absorption adjustment. + - Keep strongest targets separated by at least `min_echo_space_pulse * nech_p` samples. + - Convert to depth `D = (ir_fin - ind_range0) * delta + TD + heave` and keep `D` within + [`min_echo_depth_m`, `max_echo_depth_m`]. + + Returns: kept indices, refined `ir_fin` (1-based), refined TS/TSU, + ranges/depths, and spacing in samples. + """ + ind4 = np.asarray(c7.get("ind4", []), dtype=int) + if ind4.size == 0: + return { + "keep_ind8": np.array([], dtype=int), + "ir_fin_1b": np.array([], dtype=float), + "ts_fin": np.array([], dtype=np.float32), + "tsu_fin": np.array([], dtype=np.float32), + "targetRange_m": np.array([], dtype=float), + "d_m": np.array([], dtype=float), + "min_space_samp": float(params.get("min_echo_space_pulse", 1.0)) + * float(derived["nech_p"]), + } + + # Cond-6 survivors (0-based) + Cond-7 bounds (1-based) + iip = np.asarray(c7["iip"], dtype=int) + iid0 = np.asarray(c7["iid0"], dtype=int) + iinf_1b = np.asarray(c7["iinf_1b"], dtype=int) + isup_1b = np.asarray(c7["isup_1b"], dtype=int) + + # Subset that passed Cond-7 + ip0 = iip[ind4] + ii0 = iid0[ind4] + + plike_mat = core["plike_mat"] ts_mat = core["ts_mat"] tsu_mat = core["tsu_mat"] - Plike_mat = core["Plike_mat"] - ind_range0 = core["ind_range0"] - Ping_number_full = np.arange(1, core["npings"] + 1, dtype=int) - ping_time = ds["ping_time"].values + nech_p = float(derived["nech_p"]) + ind_range0 = float(derived["ind_range0"]) + delta = float(derived["delta"]) + alpha = float(meta["alpha"]) - # candidates per ping after Cond7 - iid = c7["iid"] - iinf = c7["iinf"] - isup = c7["isup"] - ind4 = c7["ind4"] + min_space_samp = float(params.get("min_echo_space_pulse", 1.0)) * nech_p - ip0 = i1[ind2[ind4]] - ipu = np.unique(ip0) + TD_vec = np.asarray(meta["TD_vec"], dtype=float) + heave_vec = np.asarray(meta["heave_vec"], dtype=float) - out = { - k: [] - for k in [ - "TS_comp", - "TS_uncomp", - "Target_range", - "Target_range_disp", - "Target_range_min", - "Target_range_max", - "idx_r", - "StandDev_Angles_Minor_Axis", - "StandDev_Angles_Major_Axis", - "Angle_minor_axis", - "Angle_major_axis", - "Ping_number", - "Time", - "idx_target_lin", - "pulse_env_before_lin", - "pulse_env_after_lin", - "PulseLength_Normalized_PLDL", - "Transmitted_pulse_length", - "Heave", - "Roll", - "Pitch", - "Heading", - "Dist", - ] + min_d = float(params.get("min_echo_depth_m", 0.0)) + max_d = float(params.get("max_echo_depth_m", np.inf)) + + keep_global, ir_fin_keep, ts_fin_keep, tsu_fin_keep, tr_keep, D_keep = [], [], [], [], [], [] + + for ip in np.unique(ip0): + idx_p = np.where(ip0 == ip)[0] + if idx_p.size == 0: + continue + + ir_fin = np.empty(idx_p.size, dtype=np.float64) + ts_fin = np.empty(idx_p.size, dtype=np.float64) + tsu_fin = np.empty(idx_p.size, dtype=np.float64) + + for k_local, j in enumerate(idx_p): + j_surv = int(ind4[j]) + + lo_1b = int(iinf_1b[j_surv]) + hi_1b = int(isup_1b[j_surv]) + + il6dB_1b = np.arange(lo_1b, hi_1b + 1, dtype=np.float64) + il6dB_0b = il6dB_1b.astype(int) - 1 + + P = plike_mat[int(ip), il6dB_0b].astype(np.float64) + w = 10.0 ** (P / 10.0) + + den = np.nansum(w) + ir_fin[k_local] = (np.nansum(il6dB_1b * w) / den) if den != 0.0 else np.nan + + i_ini_1b = float(ii0[j] + 1) + base_ts = float(ts_mat[int(ip), int(ii0[j])]) + base_tsu = float(tsu_mat[int(ip), int(ii0[j])]) + + corr = 40.0 * np.log10( + (ir_fin[k_local] - ind_range0) / (i_ini_1b - ind_range0) + ) + 2.0 * alpha * delta * (ir_fin[k_local] - i_ini_1b) + ts_fin[k_local] = base_ts + corr + tsu_fin[k_local] = base_tsu + corr + + # Spacing: keep strongest TS_fin targets with min separation in ir_fin + if idx_p.size > 1: + order = np.argsort(ts_fin)[::-1] + kept = [int(order[0])] + for cand in order[1:]: + cand = int(cand) + sep = np.nanmin(np.abs(ir_fin[cand] - ir_fin[np.array(kept, dtype=int)])) + if sep >= min_space_samp: + kept.append(cand) + else: + kept = [0] + + ir_keep = ir_fin[np.array(kept, dtype=int)] + ts_keep = ts_fin[np.array(kept, dtype=int)] + tsu_keep = tsu_fin[np.array(kept, dtype=int)] + + targetRange_m = (ir_keep - ind_range0) * delta + d_m = targetRange_m + float(TD_vec[int(ip)]) + float(heave_vec[int(ip)]) + + ok = np.where((d_m >= min_d) & (d_m <= max_d))[0] + + kept_surv = ind4[idx_p[np.array(kept, dtype=int)]][ok] + + keep_global.extend(kept_surv.tolist()) + ir_fin_keep.extend(ir_keep[ok].tolist()) + ts_fin_keep.extend(ts_keep[ok].tolist()) + tsu_fin_keep.extend(tsu_keep[ok].tolist()) + tr_keep.extend(targetRange_m[ok].tolist()) + D_keep.extend(d_m[ok].tolist()) + + return { + "keep_ind8": np.asarray(keep_global, dtype=int), + "ir_fin_1b": np.asarray(ir_fin_keep, dtype=float), + "ts_fin": np.asarray(ts_fin_keep, dtype=np.float32), + "tsu_fin": np.asarray(tsu_fin_keep, dtype=np.float32), + "targetRange_m": np.asarray(tr_keep, dtype=float), + "d_m": np.asarray(D_keep, dtype=float), + "min_space_samp": float(min_space_samp), } - TSfreq_rows: List[List[float]] = [] - alpha = meta["alpha"] - TD = meta["TD"] - CH_vec = meta["CH_vec"] - bot = meta["bot"] - min_space_m = derived["min_space_m"] - for ip in ipu: - ip2 = np.where(ip0 == ip)[0] - sel = ind4[ip2] +# return steps - # compute fine range ir_fin - ir_fin = np.full((sel.size,), np.nan, dtype=float) - for j, idx in enumerate(sel): - ir_fin[j] = _refine_ir_fin_1b(Plike_mat[ip, :], int(iinf[idx]), int(isup[idx])) - i_ini = iid[sel] - i_ini_1b = i_ini + 1.0 +def _matecho_out_to_dataset(out: dict, *, channel: str) -> xr.Dataset: + n = int(out.get("nb_valid_targets", 0)) + target = np.arange(n, dtype=int) - ts_fin = ts_mat[ip, i_ini] + ( - 40.0 * np.log10((ir_fin - ind_range0) / (i_ini_1b - ind_range0)) - + 2.0 * alpha * delta * (ir_fin - i_ini_1b) - ) - tsu_fin = tsu_mat[ip, i_ini] + ( - 40.0 * np.log10((ir_fin - ind_range0) / (i_ini_1b - ind_range0)) - + 2.0 * alpha * delta * (ir_fin - i_ini_1b) + if n == 0: + empty_target = np.arange(0, dtype=int) + + return xr.Dataset( + data_vars=dict( + nb_valid_targets=((), 0), + ping_time=("target", np.asarray([], dtype="datetime64[ns]")), + range_sample=("target", np.asarray([], dtype=np.int64)), + frequency_nominal=("target", np.asarray([], dtype=np.float64)), + ), + coords=dict( + target=empty_target, + channel=channel, + ), + attrs=dict(method="matecho", channel=str(channel), nb_valid_targets=0), ) - # convert to metres once - range_fin_m = (ir_fin - ind_range0) * delta - - # spacing in metres - if sel.size > 1: - inds = np.argsort(-ts_fin) - keep = [inds[0]] - for k in inds[1:]: - if np.nanmin(np.abs(range_fin_m[k] - range_fin_m[np.array(keep)])) >= min_space_m: - keep.append(k) - ind5 = np.array(keep, dtype=int) - else: - ind5 = np.array([0], dtype=int) + def _arr(key, dtype=float): + a = np.asarray(out[key]) + if dtype is not None: + a = a.astype(dtype) + if a.size != n: + raise ValueError(f"Field {key!r} has length {a.size}, expected {n}") + return a - targetRange = range_fin_m[ind5] # metres (one-way range from transducer) - D = targetRange + TD + CH_vec[ip] # metres (depth from surface) + ds_out = xr.Dataset( + data_vars=dict( + nb_valid_targets=((), n), + # --- snake_case outputs --- + ts_comp=("target", _arr("ts_comp", float), {"units": "dB"}), + ts_uncomp=("target", _arr("ts_uncomp", float), {"units": "dB"}), + target_range=("target", _arr("target_range", float), {"units": "m"}), + target_range_disp=("target", _arr("target_range_disp", float), {"units": "m"}), + target_range_min=("target", _arr("target_range_min", float), {"units": "m"}), + target_range_max=("target", _arr("target_range_max", float), {"units": "m"}), + idx_r=("target", _arr("idx_r", int)), + idx_target_lin=("target", _arr("idx_target_lin", float)), + angle_minor_axis=("target", _arr("angle_minor_axis", float), {"units": "rad"}), + angle_major_axis=("target", _arr("angle_major_axis", float), {"units": "rad"}), + stddev_angles_minor_axis=("target", _arr("stddev_angles_minor_axis", float)), + stddev_angles_major_axis=("target", _arr("stddev_angles_major_axis", float)), + heave=("target", _arr("heave", float), {"units": "m"}), + roll=("target", _arr("roll", float), {"units": "rad"}), + pitch=("target", _arr("pitch", float), {"units": "rad"}), + heading=("target", _arr("heading", float), {"units": "rad"}), + dist=("target", _arr("dist", float), {"units": "m"}), + ping_time=("target", np.asarray(out["ping_time"], dtype="datetime64[ns]")), + range_sample=("target", _arr("range_sample", int)), + frequency_nominal=("target", _arr("frequency_nominal", float)), + ), + coords=dict( + target=target, + ping_number=("target", _arr("ping_number", int)), + channel=channel, + ), + attrs=dict(method="matecho", channel=str(channel), nb_valid_targets=n), + ) + return ds_out - ok = (D >= float(Param["MinEchoDepthM"])) & (D <= float(Param["MaxEchoDepthM"])) - if not np.any(ok): - continue - targetRange = targetRange[ok] - D = D[ok] - ind5_ok = ind5[ok] - - compensatedTS = ts_fin[ind5_ok] - unCompensatedTS = tsu_fin[ind5_ok] - - samp_idx = i_ini[ind5_ok] - AlongShipAngleRad = al_mat[ip, samp_idx] - AthwartShipAngleRad = ath_mat[ip, samp_idx] - - for pidx in range(targetRange.size): - positionWorldCoord = [np.nan, np.nan, np.nan] # to add - positionBeamCoord = [ - targetRange[pidx] * np.tan(AlongShipAngleRad[pidx]), - targetRange[pidx] * np.tan(AthwartShipAngleRad[pidx]), - targetRange[pidx], - ] - - DepthIndex = int(np.argmin(np.abs(depth - D[pidx])) + 1) - IndexPingClean = int(Ping_number_full[ip]) - DepthIndexFromTransducer = int(np.round(targetRange[pidx] / delta)) - - kf = int(params.get("Frequency_index", 1)) - - TSfreq_rows.append( - [ - float(Ping_number_full[ip]), - float(D[pidx]), - float(DepthIndex), - float(compensatedTS[pidx]), - float(unCompensatedTS[pidx]), - float(AlongShipAngleRad[pidx]), - float(AthwartShipAngleRad[pidx]), - 0.0, - float(positionWorldCoord[0]), - float(positionWorldCoord[1]), - float(positionWorldCoord[2]), - float(positionBeamCoord[0]), - float(positionBeamCoord[1]), - float(positionBeamCoord[2]), - float(Param["MinEchoDepthM"]), - float(Param["MaxEchoDepthM"]), - float(bot[ip]) if np.isfinite(bot[ip]) else np.nan, - float(IndexPingClean), - float(DepthIndexFromTransducer), - float(targetRange.size), - 1.0, - 0.0, - float(kf), - ] - ) - - out["TS_comp"].append(float(compensatedTS[pidx])) - out["TS_uncomp"].append(float(unCompensatedTS[pidx])) - - out["Target_range"].append(float(D[pidx])) - out["Target_range_disp"].append(float(D[pidx])) - out["Target_range_min"].append(float(Param["MinEchoDepthM"])) - out["Target_range_max"].append(float(Param["MaxEchoDepthM"])) - - out["idx_r"].append(int(samp_idx[pidx] + 1)) - out["idx_target_lin"].append(int(samp_idx[pidx] + 1)) - - out["pulse_env_before_lin"].append(np.nan) - out["pulse_env_after_lin"].append(np.nan) - out["PulseLength_Normalized_PLDL"].append(np.nan) - out["Transmitted_pulse_length"].append(float(meta["pd"])) - - out["Angle_minor_axis"].append(float(AlongShipAngleRad[pidx])) - out["Angle_major_axis"].append(float(AthwartShipAngleRad[pidx])) - - out["StandDev_Angles_Minor_Axis"].append(float(np.nan)) - out["StandDev_Angles_Major_Axis"].append(float(np.nan)) - - out["Heave"].append(float(CH_vec[ip])) - out["Roll"].append( - float(ds.get("roll", xr.DataArray(np.nan)).isel(ping_time=ip).values) - if "roll" in ds - else np.nan - ) - out["Pitch"].append( - float(ds.get("pitch", xr.DataArray(np.nan)).isel(ping_time=ip).values) - if "pitch" in ds - else np.nan - ) - out["Heading"].append( - float(ds.get("heading", xr.DataArray(np.nan)).isel(ping_time=ip).values) - if "heading" in ds - else np.nan - ) - out["Dist"].append(np.nan) - - out["Ping_number"].append(int(Ping_number_full[ip])) - out["Time"].append(ping_time[ip]) - - return out, TSfreq_rows - - -# --- -# main -def detect_matecho(ds: xr.Dataset, params: dict) -> xr.Dataset: - """ - Single-target detector (split-beam style) on CW Sv data. +def _build_outputs_from_refined( + ds: xr.Dataset, + params: dict, + ex: dict, + meta: dict, + c7: dict, + c8: dict, +) -> dict: + keep = np.asarray(c8["keep_ind8"], dtype=int) + n = keep.size + if n == 0: + return {"nb_valid_targets": 0} + + iip = np.asarray(c7["iip"], dtype=int) + iid0 = np.asarray(c7["iid0"], dtype=int) + + ping0 = iip[keep] + ii0 = iid0[keep] + ping_number = ping0 + 1 + + ts_comp = np.asarray(c8["ts_fin"], dtype=np.float32) + ts_uncomp = np.asarray(c8["tsu_fin"], dtype=np.float32) + + target_range_m = np.asarray(c8["d_m"], dtype=float) + idx_target_lin = np.asarray(c8["ir_fin_1b"], dtype=float) + ping_time = np.asarray(ds["ping_time"].values[ping0], dtype="datetime64[ns]") + + range_sample = np.asarray(np.rint(idx_target_lin - 1.0), dtype=int) + + # API field: frequency_nominal per target + if "frequency_nominal" in ds: + fn = ds["frequency_nominal"] + if "channel" in fn.dims: + fn_val = float(fn.sel(channel=params["channel"]).values) + else: + fn_val = float(np.asarray(fn.values).squeeze()) + frequency_nominal = np.full(n, fn_val, dtype=float) + else: + frequency_nominal = np.full(n, np.nan, dtype=float) - Overview - -------- - Detect candidate single targets ping-by-ping using a peak-based workflow: + angle_minor_axis = np.asarray(ex["al_mat"][ping0, ii0], dtype=float) + angle_major_axis = np.asarray(ex["ath_mat"][ping0, ii0], dtype=float) - 1) Build core matrices - - Convert Sv (dB) to TS-like quantities (uncompensated and compensated). - - Compute one-way range (m) and apply bottom/near-surface constraints. + heave_m = np.asarray(meta["heave_vec"][ping0], dtype=float) - 2) Phase I: candidate peak selection (Conditions 1–5) - - C1: compensated TS above `MinThreshold` - - C2: beam compensation not exceeding `MaxAngleOneWayCompression` - - C3: local maximum in power-like metric (Plike) - - C4: sample index above bottom (optionally provided) with `FloorFilterTolerance` - - C5: exclude early samples (near transducer) using a fixed guard index + out = { + "nb_valid_targets": int(n), + "ts_comp": ts_comp.astype(float).tolist(), + "ts_uncomp": ts_uncomp.astype(float).tolist(), + "target_range": target_range_m.tolist(), + "target_range_disp": target_range_m.tolist(), + "target_range_min": [float(params.get("min_echo_depth_m", np.nan))] * n, + "target_range_max": [float(params.get("max_echo_depth_m", np.nan))] * n, + "ping_number": ping_number.astype(int).tolist(), + "ping_time": ping_time, + "idx_r": (ii0 + 1).astype(int).tolist(), + "idx_target_lin": idx_target_lin.tolist(), + "angle_minor_axis": angle_minor_axis.tolist(), + "angle_major_axis": angle_major_axis.tolist(), + "stddev_angles_minor_axis": [np.nan] * n, + "stddev_angles_major_axis": [np.nan] * n, + "heave": heave_m.tolist(), + "roll": [np.nan] * n, + "pitch": [np.nan] * n, + "heading": [np.nan] * n, + "dist": [np.nan] * n, + "range_sample": range_sample.tolist(), + "frequency_nominal": frequency_nominal.tolist(), + "tsfreq_matrix": np.empty((0, 23), dtype=float), + } + return out - Note: edge samples are not evaluated (requires left/right neighbours). - 3) Condition 6: angle stability - - Reject peaks whose along/athwart angle standard deviation within ±2 samples - exceeds `MaxPhaseDeviation` (in phase-step units after sensitivity scaling). +# main - 4) Condition 7: echo-length - - Define pulse envelope at -6 dB relative to the peak in Plike. - - Keep candidates whose envelope length is within - [`MinEchoLength`, `MaxEchoLength`] expressed in pulse-length units. - 5) Phase II: refinement and deconfliction (Condition 8) - - Refine target position using a weighted centroid within the -6 dB envelope. - - Convert refined index to range (m) and enforce minimum spacing - `MinEchoSpace` (pulse-length units). - - Apply depth window [`MinEchoDepthM`, `MaxEchoDepthM`] using - `transducer_depth` and `heave_compensation`. (to check further) +def detect_matecho( + ds: xr.Dataset, + params: dict, +) -> xr.Dataset: + """ + Matecho-style CW split-beam single-target detector, ported to echopype. + + Implements the Matecho Cond-1…Cond-8 pipeline on Sv (dB) and split-beam + angles for a single channel. TSU/TS/Plike are computed internally (Sv to TS transform, + range/absorption terms, and beam compensation), then candidates are filtered by + threshold, angle compression, local Plike maxima, bottom gate, phase stability, + -6 dB echo length, and a final refinement + spacing + depth gate. + + Returns an `xr.Dataset` with `target` detections (TS_comp/TS_uncomp, refined index, + ping_time/number, range/depth, angles) and `nb_valid_targets`. Parameters ---------- ds : xr.Dataset - Dataset containing Sv (dB) and optional angle/metadata variables. + Dataset containing the following variables/coords: + • `Sv` in dB with dims (`channel`, `ping_time`, `range_sample`) + • `angle_alongship`, `angle_athwartship` with the same dims + • `echo_range` aligned with `range_sample` (and optionally `ping_time`) + • `start_depth_m` (per ping; used for Matecho SHIFT alignment) + • Scalars / vectors required by Matecho physics: + - `sound_speed` + - `transmit_duration_nominal` + - `sample_interval` + - `sound_absorption` + - `equivalent_beam_angle` + - `sa_correction` + - `tau_effective` (recommended; controls `ind_range0`) + - EK80 nominal Sa term (`Sa_EK80_nominal` or equivalent if present) + - split-beam geometry terms: + `beamwidth_alongship`, `beamwidth_athwartship`, + `angle_offset_alongship`, `angle_offset_athwartship`, + `angle_sensitivity_alongship`, `angle_sensitivity_athwartship` + - `transducer_depth` + - `heave_compensation` + + Notes + ----- + * This implementation assumes **split-beam** data (valid angle fields). + * Units must be internally consistent (e.g., absorption in dB/m). + params : dict - Detection parameters (see PARAM_DEFAULTS). Must include: - - channel, pulse_length, timeSampleInterval_usec - Optional: sound_speed, alpha, bottom_da, transducer_depth, - heave_compensation, etc. (add all) + Detection configuration. Required keys: + • `channel` : str + Channel identifier to process (must match `ds["channel"]` entries). + • `seafloor` : array-like or xr.DataArray + Bottom depth/range per ping_time (used by Cond-4 bottom gate). + • `min_threshold_db` : float + Cond-1 TS threshold. + • `max_angle_oneway_compression_db` : float + Cond-2 beam compression threshold (one-way, used as 2× in the test). + • `max_phase_deviation_deg` : float + Cond-6 phase stability threshold (degrees, after Matecho scaling). + • `min_echo_length_pulse` : float + Cond-7 minimum echo length in units of pulse length (× nech_p). + • `max_echo_length_pulse` : float + Cond-7 maximum echo length in units of pulse length (× nech_p). + + Optional keys: + • `dec_tir` : int, default 8 + Cond-5 minimum range index (1-based in Matecho; implemented accordingly). + • `floor_filter_tolerance_m` : float, default 0.0 + Extra margin in the bottom gate (subtracted before converting to samples). + • `min_echo_space_pulse` : float, default 1.0 + Cond-8 minimum within-ping separation, expressed in pulse lengths. + • `min_echo_depth_m` : float, default 0.0 + Minimum accepted displayed depth (after TD + CH). + • `max_echo_depth_m` : float, default +inf + Maximum accepted displayed depth (after TD + CH). + • `DBG` : dict, optional + Debug switches / payload (implementation-specific). Returns ------- xr.Dataset - Per-target outputs (TS_comp, TS_uncomp, angles, indices, time/ping number, etc.). - Empty dataset is returned when no targets are detected. - """ + Dataset with dimension `target` (length = number of valid detections), + including (at least): + • `TS_comp`, `TS_uncomp` (dB) + • `idx_r` (coarse peak sample index, 1-based) + • `idx_target_lin` (refined index, float, 1-based) + • `Ping_number` (1-based) + • `time` (datetime64) + • angle fields at the detected position + • `nb_valid_targets` scalar + and attrs: `method="matecho"`, `channel=`, `nb_valid_targets=`. - Param = _merge_params(params) - DBG = _make_dbg() + """ - channel, Sv = _validate_inputs(ds, Param) + channel, Sv_da = validate_matecho_inputs_strict(ds, params) + params = dict(params) + params["channel"] = channel - # extract variables - ex = _extract_matrices(ds, Sv, channel) - sv_mat = ex["sv_mat"] - al_mat = ex["al_mat"] - ath_mat = ex["ath_mat"] - depth = ex["depth"] - npings = ex["npings"] + ex = _extract_matrices(ds, Sv_da, channel) + meta = _extract_metadata(ds, params) - # extract Delta + metadata + angle normalisation - delta = _compute_delta(depth) - meta = _extract_metadata(ds, Param, params, channel, npings, DBG) - al_mat, ath_mat, meta = _normalise_angles(al_mat, ath_mat, meta, DBG) + # adaptations to match matecho algorithm - # compute derived pulse quantities - derived = _compute_pulse_derived(meta, Param, delta) + # convert angles degrees to radians + deg2rad = np.pi / 180.0 + ex["al_mat"] = ex["al_mat"] * deg2rad + ex["ath_mat"] = ex["ath_mat"] * deg2rad + meta["ouv_al"] *= deg2rad + meta["ouv_at"] *= deg2rad + # Echopype angles are already offset-corrected (Matecho not), so prevent double-subtraction + meta["dec_al"] = 0.0 + meta["dec_at"] = 0.0 + #### - # get core matrices (ts_range, tsu, Plike, ts) + constraints + derived = _compute_derived_for_parity(meta, ex, params) core = _compute_core_matrices( - sv_mat=sv_mat, - al_mat=al_mat, - ath_mat=ath_mat, - depth=depth, - delta=delta, + sv_mat=ex["sv_mat"], + al_mat=ex["al_mat"], + ath_mat=ex["ath_mat"], meta=meta, - Param=Param, - channel=channel, + params=params, + derived=derived, ) - # pass conditions - i1, ind = _cond_1to5(core, Param, DBG) - if ind.size == 0: - return _empty_out(channel) - - ind2 = _cond_6_phase(i1, ind, al_mat, ath_mat, core, meta, Param, DBG) - if ind2.size == 0: - return _empty_out(channel) - - c7 = _cond_7_echo_length(i1, ind, ind2, core, derived, Param, DBG) - if c7["ind4"].size == 0: - return _empty_out(channel) + # --- Cond 1–5: candidate pixels (joint logical mask) --- + # ip_15 / ir_15 are aligned 0-based indices into (ping, sample). + ip_15, ir_15 = _cond_1to5(core, params) + + # --- Cond 6: phase stability (±2 samples window) --- + # ind6 indexes into (ip_15, ir_15); survivors are ip_15[ind6], ir_15[ind6]. + ind6, std_ph_al, std_ph_ath = _cond_6_phase( + i1=ip_15, + ind=ir_15, + al_mat=ex["al_mat"], + ath_mat=ex["ath_mat"], + core=core, + meta=meta, + params=params, + ) - out, TSfreq_rows = _cond_8_build_outputs( - ds=ds, + # --- Cond 7: echo length (-6 dB width in Plike) --- + # c7 bundles Cond-6 survivors + 6 dB bounds + pass subset (ind4). + c7 = _cond_7_echo_length( + i1=ip_15, + ind=ir_15, + ind2=ind6, + core=core, + derived=derived, params=params, - i1=i1, - ind2=ind2, + ) + + # --- Cond 8: refine (centroid), spacing, depth gate --- + # Operates on c7["ind4"]; returns final kept targets. + c8 = _cond_8_refine_and_spacing_matecho( c7=c7, core=core, meta=meta, derived=derived, - al_mat=al_mat, - ath_mat=ath_mat, - depth=depth, - Param=Param, - DBG=DBG, - delta=delta, + params=params, ) - out["nb_valid_targets"] = len(out["TS_comp"]) - out["TSfreq_matrix"] = ( - np.asarray(TSfreq_rows, dtype=float) if TSfreq_rows else np.empty((0, 23)) - ) - return _matecho_struct_to_dataset(out, channel=channel) + out = _build_outputs_from_refined(ds=ds, params=params, ex=ex, meta=meta, c7=c7, c8=c8) + return _matecho_out_to_dataset(out, channel=channel) diff --git a/echopype/tests/mask/test_mask.py b/echopype/tests/mask/test_mask.py index a57c42b10..6ad95c99d 100644 --- a/echopype/tests/mask/test_mask.py +++ b/echopype/tests/mask/test_mask.py @@ -2013,177 +2013,417 @@ def test_echoview_mincan_no_linking(): # test for single target detections -REQ = { - "channel": "chan1", - "pulse_length": 1e-3, # 1 ms dummy - "timeSampleInterval_usec": 40.0, # dummy -} +# Helpers: base coords +def _coords_ping_range(n_ping=5, n_range=10): + return { + "ping_time": np.arange(n_ping), + "range_sample": np.arange(n_range), + } -def _make_ds_single_target_stub( - n_ping=5, n_range=10, channels=("chan1",), var_name="Sv" +# Stub: Matecho-style ds (3D with channel) +def _make_ds_matecho_minimal( + n_ping=5, + n_range=10, + channels=("chan1",), ) -> xr.Dataset: - coords = {"ping_time": np.arange(n_ping), "range_sample": np.arange(n_range)} - da = xr.DataArray( - np.full((len(channels), n_ping, n_range), -90.0, dtype=float), + coords = _coords_ping_range(n_ping, n_range) + ping_time = coords["ping_time"] + range_sample = coords["range_sample"] + channel = np.array(list(channels)) + + # Core 3D fields + Sv = xr.DataArray( + np.full((len(channel), n_ping, n_range), -90.0, dtype=float), dims=("channel", "ping_time", "range_sample"), - coords={**coords, "channel": list(channels)}, - name=var_name, + coords={"channel": channel, "ping_time": ping_time, "range_sample": range_sample}, + name="Sv", ) - return da.to_dataset() -def _make_ds_single_target_wrong_dims( - n_ping=5, n_range=10, channels=("chan1",), var_name="Sv" -) -> xr.Dataset: - # Missing 'range_sample' on purpose - coords = {"ping_time": np.arange(n_ping)} + al = xr.DataArray( + np.zeros((len(channel), n_ping, n_range), dtype=float), + dims=("channel", "ping_time", "range_sample"), + coords={"channel": channel, "ping_time": ping_time, "range_sample": range_sample}, + name="angle_alongship", + ) + + ath = xr.DataArray( + np.zeros((len(channel), n_ping, n_range), dtype=float), + dims=("channel", "ping_time", "range_sample"), + coords={"channel": channel, "ping_time": ping_time, "range_sample": range_sample}, + name="angle_athwartship", + ) + + # Range axis (2D accepted) + echo_range = xr.DataArray( + np.tile(np.linspace(1.0, float(n_range), n_range)[None, :], (n_ping, 1)), + dims=("ping_time", "range_sample"), + coords={"ping_time": ping_time, "range_sample": range_sample}, + name="echo_range", + ) + + # SHIFT alignment requirement + start_depth_m = xr.DataArray( + np.zeros((n_ping, len(channel)), dtype=float), + dims=("ping_time", "channel"), + coords={"ping_time": ping_time, "channel": channel}, + name="start_depth_m", + ) + + # scalar-ish + sound_speed = xr.DataArray(1500.0, name="sound_speed") + equivalent_beam_angle = xr.DataArray(-20.0, name="equivalent_beam_angle") + sa_correction = xr.DataArray(0.0, name="sa_correction") + sample_interval = xr.DataArray(4e-5, name="sample_interval") # 40 us + + # tau_effective: channel-only + tau_effective = xr.DataArray( + np.full((len(channel),), 8e-4, dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="tau_effective", + ) + + transmit_duration_nominal = xr.DataArray( + np.full((n_ping, len(channel)), 1e-3, dtype=float), # 1 ms + dims=("ping_time", "channel"), + coords={"ping_time": ping_time, "channel": channel}, + name="transmit_duration_nominal", + ) + + # ping_time x channel required + sound_absorption = xr.DataArray( + np.full((n_ping, len(channel)), 0.003, dtype=float), + dims=("ping_time", "channel"), + coords={"ping_time": ping_time, "channel": channel}, + name="sound_absorption", + ) + transducer_depth = xr.DataArray( + np.full((n_ping, len(channel)), 2.0, dtype=float), + dims=("ping_time", "channel"), + coords={"ping_time": ping_time, "channel": channel}, + name="transducer_depth", + ) + heave_compensation = xr.DataArray( + np.zeros((n_ping, len(channel)), dtype=float), + dims=("ping_time", "channel"), + coords={"ping_time": ping_time, "channel": channel}, + name="heave_compensation", + ) + + # channel-only beam geometry + beamwidth_alongship = xr.DataArray( + np.full((len(channel),), 7.0, dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="beamwidth_alongship", + ) + beamwidth_athwartship = xr.DataArray( + np.full((len(channel),), 7.0, dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="beamwidth_athwartship", + ) + angle_offset_alongship = xr.DataArray( + np.zeros((len(channel),), dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="angle_offset_alongship", + ) + angle_offset_athwartship = xr.DataArray( + np.zeros((len(channel),), dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="angle_offset_athwartship", + ) + angle_sensitivity_alongship = xr.DataArray( + np.ones((len(channel),), dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="angle_sensitivity_alongship", + ) + angle_sensitivity_athwartship = xr.DataArray( + np.ones((len(channel),), dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="angle_sensitivity_athwartship", + ) + + beam_type = xr.DataArray( + np.full((len(channel),), 1, dtype=np.int16), + dims=("channel",), + coords={"channel": channel}, + name="beam_type", + ) + + frequency_nominal = xr.DataArray( + np.full((len(channel),), 38000.0, dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="frequency_nominal", + ) + + ds = xr.Dataset( + data_vars=dict( + Sv=Sv, + angle_alongship=al, + angle_athwartship=ath, + echo_range=echo_range, + start_depth_m=start_depth_m, + sound_speed=sound_speed, + transmit_duration_nominal=transmit_duration_nominal, + tau_effective=tau_effective, + sample_interval=sample_interval, + sound_absorption=sound_absorption, + equivalent_beam_angle=equivalent_beam_angle, + sa_correction=sa_correction, + beamwidth_alongship=beamwidth_alongship, + beamwidth_athwartship=beamwidth_athwartship, + angle_offset_alongship=angle_offset_alongship, + angle_offset_athwartship=angle_offset_athwartship, + angle_sensitivity_alongship=angle_sensitivity_alongship, + angle_sensitivity_athwartship=angle_sensitivity_athwartship, + transducer_depth=transducer_depth, + heave_compensation=heave_compensation, + beam_type=beam_type, + frequency_nominal=frequency_nominal, + ), + coords=dict(channel=channel, ping_time=ping_time, range_sample=range_sample), + ) + return ds + + +def _make_ds_matecho_wrong_dims_missing_range(channels=("chan1",), n_ping=5) -> xr.Dataset: + # Sv missing range_sample dim intentionally (and also missing range_sample coord/dim) da = xr.DataArray( np.full((len(channels), n_ping), -90.0, dtype=float), dims=("channel", "ping_time"), - coords={**coords, "channel": list(channels)}, - name=var_name, + coords={"channel": list(channels), "ping_time": np.arange(n_ping)}, + name="Sv", ) - return da.to_dataset() + ds = da.to_dataset() + ds = ds.assign( + beam_type=xr.DataArray( + np.ones((len(channels),), dtype=np.int16), + dims=("channel",), + coords={"channel": list(channels)}, + ) + ) + return ds + + +# Stub: EV split method2 ds_m2 (2D, no channel dim) +def _make_ds_ev_m2_minimal(n_ping=5, n_range=10) -> xr.Dataset: + coords = _coords_ping_range(n_ping, n_range) + + TS = xr.DataArray( + np.full((n_ping, n_range), -90.0, dtype=float), + dims=("ping_time", "range_sample"), + coords=coords, + name="TS", + ) + echo_range = xr.DataArray( + np.tile(np.linspace(1.0, float(n_range), n_range)[None, :], (n_ping, 1)), + dims=("ping_time", "range_sample"), + coords=coords, + name="echo_range", + ) + angle_al = xr.DataArray( + np.zeros((n_ping, n_range), dtype=float), + dims=("ping_time", "range_sample"), + coords=coords, + name="angle_alongship", + ) + angle_ath = xr.DataArray( + np.zeros((n_ping, n_range), dtype=float), + dims=("ping_time", "range_sample"), + coords=coords, + name="angle_athwartship", + ) + + sound_absorption = xr.DataArray(0.003, name="sound_absorption") + timeSampleInterval_usec = xr.DataArray(40.0, name="timeSampleInterval_usec") + effective_pulse_duration = xr.DataArray(8e-4, name="effective_pulse_duration") + + frequency_nominal = xr.DataArray(38000.0, name="frequency_nominal") + + # beam compensation inputs for simrad_lobe + beamwidth_al = xr.DataArray(7.0, name="beamwidth_alongship") + beamwidth_at = xr.DataArray(7.0, name="beamwidth_athwartship") + angle_off_al = xr.DataArray(0.0, name="angle_offset_alongship") + angle_off_at = xr.DataArray(0.0, name="angle_offset_athwartship") + + beam_type = xr.DataArray(np.int16(1), name="beam_type") + + ds = xr.Dataset( + data_vars=dict( + TS=TS, + echo_range=echo_range, + angle_alongship=angle_al, + angle_athwartship=angle_ath, + sound_absorption=sound_absorption, + timeSampleInterval_usec=timeSampleInterval_usec, + effective_pulse_duration=effective_pulse_duration, + frequency_nominal=frequency_nominal, + beamwidth_alongship=beamwidth_al, + beamwidth_athwartship=beamwidth_at, + angle_offset_alongship=angle_off_al, + angle_offset_athwartship=angle_off_at, + beam_type=beam_type, + ), + coords=coords, + ) + return ds + +# Tests: dispatcher generic @pytest.mark.unit -def test_detect_single_targets_var_name_missing_in_ds_raises(): - ds = _make_ds_single_target_stub(var_name="Sv") # dataset has Sv - with pytest.raises(ValueError, match=r"var_name 'Sv_corrected' not found"): +def test_detect_single_targets_unknown_method_raises(): + ds = _make_ds_matecho_minimal() + with pytest.raises(ValueError, match="Unsupported single-target method"): detect_single_targets( ds, - method="matecho", - params={**REQ, "var_name": "Sv_corrected"}, + method="__bad__", + params={"channel": "chan1", "seafloor": np.zeros(ds.sizes["ping_time"])}, ) @pytest.mark.unit -def test_detect_single_targets_default_var_name_Sv_ok(): - ds = _make_ds_single_target_stub(var_name="Sv") - out = detect_single_targets( - ds, - method="matecho", - params=REQ, # no var_name -> default Sv - ) - assert isinstance(out, xr.Dataset) - assert "target" in out.dims - assert out.sizes["target"] == 0 +def test_detect_single_targets_params_required_raises(): + ds = _make_ds_matecho_minimal() + with pytest.raises(ValueError, match="No parameters given"): + detect_single_targets(ds, method="matecho", params=None) @pytest.mark.unit -def test_detect_single_targets_var_name_wrong_dims_raises(): - ds = _make_ds_single_target_wrong_dims(var_name="Sv") - with pytest.raises(ValueError, match=r"Sv must have dims"): +def test_detect_single_targets_beam_type_missing_raises(): + ds = _make_ds_matecho_minimal().drop_vars("beam_type") + with pytest.raises(ValueError, match="beam_type variable is missing"): detect_single_targets( ds, method="matecho", - params={**REQ, "var_name": "Sv"}, + params={"channel": "chan1", "seafloor": np.zeros(ds.sizes["ping_time"])}, ) @pytest.mark.unit -def test_detect_single_targets_channel_not_in_ds_raises(): - ds = _make_ds_single_target_stub(channels=("chan1",)) - with pytest.raises(ValueError, match=r"Channel 'chanX' not found"): +def test_detect_single_targets_beam_type_not_split_raises(): + ds = _make_ds_matecho_minimal() + ds["beam_type"] = xr.DataArray( + np.full(ds.sizes["channel"], 2, dtype=np.int16), + dims=("channel",), + coords={"channel": ds["channel"]}, + ) + with pytest.raises(ValueError, match="Only split-beam data supported"): detect_single_targets( ds, method="matecho", - params={**REQ, "channel": "chanX", "var_name": "Sv"}, + params={"channel": "chan1", "seafloor": np.zeros(ds.sizes["ping_time"])}, ) +# Tests: matecho-specific validation + schema +MATECHO_REQ = { + "channel": "chan1", + "seafloor": None, + "min_threshold_db": -58.0, + "max_angle_oneway_compression_db": 2.0, + "max_phase_deviation_deg": 8.0, + "min_echo_length_pulse": 0.7, + "max_echo_length_pulse": 1.5, + # optional: + "min_echo_space_pulse": 1.0, + "min_echo_depth_m": 0.0, + "max_echo_depth_m": 1e9, + "floor_filter_tolerance_m": 0.0, +} + + @pytest.mark.unit -def test_detect_single_targets_unknown_method_raises(): - ds = _make_ds_single_target_stub() - with pytest.raises(ValueError, match="Unsupported single-target method"): - detect_single_targets(ds, method="__bad__", params={**REQ}) +def test_detect_single_targets_matecho_missing_channel_param_raises(): + ds = _make_ds_matecho_minimal() + p = dict(MATECHO_REQ) + p.pop("channel") + p["seafloor"] = np.zeros(ds.sizes["ping_time"]) + with pytest.raises(ValueError, match=r"params\['channel'\]"): + detect_single_targets(ds, method="matecho", params=p) @pytest.mark.unit -def test_detect_single_targets_params_required_raises(): - ds = _make_ds_single_target_stub() - with pytest.raises(ValueError, match="No parameters given"): - detect_single_targets(ds, method="matecho", params=None) +def test_detect_single_targets_matecho_channel_not_in_ds_raises(): + ds = _make_ds_matecho_minimal(channels=("chan1",)) + p = dict(MATECHO_REQ) + p["channel"] = "chanX" + p["seafloor"] = np.zeros(ds.sizes["ping_time"]) + with pytest.raises(ValueError, match=r"Channel 'chanX' not found"): + detect_single_targets(ds, method="matecho", params=p) @pytest.mark.unit -def test_detect_single_targets_missing_channel_raises(): - ds = _make_ds_single_target_stub() - with pytest.raises(ValueError, match=r"params\['channel'\] is required"): - p = dict(REQ) - p.pop("channel", None) +def test_detect_single_targets_matecho_sv_wrong_dims_raises(): + ds = _make_ds_matecho_wrong_dims_missing_range() + p = dict(MATECHO_REQ) + p["seafloor"] = np.zeros(ds.sizes["ping_time"]) + + with pytest.raises(ValueError, match=r"Missing required coordinate/dimension: 'range_sample'"): detect_single_targets(ds, method="matecho", params=p) @pytest.mark.unit -def test_detect_single_targets_returns_dataset_with_target_dim_zero(): - ds = _make_ds_single_target_stub(channels=("chan1",)) - out = detect_single_targets(ds, method="matecho", params=REQ) +def test_detect_single_targets_matecho_returns_empty_with_required_fields(): + ds = _make_ds_matecho_minimal() + p = dict(MATECHO_REQ) + p["seafloor"] = np.zeros(ds.sizes["ping_time"]) + + out = detect_single_targets(ds, method="matecho", params=p) assert isinstance(out, xr.Dataset) assert "target" in out.dims assert out.sizes["target"] == 0 - assert "Time" in out.coords - assert "Ping_number" in out.coords - assert out["Time"].dims == ("target",) - assert out["Ping_number"].dims == ("target",) - - assert "channel" in out.coords - assert out["channel"].ndim == 0 # scalar coordinate - - -@pytest.mark.unit -def test_detect_single_targets_schema_variables_present_even_if_empty(): - ds = _make_ds_single_target_stub() - out = detect_single_targets(ds, method="matecho", params=REQ) - - expected_vars = [ - "TS_comp", - "TS_uncomp", - "Target_range", - "Target_range_disp", - "Target_range_min", - "Target_range_max", - "idx_r", - "idx_target_lin", - "pulse_env_before_lin", - "pulse_env_after_lin", - "PulseLength_Normalized_PLDL", - "Transmitted_pulse_length", - "Angle_minor_axis", - "Angle_major_axis", - "StandDev_Angles_Minor_Axis", - "StandDev_Angles_Major_Axis", - "Heave", - "Roll", - "Pitch", - "Heading", - "Dist", - "TSfreq_matrix", - ] - for v in expected_vars: - assert v in out.data_vars - # TSfreq_matrix is 2D even if empty; the rest are 1D on target - if v == "TSfreq_matrix": - assert out[v].dims[0] == "target" - assert out[v].sizes["target"] == 0 - else: - assert out[v].dims == ("target",) - assert out[v].sizes["target"] == 0 + for v in ("ping_time", "range_sample", "frequency_nominal"): + assert v in out + assert out[v].dims == ("target",) + assert out[v].sizes["target"] == 0 @pytest.mark.unit def test_detect_single_targets_concat_empty_is_stable(): - ds = _make_ds_single_target_stub() - out1 = detect_single_targets(ds, method="matecho", params=REQ) - out2 = detect_single_targets(ds, method="matecho", params=REQ) + ds = _make_ds_matecho_minimal() + p = dict(MATECHO_REQ) + p["seafloor"] = np.zeros(ds.sizes["ping_time"]) + + out1 = detect_single_targets(ds, method="matecho", params=p) + out2 = detect_single_targets(ds, method="matecho", params=p) outc = xr.concat([out1, out2], dim="target") assert "target" in outc.dims assert outc.sizes["target"] == 0 +# Tests: EV split method2 +EV_REQ = { + "ts_threshold_db": -58.0, + "pldl_db": 6.0, + "min_norm_pulse": 0.7, + "max_norm_pulse": 1.5, + "beam_comp_model": "simrad_lobe", + "max_beam_comp_db": 4.0, + "max_sd_minor_deg": 0.6, + "max_sd_major_deg": 0.6, +} + + +@pytest.mark.unit +def test_detect_single_targets_echoview_m2_missing_required_param_raises(): + ds = _make_ds_ev_m2_minimal() + p = dict(EV_REQ) + p.pop("pldl_db") + with pytest.raises(ValueError, match=r"Missing required parameters"): + detect_single_targets(ds, method="echoview_split_method2", params=p) + + @pytest.mark.unit -def test_detect_single_targets_missing_optional_metadata_warns(): - ds = _make_ds_single_target_stub() - with pytest.warns(UserWarning, match=r"Missing beam metadata"): - detect_single_targets(ds, method="matecho", params=REQ) \ No newline at end of file +def test_detect_single_targets_echoview_m2_missing_required_ds_var_raises(): + ds = _make_ds_ev_m2_minimal().drop_vars("effective_pulse_duration") + with pytest.raises(ValueError, match=r"ds_m2 missing required variable"): + detect_single_targets(ds, method="echoview_split_method2", params=dict(EV_REQ)) \ No newline at end of file