From 9581730e4161296759f765f64066b0414a1c13b1 Mon Sep 17 00:00:00 2001 From: Ayzen Date: Thu, 19 Mar 2026 19:28:30 +0300 Subject: [PATCH] added GPR --- GPR_clean_v1.py | 653 ++++++++++ Makefile | 1 + .../common_cpp/config/include/run_config.hpp | 21 + .../common_cpp/config/src/run_config.cpp | 71 + .../common_cpp/ipc/include/shared_types.hpp | 10 + .../common_cpp/ipc/src/shared_types.cpp | 91 ++ .../data_processor/include/data_processor.hpp | 3 + .../include/processing_live_config.hpp | 13 + .../data_processor/src/data_processor.cpp | 60 +- .../src/processing_live_config.cpp | 84 ++ .../processors/include/bscan_processor.hpp | 8 +- .../processors/include/gpr_processor.hpp | 18 + .../include/passthrough_processor.hpp | 8 +- .../include/processor_interface.hpp | 10 +- .../processors/src/bscan_processor.cpp | 44 +- .../processors/src/gpr_processor.cpp | 1149 +++++++++++++++++ .../processors/src/passthrough_processor.cpp | 50 +- python_app/gui/app_window.py | 12 + .../controllers/app_window_config_mixin.py | 190 ++- .../controllers/app_window_pipeline_mixin.py | 8 +- .../gui/controllers/app_window_plot_mixin.py | 310 ++++- .../controllers/app_window_snapshot_mixin.py | 2 +- .../gui/controllers/app_window_ui_mixin.py | 18 +- .../gui/controllers/sections/__init__.py | 2 + .../sections/gpr_config_section.py | 53 + .../sections/processing_section.py | 131 +- python_app/gui/runtime/constraints.py | 48 +- python_app/models/dataset_model.py | 24 +- python_app/models/run_config_codec.py | 60 +- python_app/models/run_config_model.py | 6 + python_app/models/run_config_schema.py | 27 + python_app/models/run_config_validation.py | 33 +- .../orchestration/live_processing_config.py | 37 +- python_app/orchestration/shm/decoder.py | 99 +- python_app/runtime/run_config_smoke.json | 64 +- ...onvert_prog_libre_manual_to_vna_history.py | 384 ++++++ python_app/storage/npz/serialize.py | 90 +- python_app/storage/npz/snapshot_numpy.py | 93 ++ run_config.json | 46 +- 39 files changed, 3830 insertions(+), 201 deletions(-) create mode 100644 GPR_clean_v1.py create mode 100644 data_acq_and_processing/processing/processors/include/gpr_processor.hpp create mode 100644 data_acq_and_processing/processing/processors/src/gpr_processor.cpp create mode 100644 python_app/gui/controllers/sections/gpr_config_section.py create mode 100644 python_app/scripts/convert_prog_libre_manual_to_vna_history.py diff --git a/GPR_clean_v1.py b/GPR_clean_v1.py new file mode 100644 index 0000000..877c1d3 --- /dev/null +++ b/GPR_clean_v1.py @@ -0,0 +1,653 @@ +""" +MIMO GPR — локализация через пересечение эллипсов +================================================== + +Физика в двух словах: + Пик A-скана пары (Tx_i, Rx_j) на задержке τ означает: + |Tx → объект| + |объект → Rx| = v · τ + Это уравнение эллипса. Истинный отражатель лежит на + пересечении всех 16 эллипсов (по одному на пару). + +Алгоритм: + 1. S(f) → IFFT → 16 A-сканов + 2. Поиск пиков: SNR = пик / медиана > порог + 3. Для каждого пика → мягкий эллипс в аккумуляторе + (с компенсацией геометрического и углового затухания) + 4. CLEAN: найти максимум → убрать его эллипсы → повторить + +О параметре SHELL_SIGMA: + Аккумулятор — это «мягкое голосование». Каждый эллипс добавляет + не единицу, а гауссово-взвешенный вклад: + w = exp(−δ²/2σ²), где δ = |R_Tx + R_Rx − v·τ| + SHELL_SIGMA — ширина этой гауссовой оболочки. + Слишком широко → ghost-цели не подавляются. + Слишком узко → вклад падает до нуля из-за дискретности сетки. + Оптимум: ~ 0.4 × δZ, где δZ = v/(2B) — разрешение по глубине. + +О score: + score = количество пар (из 16), чей эллипс проходит + через данную точку с невязкой δ < 3σ. + Принимает целые значения от 0 до 16. + Максимальный score у истинного объекта = 16 (все пары согласны). + Ghost-цели имеют меньший score, т.к. согласуются только + с частью пар. + +О компенсации затухания: + При генерации S(f) сигнал ослаблен: + geo(i,j) = 1/(R_Tx · R_Rx) — геометрическое ослабление + pat(i,j) = cos²(θ_Tx)·cos²(θ_Rx) — диаграмма направленности + Без компенсации глубокий/угловой отражатель будет недооценён. + Компенсация: делим вес каждого пика на ожидаемое затухание + в точке z_apparent, вычисленное для данной пары антенн. + +Геометрия: плоскость XZ (X — вдоль антенн, Z — глубина). +""" + +""" +MIMO GPR — локализация через пересечение эллипсов +================================================== +Версия для реальных данных +""" + +import numpy as np +import matplotlib.pyplot as plt +from scipy.signal import find_peaks +from scipy.ndimage import gaussian_filter, label +from pathlib import Path + +### Изменяемые параметры +INPUT_IDX = [0,1,2,3] +OUTPUT_IDX = [0,3] + +MIN_DEPTH = 2.0 # [м] пропустить прямую волну +MAX_DEPTH = 14.0 +COMP_POWER = 0.2 # степень компенсации затухания + +F_START = 30*1e8 # Нижняя частота +F_STOP = 6*1e9 # Верхняя частота + +# Вычитание среднего фона (background removal) + +BG_SUBTRACT = True #True/False +BG_PATH = Path('/Users/ivan_root/Downloads/Telegram_dwnld/03-13-measure/2_cylinder_dif_side_and_mushrooms/preprocessed') + +DATA_PATH = Path('/Users/ivan_root/Downloads/Telegram_dwnld/03-13-measure/2_cylinder_dif_side_and_mushrooms/preprocessed/0006_id1_ns1875848749103') + + +###Конфиги + +MODE = 'point' # Один из 2х режимов `point` or `extended` + +eps_r = 1.0 # <-- Диэлектрическая проницаемость среды +v = 3e8 / np.sqrt(eps_r) # скорость света в среде + + +# КООРДИНАТЫ АНТЕНН +# Координаты вдоль оси X на поверхности (z = 0), в метрах + +# Физические координаты антенн по их реальным индексам +TX_POSITIONS = {0: 90.5*0.01, # Tx с индексом 0 → x = +0.9 м + 3: -90.5*0.01} # Tx с индексом 3 → x = -0.9 м + +RX_POSITIONS = {0: -18*0.01, + 1: 48.5*0.01, + 2: -49.0*0.01, + 3: 18.5*0.01} + +# Массивы позиций в порядке возрастания физических индексов +x_tx = np.array([TX_POSITIONS[i] for i in sorted(TX_POSITIONS)]) +x_rx = np.array([RX_POSITIONS[j] for j in sorted(RX_POSITIONS)]) + +# Параметры алгоритма +SNR_THRESH = 3.0 # минимальный SNR пика +SNR_COMP_MAX = 20.0 # Верхний порог для компенсированного значения SNR + +# Параметр: максимальное число объектов для поиска +MAX_OBJECTS = 15 # <-- настройте под вашу задачу + +# Границы сетки аккумулятора +x_min, x_max = x_tx.min() - 2.0, x_tx.max() + 2.0 # [м] +z_min, z_max = 0.10, MAX_DEPTH # <-- глубина [м] + + + +# ══════════════════════════════════════════════════════ +# 1. ЗАГРУЗКА РЕАЛЬНЫХ ДАННЫХ +# ══════════════════════════════════════════════════════ + +def load_mimo_data(data_path, input_idx, output_idx): + data_path = Path(data_path) + s21_data = {} + freq_data = {} + + for f in data_path.glob("i*_o*_s21.npy"): + name = f.stem + parts = name.split('_') + i_tx_phys = int(parts[1][1:]) + i_rx_phys = int(parts[0][1:]) + + if i_tx_phys not in output_idx or i_rx_phys not in input_idx: + continue + + # Переводим физический индекс → порядковый (0,1,2,...) + i_tx = sorted(output_idx).index(i_tx_phys) + i_rx = sorted(input_idx).index(i_rx_phys) + + s21_data[(i_tx, i_rx)] = np.load(f) + freq_file = data_path / f"i{i_rx_phys}_o{i_tx_phys}_freq.npy" + freq_data[(i_tx, i_rx)] = np.load(freq_file) + + tx_indices = sorted(set(k[0] for k in s21_data)) + rx_indices = sorted(set(k[1] for k in s21_data)) + n_tx = len(tx_indices) + n_rx = len(rx_indices) + + return s21_data, freq_data, n_tx, n_rx + + +# Загрузка данных +s21_data, freq_data, N_tx, N_rx = load_mimo_data(DATA_PATH, INPUT_IDX, OUTPUT_IDX) +N_pairs = len(s21_data) + +# ══════════════════════════════════════════════════════ +# 1б. ВЫЧИСЛЕНИЕ СРЕДНЕГО ФОНА ПО ВСЕМ СНИМКАМ +# ══════════════════════════════════════════════════════ + +def compute_background(bg_path, input_idx, output_idx): + """ + Для каждой пары (i_tx, i_rx) усредняем S21 по всем снимкам в папке. + + Возвращает: + bg : dict[(i_tx, i_rx)] → np.array (complex), усреднённый S21 + """ + bg_path = Path(bg_path) + snapshots = sorted(bg_path.glob("*/")) # каждый подкаталог — один снимок + snapshots = [s for s in snapshots if s.is_dir()] + + if len(snapshots) == 0: + print("⚠️ Снимков для фона не найдено, BG_SUBTRACT отключён.") + return None + + + # Накопитель: для каждой пары суммируем S21 + bg_sum = {} + bg_count = {} + + for snap_dir in snapshots: + for f in snap_dir.glob("i*_o*_s21.npy"): + name = f.stem + parts = name.split('_') + i_tx_phys = int(parts[1][1:]) + i_rx_phys = int(parts[0][1:]) + + if i_tx_phys not in output_idx or i_rx_phys not in input_idx: + continue + + i_tx = sorted(output_idx).index(i_tx_phys) + i_rx = sorted(input_idx).index(i_rx_phys) + key = (i_tx, i_rx) + + s21 = np.load(f) + if key not in bg_sum: + bg_sum[key] = np.zeros_like(s21, dtype=complex) + bg_count[key] = 0 + bg_sum[key] += s21 + bg_count[key] += 1 + + bg = {key: bg_sum[key] / bg_count[key] for key in bg_sum} + # print(f"готово. Пар: {len(bg)}, снимков на пару: " + # f"{list(bg_count.values())[0] if bg_count else 0}") + return bg + + +if BG_SUBTRACT: + background = compute_background(BG_PATH, INPUT_IDX, OUTPUT_IDX) + if background is None: + BG_SUBTRACT = False # автоматически выключаем если нет данных +else: + background = None + # print("BG_SUBTRACT = False, вычитание фона отключено.") + + +# Проверка частот (берём из первой пары как референс) +first_key = list(freq_data.keys())[0] +freqs = freq_data[first_key] + + + +# Проверим что частоты одинаковые для всех пар +for key, freq in freq_data.items(): + if not np.allclose(freq, freqs): + print(f"⚠️ Частоты для пары {key} отличаются!") + + + +mask_freq = (freqs >= F_START) & (freqs <= F_STOP) + +freqs = freqs[mask_freq] + +f_min, f_max = freqs[0], freqs[-1] +BW = f_max - f_min +N_f = len(freqs) + + +# ══════════════════════════════════════════════════════ +# 2. ПАРАМЕТРЫ СИСТЕМЫ +# ══════════════════════════════════════════════════════ + +# Проверка соответствия координатов антенн +assert len(x_tx) == N_tx, f"x_tx должен содержать {N_tx} элементов" +assert len(x_rx) == N_rx, f"x_rx должен содержать {N_rx} элементов" + + + +# SHELL_SIGMA — ширина гауссовой оболочки +SHELL_SIGMA = v / BW * 0.5 # [м] + + +# Сетка аккумулятора + +x_grid = np.linspace(x_min, x_max, 300) +z_grid = np.linspace(z_min, z_max, 300) +XX, ZZ = np.meshgrid(x_grid, z_grid) + +# Расстояния от сетки до каждой антенны +R_tx_grid = {i: np.sqrt((XX - x_tx[i])**2 + ZZ**2) for i in range(N_tx)} +R_rx_grid = {j: np.sqrt((XX - x_rx[j])**2 + ZZ**2) for j in range(N_rx)} + + +# ══════════════════════════════════════════════════════ +# 3. ВЫЧИСЛЕНИЕ A-СКАНОВ ИЗ РЕАЛЬНЫХ S21 +# ══════════════════════════════════════════════════════ + +def compute_ascan(s21, freq, f_start, f_stop, window=True): + """ + S21(f) → IFFT → A-скан с правильным частотным сдвигом. + + Проблема наивного подхода (buf[:n] = s21): + IFFT считает, что спектр начинается с 0 Гц. + Реальные данные начинаются с f[0] > 0, поэтому + нулевая задержка смещается и в A-скане появляются биения. + + Правильный подход — сдвиг спектра: + Шаг частотной сетки df вычисляется из данных. + Индекс первой частоты: k0 = round(f[0] / df). + Данные кладутся в H[k0 : k0+n], а не в H[0 : n]. + Тогда IFFT корректно восстанавливает временной сигнал + с нулевой задержкой в t=0. + + Размер FFT: + Минимум для покрытия всего диапазона [0, f[-1]]: + min_len = 2 * (k0 + n - 1) + Округляем вверх до степени двойки для скорости FFT. + """ + mask_freq_ = (freq >= f_start) & (freq <= f_stop) + freq = freq[mask_freq_] + s21 = s21[mask_freq_] + + n = len(freq) + if n < 2: + raise ValueError("Слишком мало частотных точек") + + # Шаг частотной сетки + df = (freq[-1] - freq[0]) / (n - 1) + if df <= 0: + raise ValueError("Частоты не возрастают") + + # Индекс первой частоты в полной сетке от 0 до f[-1] + k0 = int(np.round(freq[0] / df)) + + # Минимальный размер FFT, округлённый до степени двойки + min_len = 2 * (k0 + n - 1) + n_fft = 1 << int(np.ceil(np.log2(min_len))) + + + # Временна́я ось — пересчитываем из нового n_fft + dt = 1.0 / (n_fft * df) + t_sec = np.arange(n_fft, dtype=float) * dt + + # Оконная функция (подавление боковых лепестков IFFT) + s = s21 * np.hanning(n) if window else s21.copy() + + # Спектр со сдвигом: данные на своём месте в частотной сетке + H = np.zeros(n_fft, dtype=np.complex128) + H[k0 : k0 + n] = s + + y = np.abs(np.fft.ifft(H)) + + + return t_sec[:y.size], y[:y.size] + + + +A = {} # A[(i,j)] — амплитудный A-скан +T_h = {} # T_h[(i,j)] — временна́я ось для этой пары [с] +Z_h = {} # Z_h[(i,j)] — ось глубины [м] + +for (i, j), s21 in s21_data.items(): + s21_proc = s21.copy() + + # Вычитание фона в частотной области + if BG_SUBTRACT and background is not None and (i, j) in background: + s21_proc = s21_proc - background[(i, j)] + # Примечание: вычитаем до обрезки по частоте и до окна — + # фон вычисляется из полных (необрезанных) данных, + # поэтому вычитание корректно в полном частотном диапазоне. + + t_pair, a_pair = compute_ascan(s21_proc, freq_data[(i, j)], + f_start=F_START, f_stop=F_STOP) + T_h[(i, j)] = t_pair + Z_h[(i, j)] = t_pair * v / 2 + A[(i, j)] = a_pair + +bg_label = "с вычитанием фона" if BG_SUBTRACT else "без вычитания фона" + + + +# Общая ось z для визуализации и поиска пиков +# (берём максимальный диапазон по всем парам) +z_h = Z_h[list(Z_h.keys())[0]] # все пары дают одинаковую ось, если freq совпадают +t_h = T_h[list(T_h.keys())[0]] + + +# 4. ВИЗУАЛИЗАЦИЯ A-СКАНОВ + + +def plot_ascans(A, z_h, n_tx, n_rx): + """Отображение всех A-сканов""" + fig, axes = plt.subplots(n_tx, n_rx, figsize=(3*n_rx, 3*n_tx), + sharex=True, sharey=True) + + if n_tx == 1: + axes = axes.reshape(1, -1) + if n_rx == 1: + axes = axes.reshape(-1, 1) + + for i in range(n_tx): + for j in range(n_rx): + ax = axes[i, j] + if (i, j) in A: + ax.plot(z_h, A[(i, j)], 'b-', lw=0.8) + ax.set_title(f'Tx{i} → Rx{j}', fontsize=10) + ax.grid(True, alpha=0.3) + else: + ax.set_visible(False) + + axes[-1, 0].set_xlabel('Глубина z [м]') + axes[0, 0].set_ylabel('Амплитуда') + fig.suptitle('A-сканы всех пар Tx-Rx', fontsize=12) + plt.tight_layout() + plt.show() + + +def plot_bscan(A, z_h, x_tx, x_rx): + """B-скан: все A-сканы рядом, отсортированные по виртуальной позиции""" + pairs = sorted(A.keys(), key=lambda p: (x_tx[p[0]] + x_rx[p[1]]) / 2) + + bscan = np.array([A[p] for p in pairs]).T + x_virt = [(x_tx[p[0]] + x_rx[p[1]]) / 2 for p in pairs] + + plt.figure(figsize=(10, 6)) + plt.imshow(bscan, aspect='auto', origin='lower', + extent=[min(x_virt), max(x_virt), z_h[0], z_h[-1]], + cmap='jet') + plt.colorbar(label='Амплитуда') + plt.xlabel('Виртуальная позиция X [м]') + plt.ylabel('Глубина Z [м]') + plt.title('B-скан (все пары)') + plt.show() + + + +# ══════════════════════════════════════════════════════ +# 5. ДЕТЕКТИРОВАНИЕ ПИКОВ +# ══════════════════════════════════════════════════════ + +def attenuation_at_depth(i_tx, i_rx, z_app): + """ + Ожидаемое ослабление geo·pattern для точки прямо под виртуальным + центром пары на глубине z_app. + + Используется для компенсации: реальный SNR пика делится на это + значение, чтобы вес глубокого/углового объекта не занижался. + """ + xc = (x_tx[i_tx] + x_rx[i_rx]) / 2.0 # виртуальный центр + Rtx = np.sqrt((xc - x_tx[i_tx])**2 + z_app**2) + Rrx = np.sqrt((xc - x_rx[i_rx])**2 + z_app**2) + geo = 1.0 / (Rtx * Rrx + 1e-12) + pat = (z_app / (Rtx + 1e-12))**2 * (z_app / (Rrx + 1e-12))**2 + return geo * pat + 1e-30 # +ε чтобы не делить на ноль + + +def find_peaks_snr(i_tx, i_rx, SNR_COMP_MAX = SNR_COMP_MAX): + """ + Поиск пиков A-скана. + Возвращает список dict: + z_app — кажущаяся глубина [м] + tau — задержка [с] + snr_raw — SNR без компенсации = пик / медиана + snr_comp— SNR с компенсацией ослабления (используется в аккумуляторе) + """ + ascan = A[(i_tx, i_rx)] + z_h_ij = Z_h[(i_tx, i_rx)] + t_h_ij = T_h[(i_tx, i_rx)] + noise = np.median(ascan) + i_min = np.searchsorted(z_h_ij, MIN_DEPTH) + i_max = np.searchsorted(z_h_ij, MAX_DEPTH) + min_dist = max(4, int(v / (2*BW) / (z_h_ij[1] - z_h_ij[0]) * 0.7)) + idx, _ = find_peaks(ascan[i_min:i_max], + height=noise * SNR_THRESH, + distance=min_dist) + idx += i_min + result = [] + for p in idx: + z_app = float(z_h_ij[p]) + snr_raw = float(ascan[p] / noise) + atten = attenuation_at_depth(i_tx, i_rx, z_app) + atten_norm = atten / attenuation_at_depth(i_tx, i_rx, 2.0) # Референсная глубина - 2м + snr_comp = snr_raw / (atten_norm ** COMP_POWER + 1e-12) + snr_comp = min(snr_comp, SNR_COMP_MAX) # ← clipping + result.append({'z_app': z_app, + 'tau': float(t_h_ij[p]), + 'snr_raw': snr_raw, + 'snr_comp': snr_comp}) + return result + + +peaks = {(i, j): find_peaks_snr(i, j) + for i in range(N_tx) for j in range(N_rx)} +n_total = sum(len(v2) for v2 in peaks.values()) + +# ══════════════════════════════════════════════════════ +# 5. АККУМУЛЯТОР +# ══════════════════════════════════════════════════════ + +def build_accumulator(exclude_z_ranges): + """ + Для каждого пика строим гауссову оболочку вокруг эллипса. + """ + acc = np.zeros_like(XX) + for i in range(N_tx): + for j in range(N_rx): + for pk in peaks[(i, j)]: + if any(lo <= pk['z_app'] <= hi for lo, hi in exclude_z_ranges): + continue + R_total = v * pk['tau'] + residual = R_tx_grid[i] + R_rx_grid[j] - R_total + shell = np.exp(-0.5 * (residual / SHELL_SIGMA)**2) + acc += shell * pk['snr_comp'] # Это и есть метрика на карте накопления + return acc + + +def count_agreeing_ellipses(x_est, z_est, exclude_z_ranges): + """ + Score = количество пар, чей эллипс проходит + через точку (x_est, z_est) с невязкой δ < 3σ. + """ + count = 0 + for i in range(N_tx): + for j in range(N_rx): + for pk in peaks[(i, j)]: + if any(lo <= pk['z_app'] <= hi for lo, hi in exclude_z_ranges): + continue + Rt = np.sqrt((x_est - x_tx[i])**2 + z_est**2) + Rr = np.sqrt((x_est - x_rx[j])**2 + z_est**2) + if abs(Rt + Rr - v*pk['tau']) < SHELL_SIGMA * 6: + count += 1 + break # одна пара — один голос + return count + +# ══════════════════════════════════════════════════════ +# 6. ПОИСК ОБЪЕКТОВ +# ══════════════════════════════════════════════════════ + +def find_centroid(acc_s, iz, ix, rpz, rpx): + """ + Взвешенный центроид аккумулятора в окрестности (iz, ix). + """ + NZ, NX = acc_s.shape + iz0 = max(0, iz - rpz); iz1 = min(NZ, iz + rpz) + ix0 = max(0, ix - rpx); ix1 = min(NX, ix + rpx) + patch = acc_s[iz0:iz1, ix0:ix1].copy() + W = patch.sum() + if W <= 0: + return x_grid[ix], z_grid[iz] + rows = np.arange(iz0, iz1)[:, None] * np.ones(patch.shape) + cols = np.ones(patch.shape) * np.arange(ix0, ix1)[None, :] + iz_c = int(round(np.clip((rows * patch).sum() / W, 0, NZ-1))) + ix_c = int(round(np.clip((cols * patch).sum() / W, 0, NX-1))) + return x_grid[ix_c], z_grid[iz_c] + + +def clean_find(n_search=10, suppress_r_cm=7, thresh_frac=0.05): + """ + CLEAN-итерация для точечных объектов. + """ + dx = x_grid[1] - x_grid[0] + dz = z_grid[1] - z_grid[0] + rpx = int(suppress_r_cm / 100 / dx) + rpz = int(suppress_r_cm / 100 / dz) + + excl_z = [] + found = [] + acc_initial = build_accumulator([]) + + for step in range(n_search): + acc = build_accumulator(excl_z) + acc_s = gaussian_filter(acc, sigma=3) + + if acc_s.max() < thresh_frac * acc_initial.max(): + break + + iz, ix = np.unravel_index(acc_s.argmax(), acc_s.shape) + x_est, z_est = find_centroid(acc_s, iz, ix, rpz, rpx) + + score = count_agreeing_ellipses(x_est, z_est, excl_z) + found.append({'x': x_est, 'z': z_est, 'score': score}) + + # Пики, соответствующие этому объекту + matched = [pk['z_app'] + for i in range(N_tx) for j in range(N_rx) + for pk in peaks[(i, j)] + if not any(lo<=pk['z_app']<=hi for lo,hi in excl_z) + and abs(np.sqrt((x_est-x_tx[i])**2+z_est**2) + + np.sqrt((x_est-x_rx[j])**2+z_est**2) - + v*pk['tau']) < SHELL_SIGMA * 3] #Возможны изменения + + if matched: + margin = SHELL_SIGMA * 1.0 + excl_z.append((min(matched) - margin, max(matched) + margin)) + + return found, acc_initial + + +def extended_find(thresh_frac=0.75, min_area_cm2=2.0): + """ + Режим для протяжённых объектов (труба, плита и т.п.). + """ + acc = build_accumulator([]) + acc_s = gaussian_filter(acc, sigma=3) + binary = acc_s > thresh_frac * acc_s.max() + + dx = (x_grid[1]-x_grid[0])*100 + dz = (z_grid[1]-z_grid[0])*100 + min_pix = int(min_area_cm2 / (dx * dz)) + + labeled, n = label(binary) + regions = [] + for k in range(1, n+1): + mask = labeled == k + if mask.sum() < min_pix: + continue + w = acc_s[mask] + xs = XX[mask]; zs = ZZ[mask] + xc = (xs * w).sum() / w.sum() + zc = (zs * w).sum() / w.sum() + score = count_agreeing_ellipses(xc, zc, []) + regions.append({'x': xc, 'z': zc, 'score': score, + 'mask': mask, 'n_pix': mask.sum()}) + return regions, acc + + + +if MODE == 'point': + found, accum = clean_find(n_search=MAX_OBJECTS) + # print(f"найдено {len(found)} объектов.") + # for k, obj in enumerate(found): + # print(f" [{k+1}] x={obj['x']*100:+.1f} см, " + # f"z={obj['z']*100:.1f} см, " + # f"score={obj['score']}/{N_pairs}") +else: + regions, accum = extended_find() + found = regions + # print(f"найдено {len(regions)} регионов.") + # for k, r in enumerate(regions): + # print(f" [{k+1}] центр x={r['x']*100:+.1f} см, " + # f"z={r['z']*100:.1f} см, score={r['score']}/{N_pairs}") + +# ══════════════════════════════════════════════════════ +# 7. ГРАФИКИ +# ══════════════════════════════════════════════════════ + + + +# ─── График 2: карта накопления ──────────────────── +fig, ax = plt.subplots(figsize=(12, 7)) +acc_s = gaussian_filter(accum, sigma=3) +im = ax.imshow(acc_s, + extent=[x_grid[0]*100, x_grid[-1]*100, + z_grid[-1]*100, z_grid[0]*100], + aspect='auto', origin='upper', cmap='hot', + vmin=acc_s.max()*0.35, vmax=acc_s.max()*0.95) +plt.colorbar(im, ax=ax, label='Накопленный вес (SNR_comp × гауссова оболочка)') + +# Позиции антенн +ax.plot(x_tx*100, np.zeros(N_tx), 'r^', ms=10, label='Tx', zorder=5) +ax.plot(x_rx*100, np.zeros(N_rx), 'bv', ms=10, label='Rx', zorder=5) + +# Найденные объекты +for obj in found: + lbl = f"score={obj['score']}/{N_pairs}" + ax.plot(obj['x']*100, obj['z']*100, 'wD', ms=9, zorder=11, + markeredgecolor='black', mew=1.2) + ax.annotate(lbl, (obj['x']*100, obj['z']*100), + textcoords='offset points', xytext=(6, 4), + fontsize=8, color='white', + bbox=dict(boxstyle='round,pad=0.2', fc='black', alpha=0.5)) + +if MODE == 'extended': + for r in regions: + ax.contour(x_grid*100, z_grid*100, r['mask'].astype(float), + levels=[0.5], colors=['cyan'], linewidths=[1.2]) + +ax.plot([], [], 'wD', ms=9, markeredgecolor='k', mew=1.2, label='Найденные объекты') +ax.set_xlabel("X [см]"); ax.set_ylabel("Глубина Z [см]") +ax.set_title("Карта накопления эллипсов\n" + f"score = число пар из {N_pairs}, чей эллипс проходит через точку") +ax.set_xlim(x_grid[0]*100, x_grid[-1]*100) +ax.set_ylim(z_grid[-1]*100, z_grid[0]*100) +ax.legend(loc='lower right', fontsize=9) +ax.grid(alpha=0.25) +ax.invert_yaxis() +plt.tight_layout() +plt.show() \ No newline at end of file diff --git a/Makefile b/Makefile index 75d2816..2959e0c 100644 --- a/Makefile +++ b/Makefile @@ -45,6 +45,7 @@ PREPROC_SOURCES := \ PROCESSOR_SOURCES := \ data_acq_and_processing/processing/processors/src/bscan_processor.cpp \ + data_acq_and_processing/processing/processors/src/gpr_processor.cpp \ data_acq_and_processing/processing/processors/src/passthrough_processor.cpp \ data_acq_and_processing/processing/data_processor/src/processing_live_config.cpp \ data_acq_and_processing/processing/data_processor/src/data_processor.cpp \ diff --git a/data_acq_and_processing/common_cpp/config/include/run_config.hpp b/data_acq_and_processing/common_cpp/config/include/run_config.hpp index 3728db3..eea690c 100644 --- a/data_acq_and_processing/common_cpp/config/include/run_config.hpp +++ b/data_acq_and_processing/common_cpp/config/include/run_config.hpp @@ -81,6 +81,26 @@ struct PreprocessConfig { std::string reference_bundle_path{}; }; +struct GprTxGeometry { + // Transmitter geometry keyed by output switch position. + std::uint32_t output_pos = 0; + float x_m = 0.0F; +}; + +struct GprRxGeometry { + // Receiver geometry keyed by input switch position. + std::uint32_t input_pos = 0; + float x_m = 0.0F; +}; + +struct GprConfig { + // Stable collection-level GPR processing settings. + std::string mode = "point"; + float relative_permittivity = 1.0F; + std::vector tx_geometry{}; + std::vector rx_geometry{}; +}; + struct RunConfig { // Single configuration object shared by all C++ processes. RadarConfig radar{}; @@ -89,6 +109,7 @@ struct RunConfig { RingsConfig rings{}; RuntimeConfig runtime{}; PreprocessConfig preprocess{}; + GprConfig gpr{}; std::vector run_combos{}; }; diff --git a/data_acq_and_processing/common_cpp/config/src/run_config.cpp b/data_acq_and_processing/common_cpp/config/src/run_config.cpp index 9d8cae9..1cec13d 100644 --- a/data_acq_and_processing/common_cpp/config/src/run_config.cpp +++ b/data_acq_and_processing/common_cpp/config/src/run_config.cpp @@ -7,6 +7,8 @@ #include #include #include +#include +#include #include @@ -214,6 +216,69 @@ void validate_combo_key(const ipc::ComboKey& key, const RunConfig& config) { } } +void validate_gpr_config(const GprConfig& config, const RunConfig& run_config) { + if (config.mode != "point" && config.mode != "extended") { + throw std::runtime_error("gpr.mode must be either 'point' or 'extended'"); + } + if (!(config.relative_permittivity > 0.0F)) { + throw std::runtime_error("gpr.relative_permittivity must be > 0"); + } + + std::unordered_set seen_tx_positions{}; + for (const auto& entry : config.tx_geometry) { + if (entry.output_pos >= run_config.output_switch.positions) { + throw std::runtime_error("gpr.tx_geometry output_pos is out of range"); + } + if (!seen_tx_positions.insert(entry.output_pos).second) { + throw std::runtime_error("gpr.tx_geometry contains duplicate output_pos"); + } + } + + std::unordered_set seen_rx_positions{}; + for (const auto& entry : config.rx_geometry) { + if (entry.input_pos >= run_config.input_switch.positions) { + throw std::runtime_error("gpr.rx_geometry input_pos is out of range"); + } + if (!seen_rx_positions.insert(entry.input_pos).second) { + throw std::runtime_error("gpr.rx_geometry contains duplicate input_pos"); + } + } +} + +[[nodiscard]] auto parse_gpr_config(const Json& object, const RunConfig& run_config) -> GprConfig { + const auto* gpr_obj = as_object(object, "gpr"); + GprConfig config{}; + config.mode = optional_string(*gpr_obj, "mode", "point"); + config.relative_permittivity = optional_f32(*gpr_obj, "relative_permittivity", 1.0F); + + if (const auto* tx_value = optional_field(*gpr_obj, "tx_geometry"); tx_value != nullptr) { + const auto* tx_array = as_array(*tx_value, "gpr.tx_geometry"); + config.tx_geometry.reserve(tx_array->size()); + for (const auto& entry_value : *tx_array) { + const auto* entry_obj = as_object(entry_value, "gpr.tx_geometry[]"); + GprTxGeometry entry{}; + entry.output_pos = optional_u32(*entry_obj, "output_pos", 0U); + entry.x_m = optional_f32(*entry_obj, "x_m", 0.0F); + config.tx_geometry.push_back(std::move(entry)); + } + } + + if (const auto* rx_value = optional_field(*gpr_obj, "rx_geometry"); rx_value != nullptr) { + const auto* rx_array = as_array(*rx_value, "gpr.rx_geometry"); + config.rx_geometry.reserve(rx_array->size()); + for (const auto& entry_value : *rx_array) { + const auto* entry_obj = as_object(entry_value, "gpr.rx_geometry[]"); + GprRxGeometry entry{}; + entry.input_pos = optional_u32(*entry_obj, "input_pos", 0U); + entry.x_m = optional_f32(*entry_obj, "x_m", 0.0F); + config.rx_geometry.push_back(std::move(entry)); + } + } + + validate_gpr_config(config, run_config); + return config; +} + [[nodiscard]] auto parse_switch_config( const Json& object, const std::string& default_name, @@ -379,6 +444,12 @@ auto load_run_config(const std::string& path) -> RunConfig { config.preprocess.reference_bundle_path = optional_string(*preprocess_obj, "reference_bundle_path", ""); } + if (const auto* gpr_value = optional_field(*root_obj, "gpr"); gpr_value != nullptr) { + config.gpr = parse_gpr_config(*gpr_value, config); + } else { + validate_gpr_config(config.gpr, config); + } + { const auto* rings_obj = as_object(required_field(*root_obj, "rings"), "rings"); config.rings.raw = parse_ring_endpoint(*rings_obj, "raw", config.rings.raw); diff --git a/data_acq_and_processing/common_cpp/ipc/include/shared_types.hpp b/data_acq_and_processing/common_cpp/ipc/include/shared_types.hpp index 1118951..4507b8d 100644 --- a/data_acq_and_processing/common_cpp/ipc/include/shared_types.hpp +++ b/data_acq_and_processing/common_cpp/ipc/include/shared_types.hpp @@ -49,6 +49,8 @@ using PreprocessedCollection = RawSweepCollection; enum class ResultKind : std::uint8_t { TraceComplex = 1, ScalarF32 = 2, + ImageF32 = 3, + TableF32 = 4, }; struct ResultPayload { @@ -60,6 +62,13 @@ struct ResultPayload { std::vector trace{}; // Valid for ScalarF32 payloads. float scalar_value = 0.0F; + // Valid for ImageF32 payloads. + std::vector image_x_axis{}; + std::vector image_y_axis{}; + std::vector image_values{}; + // Valid for TableF32 payloads. + std::uint32_t table_columns = 0; + std::vector table_values{}; }; struct ResultBlock { @@ -70,6 +79,7 @@ struct ResultBlock { struct ResultCollection { std::uint64_t collection_id = 0; std::uint64_t monotonic_ns = 0; + std::vector collection_payloads{}; std::vector blocks{}; }; diff --git a/data_acq_and_processing/common_cpp/ipc/src/shared_types.cpp b/data_acq_and_processing/common_cpp/ipc/src/shared_types.cpp index 45eb1c7..4b8dd58 100644 --- a/data_acq_and_processing/common_cpp/ipc/src/shared_types.cpp +++ b/data_acq_and_processing/common_cpp/ipc/src/shared_types.cpp @@ -200,6 +200,75 @@ auto read_trace_result_payload(BinaryReader& reader, ResultPayload* payload) -> } } +void write_image_result_payload(BinaryWriter& writer, const ResultPayload& payload) { + const auto x_count = checked_count_to_u32(payload.image_x_axis.size(), "Result image x axis point count"); + const auto y_count = checked_count_to_u32(payload.image_y_axis.size(), "Result image y axis point count"); + const std::size_t expected_value_count = static_cast(x_count) * static_cast(y_count); + if (payload.image_values.size() != expected_value_count) { + throw std::runtime_error("Result image payload has inconsistent axis/value sizes"); + } + + writer.write(x_count); + writer.write(y_count); + for (const auto value : payload.image_x_axis) { + writer.write(value); + } + for (const auto value : payload.image_y_axis) { + writer.write(value); + } + for (const auto value : payload.image_values) { + writer.write(value); + } +} + +auto read_image_result_payload(BinaryReader& reader, ResultPayload* payload) -> void { + const auto x_count = reader.read(); + const auto y_count = reader.read(); + + payload->image_x_axis.reserve(x_count); + payload->image_y_axis.reserve(y_count); + payload->image_values.reserve(static_cast(x_count) * static_cast(y_count)); + + for (std::uint32_t index = 0; index < x_count; ++index) { + payload->image_x_axis.push_back(reader.read()); + } + for (std::uint32_t index = 0; index < y_count; ++index) { + payload->image_y_axis.push_back(reader.read()); + } + for (std::size_t index = 0; index < static_cast(x_count) * static_cast(y_count); ++index) { + payload->image_values.push_back(reader.read()); + } +} + +void write_table_result_payload(BinaryWriter& writer, const ResultPayload& payload) { + const auto column_count = payload.table_columns; + if (column_count == 0U && !payload.table_values.empty()) { + throw std::runtime_error("Result table payload must declare table_columns when values are present"); + } + if (column_count != 0U && (payload.table_values.size() % column_count) != 0U) { + throw std::runtime_error("Result table payload has inconsistent column/value sizes"); + } + + const std::uint32_t row_count = + column_count == 0U ? 0U : checked_count_to_u32(payload.table_values.size() / column_count, "Result table row count"); + + writer.write(column_count); + writer.write(row_count); + for (const auto value : payload.table_values) { + writer.write(value); + } +} + +auto read_table_result_payload(BinaryReader& reader, ResultPayload* payload) -> void { + const auto column_count = reader.read(); + const auto row_count = reader.read(); + payload->table_columns = column_count; + payload->table_values.reserve(static_cast(column_count) * static_cast(row_count)); + for (std::size_t index = 0; index < static_cast(column_count) * static_cast(row_count); ++index) { + payload->table_values.push_back(reader.read()); + } +} + void write_result_payload(BinaryWriter& writer, const ResultPayload& payload) { writer.write(static_cast(payload.kind)); writer.write_string(payload.processing_name); @@ -211,6 +280,12 @@ void write_result_payload(BinaryWriter& writer, const ResultPayload& payload) { case ResultKind::ScalarF32: writer.write(payload.scalar_value); return; + case ResultKind::ImageF32: + write_image_result_payload(writer, payload); + return; + case ResultKind::TableF32: + write_table_result_payload(writer, payload); + return; default: throw std::runtime_error("Unsupported result payload kind"); } @@ -228,6 +303,12 @@ void write_result_payload(BinaryWriter& writer, const ResultPayload& payload) { case ResultKind::ScalarF32: payload.scalar_value = reader.read(); return payload; + case ResultKind::ImageF32: + read_image_result_payload(reader, &payload); + return payload; + case ResultKind::TableF32: + read_table_result_payload(reader, &payload); + return payload; default: throw std::runtime_error("Unsupported result payload kind in stream"); } @@ -296,8 +377,12 @@ auto serialize_result_collection(const ResultCollection& collection) -> std::vec writer.write(kResultCollectionMagic); writer.write(collection.collection_id); writer.write(collection.monotonic_ns); + writer.write(checked_count_to_u32(collection.collection_payloads.size(), "Collection payload count")); writer.write(checked_count_to_u32(collection.blocks.size(), "Result block count")); + for (const auto& payload : collection.collection_payloads) { + write_result_payload(writer, payload); + } for (const auto& block : collection.blocks) { write_result_block(writer, block); } @@ -317,7 +402,13 @@ auto deserialize_result_collection(std::span bytes) -> Resul collection.collection_id = reader.read(); collection.monotonic_ns = reader.read(); + const auto collection_payload_count = reader.read(); const auto block_count = reader.read(); + collection.collection_payloads.reserve(collection_payload_count); + for (std::uint32_t index = 0; index < collection_payload_count; ++index) { + collection.collection_payloads.push_back(read_result_payload(reader)); + } + collection.blocks.reserve(block_count); for (std::uint32_t index = 0; index < block_count; ++index) { collection.blocks.push_back(read_result_block(reader)); diff --git a/data_acq_and_processing/processing/data_processor/include/data_processor.hpp b/data_acq_and_processing/processing/data_processor/include/data_processor.hpp index 23562e4..61e2ba3 100644 --- a/data_acq_and_processing/processing/data_processor/include/data_processor.hpp +++ b/data_acq_and_processing/processing/data_processor/include/data_processor.hpp @@ -2,6 +2,7 @@ #include #include +#include #include #include #include @@ -29,11 +30,13 @@ class DataProcessor { private: [[nodiscard]] auto process_collection( const ipc::PreprocessedCollection& preprocessed, + std::span previous_collections, ProcessorInterface& processor, const ProcessingLiveConfig& live_config ) -> ipc::ResultCollection; [[nodiscard]] auto resolve_processor(const ProcessingLiveConfig& live_config) -> ProcessorInterface&; + [[nodiscard]] auto should_replay_entire_history(const ProcessingLiveConfig& live_config) const -> bool; const config::RunConfig& config_; ipc::ShmRing& preprocessed_ring_; diff --git a/data_acq_and_processing/processing/data_processor/include/processing_live_config.hpp b/data_acq_and_processing/processing/data_processor/include/processing_live_config.hpp index a15f32c..d5da5a3 100644 --- a/data_acq_and_processing/processing/data_processor/include/processing_live_config.hpp +++ b/data_acq_and_processing/processing/data_processor/include/processing_live_config.hpp @@ -4,6 +4,7 @@ #include #include #include +#include namespace radar::processing { @@ -17,12 +18,24 @@ struct ProcessingLiveConfig { std::string processor_mode = "pass_through"; float gain_db = 0.0F; float phase_deg = 0.0F; + bool pass_through_fixed_y_enabled = false; + float pass_through_y_min_db = -100.0F; + float pass_through_y_max_db = 0.0F; std::string bscan_axis = "abs"; float bscan_cut_m = 0.824F; float bscan_max_depth_m = 1.0F; float bscan_gain = 1.0F; float bscan_start_freq_mhz = 100.0F; float bscan_stop_freq_mhz = 8800.0F; + std::vector gpr_input_positions{}; + std::vector gpr_output_positions{}; + float gpr_min_depth_m = 2.0F; + float gpr_max_depth_m = 14.0F; + float gpr_comp_power = 0.2F; + float gpr_start_freq_mhz = 3000.0F; + float gpr_stop_freq_mhz = 6000.0F; + bool gpr_background_subtract_enabled = true; + std::uint32_t gpr_background_mean_count = 10U; std::uint64_t history_command_seq = 0; HistoryCommand history_command = HistoryCommand::None; }; diff --git a/data_acq_and_processing/processing/data_processor/src/data_processor.cpp b/data_acq_and_processing/processing/data_processor/src/data_processor.cpp index d39b5a0..8d261fc 100644 --- a/data_acq_and_processing/processing/data_processor/src/data_processor.cpp +++ b/data_acq_and_processing/processing/data_processor/src/data_processor.cpp @@ -3,12 +3,13 @@ #include #include #include -#include +#include #include #include #include #include "bscan_processor.hpp" +#include "gpr_processor.hpp" #include "passthrough_processor.hpp" namespace radar::processing { @@ -54,7 +55,7 @@ DataProcessor::DataProcessor( void DataProcessor::run(const std::atomic& stop_requested) { std::vector bytes{}; - std::deque preprocessed_history{}; + std::vector preprocessed_history{}; const std::size_t history_limit = replay_history_limit(config_); std::uint64_t last_replayed_revision = live_config_loader_.revision(); std::uint64_t last_applied_history_command_seq = 0; @@ -76,13 +77,23 @@ void DataProcessor::run(const std::atomic& stop_requested) { last_applied_history_command_seq = live_config.history_command_seq; } - if (live_config.processor_mode == "bscan") { - for (const auto& cached : preprocessed_history) { - const auto replay_result = process_collection(cached, processor, live_config); + if (should_replay_entire_history(live_config)) { + for (std::size_t index = 0; index < preprocessed_history.size(); ++index) { + const auto replay_result = process_collection( + preprocessed_history[index], + std::span(preprocessed_history.data(), index), + processor, + live_config + ); publish_result_collection(replay_result, results_ring_); } } else if (!preprocessed_history.empty()) { - const auto replay_result = process_collection(preprocessed_history.back(), processor, live_config); + const auto replay_result = process_collection( + preprocessed_history.back(), + std::span(preprocessed_history.data(), preprocessed_history.size() - 1U), + processor, + live_config + ); publish_result_collection(replay_result, results_ring_); } last_replayed_revision = live_revision; @@ -92,10 +103,15 @@ void DataProcessor::run(const std::atomic& stop_requested) { auto preprocessed = ipc::deserialize_preprocessed_collection(bytes); preprocessed_history.push_back(std::move(preprocessed)); while (preprocessed_history.size() > history_limit) { - preprocessed_history.pop_front(); + preprocessed_history.erase(preprocessed_history.begin()); } - const auto result_collection = process_collection(preprocessed_history.back(), processor, live_config); + const auto result_collection = process_collection( + preprocessed_history.back(), + std::span(preprocessed_history.data(), preprocessed_history.size() - 1U), + processor, + live_config + ); publish_result_collection(result_collection, results_ring_); continue; } @@ -106,25 +122,11 @@ void DataProcessor::run(const std::atomic& stop_requested) { auto DataProcessor::process_collection( const ipc::PreprocessedCollection& preprocessed, + std::span previous_collections, ProcessorInterface& processor, const ProcessingLiveConfig& live_config ) -> ipc::ResultCollection { - ipc::ResultCollection results{}; - results.collection_id = preprocessed.collection_id; - // Keep source monotonic timestamp stable across live-config replays. - results.monotonic_ns = preprocessed.monotonic_ns; - results.blocks.reserve(preprocessed.traces.size()); - - for (const auto& trace : preprocessed.traces) { - ipc::ResultBlock block{}; - block.combo = trace.combo; - block.payloads.reserve(1U); - block.payloads.push_back(processor.process(trace, live_config)); - - results.blocks.push_back(std::move(block)); - } - - return results; + return processor.process_collection(config_, preprocessed, previous_collections, live_config); } auto DataProcessor::resolve_processor(const ProcessingLiveConfig& live_config) -> ProcessorInterface& { @@ -140,6 +142,12 @@ auto DataProcessor::resolve_processor(const ProcessingLiveConfig& live_config) - return *(processors_.begin()->second); } +auto DataProcessor::should_replay_entire_history(const ProcessingLiveConfig& live_config) const -> bool { + const std::string requested_mode = + live_config.processor_mode.empty() ? default_processor_mode_ : live_config.processor_mode; + return requested_mode == "bscan"; +} + auto create_default_processors() -> ProcessorRegistry { ProcessorRegistry processors{}; { @@ -150,6 +158,10 @@ auto create_default_processors() -> ProcessorRegistry { auto processor = std::make_unique(); processors.emplace(processor->name(), std::move(processor)); } + { + auto processor = std::make_unique(); + processors.emplace(processor->name(), std::move(processor)); + } return processors; } diff --git a/data_acq_and_processing/processing/data_processor/src/processing_live_config.cpp b/data_acq_and_processing/processing/data_processor/src/processing_live_config.cpp index fc0124a..e74e8d0 100644 --- a/data_acq_and_processing/processing/data_processor/src/processing_live_config.cpp +++ b/data_acq_and_processing/processing/data_processor/src/processing_live_config.cpp @@ -45,6 +45,27 @@ using Json = nlohmann::json; return static_cast(rounded); } +[[nodiscard]] auto parse_u32_number(const Json& value, const std::string& field_name) -> std::uint32_t { + const auto parsed = parse_u64_number(value, field_name); + if (parsed > static_cast(std::numeric_limits::max())) { + throw std::runtime_error(field_name + " is out of uint32 range"); + } + return static_cast(parsed); +} + +[[nodiscard]] auto parse_u32_array(const Json& value, const std::string& field_name) -> std::vector { + if (!value.is_array()) { + throw std::runtime_error(field_name + " must be array"); + } + + std::vector result{}; + result.reserve(value.size()); + for (std::size_t index = 0; index < value.size(); ++index) { + result.push_back(parse_u32_number(value[index], field_name + "[" + std::to_string(index) + "]")); + } + return result; +} + [[nodiscard]] auto parse_live_config(const std::string& json_text, const std::string& path) -> ProcessingLiveConfig { if (json_text.empty()) { throw std::runtime_error("Processing live config is empty: " + path); @@ -80,6 +101,24 @@ using Json = nlohmann::json; } config.phase_deg = static_cast(found->get()); } + if (const auto found = root.find("pass_through_fixed_y_enabled"); found != root.end()) { + if (!found->is_boolean()) { + throw std::runtime_error("processing.pass_through_fixed_y_enabled must be bool"); + } + config.pass_through_fixed_y_enabled = found->get(); + } + if (const auto found = root.find("pass_through_y_min_db"); found != root.end()) { + if (!found->is_number()) { + throw std::runtime_error("processing.pass_through_y_min_db must be number"); + } + config.pass_through_y_min_db = static_cast(found->get()); + } + if (const auto found = root.find("pass_through_y_max_db"); found != root.end()) { + if (!found->is_number()) { + throw std::runtime_error("processing.pass_through_y_max_db must be number"); + } + config.pass_through_y_max_db = static_cast(found->get()); + } if (const auto found = root.find("bscan_axis"); found != root.end()) { if (!found->is_string()) { throw std::runtime_error("processing.bscan_axis must be string"); @@ -116,6 +155,51 @@ using Json = nlohmann::json; } config.bscan_stop_freq_mhz = static_cast(found->get()); } + if (const auto found = root.find("gpr_input_positions"); found != root.end()) { + config.gpr_input_positions = parse_u32_array(*found, "processing.gpr_input_positions"); + } + if (const auto found = root.find("gpr_output_positions"); found != root.end()) { + config.gpr_output_positions = parse_u32_array(*found, "processing.gpr_output_positions"); + } + if (const auto found = root.find("gpr_min_depth_m"); found != root.end()) { + if (!found->is_number()) { + throw std::runtime_error("processing.gpr_min_depth_m must be number"); + } + config.gpr_min_depth_m = static_cast(found->get()); + } + if (const auto found = root.find("gpr_max_depth_m"); found != root.end()) { + if (!found->is_number()) { + throw std::runtime_error("processing.gpr_max_depth_m must be number"); + } + config.gpr_max_depth_m = static_cast(found->get()); + } + if (const auto found = root.find("gpr_comp_power"); found != root.end()) { + if (!found->is_number()) { + throw std::runtime_error("processing.gpr_comp_power must be number"); + } + config.gpr_comp_power = static_cast(found->get()); + } + if (const auto found = root.find("gpr_start_freq_mhz"); found != root.end()) { + if (!found->is_number()) { + throw std::runtime_error("processing.gpr_start_freq_mhz must be number"); + } + config.gpr_start_freq_mhz = static_cast(found->get()); + } + if (const auto found = root.find("gpr_stop_freq_mhz"); found != root.end()) { + if (!found->is_number()) { + throw std::runtime_error("processing.gpr_stop_freq_mhz must be number"); + } + config.gpr_stop_freq_mhz = static_cast(found->get()); + } + if (const auto found = root.find("gpr_background_subtract_enabled"); found != root.end()) { + if (!found->is_boolean()) { + throw std::runtime_error("processing.gpr_background_subtract_enabled must be bool"); + } + config.gpr_background_subtract_enabled = found->get(); + } + if (const auto found = root.find("gpr_background_mean_count"); found != root.end()) { + config.gpr_background_mean_count = parse_u32_number(*found, "processing.gpr_background_mean_count"); + } if (const auto found = root.find("history_command_seq"); found != root.end()) { config.history_command_seq = parse_u64_number(*found, "processing.history_command_seq"); } diff --git a/data_acq_and_processing/processing/processors/include/bscan_processor.hpp b/data_acq_and_processing/processing/processors/include/bscan_processor.hpp index fc9b931..4c0f0bb 100644 --- a/data_acq_and_processing/processing/processors/include/bscan_processor.hpp +++ b/data_acq_and_processing/processing/processors/include/bscan_processor.hpp @@ -7,10 +7,12 @@ namespace radar::processing { class BScanProcessor final : public ProcessorInterface { public: [[nodiscard]] auto name() const -> std::string override; - [[nodiscard]] auto process( - const ipc::SweepTraceBlock& trace, + [[nodiscard]] auto process_collection( + const config::RunConfig& run_config, + const ipc::PreprocessedCollection& collection, + std::span previous_collections, const ProcessingLiveConfig& live_config - ) -> ipc::ResultPayload override; + ) -> ipc::ResultCollection override; }; } // namespace radar::processing diff --git a/data_acq_and_processing/processing/processors/include/gpr_processor.hpp b/data_acq_and_processing/processing/processors/include/gpr_processor.hpp new file mode 100644 index 0000000..bc06471 --- /dev/null +++ b/data_acq_and_processing/processing/processors/include/gpr_processor.hpp @@ -0,0 +1,18 @@ +#pragma once + +#include "processor_interface.hpp" + +namespace radar::processing { + +class GprProcessor final : public ProcessorInterface { + public: + [[nodiscard]] auto name() const -> std::string override; + [[nodiscard]] auto process_collection( + const config::RunConfig& run_config, + const ipc::PreprocessedCollection& collection, + std::span previous_collections, + const ProcessingLiveConfig& live_config + ) -> ipc::ResultCollection override; +}; + +} // namespace radar::processing diff --git a/data_acq_and_processing/processing/processors/include/passthrough_processor.hpp b/data_acq_and_processing/processing/processors/include/passthrough_processor.hpp index c849913..51b10da 100644 --- a/data_acq_and_processing/processing/processors/include/passthrough_processor.hpp +++ b/data_acq_and_processing/processing/processors/include/passthrough_processor.hpp @@ -7,10 +7,12 @@ namespace radar::processing { class PassThroughProcessor final : public ProcessorInterface { public: [[nodiscard]] auto name() const -> std::string override; - [[nodiscard]] auto process( - const ipc::SweepTraceBlock& trace, + [[nodiscard]] auto process_collection( + const config::RunConfig& run_config, + const ipc::PreprocessedCollection& collection, + std::span previous_collections, const ProcessingLiveConfig& live_config - ) -> ipc::ResultPayload override; + ) -> ipc::ResultCollection override; }; } // namespace radar::processing diff --git a/data_acq_and_processing/processing/processors/include/processor_interface.hpp b/data_acq_and_processing/processing/processors/include/processor_interface.hpp index 78990f4..93cd285 100644 --- a/data_acq_and_processing/processing/processors/include/processor_interface.hpp +++ b/data_acq_and_processing/processing/processors/include/processor_interface.hpp @@ -1,8 +1,10 @@ #pragma once +#include #include #include "processing_live_config.hpp" +#include "run_config.hpp" #include "shared_types.hpp" namespace radar::processing { @@ -12,10 +14,12 @@ class ProcessorInterface { virtual ~ProcessorInterface() = default; [[nodiscard]] virtual auto name() const -> std::string = 0; - [[nodiscard]] virtual auto process( - const ipc::SweepTraceBlock& trace, + [[nodiscard]] virtual auto process_collection( + const config::RunConfig& run_config, + const ipc::PreprocessedCollection& collection, + std::span previous_collections, const ProcessingLiveConfig& live_config - ) -> ipc::ResultPayload = 0; + ) -> ipc::ResultCollection = 0; }; } // namespace radar::processing diff --git a/data_acq_and_processing/processing/processors/src/bscan_processor.cpp b/data_acq_and_processing/processing/processors/src/bscan_processor.cpp index a20b01d..b78f932 100644 --- a/data_acq_and_processing/processing/processors/src/bscan_processor.cpp +++ b/data_acq_and_processing/processing/processors/src/bscan_processor.cpp @@ -215,23 +215,39 @@ auto BScanProcessor::name() const -> std::string { return "bscan"; } -auto BScanProcessor::process(const ipc::SweepTraceBlock& trace, const ProcessingLiveConfig& live_config) - -> ipc::ResultPayload { - ipc::ResultPayload payload{}; - payload.processing_name = name(); - payload.kind = ipc::ResultKind::TraceComplex; +auto BScanProcessor::process_collection( + const config::RunConfig& /*run_config*/, + const ipc::PreprocessedCollection& collection, + std::span /*previous_collections*/, + const ProcessingLiveConfig& live_config +) -> ipc::ResultCollection { + ipc::ResultCollection results{}; + results.collection_id = collection.collection_id; + results.monotonic_ns = collection.monotonic_ns; + results.blocks.reserve(collection.traces.size()); - auto profile = compute_bscan_profile(trace, live_config); - payload.frequency_hz = std::move(profile.depth_m); - payload.trace.reserve(profile.response.size()); - for (const auto value : profile.response) { - payload.trace.push_back(ipc::Complex32{ - .re = value, - .im = 0.0F, - }); + for (const auto& trace : collection.traces) { + ipc::ResultPayload payload{}; + payload.processing_name = name(); + payload.kind = ipc::ResultKind::TraceComplex; + + auto profile = compute_bscan_profile(trace, live_config); + payload.frequency_hz = std::move(profile.depth_m); + payload.trace.reserve(profile.response.size()); + for (const auto value : profile.response) { + payload.trace.push_back(ipc::Complex32{ + .re = value, + .im = 0.0F, + }); + } + + ipc::ResultBlock block{}; + block.combo = trace.combo; + block.payloads.push_back(std::move(payload)); + results.blocks.push_back(std::move(block)); } - return payload; + return results; } } // namespace radar::processing diff --git a/data_acq_and_processing/processing/processors/src/gpr_processor.cpp b/data_acq_and_processing/processing/processors/src/gpr_processor.cpp new file mode 100644 index 0000000..7bf0e3f --- /dev/null +++ b/data_acq_and_processing/processing/processors/src/gpr_processor.cpp @@ -0,0 +1,1149 @@ +#include "gpr_processor.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace radar::processing { +namespace { + +constexpr double kPi = 3.14159265358979323846; +constexpr double kSpeedOfLightMetersPerSec = 299'792'458.0; +constexpr double kAccumulatorXMarginM = 2.0; +constexpr double kAccumulatorZMinM = 0.10; +constexpr double kSnrThresh = 3.0; +constexpr double kSnrCompMax = 20.0; +constexpr std::size_t kGridWidth = 300U; +constexpr std::size_t kGridHeight = 300U; +constexpr double kGaussianSigma = 3.0; +constexpr double kCleanSuppressRadiusM = 0.07; +constexpr double kCleanThresholdFrac = 0.05; +constexpr std::size_t kMaxObjects = 15U; +constexpr double kExtendedThresholdFrac = 0.75; +constexpr double kExtendedMinAreaCm2 = 2.0; + +using PairKey = std::uint64_t; + +struct SelectedTrace { + std::vector frequency_hz{}; + std::vector> s21{}; +}; + +struct AscanResult { + std::vector time_s{}; + std::vector depth_m{}; + std::vector amplitude{}; + double bandwidth_hz = 0.0; +}; + +struct PeakRecord { + double z_app = 0.0; + double tau = 0.0; + double snr_raw = 0.0; + double snr_comp = 0.0; +}; + +struct PointRecord { + double x_m = 0.0; + double z_m = 0.0; + double score = 0.0; +}; + +struct RegionRecord { + double x_m = 0.0; + double z_m = 0.0; + double score = 0.0; + double pixel_count = 0.0; + std::vector mask{}; +}; + +struct GeometrySelection { + std::vector input_positions{}; + std::vector output_positions{}; + std::vector x_tx{}; + std::vector x_rx{}; + std::unordered_map input_local_by_pos{}; + std::unordered_map output_local_by_pos{}; +}; + +struct BackgroundAccumulator { + std::vector> sum{}; + std::size_t count = 0U; +}; + +struct GridDefinition { + std::vector x_grid{}; + std::vector z_grid{}; + std::vector> tx_distance_grids{}; + std::vector> rx_distance_grids{}; +}; + +[[nodiscard]] auto make_pair_key(std::uint32_t tx_local_index, std::uint32_t rx_local_index) -> PairKey { + return (static_cast(tx_local_index) << 32U) | static_cast(rx_local_index); +} + +[[nodiscard]] auto next_power_of_two(std::size_t value) -> std::size_t { + if (value <= 1U) { + return 1U; + } + + std::size_t power = 1U; + while (power < value) { + if (power > (std::numeric_limits::max() >> 1U)) { + return value; + } + power <<= 1U; + } + return power; +} + +void fft_inplace(std::vector>& values, bool inverse) { + const std::size_t size = values.size(); + if (size <= 1U) { + return; + } + + for (std::size_t index = 1U, bit_reversed = 0U; index < size; ++index) { + std::size_t bit = size >> 1U; + while (bit_reversed & bit) { + bit_reversed ^= bit; + bit >>= 1U; + } + bit_reversed ^= bit; + if (index < bit_reversed) { + std::swap(values[index], values[bit_reversed]); + } + } + + for (std::size_t len = 2U; len <= size; len <<= 1U) { + const double angle = 2.0 * kPi * (inverse ? 1.0 : -1.0) / static_cast(len); + const std::complex twiddle_step(std::cos(angle), std::sin(angle)); + const std::size_t half_len = len >> 1U; + + for (std::size_t offset = 0U; offset < size; offset += len) { + std::complex twiddle(1.0, 0.0); + for (std::size_t inner = 0U; inner < half_len; ++inner) { + const auto even = values[offset + inner]; + const auto odd = values[offset + inner + half_len] * twiddle; + values[offset + inner] = even + odd; + values[offset + inner + half_len] = even - odd; + twiddle *= twiddle_step; + } + } + } + + if (!inverse) { + return; + } + + const double scale = 1.0 / static_cast(size); + for (auto& value : values) { + value *= scale; + } +} + +[[nodiscard]] auto build_axis(double min_value, double max_value, std::size_t count) -> std::vector { + std::vector axis{}; + if (count == 0U) { + return axis; + } + + axis.resize(count, min_value); + if (count == 1U) { + return axis; + } + + const double step = (max_value - min_value) / static_cast(count - 1U); + for (std::size_t index = 0U; index < count; ++index) { + axis[index] = min_value + (step * static_cast(index)); + } + return axis; +} + +[[nodiscard]] auto clamp_index(std::ptrdiff_t value, std::size_t limit) -> std::size_t { + if (limit == 0U) { + return 0U; + } + if (value < 0) { + return 0U; + } + const auto max_index = static_cast(limit - 1U); + if (value > max_index) { + return limit - 1U; + } + return static_cast(value); +} + +[[nodiscard]] auto build_gaussian_kernel(double sigma) -> std::vector { + if (!(sigma > 0.0)) { + return {1.0}; + } + + const auto radius = static_cast(std::ceil(sigma * 3.0)); + std::vector kernel(static_cast((radius * 2) + 1), 0.0); + double sum = 0.0; + for (std::ptrdiff_t offset = -radius; offset <= radius; ++offset) { + const double value = std::exp(-0.5 * std::pow(static_cast(offset) / sigma, 2.0)); + kernel[static_cast(offset + radius)] = value; + sum += value; + } + if (sum > 0.0) { + for (auto& value : kernel) { + value /= sum; + } + } + return kernel; +} + +[[nodiscard]] auto gaussian_filter_2d( + const std::vector& values, + std::size_t width, + std::size_t height, + double sigma +) -> std::vector { + if (values.empty() || width == 0U || height == 0U) { + return {}; + } + + const auto kernel = build_gaussian_kernel(sigma); + const auto radius = static_cast((kernel.size() - 1U) / 2U); + std::vector temp(values.size(), 0.0); + std::vector output(values.size(), 0.0); + + for (std::size_t row = 0U; row < height; ++row) { + for (std::size_t col = 0U; col < width; ++col) { + double sum = 0.0; + for (std::ptrdiff_t offset = -radius; offset <= radius; ++offset) { + const auto sample_col = clamp_index(static_cast(col) + offset, width); + sum += values[(row * width) + sample_col] * kernel[static_cast(offset + radius)]; + } + temp[(row * width) + col] = sum; + } + } + + for (std::size_t row = 0U; row < height; ++row) { + for (std::size_t col = 0U; col < width; ++col) { + double sum = 0.0; + for (std::ptrdiff_t offset = -radius; offset <= radius; ++offset) { + const auto sample_row = clamp_index(static_cast(row) + offset, height); + sum += temp[(sample_row * width) + col] * kernel[static_cast(offset + radius)]; + } + output[(row * width) + col] = sum; + } + } + + return output; +} + +[[nodiscard]] auto median_copy(std::vector values) -> double { + if (values.empty()) { + return 0.0; + } + + const auto middle = values.begin() + static_cast(values.size() / 2U); + std::nth_element(values.begin(), middle, values.end()); + double median = *middle; + if ((values.size() % 2U) == 0U) { + const auto lower_middle = values.begin() + static_cast((values.size() / 2U) - 1U); + std::nth_element(values.begin(), lower_middle, values.end()); + median = 0.5 * (median + *lower_middle); + } + return median; +} + +[[nodiscard]] auto find_peak_indices( + const std::vector& values, + std::size_t start_index, + std::size_t stop_index, + double threshold, + std::size_t min_distance +) -> std::vector { + if (stop_index <= start_index + 2U) { + return {}; + } + + std::vector candidates{}; + for (std::size_t index = start_index + 1U; index + 1U < stop_index; ++index) { + if (values[index] < threshold) { + continue; + } + if (values[index] <= values[index - 1U]) { + continue; + } + if (values[index] < values[index + 1U]) { + continue; + } + candidates.push_back(index); + } + + std::sort(candidates.begin(), candidates.end(), [&](const std::size_t left, const std::size_t right) { + return values[left] > values[right]; + }); + + std::vector selected{}; + for (const auto index : candidates) { + bool keep = true; + for (const auto accepted : selected) { + const auto distance = accepted > index ? accepted - index : index - accepted; + if (distance < min_distance) { + keep = false; + break; + } + } + if (keep) { + selected.push_back(index); + } + } + + std::sort(selected.begin(), selected.end()); + return selected; +} + +[[nodiscard]] auto attenuation_at_depth( + std::size_t tx_index, + std::size_t rx_index, + double z_app, + const std::vector& x_tx, + const std::vector& x_rx +) -> double { + const double xc = (x_tx[tx_index] + x_rx[rx_index]) / 2.0; + const double r_tx = std::sqrt(std::pow(xc - x_tx[tx_index], 2.0) + std::pow(z_app, 2.0)); + const double r_rx = std::sqrt(std::pow(xc - x_rx[rx_index], 2.0) + std::pow(z_app, 2.0)); + const double geo = 1.0 / ((r_tx * r_rx) + 1e-12); + const double pat = + std::pow(z_app / (r_tx + 1e-12), 2.0) * std::pow(z_app / (r_rx + 1e-12), 2.0); + return geo * pat + 1e-30; +} + +[[nodiscard]] auto lower_bound_index(const std::vector& axis, double value) -> std::size_t { + return static_cast(std::distance(axis.begin(), std::lower_bound(axis.begin(), axis.end(), value))); +} + +[[nodiscard]] auto build_geometry_selection( + const config::RunConfig& run_config, + const ProcessingLiveConfig& live_config +) -> GeometrySelection { + GeometrySelection selection{}; + + std::unordered_map tx_x_by_pos{}; + for (const auto& entry : run_config.gpr.tx_geometry) { + tx_x_by_pos[entry.output_pos] = static_cast(entry.x_m); + } + + std::unordered_map rx_x_by_pos{}; + for (const auto& entry : run_config.gpr.rx_geometry) { + rx_x_by_pos[entry.input_pos] = static_cast(entry.x_m); + } + + auto select_positions = [](const std::vector& requested, std::vector available) { + std::sort(available.begin(), available.end()); + available.erase(std::unique(available.begin(), available.end()), available.end()); + if (requested.empty()) { + return available; + } + + std::vector result{}; + for (const auto value : requested) { + if (std::find(available.begin(), available.end(), value) != available.end()) { + result.push_back(value); + } + } + std::sort(result.begin(), result.end()); + result.erase(std::unique(result.begin(), result.end()), result.end()); + return result; + }; + + std::vector available_outputs{}; + available_outputs.reserve(tx_x_by_pos.size()); + for (const auto& [position, _] : tx_x_by_pos) { + available_outputs.push_back(position); + } + + std::vector available_inputs{}; + available_inputs.reserve(rx_x_by_pos.size()); + for (const auto& [position, _] : rx_x_by_pos) { + available_inputs.push_back(position); + } + + selection.output_positions = select_positions(live_config.gpr_output_positions, std::move(available_outputs)); + selection.input_positions = select_positions(live_config.gpr_input_positions, std::move(available_inputs)); + + selection.x_tx.reserve(selection.output_positions.size()); + for (std::size_t index = 0U; index < selection.output_positions.size(); ++index) { + const auto position = selection.output_positions[index]; + selection.output_local_by_pos[position] = static_cast(index); + selection.x_tx.push_back(tx_x_by_pos[position]); + } + + selection.x_rx.reserve(selection.input_positions.size()); + for (std::size_t index = 0U; index < selection.input_positions.size(); ++index) { + const auto position = selection.input_positions[index]; + selection.input_local_by_pos[position] = static_cast(index); + selection.x_rx.push_back(rx_x_by_pos[position]); + } + + return selection; +} + +[[nodiscard]] auto build_background_mean( + std::span previous_collections, + const GeometrySelection& selection, + const ProcessingLiveConfig& live_config +) -> std::unordered_map>> { + std::unordered_map>> result{}; + if (!live_config.gpr_background_subtract_enabled || live_config.gpr_background_mean_count == 0U) { + return result; + } + + const std::size_t mean_count = static_cast(live_config.gpr_background_mean_count); + const std::size_t start_index = + previous_collections.size() > mean_count ? previous_collections.size() - mean_count : 0U; + + std::unordered_map accumulators{}; + for (std::size_t collection_index = start_index; collection_index < previous_collections.size(); ++collection_index) { + const auto& collection = previous_collections[collection_index]; + for (const auto& trace : collection.traces) { + const auto output_it = selection.output_local_by_pos.find(trace.combo.output_pos); + const auto input_it = selection.input_local_by_pos.find(trace.combo.input_pos); + if (output_it == selection.output_local_by_pos.end() || input_it == selection.input_local_by_pos.end()) { + continue; + } + + const auto key = make_pair_key(output_it->second, input_it->second); + auto& accumulator = accumulators[key]; + if (accumulator.sum.empty()) { + accumulator.sum.assign(trace.s21.size(), std::complex(0.0, 0.0)); + } + if (accumulator.sum.size() != trace.s21.size()) { + continue; + } + + for (std::size_t sample_index = 0U; sample_index < trace.s21.size(); ++sample_index) { + const auto& sample = trace.s21[sample_index]; + accumulator.sum[sample_index] += std::complex(sample.re, sample.im); + } + accumulator.count += 1U; + } + } + + for (auto& [key, accumulator] : accumulators) { + if (accumulator.count == 0U) { + continue; + } + + auto& mean_trace = result[key]; + mean_trace = std::move(accumulator.sum); + const double inverse_count = 1.0 / static_cast(accumulator.count); + for (auto& sample : mean_trace) { + sample *= inverse_count; + } + } + + return result; +} + +[[nodiscard]] auto collect_selected_traces( + const ipc::PreprocessedCollection& collection, + const GeometrySelection& selection, + const std::unordered_map>>& background_mean +) -> std::unordered_map { + std::unordered_map traces{}; + + for (const auto& trace : collection.traces) { + const auto output_it = selection.output_local_by_pos.find(trace.combo.output_pos); + const auto input_it = selection.input_local_by_pos.find(trace.combo.input_pos); + if (output_it == selection.output_local_by_pos.end() || input_it == selection.input_local_by_pos.end()) { + continue; + } + + SelectedTrace selected{}; + selected.frequency_hz.reserve(trace.frequency_hz.size()); + for (const auto value : trace.frequency_hz) { + selected.frequency_hz.push_back(static_cast(value)); + } + + selected.s21.reserve(trace.s21.size()); + const auto key = make_pair_key(output_it->second, input_it->second); + const auto background_it = background_mean.find(key); + for (std::size_t sample_index = 0U; sample_index < trace.s21.size(); ++sample_index) { + std::complex sample(trace.s21[sample_index].re, trace.s21[sample_index].im); + if (background_it != background_mean.end() && background_it->second.size() == trace.s21.size()) { + sample -= background_it->second[sample_index]; + } + selected.s21.push_back(sample); + } + + traces[key] = std::move(selected); + } + + return traces; +} + +[[nodiscard]] auto compute_ascan( + const SelectedTrace& trace, + double start_hz, + double stop_hz, + double velocity_mps +) -> AscanResult { + AscanResult result{}; + if (trace.frequency_hz.size() != trace.s21.size()) { + return result; + } + + std::vector frequency_hz{}; + std::vector> s21{}; + frequency_hz.reserve(trace.frequency_hz.size()); + s21.reserve(trace.s21.size()); + + const double low_hz = std::min(start_hz, stop_hz); + const double high_hz = std::max(start_hz, stop_hz); + for (std::size_t index = 0U; index < trace.frequency_hz.size(); ++index) { + const double frequency_value = trace.frequency_hz[index]; + if (frequency_value < low_hz || frequency_value > high_hz) { + continue; + } + frequency_hz.push_back(frequency_value); + s21.push_back(trace.s21[index]); + } + + if (frequency_hz.size() < 2U) { + return result; + } + + const std::size_t point_count = frequency_hz.size(); + const double df = (frequency_hz.back() - frequency_hz.front()) / static_cast(point_count - 1U); + if (!(df > 0.0)) { + return result; + } + + const auto start_bin = static_cast(std::llround(frequency_hz.front() / df)); + if (start_bin < 0) { + return result; + } + + const auto start_index = static_cast(start_bin); + const std::size_t min_fft_len = 2U * (start_index + point_count - 1U); + const std::size_t fft_len = next_power_of_two(min_fft_len); + if (fft_len < min_fft_len || start_index > fft_len || point_count > (fft_len - start_index)) { + return result; + } + + std::vector> spectrum(fft_len, std::complex(0.0, 0.0)); + for (std::size_t index = 0U; index < point_count; ++index) { + const double window = point_count > 1U + ? 0.5 - (0.5 * std::cos((2.0 * kPi * static_cast(index)) / static_cast(point_count - 1U))) + : 1.0; + spectrum[start_index + index] = s21[index] * window; + } + + fft_inplace(spectrum, true); + + result.bandwidth_hz = frequency_hz.back() - frequency_hz.front(); + const double dt = 1.0 / (static_cast(fft_len) * df); + result.time_s.resize(fft_len, 0.0); + result.depth_m.resize(fft_len, 0.0); + result.amplitude.resize(fft_len, 0.0); + for (std::size_t index = 0U; index < fft_len; ++index) { + result.time_s[index] = static_cast(index) * dt; + result.depth_m[index] = result.time_s[index] * velocity_mps / 2.0; + result.amplitude[index] = std::abs(spectrum[index]); + } + + return result; +} + +[[nodiscard]] auto build_grid( + const std::vector& x_tx, + const std::vector& x_rx, + double max_depth_m +) -> GridDefinition { + GridDefinition grid{}; + if (x_tx.empty() || x_rx.empty() || !(max_depth_m > kAccumulatorZMinM)) { + return grid; + } + + const auto [tx_min_it, tx_max_it] = std::minmax_element(x_tx.begin(), x_tx.end()); + const auto [rx_min_it, rx_max_it] = std::minmax_element(x_rx.begin(), x_rx.end()); + const double x_min = std::min(*tx_min_it, *rx_min_it) - kAccumulatorXMarginM; + const double x_max = std::max(*tx_max_it, *rx_max_it) + kAccumulatorXMarginM; + + grid.x_grid = build_axis(x_min, x_max, kGridWidth); + grid.z_grid = build_axis(kAccumulatorZMinM, max_depth_m, kGridHeight); + + const std::size_t cell_count = grid.x_grid.size() * grid.z_grid.size(); + grid.tx_distance_grids.assign(x_tx.size(), std::vector(cell_count, 0.0)); + grid.rx_distance_grids.assign(x_rx.size(), std::vector(cell_count, 0.0)); + + for (std::size_t tx_index = 0U; tx_index < x_tx.size(); ++tx_index) { + for (std::size_t row = 0U; row < grid.z_grid.size(); ++row) { + const double z_value = grid.z_grid[row]; + for (std::size_t col = 0U; col < grid.x_grid.size(); ++col) { + const double x_value = grid.x_grid[col]; + const auto cell_index = (row * grid.x_grid.size()) + col; + grid.tx_distance_grids[tx_index][cell_index] = + std::sqrt(std::pow(x_value - x_tx[tx_index], 2.0) + std::pow(z_value, 2.0)); + } + } + } + + for (std::size_t rx_index = 0U; rx_index < x_rx.size(); ++rx_index) { + for (std::size_t row = 0U; row < grid.z_grid.size(); ++row) { + const double z_value = grid.z_grid[row]; + for (std::size_t col = 0U; col < grid.x_grid.size(); ++col) { + const double x_value = grid.x_grid[col]; + const auto cell_index = (row * grid.x_grid.size()) + col; + grid.rx_distance_grids[rx_index][cell_index] = + std::sqrt(std::pow(x_value - x_rx[rx_index], 2.0) + std::pow(z_value, 2.0)); + } + } + } + + return grid; +} + +[[nodiscard]] auto is_excluded(double z_value, const std::vector>& ranges) -> bool { + for (const auto& [low, high] : ranges) { + if (z_value >= low && z_value <= high) { + return true; + } + } + return false; +} + +[[nodiscard]] auto build_accumulator( + const GridDefinition& grid, + const std::unordered_map>& peaks_by_pair, + const std::vector>& exclude_ranges, + double velocity_mps, + double shell_sigma_m, + const std::vector& x_tx, + const std::vector& x_rx +) -> std::vector { + const std::size_t width = grid.x_grid.size(); + const std::size_t height = grid.z_grid.size(); + std::vector accumulator(width * height, 0.0); + if (!(shell_sigma_m > 0.0)) { + return accumulator; + } + + for (std::size_t tx_index = 0U; tx_index < x_tx.size(); ++tx_index) { + for (std::size_t rx_index = 0U; rx_index < x_rx.size(); ++rx_index) { + const auto peak_it = peaks_by_pair.find(make_pair_key(static_cast(tx_index), static_cast(rx_index))); + if (peak_it == peaks_by_pair.end()) { + continue; + } + + const auto& tx_grid = grid.tx_distance_grids[tx_index]; + const auto& rx_grid = grid.rx_distance_grids[rx_index]; + for (const auto& peak : peak_it->second) { + if (is_excluded(peak.z_app, exclude_ranges)) { + continue; + } + + const double range_total = velocity_mps * peak.tau; + for (std::size_t cell_index = 0U; cell_index < accumulator.size(); ++cell_index) { + const double residual = tx_grid[cell_index] + rx_grid[cell_index] - range_total; + const double shell = std::exp(-0.5 * std::pow(residual / shell_sigma_m, 2.0)); + accumulator[cell_index] += shell * peak.snr_comp; + } + } + } + } + + return accumulator; +} + +[[nodiscard]] auto count_agreeing_ellipses( + double x_est, + double z_est, + const std::unordered_map>& peaks_by_pair, + const std::vector>& exclude_ranges, + const std::vector& x_tx, + const std::vector& x_rx, + double velocity_mps, + double shell_sigma_m +) -> double { + std::size_t count = 0U; + for (std::size_t tx_index = 0U; tx_index < x_tx.size(); ++tx_index) { + for (std::size_t rx_index = 0U; rx_index < x_rx.size(); ++rx_index) { + const auto peak_it = peaks_by_pair.find(make_pair_key(static_cast(tx_index), static_cast(rx_index))); + if (peak_it == peaks_by_pair.end()) { + continue; + } + + for (const auto& peak : peak_it->second) { + if (is_excluded(peak.z_app, exclude_ranges)) { + continue; + } + + const double rt = std::sqrt(std::pow(x_est - x_tx[tx_index], 2.0) + std::pow(z_est, 2.0)); + const double rr = std::sqrt(std::pow(x_est - x_rx[rx_index], 2.0) + std::pow(z_est, 2.0)); + if (std::abs((rt + rr) - (velocity_mps * peak.tau)) < shell_sigma_m * 6.0) { + count += 1U; + break; + } + } + } + } + + return static_cast(count); +} + +[[nodiscard]] auto find_centroid( + const std::vector& values, + std::size_t width, + std::size_t height, + std::size_t iz, + std::size_t ix, + std::size_t radius_z, + std::size_t radius_x, + const std::vector& x_grid, + const std::vector& z_grid +) -> std::pair { + if (values.empty()) { + return {x_grid[ix], z_grid[iz]}; + } + + const auto row_start = iz > radius_z ? iz - radius_z : 0U; + const auto row_stop = std::min(height, iz + radius_z + 1U); + const auto col_start = ix > radius_x ? ix - radius_x : 0U; + const auto col_stop = std::min(width, ix + radius_x + 1U); + + double weight_sum = 0.0; + double row_weighted_sum = 0.0; + double col_weighted_sum = 0.0; + for (std::size_t row = row_start; row < row_stop; ++row) { + for (std::size_t col = col_start; col < col_stop; ++col) { + const double weight = values[(row * width) + col]; + weight_sum += weight; + row_weighted_sum += static_cast(row) * weight; + col_weighted_sum += static_cast(col) * weight; + } + } + + if (!(weight_sum > 0.0)) { + return {x_grid[ix], z_grid[iz]}; + } + + const auto centroid_row = clamp_index(static_cast(std::llround(row_weighted_sum / weight_sum)), height); + const auto centroid_col = clamp_index(static_cast(std::llround(col_weighted_sum / weight_sum)), width); + return {x_grid[centroid_col], z_grid[centroid_row]}; +} + +[[nodiscard]] auto flatten_table( + const std::vector>& rows, + std::uint32_t column_count +) -> std::vector { + std::vector values{}; + values.reserve(rows.size() * column_count); + for (const auto& row : rows) { + values.insert(values.end(), row.begin(), row.end()); + } + return values; +} + +[[nodiscard]] auto build_table_payload( + const std::string& processing_name, + const std::vector>& rows, + std::uint32_t column_count +) -> ipc::ResultPayload { + ipc::ResultPayload payload{}; + payload.processing_name = processing_name; + payload.kind = ipc::ResultKind::TableF32; + payload.table_columns = column_count; + payload.table_values = flatten_table(rows, column_count); + return payload; +} + +[[nodiscard]] auto build_image_payload( + const std::string& processing_name, + const std::vector& x_axis, + const std::vector& y_axis, + const std::vector& values +) -> ipc::ResultPayload { + ipc::ResultPayload payload{}; + payload.processing_name = processing_name; + payload.kind = ipc::ResultKind::ImageF32; + payload.image_x_axis.reserve(x_axis.size()); + for (const auto value : x_axis) { + payload.image_x_axis.push_back(static_cast(value)); + } + payload.image_y_axis.reserve(y_axis.size()); + for (const auto value : y_axis) { + payload.image_y_axis.push_back(static_cast(value)); + } + payload.image_values.reserve(values.size()); + for (const auto value : values) { + payload.image_values.push_back(static_cast(value)); + } + return payload; +} + +[[nodiscard]] auto accumulator_max(const std::vector& values) -> double { + if (values.empty()) { + return 0.0; + } + return *std::max_element(values.begin(), values.end()); +} + +[[nodiscard]] auto clean_find_points( + const GridDefinition& grid, + const std::unordered_map>& peaks_by_pair, + const std::vector& x_tx, + const std::vector& x_rx, + double velocity_mps, + double shell_sigma_m +) -> std::pair, std::vector> { + std::vector found{}; + const auto accumulator = build_accumulator(grid, peaks_by_pair, {}, velocity_mps, shell_sigma_m, x_tx, x_rx); + const double initial_max = accumulator_max(accumulator); + if (!(initial_max > 0.0) || grid.x_grid.size() < 2U || grid.z_grid.size() < 2U) { + return {found, gaussian_filter_2d(accumulator, grid.x_grid.size(), grid.z_grid.size(), kGaussianSigma)}; + } + + const double dx = grid.x_grid[1] - grid.x_grid[0]; + const double dz = grid.z_grid[1] - grid.z_grid[0]; + const auto radius_x = static_cast(std::max(1.0, std::round(kCleanSuppressRadiusM / std::max(dx, 1e-6)))); + const auto radius_z = static_cast(std::max(1.0, std::round(kCleanSuppressRadiusM / std::max(dz, 1e-6)))); + + std::vector> excluded_ranges{}; + for (std::size_t step = 0U; step < kMaxObjects; ++step) { + const auto current = build_accumulator(grid, peaks_by_pair, excluded_ranges, velocity_mps, shell_sigma_m, x_tx, x_rx); + const auto smoothed = gaussian_filter_2d(current, grid.x_grid.size(), grid.z_grid.size(), kGaussianSigma); + const double smoothed_max = accumulator_max(smoothed); + if (!(smoothed_max > (kCleanThresholdFrac * initial_max))) { + break; + } + + const auto max_it = std::max_element(smoothed.begin(), smoothed.end()); + const auto max_index = static_cast(std::distance(smoothed.begin(), max_it)); + const auto iz = max_index / grid.x_grid.size(); + const auto ix = max_index % grid.x_grid.size(); + const auto [x_est, z_est] = find_centroid( + smoothed, + grid.x_grid.size(), + grid.z_grid.size(), + iz, + ix, + radius_z, + radius_x, + grid.x_grid, + grid.z_grid + ); + + PointRecord point{}; + point.x_m = x_est; + point.z_m = z_est; + point.score = count_agreeing_ellipses(x_est, z_est, peaks_by_pair, excluded_ranges, x_tx, x_rx, velocity_mps, shell_sigma_m); + found.push_back(point); + + std::vector matched_depths{}; + for (std::size_t tx_index = 0U; tx_index < x_tx.size(); ++tx_index) { + for (std::size_t rx_index = 0U; rx_index < x_rx.size(); ++rx_index) { + const auto peak_it = peaks_by_pair.find(make_pair_key(static_cast(tx_index), static_cast(rx_index))); + if (peak_it == peaks_by_pair.end()) { + continue; + } + for (const auto& peak : peak_it->second) { + if (is_excluded(peak.z_app, excluded_ranges)) { + continue; + } + const double rt = std::sqrt(std::pow(x_est - x_tx[tx_index], 2.0) + std::pow(z_est, 2.0)); + const double rr = std::sqrt(std::pow(x_est - x_rx[rx_index], 2.0) + std::pow(z_est, 2.0)); + if (std::abs((rt + rr) - (velocity_mps * peak.tau)) < shell_sigma_m * 3.0) { + matched_depths.push_back(peak.z_app); + } + } + } + } + + if (!matched_depths.empty()) { + const auto [min_it, max_it_depth] = std::minmax_element(matched_depths.begin(), matched_depths.end()); + excluded_ranges.emplace_back(*min_it - shell_sigma_m, *max_it_depth + shell_sigma_m); + } + } + + return {found, gaussian_filter_2d(accumulator, grid.x_grid.size(), grid.z_grid.size(), kGaussianSigma)}; +} + +[[nodiscard]] auto extended_find_regions( + const GridDefinition& grid, + const std::unordered_map>& peaks_by_pair, + const std::vector& x_tx, + const std::vector& x_rx, + double velocity_mps, + double shell_sigma_m +) -> std::pair, std::vector> { + std::vector regions{}; + const auto accumulator = build_accumulator(grid, peaks_by_pair, {}, velocity_mps, shell_sigma_m, x_tx, x_rx); + const auto smoothed = gaussian_filter_2d(accumulator, grid.x_grid.size(), grid.z_grid.size(), kGaussianSigma); + const double smoothed_max = accumulator_max(smoothed); + if (!(smoothed_max > 0.0) || grid.x_grid.size() < 2U || grid.z_grid.size() < 2U) { + return {regions, smoothed}; + } + + const double dx_cm = (grid.x_grid[1] - grid.x_grid[0]) * 100.0; + const double dz_cm = (grid.z_grid[1] - grid.z_grid[0]) * 100.0; + const auto min_pixels = static_cast(std::max(1.0, std::floor(kExtendedMinAreaCm2 / std::max(dx_cm * dz_cm, 1e-6)))); + + const std::size_t width = grid.x_grid.size(); + const std::size_t height = grid.z_grid.size(); + const double threshold = kExtendedThresholdFrac * smoothed_max; + std::vector visited(width * height, 0U); + + for (std::size_t row = 0U; row < height; ++row) { + for (std::size_t col = 0U; col < width; ++col) { + const auto start_index = (row * width) + col; + if (visited[start_index] != 0U || smoothed[start_index] <= threshold) { + continue; + } + + std::vector stack{start_index}; + std::vector component{}; + visited[start_index] = 1U; + + while (!stack.empty()) { + const auto cell_index = stack.back(); + stack.pop_back(); + component.push_back(cell_index); + + const auto cell_row = cell_index / width; + const auto cell_col = cell_index % width; + const std::pair offsets[] = { + {-1, 0}, + {1, 0}, + {0, -1}, + {0, 1}, + }; + + for (const auto& [row_offset, col_offset] : offsets) { + const auto next_row = static_cast(cell_row) + row_offset; + const auto next_col = static_cast(cell_col) + col_offset; + if (next_row < 0 || next_col < 0) { + continue; + } + if (next_row >= static_cast(height) || next_col >= static_cast(width)) { + continue; + } + const auto next_index = + (static_cast(next_row) * width) + static_cast(next_col); + if (visited[next_index] != 0U || smoothed[next_index] <= threshold) { + continue; + } + visited[next_index] = 1U; + stack.push_back(next_index); + } + } + + if (component.size() < min_pixels) { + continue; + } + + RegionRecord region{}; + region.mask.assign(width * height, 0.0F); + double weight_sum = 0.0; + double x_weight_sum = 0.0; + double z_weight_sum = 0.0; + for (const auto cell_index : component) { + const auto cell_row = cell_index / width; + const auto cell_col = cell_index % width; + const double weight = smoothed[cell_index]; + weight_sum += weight; + x_weight_sum += grid.x_grid[cell_col] * weight; + z_weight_sum += grid.z_grid[cell_row] * weight; + region.mask[cell_index] = 1.0F; + } + + if (!(weight_sum > 0.0)) { + continue; + } + + region.x_m = x_weight_sum / weight_sum; + region.z_m = z_weight_sum / weight_sum; + region.score = count_agreeing_ellipses(region.x_m, region.z_m, peaks_by_pair, {}, x_tx, x_rx, velocity_mps, shell_sigma_m); + region.pixel_count = static_cast(component.size()); + regions.push_back(std::move(region)); + } + } + + return {regions, smoothed}; +} + +} // namespace + +auto GprProcessor::name() const -> std::string { + return "gpr"; +} + +auto GprProcessor::process_collection( + const config::RunConfig& run_config, + const ipc::PreprocessedCollection& collection, + std::span previous_collections, + const ProcessingLiveConfig& live_config +) -> ipc::ResultCollection { + ipc::ResultCollection results{}; + results.collection_id = collection.collection_id; + results.monotonic_ns = collection.monotonic_ns; + + const auto selection = build_geometry_selection(run_config, live_config); + if (selection.input_positions.empty() || selection.output_positions.empty()) { + return results; + } + + const auto background_mean = build_background_mean(previous_collections, selection, live_config); + const auto traces_by_pair = collect_selected_traces(collection, selection, background_mean); + if (traces_by_pair.empty()) { + return results; + } + + const double velocity_mps = + kSpeedOfLightMetersPerSec / std::sqrt(std::max(1e-6, static_cast(run_config.gpr.relative_permittivity))); + const double start_hz = static_cast(live_config.gpr_start_freq_mhz) * 1'000'000.0; + const double stop_hz = static_cast(live_config.gpr_stop_freq_mhz) * 1'000'000.0; + + std::unordered_map ascans_by_pair{}; + double bandwidth_hz = 0.0; + for (const auto& [key, trace] : traces_by_pair) { + auto ascan = compute_ascan(trace, start_hz, stop_hz, velocity_mps); + if (ascan.amplitude.empty() || !(ascan.bandwidth_hz > 0.0)) { + continue; + } + bandwidth_hz = std::max(bandwidth_hz, ascan.bandwidth_hz); + ascans_by_pair.emplace(key, std::move(ascan)); + } + if (ascans_by_pair.empty() || !(bandwidth_hz > 0.0)) { + return results; + } + + const double max_depth_m = static_cast(live_config.gpr_max_depth_m); + const auto grid = build_grid(selection.x_tx, selection.x_rx, max_depth_m); + if (grid.x_grid.empty() || grid.z_grid.empty()) { + return results; + } + + const double shell_sigma_m = velocity_mps / bandwidth_hz * 0.5; + std::unordered_map> peaks_by_pair{}; + for (std::size_t tx_index = 0U; tx_index < selection.x_tx.size(); ++tx_index) { + for (std::size_t rx_index = 0U; rx_index < selection.x_rx.size(); ++rx_index) { + const auto key = make_pair_key(static_cast(tx_index), static_cast(rx_index)); + const auto ascan_it = ascans_by_pair.find(key); + if (ascan_it == ascans_by_pair.end()) { + continue; + } + + const auto& ascan = ascan_it->second; + if (ascan.depth_m.size() < 3U || ascan.amplitude.size() < 3U) { + continue; + } + + const double noise = median_copy(ascan.amplitude); + const auto min_index = lower_bound_index(ascan.depth_m, static_cast(live_config.gpr_min_depth_m)); + const auto max_index = lower_bound_index(ascan.depth_m, static_cast(live_config.gpr_max_depth_m)); + if (max_index <= min_index + 2U) { + continue; + } + + const double z_step = std::max(ascan.depth_m[1] - ascan.depth_m[0], 1e-6); + const std::size_t min_distance = static_cast(std::max( + 4.0, + std::floor(((velocity_mps / (2.0 * bandwidth_hz)) / z_step) * 0.7) + )); + const auto peak_indices = find_peak_indices(ascan.amplitude, min_index, max_index, noise * kSnrThresh, min_distance); + auto& peaks = peaks_by_pair[key]; + peaks.reserve(peak_indices.size()); + for (const auto peak_index : peak_indices) { + const double z_app = ascan.depth_m[peak_index]; + const double snr_raw = ascan.amplitude[peak_index] / std::max(noise, 1e-12); + const double attenuation = attenuation_at_depth(tx_index, rx_index, z_app, selection.x_tx, selection.x_rx); + const double attenuation_norm = + attenuation / attenuation_at_depth(tx_index, rx_index, 2.0, selection.x_tx, selection.x_rx); + const double snr_comp = std::min( + snr_raw / (std::pow(attenuation_norm, static_cast(live_config.gpr_comp_power)) + 1e-12), + kSnrCompMax + ); + peaks.push_back(PeakRecord{ + .z_app = z_app, + .tau = ascan.time_s[peak_index], + .snr_raw = snr_raw, + .snr_comp = snr_comp, + }); + } + } + } + + if (peaks_by_pair.empty()) { + return results; + } + + if (run_config.gpr.mode == "extended") { + const auto [regions, smoothed_accumulator] = extended_find_regions( + grid, + peaks_by_pair, + selection.x_tx, + selection.x_rx, + velocity_mps, + shell_sigma_m + ); + results.collection_payloads.push_back( + build_image_payload("gpr_accumulator", grid.x_grid, grid.z_grid, smoothed_accumulator) + ); + + std::vector> region_rows{}; + region_rows.reserve(regions.size()); + for (std::size_t index = 0U; index < regions.size(); ++index) { + const auto& region = regions[index]; + region_rows.push_back( + { + static_cast(region.x_m), + static_cast(region.z_m), + static_cast(region.score), + static_cast(region.pixel_count), + } + ); + results.collection_payloads.push_back( + build_image_payload( + "gpr_region_mask_" + std::to_string(index), + grid.x_grid, + grid.z_grid, + std::vector(region.mask.begin(), region.mask.end()) + ) + ); + } + results.collection_payloads.push_back(build_table_payload("gpr_region_centers", region_rows, 4U)); + return results; + } + + const auto [points, smoothed_accumulator] = clean_find_points( + grid, + peaks_by_pair, + selection.x_tx, + selection.x_rx, + velocity_mps, + shell_sigma_m + ); + results.collection_payloads.push_back(build_image_payload("gpr_accumulator", grid.x_grid, grid.z_grid, smoothed_accumulator)); + + std::vector> point_rows{}; + point_rows.reserve(points.size()); + for (const auto& point : points) { + point_rows.push_back( + { + static_cast(point.x_m), + static_cast(point.z_m), + static_cast(point.score), + } + ); + } + results.collection_payloads.push_back(build_table_payload("gpr_points", point_rows, 3U)); + return results; +} + +} // namespace radar::processing diff --git a/data_acq_and_processing/processing/processors/src/passthrough_processor.cpp b/data_acq_and_processing/processing/processors/src/passthrough_processor.cpp index cba60ed..d55244a 100644 --- a/data_acq_and_processing/processing/processors/src/passthrough_processor.cpp +++ b/data_acq_and_processing/processing/processors/src/passthrough_processor.cpp @@ -13,27 +13,43 @@ auto PassThroughProcessor::name() const -> std::string { return "pass_through"; } -auto PassThroughProcessor::process(const ipc::SweepTraceBlock& trace, const ProcessingLiveConfig& live_config) - -> ipc::ResultPayload { - ipc::ResultPayload payload{}; - payload.processing_name = name(); - payload.kind = ipc::ResultKind::TraceComplex; - payload.frequency_hz = trace.frequency_hz; - payload.trace = trace.s21; +auto PassThroughProcessor::process_collection( + const config::RunConfig& /*run_config*/, + const ipc::PreprocessedCollection& collection, + std::span /*previous_collections*/, + const ProcessingLiveConfig& live_config +) -> ipc::ResultCollection { + ipc::ResultCollection results{}; + results.collection_id = collection.collection_id; + results.monotonic_ns = collection.monotonic_ns; + results.blocks.reserve(collection.traces.size()); - const float linear_gain = std::pow(10.0F, live_config.gain_db / 20.0F); - const float phase_rad = live_config.phase_deg * (kPi / 180.0F); - const float cos_phase = std::cos(phase_rad); - const float sin_phase = std::sin(phase_rad); + for (const auto& trace : collection.traces) { + ipc::ResultPayload payload{}; + payload.processing_name = name(); + payload.kind = ipc::ResultKind::TraceComplex; + payload.frequency_hz = trace.frequency_hz; + payload.trace = trace.s21; - for (auto& sample : payload.trace) { - const float re = sample.re; - const float im = sample.im; - sample.re = linear_gain * ((re * cos_phase) - (im * sin_phase)); - sample.im = linear_gain * ((re * sin_phase) + (im * cos_phase)); + const float linear_gain = std::pow(10.0F, live_config.gain_db / 20.0F); + const float phase_rad = live_config.phase_deg * (kPi / 180.0F); + const float cos_phase = std::cos(phase_rad); + const float sin_phase = std::sin(phase_rad); + + for (auto& sample : payload.trace) { + const float re = sample.re; + const float im = sample.im; + sample.re = linear_gain * ((re * cos_phase) - (im * sin_phase)); + sample.im = linear_gain * ((re * sin_phase) + (im * cos_phase)); + } + + ipc::ResultBlock block{}; + block.combo = trace.combo; + block.payloads.push_back(std::move(payload)); + results.blocks.push_back(std::move(block)); } - return payload; + return results; } } // namespace radar::processing diff --git a/python_app/gui/app_window.py b/python_app/gui/app_window.py index bdd2a8b..c4a188b 100644 --- a/python_app/gui/app_window.py +++ b/python_app/gui/app_window.py @@ -104,6 +104,18 @@ class AppWindow( self._bscan_depth_axis_by_combo = {} self._bscan_history_floor_collection_id = 0 self._bscan_render_signature = None + self._gpr_lookup_table = None + self._gpr_image_item = None + self._gpr_tx_item = None + self._gpr_rx_item = None + self._gpr_points_item = None + self._gpr_region_centers_item = None + self._gpr_point_labels = [] + self._gpr_region_center_labels = [] + self._gpr_region_mask_items = [] + self._gpr_region_contours = [] + self._gpr_geometry_signature = None + self._gpr_selected_geometry = None self._phase_viewbox = None self._history_run_signature = None self._radar_limits: dict[str, float | int] | None = None diff --git a/python_app/gui/controllers/app_window_config_mixin.py b/python_app/gui/controllers/app_window_config_mixin.py index 36a3991..9a87131 100644 --- a/python_app/gui/controllers/app_window_config_mixin.py +++ b/python_app/gui/controllers/app_window_config_mixin.py @@ -4,7 +4,8 @@ from __future__ import annotations from python_app.gui.runtime.history import remove_last_aligned_histories from python_app.hardware_full.librevna_service import LibreVnaService -from python_app.models.run_config_model import ComboModel, RunConfigModel +from python_app.models.run_config_model import ComboModel, GprRxGeometryModel, GprTxGeometryModel, RunConfigModel +from python_app.models.run_config_validation import validate_gpr_model from python_app.orchestration.config_writer import parse_combos_from_text from python_app.orchestration.live_processing_config import ProcessingLiveConfig from python_app.storage.npz_store import radar_key_from_config @@ -13,6 +14,58 @@ from python_app.storage.npz_store import radar_key_from_config class AppWindowConfigMixin: """Builds runtime config models from current UI state.""" + @staticmethod + def _parse_csv_int_list(text: str) -> list[int]: + """Parse comma-separated integer selection list.""" + cleaned = text.strip() + if not cleaned: + return [] + values: list[int] = [] + for part in cleaned.split(","): + token = part.strip() + if not token: + continue + values.append(int(token)) + return values + + @staticmethod + def _parse_gpr_tx_geometry_text(text: str) -> list[GprTxGeometryModel]: + """Parse line-based Tx geometry editor text.""" + entries: list[GprTxGeometryModel] = [] + for line_number, raw_line in enumerate(text.splitlines(), start=1): + line = raw_line.strip() + if not line: + continue + parts = line.split() + if len(parts) != 2: + raise ValueError(f"Invalid Tx geometry line {line_number}: expected `output_pos x_m`") + entries.append( + GprTxGeometryModel( + output_pos=int(parts[0]), + x_m=float(parts[1]), + ) + ) + return entries + + @staticmethod + def _parse_gpr_rx_geometry_text(text: str) -> list[GprRxGeometryModel]: + """Parse line-based Rx geometry editor text.""" + entries: list[GprRxGeometryModel] = [] + for line_number, raw_line in enumerate(text.splitlines(), start=1): + line = raw_line.strip() + if not line: + continue + parts = line.split() + if len(parts) != 2: + raise ValueError(f"Invalid Rx geometry line {line_number}: expected `input_pos x_m`") + entries.append( + GprRxGeometryModel( + input_pos=int(parts[0]), + x_m=float(parts[1]), + ) + ) + return entries + def _save_current_config(self) -> None: """Persist currently selected GUI settings into root run_config.json.""" try: @@ -68,6 +121,15 @@ class AppWindowConfigMixin: config.preprocess.calibration_set = self._selected_calibration_set config.preprocess.reference_set = self._selected_reference_set + config.gpr.mode = self._gpr_config_mode.currentText() + config.gpr.relative_permittivity = float(self._gpr_relative_permittivity.value()) + config.gpr.tx_geometry = self._parse_gpr_tx_geometry_text(self._gpr_tx_geometry_input.toPlainText()) + config.gpr.rx_geometry = self._parse_gpr_rx_geometry_text(self._gpr_rx_geometry_input.toPlainText()) + validate_gpr_model( + config.gpr, + input_switch_positions=config.input_switch.positions, + output_switch_positions=config.output_switch.positions, + ) return config def _radar_key(self, config: RunConfigModel) -> str: @@ -85,16 +147,31 @@ class AppWindowConfigMixin: def _live_processing_config(self, *, history_command: str = "none") -> ProcessingLiveConfig: """Build live processing config from current processing widgets.""" self._sync_bscan_frequency_limits_with_radar() + self._sync_gpr_frequency_limits_with_radar() + y_min_db = float(self._pass_through_y_min_db.value()) + y_max_db = float(self._pass_through_y_max_db.value()) return ProcessingLiveConfig( processor_mode=self._processing_mode.currentText(), gain_db=float(self._processing_gain_db.value()), phase_deg=float(self._processing_phase_deg.value()), + pass_through_fixed_y_enabled=bool(self._pass_through_fixed_y_enabled.isChecked()), + pass_through_y_min_db=min(y_min_db, y_max_db), + pass_through_y_max_db=max(y_min_db, y_max_db), bscan_axis=self._bscan_axis.currentText(), bscan_cut_m=float(self._bscan_cut_m.value()), bscan_max_depth_m=float(self._bscan_max_depth_m.value()), bscan_gain=float(self._bscan_gain.value()), bscan_start_freq_mhz=float(self._bscan_start_freq_mhz.value()), bscan_stop_freq_mhz=float(self._bscan_stop_freq_mhz.value()), + gpr_input_positions=self._parse_csv_int_list(self._gpr_input_positions_input.text()), + gpr_output_positions=self._parse_csv_int_list(self._gpr_output_positions_input.text()), + gpr_min_depth_m=float(self._gpr_min_depth_m.value()), + gpr_max_depth_m=float(self._gpr_max_depth_m.value()), + gpr_comp_power=float(self._gpr_comp_power.value()), + gpr_start_freq_mhz=float(self._gpr_start_freq_mhz.value()), + gpr_stop_freq_mhz=float(self._gpr_stop_freq_mhz.value()), + gpr_background_subtract_enabled=bool(self._gpr_background_subtract_enabled.isChecked()), + gpr_background_mean_count=int(self._gpr_background_mean_count.value()), history_command_seq=int(self._history_command_seq), history_command=str(history_command), ) @@ -109,10 +186,18 @@ class AppWindowConfigMixin: """Handle live-processing setting changes and trigger redraw when needed.""" try: self._write_live_processing_config() - if self._processing_mode.currentText() == "bscan": + current_mode = self._processing_mode.currentText() + if current_mode == "bscan": self._drain_results_until_quiet(timeout_s=0.25, poll_s=0.01) self._sync_bscan_history_from_results() self._draw_bscan_heatmap_from_history() + elif current_mode == "gpr": + self._drain_results_until_quiet(timeout_s=0.35, poll_s=0.01) + if self._result_history: + if not self._draw_results(self._result_history[-1]): + self._clear_gpr_plot() + else: + self._clear_gpr_plot() elif self._result_history: self._draw_results(self._result_history[-1]) except Exception as exc: # noqa: BLE001 @@ -123,6 +208,7 @@ class AppWindowConfigMixin: mode_to_page = { "pass_through": 0, "bscan": 1, + "gpr": 2, } self._set_plot_mode(mode) self._processing_mode_pages.setCurrentIndex(mode_to_page.get(mode, 0)) @@ -134,17 +220,31 @@ class AppWindowConfigMixin: def _on_bscan_clear_history_clicked(self) -> None: """Permanently clear all runtime histories, ring backlogs, and B-scan cache.""" - self._apply_bscan_history_deletion(remove_last_only=False) + self._apply_history_mode_deletion(mode_label="B-scan", remove_last_only=False) def _on_bscan_remove_last_sweep_clicked(self) -> None: """Permanently delete the latest sweep from runtime histories and rings.""" - self._apply_bscan_history_deletion(remove_last_only=True) + self._apply_history_mode_deletion(mode_label="B-scan", remove_last_only=True) - def _apply_bscan_history_deletion(self, *, remove_last_only: bool) -> None: - """Apply destructive B-scan history deletion via C++ processor history commands.""" + def _on_gpr_clear_history_clicked(self) -> None: + """Permanently clear all runtime histories, ring backlogs, and GPR cache.""" + self._apply_history_mode_deletion(mode_label="GPR", remove_last_only=False) + + def _on_gpr_remove_last_measurement_clicked(self) -> None: + """Permanently delete the latest measurement from runtime histories and rings.""" + self._apply_history_mode_deletion(mode_label="GPR", remove_last_only=True) + + def _clear_history_mode_caches(self) -> None: + """Drop mode-specific cached render state.""" + self._clear_bscan_plot_history() + if hasattr(self, "_gpr_plot"): + self._clear_gpr_plot() + + def _apply_history_mode_deletion(self, *, mode_label: str, remove_last_only: bool) -> None: + """Apply destructive history deletion via C++ processor history commands.""" resume_acquisition = self._supervisor.is_running() history_command = "remove_last" if remove_last_only else "clear_all" - action = "last sweep removed" if remove_last_only else "history fully cleared" + action = "last measurement removed" if remove_last_only else "history fully cleared" dropped_results = 0 try: @@ -169,7 +269,7 @@ class AppWindowConfigMixin: retained_pre=retained_pre, retained_result=retained_result, ) - self._clear_bscan_plot_history() + self._clear_history_mode_caches() self._write_live_processing_config(history_command=history_command, bump_history_seq=True) @@ -180,9 +280,9 @@ class AppWindowConfigMixin: self._redraw_after_history_deletion() if resume_acquisition: self._start_run() - self._log(f"B-scan {action}; dropped pending results={dropped_results}") + self._log(f"{mode_label} {action}; dropped pending results={dropped_results}") except Exception as exc: # noqa: BLE001 - self._show_error(f"Failed to delete B-scan history: {exc}") + self._show_error(f"Failed to delete {mode_label} history: {exc}") def _redraw_after_history_deletion(self) -> None: """Refresh plot immediately after destructive history deletion.""" @@ -192,6 +292,11 @@ class AppWindowConfigMixin: if not self._draw_bscan_heatmap_from_history(): self._bscan_plot.clear() return + if self._processing_mode.currentText() == "gpr": + if self._result_history and self._draw_results(self._result_history[-1]): + return + self._clear_gpr_plot() + return if self._result_history: self._draw_results(self._result_history[-1]) return @@ -208,8 +313,8 @@ class AppWindowConfigMixin: self._on_processing_live_settings_changed() def _on_radar_sweep_limits_changed(self) -> None: - """Clamp B-scan frequency bounds after sweep start/stop edits.""" - if self._sync_bscan_frequency_limits_with_radar(): + """Clamp processing frequency bounds after sweep start/stop edits.""" + if self._sync_processing_frequency_limits_with_radar(): self._on_processing_live_settings_changed() def _refresh_radar_limits_from_device(self) -> bool: @@ -301,7 +406,7 @@ class AppWindowConfigMixin: or prev_power != self._power_input.text().strip() ) - self._sync_bscan_frequency_limits_with_radar() + self._sync_processing_frequency_limits_with_radar() return changed @staticmethod @@ -326,16 +431,16 @@ class AppWindowConfigMixin: widget.setText(str(value)) return value - def _sync_bscan_frequency_limits_with_radar(self) -> bool: - """Synchronize B-scan start/stop MHz widget ranges with radar sweep bounds.""" - required_widgets = ( - "_start_hz_input", - "_stop_hz_input", - "_bscan_start_freq_mhz", - "_bscan_stop_freq_mhz", - ) + def _sync_processing_frequency_limits_with_radar(self) -> bool: + """Synchronize all processing frequency widgets with radar sweep bounds.""" + bscan_changed = self._sync_bscan_frequency_limits_with_radar() + gpr_changed = self._sync_gpr_frequency_limits_with_radar() + return bscan_changed or gpr_changed + + def _sync_frequency_spinboxes_with_radar(self, widget_names: tuple[str, ...]) -> bool: + """Synchronize one or more MHz spin boxes with current radar sweep bounds.""" + required_widgets = ("_start_hz_input", "_stop_hz_input", *widget_names) if not all(hasattr(self, widget_name) for widget_name in required_widgets): - # Processing callbacks can fire while UI groups are still being built. return False try: @@ -348,28 +453,41 @@ class AppWindowConfigMixin: radar_max_mhz = max(radar_start_hz, radar_stop_hz) / 1_000_000.0 changed = False - for widget in (self._bscan_start_freq_mhz, self._bscan_stop_freq_mhz): + widgets = [getattr(self, widget_name) for widget_name in widget_names] + for widget in widgets: if widget.minimum() != radar_min_mhz or widget.maximum() != radar_max_mhz: changed = True widget.blockSignals(True) widget.setRange(radar_min_mhz, radar_max_mhz) widget.blockSignals(False) - clamped_start_mhz = min(max(self._bscan_start_freq_mhz.value(), radar_min_mhz), radar_max_mhz) - clamped_stop_mhz = min(max(self._bscan_stop_freq_mhz.value(), radar_min_mhz), radar_max_mhz) - if clamped_start_mhz != self._bscan_start_freq_mhz.value(): - changed = True - self._bscan_start_freq_mhz.blockSignals(True) - self._bscan_start_freq_mhz.setValue(clamped_start_mhz) - self._bscan_start_freq_mhz.blockSignals(False) - if clamped_stop_mhz != self._bscan_stop_freq_mhz.value(): - changed = True - self._bscan_stop_freq_mhz.blockSignals(True) - self._bscan_stop_freq_mhz.setValue(clamped_stop_mhz) - self._bscan_stop_freq_mhz.blockSignals(False) - + for widget in widgets: + clamped_value = min(max(widget.value(), radar_min_mhz), radar_max_mhz) + if clamped_value != widget.value(): + changed = True + widget.blockSignals(True) + widget.setValue(clamped_value) + widget.blockSignals(False) return changed + def _sync_bscan_frequency_limits_with_radar(self) -> bool: + """Synchronize B-scan start/stop MHz widget ranges with radar sweep bounds.""" + return self._sync_frequency_spinboxes_with_radar( + ( + "_bscan_start_freq_mhz", + "_bscan_stop_freq_mhz", + ) + ) + + def _sync_gpr_frequency_limits_with_radar(self) -> bool: + """Synchronize GPR start/stop MHz widget ranges with radar sweep bounds.""" + return self._sync_frequency_spinboxes_with_radar( + ( + "_gpr_start_freq_mhz", + "_gpr_stop_freq_mhz", + ) + ) + @staticmethod def _switches_are_effectively_static(config: RunConfigModel) -> bool: """Return `True` when switch setup effectively yields one fixed combo.""" diff --git a/python_app/gui/controllers/app_window_pipeline_mixin.py b/python_app/gui/controllers/app_window_pipeline_mixin.py index 5e3a0e1..3628249 100644 --- a/python_app/gui/controllers/app_window_pipeline_mixin.py +++ b/python_app/gui/controllers/app_window_pipeline_mixin.py @@ -353,7 +353,7 @@ class AppWindowPipelineMixin: """Reset runtime history and B-scan caches.""" self._replace_runtime_history(retained_raw=[], retained_pre=[], retained_result=[]) self._bscan_history_floor_collection_id = 0 - self._clear_bscan_plot_history() + self._clear_history_mode_caches() self._update_history_indicator() def _replace_runtime_history( @@ -385,4 +385,8 @@ class AppWindowPipelineMixin: def _validate_processing_mode_constraints(self, config: RunConfigModel) -> None: """Validate processing-mode constraints for run start.""" - validate_processing_mode_constraints(self._processing_mode.currentText(), config) + validate_processing_mode_constraints( + self._processing_mode.currentText(), + config, + self._live_processing_config(), + ) diff --git a/python_app/gui/controllers/app_window_plot_mixin.py b/python_app/gui/controllers/app_window_plot_mixin.py index cda4164..7224736 100644 --- a/python_app/gui/controllers/app_window_plot_mixin.py +++ b/python_app/gui/controllers/app_window_plot_mixin.py @@ -36,6 +36,8 @@ class AppWindowPlotMixin: """Draw collection based on currently selected processing mode.""" if self._processing_mode.currentText() == "bscan": return self._draw_bscan_heatmap(collection) + if self._processing_mode.currentText() == "gpr": + return self._draw_gpr_map(collection) return self._draw_trace_lines(collection) def _show_magnitude_curves(self) -> bool: @@ -46,9 +48,24 @@ class AppWindowPlotMixin: """Return whether phase curves should be rendered.""" return self._show_phase_checkbox.isChecked() + def _pass_through_fixed_y_range(self) -> tuple[bool, float, float]: + """Return normalized magnitude Y-range override for pass-through mode.""" + y_min = float(self._pass_through_y_min_db.value()) + y_max = float(self._pass_through_y_max_db.value()) + return bool(self._pass_through_fixed_y_enabled.isChecked()), min(y_min, y_max), max(y_min, y_max) + + def _configure_pass_through_magnitude_axis(self, plot: pg.PlotWidget) -> None: + """Apply pass-through magnitude-axis autorange or fixed Y window.""" + fixed_y_enabled, y_min, y_max = self._pass_through_fixed_y_range() + view_box = plot.getViewBox() + view_box.invertY(False) + view_box.enableAutoRange(x=True, y=not fixed_y_enabled) + if fixed_y_enabled: + plot.setYRange(y_min, y_max, padding=0.0) + def _on_trace_visibility_changed(self, *_args) -> None: """Redraw pass-through traces when magnitude/phase toggles changed.""" - if self._processing_mode.currentText() == "bscan": + if self._processing_mode.currentText() in {"bscan", "gpr"}: return if self._result_history: self._draw_results(self._result_history[-1]) @@ -98,8 +115,7 @@ class AppWindowPlotMixin: if show_magnitude: mag_item = magnitude_plot.getPlotItem() - magnitude_plot.getViewBox().invertY(False) - magnitude_plot.getViewBox().enableAutoRange(x=True, y=True) + self._configure_pass_through_magnitude_axis(magnitude_plot) mag_item.showAxis("left", show=True) mag_item.showAxis("bottom", show=not show_phase) magnitude_plot.setLabel("left", "Magnitude", units="dB") @@ -219,6 +235,8 @@ class AppWindowPlotMixin: phase_plot.setXRange(x_min, x_max, padding=0.02) if show_phase: phase_plot.setYRange(-180.0, 180.0, padding=0.02) + if show_magnitude: + self._configure_pass_through_magnitude_axis(magnitude_plot) return has_data @staticmethod @@ -481,8 +499,290 @@ class AppWindowPlotMixin: """Hide right axis and clear phase overlay when phase is not rendered.""" self._clear_phase_overlay() + def _clear_gpr_plot(self) -> None: + """Clear latest GPR plot surface.""" + if not hasattr(self, "_gpr_plot"): + return + self._clear_gpr_point_labels() + self._clear_gpr_region_labels() + self._clear_gpr_region_masks() + if self._gpr_image_item is not None: + self._gpr_image_item.hide() + if self._gpr_tx_item is not None: + self._gpr_tx_item.setData(x=[], y=[]) + self._gpr_tx_item.hide() + if self._gpr_rx_item is not None: + self._gpr_rx_item.setData(x=[], y=[]) + self._gpr_rx_item.hide() + if self._gpr_points_item is not None: + self._gpr_points_item.setData(x=[], y=[]) + self._gpr_points_item.hide() + if self._gpr_region_centers_item is not None: + self._gpr_region_centers_item.setData(x=[], y=[]) + self._gpr_region_centers_item.hide() + self._gpr_plot.setTitle(f"GPR {self._gpr_config_mode.currentText()}") + + def _ensure_gpr_plot_items(self) -> None: + """Create persistent GPR plot items once and reuse them on redraw.""" + if self._gpr_image_item is not None: + return + + plot = self._gpr_plot + plot_item = plot.getPlotItem() + plot_item.showAxis("left", show=True) + plot_item.showAxis("bottom", show=True) + plot_item.setClipToView(True) + plot.setLabel("bottom", "X", units="m") + plot.setLabel("left", "Depth", units="m") + + view_box = plot.getViewBox() + view_box.invertY(True) + view_box.enableAutoRange(x=False, y=False) + + if self._gpr_lookup_table is None: + self._gpr_lookup_table = self._build_lut(["#081c15", "#1b4332", "#ffd166", "#f94144"]) + + self._gpr_image_item = pg.ImageItem(axisOrder="row-major") + self._gpr_image_item.setZValue(0) + self._gpr_image_item.hide() + plot.addItem(self._gpr_image_item) + + self._gpr_tx_item = pg.ScatterPlotItem() + self._gpr_tx_item.setZValue(20) + self._gpr_tx_item.hide() + plot.addItem(self._gpr_tx_item) + + self._gpr_rx_item = pg.ScatterPlotItem() + self._gpr_rx_item.setZValue(20) + self._gpr_rx_item.hide() + plot.addItem(self._gpr_rx_item) + + self._gpr_points_item = pg.ScatterPlotItem() + self._gpr_points_item.setZValue(30) + self._gpr_points_item.hide() + plot.addItem(self._gpr_points_item) + + self._gpr_region_centers_item = pg.ScatterPlotItem() + self._gpr_region_centers_item.setZValue(30) + self._gpr_region_centers_item.hide() + plot.addItem(self._gpr_region_centers_item) + + def _clear_gpr_point_labels(self) -> None: + """Remove dynamic point-score labels from GPR plot.""" + for item in self._gpr_point_labels: + try: + self._gpr_plot.removeItem(item) + except Exception: # noqa: BLE001 + pass + self._gpr_point_labels.clear() + + def _clear_gpr_region_labels(self) -> None: + """Remove dynamic region labels from GPR plot.""" + for item in self._gpr_region_center_labels: + try: + self._gpr_plot.removeItem(item) + except Exception: # noqa: BLE001 + pass + self._gpr_region_center_labels.clear() + + def _clear_gpr_region_masks(self) -> None: + """Remove dynamic region contour carriers from GPR plot.""" + for item in self._gpr_region_mask_items: + try: + self._gpr_plot.removeItem(item) + except Exception: # noqa: BLE001 + pass + self._gpr_region_mask_items.clear() + self._gpr_region_contours.clear() + + @staticmethod + def _collection_payload_by_name(collection: ResultCollection, name: str, kind: int | None = None): + """Return first collection payload matching name and optional kind.""" + for payload in collection.collection_payloads: + if payload.processing_name != name: + continue + if kind is not None and int(payload.kind) != int(kind): + continue + return payload + return None + + @staticmethod + def _collection_payloads_by_prefix(collection: ResultCollection, prefix: str, kind: int | None = None): + """Return collection payloads matching processing-name prefix.""" + payloads = [] + for payload in collection.collection_payloads: + if not str(payload.processing_name).startswith(prefix): + continue + if kind is not None and int(payload.kind) != int(kind): + continue + payloads.append(payload) + return payloads + + def _selected_gpr_geometry(self) -> tuple[np.ndarray, np.ndarray]: + """Resolve selected Tx/Rx geometry arrays for current GPR selection.""" + requested_inputs = tuple(self._parse_csv_int_list(self._gpr_input_positions_input.text())) + requested_outputs = tuple(self._parse_csv_int_list(self._gpr_output_positions_input.text())) + signature = ( + self._gpr_tx_geometry_input.toPlainText(), + self._gpr_rx_geometry_input.toPlainText(), + requested_inputs, + requested_outputs, + ) + if signature == self._gpr_geometry_signature and self._gpr_selected_geometry is not None: + return self._gpr_selected_geometry + + tx_entries = self._parse_gpr_tx_geometry_text(signature[0]) + rx_entries = self._parse_gpr_rx_geometry_text(signature[1]) + requested_input_set = set(requested_inputs) + requested_output_set = set(requested_outputs) + + rx_entries = sorted(rx_entries, key=lambda entry: int(entry.input_pos)) + tx_entries = sorted(tx_entries, key=lambda entry: int(entry.output_pos)) + if requested_input_set: + rx_entries = [entry for entry in rx_entries if int(entry.input_pos) in requested_input_set] + if requested_output_set: + tx_entries = [entry for entry in tx_entries if int(entry.output_pos) in requested_output_set] + + x_tx = np.asarray([float(entry.x_m) for entry in tx_entries], dtype=np.float32) + x_rx = np.asarray([float(entry.x_m) for entry in rx_entries], dtype=np.float32) + self._gpr_geometry_signature = signature + self._gpr_selected_geometry = (x_tx, x_rx) + return self._gpr_selected_geometry + + def _draw_gpr_map(self, collection: ResultCollection) -> bool: + """Draw latest collection-level GPR accumulator and annotations.""" + accumulator_payload = self._collection_payload_by_name(collection, "gpr_accumulator", kind=3) + if accumulator_payload is None: + self._clear_gpr_plot() + return False + + image = np.asarray(accumulator_payload.image, dtype=np.float32) + x_axis = np.asarray(accumulator_payload.image_x_axis, dtype=np.float32) + y_axis = np.asarray(accumulator_payload.image_y_axis, dtype=np.float32) + if image.ndim != 2 or image.size == 0 or x_axis.size == 0 or y_axis.size == 0: + self._clear_gpr_plot() + return False + + x_min = float(x_axis[0]) + x_max = float(x_axis[-1]) + y_min = float(y_axis[0]) + y_max = float(y_axis[-1]) + rect = QRectF(x_min, y_min, max(x_max - x_min, 1e-6), max(y_max - y_min, 1e-6)) + + plot = self._gpr_plot + plot.setUpdatesEnabled(False) + try: + self._ensure_gpr_plot_items() + self._clear_gpr_point_labels() + self._clear_gpr_region_labels() + self._clear_gpr_region_masks() + + self._gpr_image_item.setImage(image, autoLevels=False) + self._gpr_image_item.setRect(rect) + self._gpr_image_item.setLookupTable(self._gpr_lookup_table) + self._gpr_image_item.setLevels((float(np.min(image)), float(np.max(image) + 1e-6))) + self._gpr_image_item.show() + + plot.setXRange(x_min, x_max, padding=0.02) + plot.setYRange(y_min, y_max, padding=0.02) + + x_tx, x_rx = self._selected_gpr_geometry() + if x_tx.size > 0: + self._gpr_tx_item.setData( + x=x_tx, + y=np.zeros_like(x_tx), + symbol="t", + size=13, + brush=pg.mkBrush("#ff595e"), + pen=pg.mkPen("#ffca3a", width=1.0), + ) + self._gpr_tx_item.show() + else: + self._gpr_tx_item.setData(x=[], y=[]) + self._gpr_tx_item.hide() + + if x_rx.size > 0: + self._gpr_rx_item.setData( + x=x_rx, + y=np.zeros_like(x_rx), + symbol="t1", + size=13, + brush=pg.mkBrush("#4cc9f0"), + pen=pg.mkPen("#e0fbfc", width=1.0), + ) + self._gpr_rx_item.show() + else: + self._gpr_rx_item.setData(x=[], y=[]) + self._gpr_rx_item.hide() + + points_payload = self._collection_payload_by_name(collection, "gpr_points", kind=4) + if points_payload is not None and np.asarray(points_payload.table).size > 0: + points = np.asarray(points_payload.table, dtype=np.float32) + self._gpr_points_item.setData( + x=points[:, 0], + y=points[:, 1], + symbol="d", + size=11, + brush=pg.mkBrush("#ffffff"), + pen=pg.mkPen("#111111", width=1.1), + ) + self._gpr_points_item.show() + for x_value, y_value, score in points: + label = pg.TextItem(text=f"{float(score):.0f}", color="#ffffff", anchor=(0.0, 1.0)) + label.setZValue(40) + label.setPos(float(x_value), float(y_value)) + plot.addItem(label) + self._gpr_point_labels.append(label) + else: + self._gpr_points_item.setData(x=[], y=[]) + self._gpr_points_item.hide() + + region_centers_payload = self._collection_payload_by_name(collection, "gpr_region_centers", kind=4) + if region_centers_payload is not None and np.asarray(region_centers_payload.table).size > 0: + centers = np.asarray(region_centers_payload.table, dtype=np.float32) + self._gpr_region_centers_item.setData( + x=centers[:, 0], + y=centers[:, 1], + symbol="o", + size=10, + brush=pg.mkBrush("#80ed99"), + pen=pg.mkPen("#081c15", width=1.1), + ) + self._gpr_region_centers_item.show() + for row in centers: + label = pg.TextItem(text=f"{float(row[2]):.0f}", color="#d8f3dc", anchor=(0.0, 1.0)) + label.setZValue(40) + label.setPos(float(row[0]), float(row[1])) + plot.addItem(label) + self._gpr_region_center_labels.append(label) + else: + self._gpr_region_centers_item.setData(x=[], y=[]) + self._gpr_region_centers_item.hide() + + for payload in self._collection_payloads_by_prefix(collection, "gpr_region_mask_", kind=3): + mask = np.asarray(payload.image, dtype=np.float32) + if mask.ndim != 2 or mask.size == 0: + continue + mask_image = pg.ImageItem(axisOrder="row-major") + mask_image.setZValue(5) + mask_image.setImage(mask, autoLevels=False) + mask_image.setRect(rect) + mask_image.setOpacity(0.0) + plot.addItem(mask_image) + contour = pg.IsocurveItem(data=mask, level=0.5, pen=pg.mkPen("#4cc9f0", width=1.3)) + contour.setParentItem(mask_image) + self._gpr_region_mask_items.append(mask_image) + self._gpr_region_contours.append(contour) + + plot.setTitle(f"GPR {self._gpr_config_mode.currentText()}") + finally: + plot.setUpdatesEnabled(True) + return True + def _result_collection_has_trace(self, collection: ResultCollection) -> bool: """Return `True` when collection contains at least one trace payload.""" + if collection.collection_payloads: + return True for block in collection.blocks: for payload in block.payloads: if payload.kind == 1 and payload.trace.size > 0: @@ -503,8 +803,7 @@ class AppWindowPlotMixin: return if show_magnitude: - magnitude_plot.getViewBox().invertY(False) - magnitude_plot.getViewBox().enableAutoRange(x=True, y=True) + self._configure_pass_through_magnitude_axis(magnitude_plot) magnitude_plot.getPlotItem().showAxis("bottom", show=not show_phase) magnitude_plot.setLabel("left", "Magnitude", units="dB") magnitude_plot.setTitle(title) @@ -548,5 +847,6 @@ class AppWindowPlotMixin: x_max = float(np.max(trace.frequency_hz)) if show_magnitude: magnitude_plot.setXRange(x_min, x_max, padding=0.02) + self._configure_pass_through_magnitude_axis(magnitude_plot) if show_phase: phase_plot.setXRange(x_min, x_max, padding=0.02) diff --git a/python_app/gui/controllers/app_window_snapshot_mixin.py b/python_app/gui/controllers/app_window_snapshot_mixin.py index 85d8b65..c6038d8 100644 --- a/python_app/gui/controllers/app_window_snapshot_mixin.py +++ b/python_app/gui/controllers/app_window_snapshot_mixin.py @@ -111,7 +111,7 @@ class AppWindowSnapshotMixin: self._replace_runtime_history(retained_raw=[], retained_pre=[], retained_result=[]) self._bscan_history_floor_collection_id = 0 - self._clear_bscan_plot_history() + self._clear_history_mode_caches() # Clear processor-side replay cache so newly rendered B-scan starts clean. self._write_live_processing_config(history_command="clear_all", bump_history_seq=True) diff --git a/python_app/gui/controllers/app_window_ui_mixin.py b/python_app/gui/controllers/app_window_ui_mixin.py index 720ea7f..f500fa3 100644 --- a/python_app/gui/controllers/app_window_ui_mixin.py +++ b/python_app/gui/controllers/app_window_ui_mixin.py @@ -2,6 +2,7 @@ Layout is intentionally split into two independent plot surfaces: - single `PlotWidget` for B-scan heatmap rendering; +- single `PlotWidget` for GPR accumulator/annotation rendering; - stacked magnitude/phase `PlotWidget`s for pass-through traces. `_set_plot_mode()` switches between these surfaces via `QStackedWidget`. @@ -26,6 +27,7 @@ import pyqtgraph as pg from python_app.gui.controllers.sections import ( build_data_actions_group, + build_gpr_config_group, build_hardware_actions_group, build_pipeline_group, build_preprocess_summary_group, @@ -72,6 +74,7 @@ class AppWindowUiMixin: # We create both upfront and only switch active page at runtime. self._plot_stack = QStackedWidget(root) self._build_bscan_plot_page() + self._build_gpr_plot_page() self._build_trace_plot_page() # Default view on startup is pass-through traces. @@ -121,6 +124,12 @@ class AppWindowUiMixin: self._plot_stack.addWidget(self._trace_plots_container) + def _build_gpr_plot_page(self) -> None: + """Create GPR page in plot stack.""" + self._gpr_plot = pg.PlotWidget(background="#0f141c") + self._gpr_plot.showGrid(x=True, y=True, alpha=0.2) + self._plot_stack.addWidget(self._gpr_plot) + def _build_settings_toggle(self, root_layout: QHBoxLayout) -> None: """Create narrow button used to collapse or show settings panel.""" self._settings_toggle_button = QPushButton("<") @@ -182,6 +191,7 @@ class AppWindowUiMixin: build_data_actions_group(self), build_preprocess_summary_group(self), build_processing_group(self), + build_gpr_config_group(self), build_radar_group(self), build_switch_group(self), ] @@ -205,12 +215,16 @@ class AppWindowUiMixin: def _set_plot_mode(self, mode: str) -> None: """Switch visible plot page according to processing mode. - `bscan` -> show `self._bscan_plot` (single heatmap surface) - otherwise -> show `self._trace_plots_container` (magnitude + phase) + `bscan` -> show `self._bscan_plot` + `gpr` -> show `self._gpr_plot` + otherwise -> show `self._trace_plots_container` """ if mode == "bscan": self._plot_stack.setCurrentWidget(self._bscan_plot) return + if mode == "gpr": + self._plot_stack.setCurrentWidget(self._gpr_plot) + return self._plot_stack.setCurrentWidget(self._trace_plots_container) @staticmethod diff --git a/python_app/gui/controllers/sections/__init__.py b/python_app/gui/controllers/sections/__init__.py index 62d0256..72bbf3f 100644 --- a/python_app/gui/controllers/sections/__init__.py +++ b/python_app/gui/controllers/sections/__init__.py @@ -1,6 +1,7 @@ """Composable UI section builders used by AppWindow UI mixin.""" from python_app.gui.controllers.sections.data_actions_section import build_data_actions_group +from python_app.gui.controllers.sections.gpr_config_section import build_gpr_config_group from python_app.gui.controllers.sections.hardware_actions_section import build_hardware_actions_group from python_app.gui.controllers.sections.pipeline_section import build_pipeline_group from python_app.gui.controllers.sections.preprocess_summary_section import build_preprocess_summary_group @@ -10,6 +11,7 @@ from python_app.gui.controllers.sections.switch_section import build_switch_grou __all__ = [ "build_data_actions_group", + "build_gpr_config_group", "build_hardware_actions_group", "build_pipeline_group", "build_preprocess_summary_group", diff --git a/python_app/gui/controllers/sections/gpr_config_section.py b/python_app/gui/controllers/sections/gpr_config_section.py new file mode 100644 index 0000000..56d6b7a --- /dev/null +++ b/python_app/gui/controllers/sections/gpr_config_section.py @@ -0,0 +1,53 @@ +"""Builder for stable GPR configuration section.""" + +from __future__ import annotations + +from PyQt6.QtWidgets import QComboBox, QDoubleSpinBox, QFormLayout, QGroupBox, QPlainTextEdit + + +def _format_tx_geometry(owner) -> str: + """Render Tx geometry defaults into editable line-based text.""" + return "\n".join( + f"{int(entry.output_pos)} {float(entry.x_m):g}" + for entry in owner._defaults_config.gpr.tx_geometry + ) + + +def _format_rx_geometry(owner) -> str: + """Render Rx geometry defaults into editable line-based text.""" + return "\n".join( + f"{int(entry.input_pos)} {float(entry.x_m):g}" + for entry in owner._defaults_config.gpr.rx_geometry + ) + + +def build_gpr_config_group(owner) -> QGroupBox: + """Create stable GPR config controls backed by run_config.json.""" + group = QGroupBox("GPR Config") + form = QFormLayout(group) + form.setFieldGrowthPolicy(QFormLayout.FieldGrowthPolicy.AllNonFixedFieldsGrow) + defaults = owner._defaults_config.gpr + + owner._gpr_config_mode = QComboBox() + owner._gpr_config_mode.addItems(["point", "extended"]) + owner._set_combo_current_text(owner._gpr_config_mode, defaults.mode) + + owner._gpr_relative_permittivity = QDoubleSpinBox() + owner._gpr_relative_permittivity.setDecimals(4) + owner._gpr_relative_permittivity.setRange(0.0001, 1000.0) + owner._gpr_relative_permittivity.setSingleStep(0.05) + owner._gpr_relative_permittivity.setValue(float(defaults.relative_permittivity)) + + owner._gpr_tx_geometry_input = QPlainTextEdit(_format_tx_geometry(owner)) + owner._gpr_tx_geometry_input.setPlaceholderText("output_pos x_m") + owner._gpr_tx_geometry_input.setMinimumHeight(88) + + owner._gpr_rx_geometry_input = QPlainTextEdit(_format_rx_geometry(owner)) + owner._gpr_rx_geometry_input.setPlaceholderText("input_pos x_m") + owner._gpr_rx_geometry_input.setMinimumHeight(120) + + form.addRow("Mode", owner._gpr_config_mode) + form.addRow("Relative Permittivity", owner._gpr_relative_permittivity) + form.addRow("Tx Geometry", owner._gpr_tx_geometry_input) + form.addRow("Rx Geometry", owner._gpr_rx_geometry_input) + return group diff --git a/python_app/gui/controllers/sections/processing_section.py b/python_app/gui/controllers/sections/processing_section.py index 984d94b..4ae13c9 100644 --- a/python_app/gui/controllers/sections/processing_section.py +++ b/python_app/gui/controllers/sections/processing_section.py @@ -9,21 +9,39 @@ from PyQt6.QtWidgets import ( QFormLayout, QGroupBox, QHBoxLayout, + QLineEdit, QPushButton, QSizePolicy, + QSpinBox, QStackedWidget, QWidget, ) +def _default_gpr_input_positions(owner) -> str: + """Build default live input-position selection from stable GPR config.""" + geometry_values = {int(entry.input_pos) for entry in owner._defaults_config.gpr.rx_geometry} + combo_values = {int(combo.input) for combo in owner._defaults_config.combos} + values = sorted(geometry_values & combo_values) or sorted(geometry_values) + return ",".join(str(value) for value in values) + + +def _default_gpr_output_positions(owner) -> str: + """Build default live output-position selection from stable GPR config.""" + geometry_values = {int(entry.output_pos) for entry in owner._defaults_config.gpr.tx_geometry} + combo_values = {int(combo.output) for combo in owner._defaults_config.combos} + values = sorted(geometry_values & combo_values) or sorted(geometry_values) + return ",".join(str(value) for value in values) + + def build_processing_group(owner) -> QGroupBox: - """Create processing mode section with pass-through and B-scan pages.""" + """Create processing mode section with pass-through, B-scan, and GPR pages.""" group = QGroupBox("Processing") form = QFormLayout(group) form.setFieldGrowthPolicy(QFormLayout.FieldGrowthPolicy.AllNonFixedFieldsGrow) owner._processing_mode = QComboBox() - owner._processing_mode.addItems(["pass_through", "bscan"]) + owner._processing_mode.addItems(["pass_through", "bscan", "gpr"]) owner._processing_mode_pages = QStackedWidget(group) owner._processing_mode_pages.setSizePolicy(QSizePolicy.Policy.Preferred, QSizePolicy.Policy.Fixed) @@ -51,10 +69,35 @@ def build_processing_group(owner) -> QGroupBox: owner._show_phase_checkbox = QCheckBox("Show phase") owner._show_phase_checkbox.setChecked(True) + owner._pass_through_fixed_y_enabled = QCheckBox("Fix magnitude Y range") + owner._pass_through_fixed_y_enabled.setChecked(False) + + owner._pass_through_y_min_db = QDoubleSpinBox() + owner._pass_through_y_min_db.setDecimals(1) + owner._pass_through_y_min_db.setRange(-240.0, 240.0) + owner._pass_through_y_min_db.setSingleStep(1.0) + owner._pass_through_y_min_db.setValue(-100.0) + + owner._pass_through_y_max_db = QDoubleSpinBox() + owner._pass_through_y_max_db.setDecimals(1) + owner._pass_through_y_max_db.setRange(-240.0, 240.0) + owner._pass_through_y_max_db.setSingleStep(1.0) + owner._pass_through_y_max_db.setValue(0.0) + + def sync_pass_through_y_controls() -> None: + enabled = owner._pass_through_fixed_y_enabled.isChecked() + owner._pass_through_y_min_db.setEnabled(enabled) + owner._pass_through_y_max_db.setEnabled(enabled) + + sync_pass_through_y_controls() + pass_through_form.addRow("Gain dB (live)", owner._processing_gain_db) pass_through_form.addRow("Phase deg (live)", owner._processing_phase_deg) pass_through_form.addRow(owner._show_magnitude_checkbox) pass_through_form.addRow(owner._show_phase_checkbox) + pass_through_form.addRow(owner._pass_through_fixed_y_enabled) + pass_through_form.addRow("Y min dB", owner._pass_through_y_min_db) + pass_through_form.addRow("Y max dB", owner._pass_through_y_max_db) owner._processing_mode_pages.addWidget(pass_through_page) bscan_page = QWidget(owner._processing_mode_pages) @@ -115,11 +158,86 @@ def build_processing_group(owner) -> QGroupBox: bscan_form.addRow(bscan_actions) owner._processing_mode_pages.addWidget(bscan_page) + gpr_page = QWidget(owner._processing_mode_pages) + gpr_page.setSizePolicy(QSizePolicy.Policy.Preferred, QSizePolicy.Policy.Fixed) + gpr_form = QFormLayout(gpr_page) + gpr_form.setFieldGrowthPolicy(QFormLayout.FieldGrowthPolicy.AllNonFixedFieldsGrow) + + owner._gpr_input_positions_input = QLineEdit(_default_gpr_input_positions(owner)) + owner._gpr_input_positions_input.setPlaceholderText("0,1,2") + + owner._gpr_output_positions_input = QLineEdit(_default_gpr_output_positions(owner)) + owner._gpr_output_positions_input.setPlaceholderText("0,1") + + owner._gpr_min_depth_m = QDoubleSpinBox() + owner._gpr_min_depth_m.setDecimals(2) + owner._gpr_min_depth_m.setRange(0.0, 50.0) + owner._gpr_min_depth_m.setSingleStep(0.1) + owner._gpr_min_depth_m.setValue(2.0) + + owner._gpr_max_depth_m = QDoubleSpinBox() + owner._gpr_max_depth_m.setDecimals(2) + owner._gpr_max_depth_m.setRange(0.1, 50.0) + owner._gpr_max_depth_m.setSingleStep(0.1) + owner._gpr_max_depth_m.setValue(14.0) + + owner._gpr_comp_power = QDoubleSpinBox() + owner._gpr_comp_power.setDecimals(3) + owner._gpr_comp_power.setRange(0.0, 5.0) + owner._gpr_comp_power.setSingleStep(0.05) + owner._gpr_comp_power.setValue(0.2) + + owner._gpr_start_freq_mhz = QDoubleSpinBox() + owner._gpr_start_freq_mhz.setDecimals(1) + owner._gpr_start_freq_mhz.setRange(100.0, 8800.0) + owner._gpr_start_freq_mhz.setSingleStep(10.0) + owner._gpr_start_freq_mhz.setValue(3000.0) + + owner._gpr_stop_freq_mhz = QDoubleSpinBox() + owner._gpr_stop_freq_mhz.setDecimals(1) + owner._gpr_stop_freq_mhz.setRange(100.0, 8800.0) + owner._gpr_stop_freq_mhz.setSingleStep(10.0) + owner._gpr_stop_freq_mhz.setValue(6000.0) + + owner._gpr_background_subtract_enabled = QCheckBox("Subtract mean of previous collections") + owner._gpr_background_subtract_enabled.setChecked(True) + + owner._gpr_background_mean_count = QSpinBox() + owner._gpr_background_mean_count.setRange(0, 10_000) + owner._gpr_background_mean_count.setValue(10) + + owner._gpr_clear_history_button = QPushButton("Clear GPR History") + owner._gpr_clear_history_button.clicked.connect(owner._on_gpr_clear_history_clicked) + owner._gpr_remove_last_button = QPushButton("Remove Last Measurement") + owner._gpr_remove_last_button.clicked.connect(owner._on_gpr_remove_last_measurement_clicked) + gpr_actions = QWidget(owner._processing_mode_pages) + gpr_actions_layout = QHBoxLayout(gpr_actions) + gpr_actions_layout.setContentsMargins(0, 0, 0, 0) + gpr_actions_layout.setSpacing(8) + gpr_actions_layout.addWidget(owner._gpr_remove_last_button) + gpr_actions_layout.addWidget(owner._gpr_clear_history_button) + + gpr_form.addRow("Input positions", owner._gpr_input_positions_input) + gpr_form.addRow("Output positions", owner._gpr_output_positions_input) + gpr_form.addRow("Min depth m", owner._gpr_min_depth_m) + gpr_form.addRow("Max depth m", owner._gpr_max_depth_m) + gpr_form.addRow("Comp power", owner._gpr_comp_power) + gpr_form.addRow("Start MHz", owner._gpr_start_freq_mhz) + gpr_form.addRow("Stop MHz", owner._gpr_stop_freq_mhz) + gpr_form.addRow(owner._gpr_background_subtract_enabled) + gpr_form.addRow("Mean count", owner._gpr_background_mean_count) + gpr_form.addRow(gpr_actions) + owner._processing_mode_pages.addWidget(gpr_page) + owner._processing_mode.currentTextChanged.connect(owner._on_processing_mode_changed) owner._processing_gain_db.valueChanged.connect(owner._on_processing_live_settings_changed) owner._processing_phase_deg.valueChanged.connect(owner._on_processing_live_settings_changed) owner._show_magnitude_checkbox.toggled.connect(owner._on_trace_visibility_changed) owner._show_phase_checkbox.toggled.connect(owner._on_trace_visibility_changed) + owner._pass_through_fixed_y_enabled.toggled.connect(sync_pass_through_y_controls) + owner._pass_through_fixed_y_enabled.toggled.connect(owner._on_processing_live_settings_changed) + owner._pass_through_y_min_db.valueChanged.connect(owner._on_processing_live_settings_changed) + owner._pass_through_y_max_db.valueChanged.connect(owner._on_processing_live_settings_changed) form.addRow("Mode", owner._processing_mode) form.addRow(owner._processing_mode_pages) @@ -129,6 +247,15 @@ def build_processing_group(owner) -> QGroupBox: owner._bscan_gain.valueChanged.connect(owner._on_processing_live_settings_changed) owner._bscan_start_freq_mhz.valueChanged.connect(owner._on_processing_live_settings_changed) owner._bscan_stop_freq_mhz.valueChanged.connect(owner._on_processing_live_settings_changed) + owner._gpr_input_positions_input.editingFinished.connect(owner._on_processing_live_settings_changed) + owner._gpr_output_positions_input.editingFinished.connect(owner._on_processing_live_settings_changed) + owner._gpr_min_depth_m.valueChanged.connect(owner._on_processing_live_settings_changed) + owner._gpr_max_depth_m.valueChanged.connect(owner._on_processing_live_settings_changed) + owner._gpr_comp_power.valueChanged.connect(owner._on_processing_live_settings_changed) + owner._gpr_start_freq_mhz.valueChanged.connect(owner._on_processing_live_settings_changed) + owner._gpr_stop_freq_mhz.valueChanged.connect(owner._on_processing_live_settings_changed) + owner._gpr_background_subtract_enabled.toggled.connect(owner._on_processing_live_settings_changed) + owner._gpr_background_mean_count.valueChanged.connect(owner._on_processing_live_settings_changed) owner._on_processing_mode_changed(owner._processing_mode.currentText()) return group diff --git a/python_app/gui/runtime/constraints.py b/python_app/gui/runtime/constraints.py index 9ef7117..13b5995 100644 --- a/python_app/gui/runtime/constraints.py +++ b/python_app/gui/runtime/constraints.py @@ -3,19 +3,53 @@ from __future__ import annotations from python_app.models.run_config_model import RunConfigModel +from python_app.orchestration.live_processing_config import ProcessingLiveConfig -def validate_processing_mode_constraints(processing_mode: str, config: RunConfigModel) -> None: +def validate_processing_mode_constraints( + processing_mode: str, + config: RunConfigModel, + live_config: ProcessingLiveConfig, +) -> None: """Validate mode-specific constraints for current run configuration.""" - if processing_mode != "bscan": + if processing_mode == "bscan": + any_native_switch = config.input_switch.driver_mode == "native" or config.output_switch.driver_mode == "native" + if not any_native_switch: + return + + combo_count = len({(int(combo.input), int(combo.output)) for combo in config.combos}) + if combo_count != 1: + raise RuntimeError( + f"B-scan with native switches requires exactly one run combo (now {combo_count})" + ) return - any_native_switch = config.input_switch.driver_mode == "native" or config.output_switch.driver_mode == "native" - if not any_native_switch: + if processing_mode != "gpr": return - combo_count = len({(int(combo.input), int(combo.output)) for combo in config.combos}) - if combo_count != 1: + available_inputs = sorted({int(entry.input_pos) for entry in config.gpr.rx_geometry}) + available_outputs = sorted({int(entry.output_pos) for entry in config.gpr.tx_geometry}) + if not available_inputs or not available_outputs: + raise RuntimeError("GPR requires non-empty Tx/Rx geometry in run_config") + + requested_inputs = sorted({int(value) for value in live_config.gpr_input_positions}) + requested_outputs = sorted({int(value) for value in live_config.gpr_output_positions}) + + selected_inputs = requested_inputs or available_inputs + selected_outputs = requested_outputs or available_outputs + + missing_inputs = [value for value in selected_inputs if value not in available_inputs] + missing_outputs = [value for value in selected_outputs if value not in available_outputs] + if missing_inputs: + raise RuntimeError(f"GPR input positions are missing from geometry config: {missing_inputs}") + if missing_outputs: + raise RuntimeError(f"GPR output positions are missing from geometry config: {missing_outputs}") + + required_combos = {(int(input_pos), int(output_pos)) for input_pos in selected_inputs for output_pos in selected_outputs} + configured_combos = {(int(combo.input), int(combo.output)) for combo in config.combos} + missing_combos = sorted(required_combos - configured_combos) + if missing_combos: raise RuntimeError( - f"B-scan with native switches requires exactly one run combo (now {combo_count})" + "GPR run combos do not cover selected input/output positions: " + f"{missing_combos}" ) diff --git a/python_app/models/dataset_model.py b/python_app/models/dataset_model.py index 0d1ebc7..3e7fe24 100644 --- a/python_app/models/dataset_model.py +++ b/python_app/models/dataset_model.py @@ -7,6 +7,21 @@ from dataclasses import dataclass, field import numpy as np +def _empty_f32_array() -> np.ndarray: + """Return empty float32 array used by payload defaults.""" + return np.array([], dtype=np.float32) + + +def _empty_c64_array() -> np.ndarray: + """Return empty complex64 array used by payload defaults.""" + return np.array([], dtype=np.complex64) + + +def _empty_f32_matrix() -> np.ndarray: + """Return empty 2D float32 matrix used by payload defaults.""" + return np.zeros((0, 0), dtype=np.float32) + + @dataclass(frozen=True, slots=True) class ComboKey: """Switch combination key: input position + output position.""" @@ -39,9 +54,13 @@ class ResultPayload: processing_name: str kind: int - frequency_hz: np.ndarray - trace: np.ndarray + frequency_hz: np.ndarray = field(default_factory=_empty_f32_array) + trace: np.ndarray = field(default_factory=_empty_c64_array) scalar_value: float = 0.0 + image_x_axis: np.ndarray = field(default_factory=_empty_f32_array) + image_y_axis: np.ndarray = field(default_factory=_empty_f32_array) + image: np.ndarray = field(default_factory=_empty_f32_matrix) + table: np.ndarray = field(default_factory=_empty_f32_matrix) @dataclass(slots=True) @@ -58,4 +77,5 @@ class ResultCollection: collection_id: int monotonic_ns: int + collection_payloads: list[ResultPayload] = field(default_factory=list) blocks: list[ResultBlock] = field(default_factory=list) diff --git a/python_app/models/run_config_codec.py b/python_app/models/run_config_codec.py index faf80df..c423d90 100644 --- a/python_app/models/run_config_codec.py +++ b/python_app/models/run_config_codec.py @@ -4,8 +4,13 @@ from __future__ import annotations from typing import Any -from python_app.models.run_config_schema import ComboModel, RunConfigModel -from python_app.models.run_config_validation import load_ring_payload, load_switch_payload +from python_app.models.run_config_schema import ( + ComboModel, + GprRxGeometryModel, + GprTxGeometryModel, + RunConfigModel, +) +from python_app.models.run_config_validation import load_ring_payload, load_switch_payload, validate_gpr_model def _as_dict(value: Any, context: str) -> dict[str, Any]: @@ -29,6 +34,7 @@ def run_config_from_dict(payload: dict[str, Any]) -> RunConfigModel: port2_payload = _as_dict(switches_payload.get("port2"), "switches.port2") run_payload = _as_dict(payload.get("run"), "run") preprocess_payload = _as_dict(payload.get("preprocess"), "preprocess") + gpr_payload = _as_dict(payload.get("gpr"), "gpr") rings_payload = _as_dict(payload.get("rings"), "rings") raw_ring_payload = _as_dict(rings_payload.get("raw"), "rings.raw") raw_tap_ring_payload = _as_dict(rings_payload.get("raw_tap"), "rings.raw_tap") @@ -68,6 +74,38 @@ def run_config_from_dict(payload: dict[str, Any]) -> RunConfigModel: preprocess_payload.get("reference_bundle_path", model.preprocess.reference_bundle_path) ) + model.gpr.mode = str(gpr_payload.get("mode", model.gpr.mode)) + model.gpr.relative_permittivity = float( + gpr_payload.get("relative_permittivity", model.gpr.relative_permittivity) + ) + model.gpr.tx_geometry = [] + tx_geometry_payload = gpr_payload.get("tx_geometry", []) + if isinstance(tx_geometry_payload, list): + for entry in tx_geometry_payload: + entry_payload = _as_dict(entry, "gpr.tx_geometry[]") + model.gpr.tx_geometry.append( + GprTxGeometryModel( + output_pos=int(entry_payload.get("output_pos", 0)), + x_m=float(entry_payload.get("x_m", 0.0)), + ) + ) + model.gpr.rx_geometry = [] + rx_geometry_payload = gpr_payload.get("rx_geometry", []) + if isinstance(rx_geometry_payload, list): + for entry in rx_geometry_payload: + entry_payload = _as_dict(entry, "gpr.rx_geometry[]") + model.gpr.rx_geometry.append( + GprRxGeometryModel( + input_pos=int(entry_payload.get("input_pos", 0)), + x_m=float(entry_payload.get("x_m", 0.0)), + ) + ) + validate_gpr_model( + model.gpr, + input_switch_positions=model.input_switch.positions, + output_switch_positions=model.output_switch.positions, + ) + load_ring_payload(raw_ring_payload, model.rings.raw) load_ring_payload(raw_tap_ring_payload, model.rings.raw_tap) load_ring_payload(pre_ring_payload, model.rings.preprocessed) @@ -145,6 +183,24 @@ def run_config_to_dict(model: RunConfigModel) -> dict[str, Any]: "calibration_bundle_path": model.preprocess.calibration_bundle_path, "reference_bundle_path": model.preprocess.reference_bundle_path, }, + "gpr": { + "mode": model.gpr.mode, + "relative_permittivity": model.gpr.relative_permittivity, + "tx_geometry": [ + { + "output_pos": entry.output_pos, + "x_m": entry.x_m, + } + for entry in model.gpr.tx_geometry + ], + "rx_geometry": [ + { + "input_pos": entry.input_pos, + "x_m": entry.x_m, + } + for entry in model.gpr.rx_geometry + ], + }, "rings": { "raw": { "name": model.rings.raw.name, diff --git a/python_app/models/run_config_model.py b/python_app/models/run_config_model.py index d7c4908..ecf8686 100644 --- a/python_app/models/run_config_model.py +++ b/python_app/models/run_config_model.py @@ -3,6 +3,9 @@ from python_app.models.run_config_codec import run_config_from_dict, run_config_to_dict from python_app.models.run_config_schema import ( ComboModel, + GprModel, + GprRxGeometryModel, + GprTxGeometryModel, PreprocessModel, RadarModel, RadarSweepModel, @@ -20,6 +23,9 @@ from python_app.models.run_config_validation import ( __all__ = [ "ComboModel", + "GprModel", + "GprRxGeometryModel", + "GprTxGeometryModel", "PreprocessModel", "RadarModel", "RadarSweepModel", diff --git a/python_app/models/run_config_schema.py b/python_app/models/run_config_schema.py index cc91ebe..995bd62 100644 --- a/python_app/models/run_config_schema.py +++ b/python_app/models/run_config_schema.py @@ -95,6 +95,32 @@ class PreprocessModel: reference_bundle_path: str = "" +@dataclass(slots=True) +class GprTxGeometryModel: + """One transmitter geometry record keyed by output switch position.""" + + output_pos: int = 0 + x_m: float = 0.0 + + +@dataclass(slots=True) +class GprRxGeometryModel: + """One receiver geometry record keyed by input switch position.""" + + input_pos: int = 0 + x_m: float = 0.0 + + +@dataclass(slots=True) +class GprModel: + """Stable GPR configuration saved in run_config.json.""" + + mode: str = "point" + relative_permittivity: float = 1.0 + tx_geometry: list[GprTxGeometryModel] = field(default_factory=list) + rx_geometry: list[GprRxGeometryModel] = field(default_factory=list) + + @dataclass(slots=True) class RunConfigModel: """Top-level runtime config model consumed by C++ processes and GUI.""" @@ -105,6 +131,7 @@ class RunConfigModel: rings: RingsModel = field(default_factory=RingsModel) runtime: RuntimeModel = field(default_factory=RuntimeModel) preprocess: PreprocessModel = field(default_factory=PreprocessModel) + gpr: GprModel = field(default_factory=GprModel) combos: list[ComboModel] = field(default_factory=list) @staticmethod diff --git a/python_app/models/run_config_validation.py b/python_app/models/run_config_validation.py index 26080d0..ae2b343 100644 --- a/python_app/models/run_config_validation.py +++ b/python_app/models/run_config_validation.py @@ -4,7 +4,7 @@ from __future__ import annotations from typing import Any -from python_app.models.run_config_schema import ComboModel, RingEndpointModel, SwitchModel +from python_app.models.run_config_schema import ComboModel, GprModel, RingEndpointModel, SwitchModel def load_switch_payload( payload: dict[str, Any], @@ -30,6 +30,37 @@ def load_ring_payload(payload: dict[str, Any], target: RingEndpointModel) -> Non target.slot_size_bytes = int(payload.get("slot_size_bytes", target.slot_size_bytes)) +def validate_gpr_model( + gpr: GprModel, + *, + input_switch_positions: int, + output_switch_positions: int, +) -> None: + """Validate stable GPR config against current switch dimensions.""" + if gpr.mode not in {"point", "extended"}: + raise ValueError("gpr.mode must be either 'point' or 'extended'") + if float(gpr.relative_permittivity) <= 0.0: + raise ValueError("gpr.relative_permittivity must be > 0") + + seen_output_positions: set[int] = set() + for entry in gpr.tx_geometry: + output_pos = int(entry.output_pos) + if output_pos < 0 or output_pos >= int(output_switch_positions): + raise ValueError("gpr.tx_geometry output_pos is out of range") + if output_pos in seen_output_positions: + raise ValueError("gpr.tx_geometry contains duplicate output_pos") + seen_output_positions.add(output_pos) + + seen_input_positions: set[int] = set() + for entry in gpr.rx_geometry: + input_pos = int(entry.input_pos) + if input_pos < 0 or input_pos >= int(input_switch_positions): + raise ValueError("gpr.rx_geometry input_pos is out of range") + if input_pos in seen_input_positions: + raise ValueError("gpr.rx_geometry contains duplicate input_pos") + seen_input_positions.add(input_pos) + + def parse_combos_from_text(text: str) -> list[ComboModel]: """Parse UI combos string in `input:output,input:output` format.""" cleaned = text.strip() diff --git a/python_app/orchestration/live_processing_config.py b/python_app/orchestration/live_processing_config.py index 155a13a..1b67761 100644 --- a/python_app/orchestration/live_processing_config.py +++ b/python_app/orchestration/live_processing_config.py @@ -14,27 +14,62 @@ class ProcessingLiveConfig: processor_mode: str = "pass_through" gain_db: float = 0.0 phase_deg: float = 0.0 + pass_through_fixed_y_enabled: bool = False + pass_through_y_min_db: float = -100.0 + pass_through_y_max_db: float = 0.0 bscan_axis: str = "abs" bscan_cut_m: float = 0.824 bscan_max_depth_m: float = 1.0 bscan_gain: float = 1.0 bscan_start_freq_mhz: float = 100.0 bscan_stop_freq_mhz: float = 8800.0 + gpr_input_positions: list[int] | None = None + gpr_output_positions: list[int] | None = None + gpr_min_depth_m: float = 2.0 + gpr_max_depth_m: float = 14.0 + gpr_comp_power: float = 0.2 + gpr_start_freq_mhz: float = 3000.0 + gpr_stop_freq_mhz: float = 6000.0 + gpr_background_subtract_enabled: bool = True + gpr_background_mean_count: int = 10 history_command_seq: int = 0 history_command: str = "none" - def to_dict(self) -> dict[str, float | str | int]: + def __post_init__(self) -> None: + """Normalize optional list fields to concrete integer lists.""" + if self.gpr_input_positions is None: + self.gpr_input_positions = [] + else: + self.gpr_input_positions = [int(value) for value in self.gpr_input_positions] + if self.gpr_output_positions is None: + self.gpr_output_positions = [] + else: + self.gpr_output_positions = [int(value) for value in self.gpr_output_positions] + + def to_dict(self) -> dict[str, object]: """Convert live config to JSON-serializable dictionary.""" return { "processor_mode": str(self.processor_mode), "gain_db": float(self.gain_db), "phase_deg": float(self.phase_deg), + "pass_through_fixed_y_enabled": bool(self.pass_through_fixed_y_enabled), + "pass_through_y_min_db": float(self.pass_through_y_min_db), + "pass_through_y_max_db": float(self.pass_through_y_max_db), "bscan_axis": str(self.bscan_axis), "bscan_cut_m": float(self.bscan_cut_m), "bscan_max_depth_m": float(self.bscan_max_depth_m), "bscan_gain": float(self.bscan_gain), "bscan_start_freq_mhz": float(self.bscan_start_freq_mhz), "bscan_stop_freq_mhz": float(self.bscan_stop_freq_mhz), + "gpr_input_positions": [int(value) for value in self.gpr_input_positions], + "gpr_output_positions": [int(value) for value in self.gpr_output_positions], + "gpr_min_depth_m": float(self.gpr_min_depth_m), + "gpr_max_depth_m": float(self.gpr_max_depth_m), + "gpr_comp_power": float(self.gpr_comp_power), + "gpr_start_freq_mhz": float(self.gpr_start_freq_mhz), + "gpr_stop_freq_mhz": float(self.gpr_stop_freq_mhz), + "gpr_background_subtract_enabled": bool(self.gpr_background_subtract_enabled), + "gpr_background_mean_count": int(self.gpr_background_mean_count), "history_command_seq": int(self.history_command_seq), "history_command": str(self.history_command), } diff --git a/python_app/orchestration/shm/decoder.py b/python_app/orchestration/shm/decoder.py index 1daac7a..0d643cd 100644 --- a/python_app/orchestration/shm/decoder.py +++ b/python_app/orchestration/shm/decoder.py @@ -56,6 +56,62 @@ def decode_trace_collection(payload: bytes, expected_magic: int) -> SweepCollect def decode_result_collection(payload: bytes) -> ResultCollection: """Decode one processed result collection from binary payload.""" + + def read_payload(cursor: ByteCursor) -> ResultPayload: + """Decode one result payload from stream.""" + kind = cursor.read_u8() + name_size = cursor.read_u16() + name = cursor.read_bytes(name_size).decode("utf-8") + + if kind == 1: + point_count = cursor.read_u32() + freq = np.frombuffer(cursor.read_bytes(point_count * 4), dtype=" 0 else np.zeros((0, 0), dtype=np.float32) + return ResultPayload( + processing_name=name, + kind=kind, + image_x_axis=image_x_axis, + image_y_axis=image_y_axis, + image=image, + ) + if kind == 4: + table_columns = cursor.read_u32() + table_rows = cursor.read_u32() + value_count = table_columns * table_rows + table_values = np.frombuffer(cursor.read_bytes(value_count * 4), dtype=" 0 and table_columns > 0 + else np.zeros((0, 0), dtype=np.float32) + ) + return ResultPayload( + processing_name=name, + kind=kind, + table=table, + ) + raise ValueError(f"Unsupported result payload kind: {kind}") + cursor = ByteCursor(payload) magic = cursor.read_u32() if magic != RESULT_MAGIC: @@ -63,8 +119,13 @@ def decode_result_collection(payload: bytes) -> ResultCollection: collection_id = cursor.read_u64() monotonic_ns = cursor.read_u64() + collection_payload_count = cursor.read_u32() block_count = cursor.read_u32() + collection_payloads: list[ResultPayload] = [] + for _ in range(collection_payload_count): + collection_payloads.append(read_payload(cursor)) + blocks: list[ResultBlock] = [] for _ in range(block_count): input_pos = cursor.read_u32() @@ -73,36 +134,7 @@ def decode_result_collection(payload: bytes) -> ResultCollection: payloads: list[ResultPayload] = [] for _ in range(payload_count): - kind = cursor.read_u8() - name_size = cursor.read_u16() - name = cursor.read_bytes(name_size).decode("utf-8") - - if kind == 1: - point_count = cursor.read_u32() - freq = np.frombuffer(cursor.read_bytes(point_count * 4), dtype=" ResultCollection: ) ) - return ResultCollection(collection_id=collection_id, monotonic_ns=monotonic_ns, blocks=blocks) + return ResultCollection( + collection_id=collection_id, + monotonic_ns=monotonic_ns, + collection_payloads=collection_payloads, + blocks=blocks, + ) diff --git a/python_app/runtime/run_config_smoke.json b/python_app/runtime/run_config_smoke.json index c7e706b..cb081c4 100644 --- a/python_app/runtime/run_config_smoke.json +++ b/python_app/runtime/run_config_smoke.json @@ -28,13 +28,13 @@ "port2": { "name": "port2", "driver_mode": "mock", - "driver": "hmc349a", + "driver": "h7992", "radar_port": 2, - "positions": 2, + "positions": 4, "default_position": 0, "gpio_chip": "/dev/gpiochip0", "pin_a": 22, - "pin_b": -1, + "pin_b": 23, "invert_logic": false } }, @@ -53,20 +53,12 @@ "output": 0 }, { - "input": 0, - "output": 1 + "input": 2, + "output": 0 }, { - "input": 1, - "output": 1 - }, - { - "input": 0, - "output": 2 - }, - { - "input": 1, - "output": 2 + "input": 3, + "output": 0 }, { "input": 0, @@ -75,6 +67,14 @@ { "input": 1, "output": 3 + }, + { + "input": 2, + "output": 3 + }, + { + "input": 3, + "output": 3 } ] }, @@ -84,6 +84,38 @@ "calibration_bundle_path": "/home/europa/Documents/radar_system/python_app/runtime/calibration_bundle.bin", "reference_bundle_path": "/home/europa/Documents/radar_system/python_app/runtime/reference_bundle.bin" }, + "gpr": { + "mode": "point", + "relative_permittivity": 1.0, + "tx_geometry": [ + { + "output_pos": 0, + "x_m": 0.905 + }, + { + "output_pos": 3, + "x_m": -0.905 + } + ], + "rx_geometry": [ + { + "input_pos": 0, + "x_m": -0.18 + }, + { + "input_pos": 1, + "x_m": 0.485 + }, + { + "input_pos": 2, + "x_m": -0.49 + }, + { + "input_pos": 3, + "x_m": 0.185 + } + ] + }, "rings": { "raw": { "name": "/radar_raw_smoke_1703912_791574940686872", @@ -111,4 +143,4 @@ "slot_size_bytes": 2097152 } } -} \ No newline at end of file +} diff --git a/python_app/scripts/convert_prog_libre_manual_to_vna_history.py b/python_app/scripts/convert_prog_libre_manual_to_vna_history.py new file mode 100644 index 0000000..57fc4f8 --- /dev/null +++ b/python_app/scripts/convert_prog_libre_manual_to_vna_history.py @@ -0,0 +1,384 @@ +"""Convert manual LibreVNA S11/S21 CSV captures in prog_libre to vna history JSON.""" + +from __future__ import annotations + +import argparse +import csv +from datetime import datetime, timezone +import json +from pathlib import Path +from typing import Any + +import numpy as np + + +def _real_imag_keys(trace_prefix: str) -> tuple[str, str]: + return f"{trace_prefix}_Real", f"{trace_prefix}_Imaginary" + + +def _load_complex_trace(csv_path: Path, trace_prefix: str) -> tuple[np.ndarray, np.ndarray]: + real_key, imag_key = _real_imag_keys(trace_prefix) + frequencies: list[float] = [] + values: list[complex] = [] + + with csv_path.open(encoding="utf-8", newline="") as handle: + reader = csv.DictReader(handle) + for row in reader: + frequencies.append(float(row["Frequency"])) + values.append(complex(float(row[real_key]), float(row[imag_key]))) + + frequency_hz = np.asarray(frequencies, dtype=np.float64) + trace = np.asarray(values, dtype=np.complex128) + if frequency_hz.size == 0 or trace.size == 0: + raise ValueError(f"CSV has no points: {csv_path}") + if frequency_hz.shape != trace.shape: + raise ValueError(f"Frequency/trace size mismatch: {csv_path}") + return frequency_hz, trace + + +def _solve_one_port_osl( + open_trace: np.ndarray, + short_trace: np.ndarray, + load_trace: np.ndarray, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Solve ideal OSL one-port calibration coefficients.""" + directivity = load_trace + open_delta = open_trace - directivity + short_delta = short_trace - directivity + denom = open_delta - short_delta + + source_match = np.zeros_like(directivity) + reflection_tracking = np.ones_like(directivity) + + stable_mask = np.abs(denom) > 1e-18 + source_match[stable_mask] = (open_delta[stable_mask] + short_delta[stable_mask]) / denom[stable_mask] + reflection_tracking[stable_mask] = open_delta[stable_mask] * (1.0 - source_match[stable_mask]) + return directivity, source_match, reflection_tracking + + +def _apply_one_port_osl( + measured_trace: np.ndarray, + directivity: np.ndarray, + source_match: np.ndarray, + reflection_tracking: np.ndarray, +) -> np.ndarray: + numerator = measured_trace - directivity + denominator = reflection_tracking + (source_match * numerator) + + corrected = np.array(numerator, copy=True) + stable_mask = np.abs(denominator) > 1e-18 + corrected[stable_mask] = numerator[stable_mask] / denominator[stable_mask] + return corrected + + +def _apply_through_calibration(measured_trace: np.ndarray, through_trace: np.ndarray) -> np.ndarray: + corrected = np.array(measured_trace, copy=True) + stable_mask = np.abs(through_trace) > 1e-18 + corrected[stable_mask] = measured_trace[stable_mask] / through_trace[stable_mask] + return corrected + + +def _complex_to_points(values: np.ndarray) -> list[list[float]]: + return [[float(value.real), float(value.imag)] for value in values] + + +def _scan_file_sort_key(csv_path: Path) -> tuple[int, str]: + stem = csv_path.stem + return (int(stem), stem) if stem.isdigit() else (10**9, stem) + + +def _load_scan_series(folder: Path, trace_prefix: str) -> tuple[np.ndarray, list[tuple[str, np.ndarray]]]: + scan_paths = [ + path + for path in sorted(folder.glob("*.csv"), key=_scan_file_sort_key) + if path.stem.isdigit() + ] + if not scan_paths: + raise FileNotFoundError(f"No numbered scan CSV files found in {folder}") + + base_frequency_hz: np.ndarray | None = None + scans: list[tuple[str, np.ndarray]] = [] + for path in scan_paths: + frequency_hz, trace = _load_complex_trace(path, trace_prefix) + if base_frequency_hz is None: + base_frequency_hz = frequency_hz + elif not np.allclose(base_frequency_hz, frequency_hz, rtol=0.0, atol=1e-6): + raise ValueError(f"Frequency axis mismatch in {path}") + scans.append((path.name, trace)) + + assert base_frequency_hz is not None + return base_frequency_hz, scans + + +def _require_matching_frequency_axis(label: str, left: np.ndarray, right: np.ndarray) -> None: + if not np.allclose(left, right, rtol=0.0, atol=1e-6): + raise ValueError(f"{label} frequency axes do not match") + + +def _build_history_payload( + *, + source_dir: Path, + mode: str, + frequency_hz: np.ndarray, + sweep_scans: list[tuple[str, np.ndarray]], + calibrated_scans: list[tuple[str, np.ndarray]], + reference_trace: np.ndarray, + primary_stage: str, + raw_record_count: int, + preprocessed_record_count: int, +) -> dict[str, Any]: + if len(sweep_scans) != len(calibrated_scans): + raise ValueError("Sweep/calibrated scan counts do not match") + + sweep_history: list[dict[str, Any]] = [] + reference_points = _complex_to_points(reference_trace) + for index, ((scan_name, sweep_trace), (cal_name, calibrated_trace)) in enumerate( + zip(sweep_scans, calibrated_scans, strict=True) + ): + if scan_name != cal_name: + raise ValueError(f"Scan ordering mismatch: {scan_name} vs {cal_name}") + + sweep_history.append( + { + "timestamp": float(index), + "sweep_points": _complex_to_points(sweep_trace), + "calibrated_points": _complex_to_points(calibrated_trace), + "reference_points": reference_points, + "vna_config": { + "mode": mode, + "start_freq": float(frequency_hz[0]), + "stop_freq": float(frequency_hz[-1]), + "points": int(frequency_hz.size), + }, + } + ) + + return { + "format": "vna-system-history-v1", + "converter": "python_app/scripts/convert_prog_libre_manual_to_vna_history.py", + "converted_at_utc": datetime.now(timezone.utc).isoformat(), + "source_snapshot_dir": str(source_dir.resolve()), + "input_index": 0, + "output_index": 0, + "primary_stage": primary_stage, + "raw_record_count": int(raw_record_count), + "preprocessed_record_count": int(preprocessed_record_count), + "sweep_history": sweep_history, + } + + +def _write_payload(output_path: Path, payload: dict[str, Any]) -> None: + output_path.parent.mkdir(parents=True, exist_ok=True) + output_path.write_text(json.dumps(payload, ensure_ascii=False, indent=2), encoding="utf-8") + + +def _rmse(left: np.ndarray, right: np.ndarray) -> float: + return float(np.sqrt(np.mean(np.abs(left - right) ** 2))) + + +def _build_parser() -> argparse.ArgumentParser: + parser = argparse.ArgumentParser( + description="Convert manual LibreVNA S11/S21 CSV captures in prog_libre to vna history JSON.", + ) + parser.add_argument( + "--calibration-dir", + type=Path, + default=Path("prog_libre/calibration"), + help="Directory with calibration CSV files.", + ) + parser.add_argument( + "--raw-dir", + type=Path, + default=Path("prog_libre/1-6000mhz_no-calibrated_libre"), + help="Directory with uncalibrated scan CSV files and ref.csv.", + ) + parser.add_argument( + "--calibrated-dir", + type=Path, + default=Path("prog_libre/1-6000mhz_calibrated_libre"), + help="Directory with already calibrated scan CSV files and ref.csv.", + ) + parser.add_argument( + "--raw-output", + "--s11-raw-output", + dest="s11_raw_output", + type=Path, + default=Path("prog_libre/1-6000mhz_no-calibrated_libre_s11_osl_p1_vna_bscan_history.json"), + help="Output JSON for uncalibrated S11 scans after applying OSL calibration.", + ) + parser.add_argument( + "--calibrated-output", + "--s11-calibrated-output", + dest="s11_calibrated_output", + type=Path, + default=Path("prog_libre/1-6000mhz_calibrated_libre_s11_passthrough_vna_bscan_history.json"), + help="Output JSON for already calibrated S11 scans without extra calibration.", + ) + parser.add_argument( + "--s21-raw-output", + dest="s21_raw_output", + type=Path, + default=Path("prog_libre/1-6000mhz_no-calibrated_libre_s21_through_vna_bscan_history.json"), + help="Output JSON for uncalibrated S21 scans after applying through calibration.", + ) + parser.add_argument( + "--s21-calibrated-output", + dest="s21_calibrated_output", + type=Path, + default=Path("prog_libre/1-6000mhz_calibrated_libre_s21_passthrough_vna_bscan_history.json"), + help="Output JSON for already calibrated S21 scans without extra calibration.", + ) + return parser + + +def main() -> None: + args = _build_parser().parse_args() + + calibration_dir = args.calibration_dir.expanduser().resolve() + raw_dir = args.raw_dir.expanduser().resolve() + calibrated_dir = args.calibrated_dir.expanduser().resolve() + + s11_raw_output = args.s11_raw_output.expanduser().resolve() + s11_calibrated_output = args.s11_calibrated_output.expanduser().resolve() + s21_raw_output = args.s21_raw_output.expanduser().resolve() + s21_calibrated_output = args.s21_calibrated_output.expanduser().resolve() + + s11_cal_frequency_hz, open_trace = _load_complex_trace(calibration_dir / "open_rfc18_p1.csv", "S11") + short_frequency_hz, short_trace = _load_complex_trace(calibration_dir / "short_rfc18_p1.csv", "S11") + load_frequency_hz, load_trace = _load_complex_trace(calibration_dir / "load_rfc18_p1.csv", "S11") + _require_matching_frequency_axis("S11 calibration", s11_cal_frequency_hz, short_frequency_hz) + _require_matching_frequency_axis("S11 calibration", s11_cal_frequency_hz, load_frequency_hz) + + directivity, source_match, reflection_tracking = _solve_one_port_osl(open_trace, short_trace, load_trace) + + s11_raw_frequency_hz, s11_raw_scans = _load_scan_series(raw_dir, "S11") + s11_calibrated_frequency_hz, s11_passthrough_scans = _load_scan_series(calibrated_dir, "S11") + _require_matching_frequency_axis("S11 raw vs calibration", s11_raw_frequency_hz, s11_cal_frequency_hz) + _require_matching_frequency_axis("S11 calibrated vs calibration", s11_calibrated_frequency_hz, s11_cal_frequency_hz) + + s11_raw_reference_frequency_hz, s11_raw_reference = _load_complex_trace(raw_dir / "ref.csv", "S11") + s11_calibrated_reference_frequency_hz, s11_calibrated_reference = _load_complex_trace(calibrated_dir / "ref.csv", "S11") + _require_matching_frequency_axis("S11 raw reference vs calibration", s11_raw_reference_frequency_hz, s11_cal_frequency_hz) + _require_matching_frequency_axis( + "S11 calibrated reference vs calibration", + s11_calibrated_reference_frequency_hz, + s11_cal_frequency_hz, + ) + + s11_corrected_scans = [ + ( + scan_name, + _apply_one_port_osl(trace, directivity, source_match, reflection_tracking), + ) + for scan_name, trace in s11_raw_scans + ] + s11_corrected_reference = _apply_one_port_osl(s11_raw_reference, directivity, source_match, reflection_tracking) + + s11_raw_payload = _build_history_payload( + source_dir=raw_dir, + mode="s11", + frequency_hz=s11_raw_frequency_hz, + sweep_scans=s11_raw_scans, + calibrated_scans=s11_corrected_scans, + reference_trace=s11_corrected_reference, + primary_stage="raw", + raw_record_count=len(s11_raw_scans), + preprocessed_record_count=len(s11_corrected_scans), + ) + s11_calibrated_payload = _build_history_payload( + source_dir=calibrated_dir, + mode="s11", + frequency_hz=s11_calibrated_frequency_hz, + sweep_scans=s11_passthrough_scans, + calibrated_scans=s11_passthrough_scans, + reference_trace=s11_calibrated_reference, + primary_stage="preprocessed", + raw_record_count=0, + preprocessed_record_count=len(s11_passthrough_scans), + ) + + s21_cal_frequency_hz, through_trace = _load_complex_trace(calibration_dir / "through21_rfc18.csv", "S21") + + s21_raw_frequency_hz, s21_raw_scans = _load_scan_series(raw_dir, "S21") + s21_calibrated_frequency_hz, s21_passthrough_scans = _load_scan_series(calibrated_dir, "S21") + _require_matching_frequency_axis("S21 raw vs calibration", s21_raw_frequency_hz, s21_cal_frequency_hz) + _require_matching_frequency_axis("S21 calibrated vs calibration", s21_calibrated_frequency_hz, s21_cal_frequency_hz) + + s21_raw_reference_frequency_hz, s21_raw_reference = _load_complex_trace(raw_dir / "ref.csv", "S21") + s21_calibrated_reference_frequency_hz, s21_calibrated_reference = _load_complex_trace(calibrated_dir / "ref.csv", "S21") + _require_matching_frequency_axis("S21 raw reference vs calibration", s21_raw_reference_frequency_hz, s21_cal_frequency_hz) + _require_matching_frequency_axis( + "S21 calibrated reference vs calibration", + s21_calibrated_reference_frequency_hz, + s21_cal_frequency_hz, + ) + + s21_corrected_scans = [ + ( + scan_name, + _apply_through_calibration(trace, through_trace), + ) + for scan_name, trace in s21_raw_scans + ] + s21_corrected_reference = _apply_through_calibration(s21_raw_reference, through_trace) + + s21_raw_payload = _build_history_payload( + source_dir=raw_dir, + mode="s21", + frequency_hz=s21_raw_frequency_hz, + sweep_scans=s21_raw_scans, + calibrated_scans=s21_corrected_scans, + reference_trace=s21_corrected_reference, + primary_stage="raw", + raw_record_count=len(s21_raw_scans), + preprocessed_record_count=len(s21_corrected_scans), + ) + s21_calibrated_payload = _build_history_payload( + source_dir=calibrated_dir, + mode="s21", + frequency_hz=s21_calibrated_frequency_hz, + sweep_scans=s21_passthrough_scans, + calibrated_scans=s21_passthrough_scans, + reference_trace=s21_calibrated_reference, + primary_stage="preprocessed", + raw_record_count=0, + preprocessed_record_count=len(s21_passthrough_scans), + ) + + _write_payload(s11_raw_output, s11_raw_payload) + _write_payload(s11_calibrated_output, s11_calibrated_payload) + _write_payload(s21_raw_output, s21_raw_payload) + _write_payload(s21_calibrated_output, s21_calibrated_payload) + + s11_rmse_values = [ + _rmse(corrected_trace, passthrough_trace) + for (_, corrected_trace), (_, passthrough_trace) in zip(s11_corrected_scans, s11_passthrough_scans, strict=True) + ] + s21_rmse_values = [ + _rmse(corrected_trace, passthrough_trace) + for (_, corrected_trace), (_, passthrough_trace) in zip(s21_corrected_scans, s21_passthrough_scans, strict=True) + ] + s11_reference_rmse = _rmse(s11_corrected_reference, s11_calibrated_reference) + s21_reference_rmse = _rmse(s21_corrected_reference, s21_calibrated_reference) + + print( + "Converted manual prog_libre captures to vna history JSON:\n" + f" S11 raw output: {s11_raw_output}\n" + f" S11 calibrated output: {s11_calibrated_output}\n" + f" S21 raw output: {s21_raw_output}\n" + f" S21 calibrated output: {s21_calibrated_output}\n" + f" sweep count: {len(s11_raw_scans)}\n" + f" points per sweep: {s11_raw_frequency_hz.size}\n" + f" S11 calibration: p1 ideal OSL\n" + f" S11 mean sweep RMSE vs provided calibrated folder: {float(np.mean(s11_rmse_values)):.6f}\n" + f" S11 max sweep RMSE vs provided calibrated folder: {float(np.max(s11_rmse_values)):.6f}\n" + f" S11 reference RMSE vs provided calibrated ref: {s11_reference_rmse:.6f}\n" + f" S21 calibration: through21 complex division\n" + f" S21 mean sweep RMSE vs provided calibrated folder: {float(np.mean(s21_rmse_values)):.6f}\n" + f" S21 max sweep RMSE vs provided calibrated folder: {float(np.max(s21_rmse_values)):.6f}\n" + f" S21 reference RMSE vs provided calibrated ref: {s21_reference_rmse:.6f}" + ) + + +if __name__ == "__main__": + main() diff --git a/python_app/storage/npz/serialize.py b/python_app/storage/npz/serialize.py index 3a3210b..aaa94d6 100644 --- a/python_app/storage/npz/serialize.py +++ b/python_app/storage/npz/serialize.py @@ -39,39 +39,79 @@ def serialize_trace_collection(collection: SweepCollection, magic: int) -> bytes def serialize_result_collection(collection: ResultCollection) -> bytes: """Serialize one processed collection with result blocks/payloads.""" + + def serialize_payload(buffer: bytearray, payload) -> None: + """Append one payload in ring-compatible result format.""" + name_bytes = payload.processing_name.encode("utf-8") + if len(name_bytes) > 0xFFFF: + raise ValueError("processing_name is too long") + + buffer.extend(struct.pack(" 0xFFFF: - raise ValueError("processing_name is too long") - - buffer.extend(struct.pack("