Compare commits

4 Commits
Author SHA1 Message Date
aweandClaude Fable 5 7329750c0c sar3d: fixed default scene window X [-2,2], Y [-0.5,2], Z [0,8] m
The grid bounds are now explicit config (sar3d_x/y/z_min/max_m) instead of
being derived from the antenna span and the depth gate. The depth mask
still zeroes voxels outside [min_depth, max_depth], so the wider Z window
just shows the masked wedge. The bound fields are double: the parity tool
recomputes the notebook-derived bounds (notebook_grid_bounds) and passes
them through the live config, and a float roundtrip shifted the axes by
~1e-7 m, flipping depth-mask knife-edge voxels against the reference npz.
Parity on frames 0034-0036 unchanged (1.3e-6, axes bit-match).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
2026-09-25 17:12:18 +03:00
aweandClaude Fable 5 60f6c9e567 sar3d: chunk the aperture kernel over voxels to survive low free VRAM
Full-volume acc/cnt device buffers (~115 MB for the 9.6M-voxel grid) plus
the CUDA context did not fit next to a desktop session that had taken the
T600's 4 GB down to ~190 MB free, so every frame died with cuMemAlloc OOM
and the GUI showed no sar3d picture. The kernel now processes
[gid_begin, gid_end) ranges and the host launches it in 2M-voxel chunks
reusing one ~24 MB buffer pair. Parity vs the reference volumes is
unchanged (max diff 1.3e-6).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
2026-09-25 17:01:22 +03:00
aweandClaude Fable 5 e1dfbec712 sar3d: python plumbing for the new parameters and real replay timing
- replay_raw_producer carries the RECORDED per-trace capture windows
  from each snapshot's meta.json (re-anchored by the same shift as the
  frame time) instead of one synthetic window for all traces — without
  this the segmented timing model degenerates in replay;
- GUI/profile/orchestration swap the removed window knobs
  (center_depth, half_width_*, grid_n*) for the notebook parameters
  (ground depth, pitch override, voxel pitch, Y range, X margin, norm
  edge margin, smooth sigma, comp power, tx per segment).

Headless replay smoke on the 20260922_mine_fast_true dataset passes:
51 sar3d result collections, transparent preprocessor.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
2026-09-25 16:42:14 +03:00
aweandClaude Fable 5 b8cb572bc4 sar3d: port the Multi_antenna_3D_clean_v1 notebook algorithm
The old sar3d core was a port of the superseded Horns_3D notebook and
produced wrong volumes. The compute core now follows the reference
notebook exactly:

- segmented_interleaved frame timing rebuilt from per-trace capture
  stamps (uniform tx_per_segment split as the no-stamps fallback), with
  the coarse inter-segment dtau and the intra-sweep dtau applied as one
  phase ramp before the IFFT;
- ground scene coordinates with the ground_depth_m lift (scene Y=0 is
  the ground), full-region grid instead of a window around a center
  depth, geo_ref compensation anchored on the ground plane, power 0.3;
- notebook finalization in C++: |acc/cnt|, depth mask, normalization
  without the Y-edge planes, separable Gaussian over (Z, X), level;
- CUDA kernel computes geometry in double: float32 loses ~3e-4 of an
  A-scan sample in tau and drifts ~1e-4 from the float64 reference.

Live-config keys follow: the center_depth/half_width/grid_n window knobs
are replaced by ground/voxel/y-range/margin/sigma/power/tx_per_segment.

Parity against the notebook npz exports (frames 0034-0036 of the
20260922_mine_fast_true dataset): max|diff| <= 1.6e-6 off the depth-mask
knife edge, identical argmax, level rel diff <= 7e-7. Verified via
make sar3d_parity --dump-dir + scripts/sar3d_volume_compare.py; the
executable spec lives in scripts/sar3d_reference_check.py (bit-tight
against the npz, including the animation shifted-grid arithmetic).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
2026-09-25 16:42:01 +03:00
26 changed files with 1777 additions and 520 deletions
+6 -3
View File
@@ -145,11 +145,14 @@ all: $(TARGETS)
kamil_adc_collector: $(BIN_DIR)/kamil_adc_collector
# On-demand parity check of the sar3d processor against the gpr_cpp
# standalone_bp example (same .npy dataset). Not part of `all`.
# On-demand parity check of the sar3d processor against the reference notebook
# volumes (scripts/sar3d_volume_compare.py diffs the dumped volume against a
# *_volume.npz export). Not part of `all`.
SAR3D_PARITY_OBJS := $(addprefix $(BUILD_DIR)/, \
data_acq_and_processing/processing/tools/sar3d_parity_main.o \
data_acq_and_processing/processing/processors/src/sar3d_processor.o)
data_acq_and_processing/processing/processors/src/sar3d_processor.o \
data_acq_and_processing/common_cpp/config/src/run_config.o \
data_acq_and_processing/common_cpp/ipc/src/shared_types.o)
sar3d_parity: $(BIN_DIR)/sar3d_parity
+28 -30
View File
@@ -1,23 +1,24 @@
# gpr_cpp — C++/CUDA счёт-ядро MIMO GPR 3D SAR
> Вендорено из репозитория VNA_AlGo (`gpr_cpp/`, коммит `544dc59`); python-биндинги
> (`bindings/`, pybind11) при переносе не взяты — сборка здесь только C++
> (`GPR_BUILD_PYTHON=OFF`). Собирается из корневого Makefile radar_system
> Инфраструктура (сборка, NVRTC, PocketFFT, npy) вендорена из репозитория
> VNA_AlGo (`gpr_cpp/`, коммит `544dc59`) без python-биндингов
> (`GPR_BUILD_PYTHON=OFF`). АЛГОРИТМ с сентября 2026 переписан: это порт
> эталонного ноутбука `Multi_antenna_3D_clean_v1.ipynb` (ячейки 0 и 4), а не
> старого `gpr_core`/Horns_3D. Паритет с ноутбуком проверяется инструментом
> `sar3d_parity` + `python_app/scripts/sar3d_volume_compare.py` по
> референсным `*_volume.npz`. Собирается из корневого Makefile radar_system
> (цель `$(GPR_CPP_LIB)` вызывает `build_cmake.sh`).
Порт вычислительного ядра из Python/CuPy (`gpr_core`) в самостоятельную
C++/CUDA библиотеку для интеграции в C++ проект.
Алгоритм — см. `ALGORITHM.md` в репозитории VNA_AlGo.
## Что внутри
## Что внутри (граница C++ ↔ Python)
**В C++ (`libgpr_core_cpp`):** геометрия, препроцессинг A-сканов (bg-subtract,
фазовая коррекция движения, k0-вставка + Ханн + IFFT через PocketFFT, нормировка),
3D SAR апертурный объём (одно CUDA-ядро на всю апертуру). Возвращает **сырые**
`acc`(complex64)/`cnt`(float32).
**Остаётся в Python:** загрузка `.npy`, финализация объёма (abs+gate+gaussian+
нормировка), CLEAN-детекция, трек, оценка скорости, plotly.
**В C++ (`libgpr_core_cpp`):** сегментированный тайминг кадра по per-trace
меткам захвата (`compute_pair_timestamps`), препроцессинг A-сканов (bg-subtract,
полная фазовая коррекция движения intra+pair до IFFT, k0-вставка + Ханн + IFFT
через PocketFFT, нормировка пар), полнообластная ground-сетка
(`build_full_region_grid`), 3D SAR апертурный объём (одно CUDA-ядро на всю
апертуру; геометрия в double, A-сканы float32) и финализация
(`finish_volume`: |acc/cnt| → маска глубины → нормировка без Y-краёв → гаусс
по (Z,X) → нормировка; возвращает объём [0,1] + level).
## Ключевое решение: CUDA без nvcc (NVRTC + Driver API)
@@ -42,14 +43,14 @@ libnvrtc прописан — бинарники/модуль грузятся
## Публичный C++ API ([include/gpr/](include/gpr/))
```cpp
#include "gpr/config.hpp" // GeometryConfig/BPConfig/MotionConfig/SARConfig (дефолты = config.py)
#include "gpr/ascan.hpp" // preprocess_snapshot(...) -> AScanSet
#include "gpr/sar3d.hpp" // SarEngine::aperture_volume(frames, grid) -> VolumeRaw{acc,cnt}
#include "gpr/geometry.hpp" // frame_offsets, frame_weights, build_grid_scene, point_to_scene
#include "gpr/config.hpp" // GeometryConfig/BPConfig/MotionConfig/SARConfig (дефолты = ноутбук)
#include "gpr/ascan.hpp" // compute_pair_timestamps, preprocess_snapshot(...) -> AScanSet
#include "gpr/sar3d.hpp" // SarEngine::aperture_volume -> VolumeRaw; finish_volume -> FinishedVolume
#include "gpr/geometry.hpp" // frame_offsets, frame_weights, build_full_region_grid
gpr::SarEngine engine(geom, bp, mot, sar); // компилит ядро под устройство
auto vol = engine.aperture_volume(frames, grid); // один запуск ядра на апертуру
// vol.acc / vol.cnt — сырые; финализацию (gaussian/normalize) делает вызывающий
auto raw = engine.aperture_volume(frames, grid); // один запуск ядра на апертуру
auto fin = gpr::finish_volume(raw, grid, mot.pitch_deg, bp, sar); // объём [0,1] + level
```
Пример без Python: [examples/standalone_bp.cpp](examples/standalone_bp.cpp) — грузит
S21 из `.npy`, строит объём, печатает максимум.
@@ -67,16 +68,13 @@ a = _gpr_cpp.preprocess_snapshot(s21, freq, background, mot, geom, bp, True) #
acc, cnt = eng.aperture_volume([(a, dx, dy, dz, w), ...], gx, gy, gz) # numpy
```
## Паритет (проверено на x86 T600)
## Паритет с ноутбуком (проверено на x86 T600, датасет 20260922_mine_fast_true)
| Этап | Порог | Факт |
|---|---|---|
| config дефолты vs dataclass | точно | OK |
| A-сканы raw H vs numpy (float64) | <1e-9 | 1.7e-12 |
| A-сканы норм. H (FFT-ULP через медиану) | <1e-4 | 2.5e-6 |
| **3D объём C++ ядро vs cupy** (float32) | <5e-3 | **3e-7, argmax идентичен** |
Тесты: [../tests/test_parity_cpp.py](../tests/test_parity_cpp.py) (C++ vs cupy).
Кадры 0034/0035/0036 против референсных `*_volume.npz` ноутбука
(`make sar3d_parity` + `python_app/scripts/sar3d_volume_compare.py`):
max|Δ| ≤ 1.6e-6 вне ножевой границы маски глубины, argmax и level совпадают
(level rel ≤ 7e-7). Python-эталон той же цепочки:
`python_app/scripts/sar3d_reference_check.py`.
## Статус
@@ -95,7 +95,20 @@ int main(int argc, char** argv) {
}
};
auto ts = compute_pair_timestamps(mot, geom);
// Демо без чтения meta.json: равномерные сегменты по интервалу между
// соседними снапшотами (реальные per-trace окна читает sar3d_parity).
double acq_s = 0.0;
{
std::vector<double> stamps;
for (auto& s : snaps) {
auto pos = s.rfind("_ns");
if (pos != std::string::npos) stamps.push_back(std::stod(s.substr(pos + 3)) * 1e-9);
}
std::sort(stamps.begin(), stamps.end());
if (stamps.size() >= 2) acq_s = stamps[1] - stamps[0];
}
auto windows = uniform_tx_windows(geom.output_idx, acq_s, mot.tx_per_segment);
auto ts = compute_pair_timestamps(windows, acq_s, geom.output_idx, geom, mot);
SarEngine engine(geom, bp, mot, sar);
// держим AScanSet живыми
@@ -122,22 +135,23 @@ int main(int argc, char** argv) {
frames.push_back(ApertureFrame{a.get(), offs[f], wts[f]});
}
// центр сетки: под центральным снапшотом на глубине 5.6 м (сцена)
auto center_scene = point_to_scene(0.0, 0.0, 5.6, mot, sar);
Grid grid = build_grid_scene(center_scene, sar);
// ноутбучные границы сцены (размах антенн ±2 м, полоса глубин) вместо
// фиксированного окна по умолчанию
notebook_grid_bounds(geom, bp, sar, mot.pitch_deg, 2.0,
sar.x_min_m, sar.x_max_m, sar.z_min_m, sar.z_max_m);
Grid grid = build_full_region_grid(sar);
printf("сетка %dx%dx%d, счёт объёма...\n", grid.NX(), grid.NY(), grid.NZ());
VolumeRaw vol = engine.aperture_volume(frames, grid);
FinishedVolume fin = finish_volume(vol, grid, mot.pitch_deg, bp, sar);
// финализация: |acc/cnt|, максимум
double vmax = 0; size_t amax = 0;
for (size_t i = 0; i < vol.acc.size(); ++i) {
double v = std::abs(vol.acc[i]) / (vol.cnt[i] + 1e-12);
if (v > vmax) { vmax = v; amax = i; }
}
for (size_t i = 0; i < fin.volume.size(); ++i)
if (fin.volume[i] > vmax) { vmax = fin.volume[i]; amax = i; }
int NX = vol.NX, NZ = vol.NZ;
int ix = amax % NX, iz = (amax / NX) % NZ, iy = amax / ((size_t)NX * NZ);
printf("МАКСИМУМ объёма: |acc/cnt|=%.4e в вокселе (ix=%d,iy=%d,iz=%d)\n", vmax, ix, iy, iz);
printf("МАКСИМУМ объёма: %.4f (level %.4e) в вокселе (ix=%d,iy=%d,iz=%d)\n",
vmax, fin.level, ix, iy, iz);
printf(" X=%.1f см, Y=%.1f см, Z=%.1f см\n",
grid.x[ix] * 100, grid.y[iy] * 100, grid.z[iz] * 100);
printf("OK: автономный C++ путь отработал без Python.\n");
@@ -1,4 +1,6 @@
// Препроцессинг S21 → комплексные A-сканы (host, float64). Порт preprocess.py.
// Препроцессинг S21 → комплексные A-сканы (host, float64).
// Порт ячейки 0 Multi_antenna_3D_clean_v1.ipynb (segmented_interleaved тайминг,
// frequency-domain коррекция движения, k0-IFFT, нормировка пар).
#pragma once
#include <map>
#include <string>
@@ -12,15 +14,34 @@ namespace gpr {
using PairKey = std::pair<int, int>;
// Тайминги пары (нужны slot/n_slots/t_center для фазовой коррекции).
// Тайминг пары в модели segmented_interleaved: пара живёт в сегменте своего Tx
// (окно [segment_start, segment_start+segment_time] от начала акквизиции),
// внутри сегмента передатчики чередуются по частоте (slot из n_slots).
struct PairTiming {
int tx_interleave_slot = 0;
int tx_interleave_slots = 1;
double t_center_s = 0.0;
double segment_start_s = 0.0;
double segment_time_s = 0.0;
double t_center_s = 0.0; // центр окна сегмента
// Грубый межсегментный сдвиг: dtau центра сегмента относительно центра
// акквизиции кадра (2*dz/v, dz = dir*speed*dt*cos(look)).
double dtau_pair_s = 0.0;
};
using PairTimestamps = std::map<PairKey, PairTiming>;
PairTimestamps compute_pair_timestamps(const MotionConfig&, const GeometryConfig&);
// Окна захвата по Tx (output_pos → [start,end] сек от начала акквизиции).
using TxWindows = std::map<int, std::pair<double, double>>;
// Fallback без per-trace меток: сегменты по tx_per_segment штук делят кадр
// [0, acq_s] поровну (медианная модель ноутбука compute_pair_timestamps).
TxWindows uniform_tx_windows(const std::vector<int>& tx_order, double acq_s, int tx_per_segment);
// Порт compute_pair_timestamps_for_snapshot: Tx с одинаковым окном — один
// сегмент (сегменты по возрастанию окна), слот — порядок в tx_order (полный
// аппаратный порядок из run.combos; выбранные каналы его не меняют).
PairTimestamps compute_pair_timestamps(const TxWindows& tx_windows_s, double acq_s,
const std::vector<int>& tx_order,
const GeometryConfig&, const MotionConfig&);
// A-сканы одного снапшота, общая ось времени у всех пар (assert).
struct AScanSet {
@@ -32,7 +53,7 @@ struct AScanSet {
const c128* pair_ptr(int p) const { return H.data() + (size_t)p * n_fft; }
};
// Порт compute_pair_ascans: bg-subtract → phase correction → k0-IFFT → норма.
// Порт sar_pair_ascans: bg-subtract → phase correction → k0-IFFT → норма.
// s21/freq/background — по ключам (i_tx,i_rx); background может быть nullptr.
AScanSet preprocess_snapshot(
const std::map<PairKey, std::vector<c128>>& s21,
@@ -42,10 +63,12 @@ AScanSet preprocess_snapshot(
const PairTimestamps&, bool normalize = true);
// Отдельные шаги (для тестов/переиспользования):
// Полная фазовая модель движения до IFFT: dtau_intra(f) + dtau_pair одной
// фазой exp(sign*j*2*pi*f*dtau_total) (порт apply_intra_sweep_phase_correction).
std::vector<c128> apply_intra_sweep_phase_correction(
const std::vector<c128>& s21, const std::vector<double>& freq_full,
const PairTiming&, const MotionConfig&, double wave_speed,
double f_start, double f_stop, int n_slots_default);
double f_start, double f_stop);
// Возвращает t (ось), h (комплексный A-скан). n_fft через out-параметр.
void compute_ascan_bp(const std::vector<c128>& s21, const std::vector<double>& freq,
@@ -1,4 +1,5 @@
// Конфигурация счёт-ядра. Дефолты 1:1 повторяют gpr_core/config.py.
// Конфигурация счёт-ядра. Дефолты 1:1 повторяют Multi_antenna_3D_clean_v1.ipynb
// (ячейки 0 и 4) — эталонный ноутбук 3D SAR.
#pragma once
#include <algorithm>
#include <array>
@@ -11,13 +12,13 @@ namespace gpr {
struct GeometryConfig {
std::vector<int> input_idx{0, 1, 2, 3};
std::vector<int> output_idx{0, 1};
// (x,y,z) метры; порядок по физическим индексам как в config.py
std::vector<std::array<double, 3>> tx_positions{{0.785, 0.0, 0.0}, {-0.785, 0.0, 0.0}};
std::vector<int> output_idx{0, 1, 2, 3};
// (x,y,z) метры; порядок соответствует input_idx/output_idx
std::vector<std::array<double, 3>> tx_positions{
{-0.88, 0.0, 0.0}, {0.89, 0.0, 0.0}, {-0.33, 0.445, 0.0}, {0.33, 0.445, 0.0}};
std::vector<std::array<double, 3>> rx_positions{
{-0.42, 0.0, 0.0}, {0.48, 0.0, 0.0}, {0.18, 0.0, 0.0}, {-0.12, 0.0, 0.0}};
{0.23, 0.0, 0.0}, {0.54, 0.15, 0.0}, {-0.55, 0.3, 0.0}, {-0.23, 0.0, 0.0}};
double eps_r = 1.0;
double bp_plane_y = 0.0;
// --- производные (заполняются finalize) ---
std::vector<double> x_tx, y_tx, z_tx; // отсортированы по phys-индексу
@@ -50,53 +51,64 @@ struct GeometryConfig {
};
struct BPConfig {
double f_start = 20 * 1e8; // 2 ГГц
double f_start = 2.0e9;
double f_stop = 6.0e9;
double min_depth = 3.2;
double max_depth = 10.0;
double grid_x_pad = 2.0;
double z_min = 0.1;
int nx = 300;
int nz = 300;
double min_depth = 2.5;
double max_depth = 8.0;
int oversample = 8;
bool window = true;
bool window = true; // окно Ханна на полосе
bool pair_normalize = true;
double pair_norm_percentile = 50.0;
double pair_norm_depth_min = 3.2;
double pair_norm_depth_max = 10.0;
double pair_norm_depth_min = 2.5;
double pair_norm_depth_max = 8.0;
double pair_norm_eps = 1e-15;
double comp_range_power = 0.25;
// Компенсация геометрического расхождения 1/(Rtx*Rrx), нормированная на
// точку под виртуальным центром пары на глубине comp_ref_depth НА ПЛОСКОСТИ
// ЗЕМЛИ (высота (-ground_depth + d*sin p)/cos p).
double comp_range_power = 0.3;
double comp_weight_max = 5.0;
double comp_ref_depth = 5.0;
double smooth_sigma = 1.5;
bool bg_subtract = true;
double cf_eps = 1e-12;
};
struct MotionConfig {
double speed_m_s = 0.854;
double look_angle_deg = 6.0;
double pitch_deg = 6.0;
double tx_sweep_time_s = 0.079;
double tx_switch_time_s = 1e-5;
// frequency_timing_mode всегда interleaved_tx_by_frequency
double speed_m_s = 1.45;
double look_angle_deg = 17.5;
// Наклон рамы для ground-координат и смещений кадров; в ноутбуке
// BP_GROUND_PITCH_DEG=None означает «равен look_angle».
double pitch_deg = 17.5;
double direction_sign = +1.0;
// Внутрисвиповая (по частоте) фазовая коррекция движения.
bool apply_freq_phase_correction = true;
// Грубая межсегментная коррекция (dtau сегмента относительно центра кадра),
// вносится той же фазовой моделью до IFFT (режим 'frequency_domain').
bool apply_pair_motion_correction = true;
double motion_phase_sign = -1.0;
double frame_time_s(int n_tx) const {
return n_tx * tx_sweep_time_s + std::max(0, n_tx - 1) * tx_switch_time_s;
}
// Fallback-разбиение кадра на сегменты, когда per-trace меток захвата нет.
int tx_per_segment = 2;
};
struct SARConfig {
int frame_radius = 6;
double edge_frame_weight = 0.75;
std::array<double, 3> half_width_m{1.8, 1.0, 1.8};
std::array<int, 3> grid{150, 120, 160}; // NX, NY, NZ
bool coord_mode_ground = true; // 'ground' vs 'radar_local'
int frame_radius = 0; // 0 = один снимок (физическая вертикальная база)
double edge_frame_weight = 0.8; // вес крайнего кадра апертуры
bool coord_mode_ground = true; // 'ground' vs 'radar_local'
bool use_compensation = true;
bool normalize_each_frame = true;
// Высота НИЖНЕГО ряда антенн над землёй; Y=0 сцены — поверхность земли.
double ground_depth_m = 1.25;
// Явные границы сцены [м]: X — поперёк, Y — высота над землёй, Z — вдаль.
// Воксели с радарной глубиной вне [min_depth, max_depth] всё равно
// обнуляются маской, так что Z-окно шире полосы глубин безопасно.
double x_min_m = -2.0;
double x_max_m = 2.0;
double y_min_m = -0.5;
double y_max_m = 2.0;
double z_min_m = 0.0;
double z_max_m = 8.0;
std::array<double, 3> voxel_m{0.025, 0.02, 0.02}; // шаг (X, Y, Z)
// Финализация объёма: нормировка без верхних/нижних norm_edge_margin_m по Y
// (там ложный гребень постоянной дальности), сглаживание по (Z,X) в каждом
// Y-плане.
double norm_edge_margin_m = 0.3;
double smooth_sigma = 1.5;
};
} // namespace gpr
@@ -18,10 +18,23 @@ std::vector<FrameOffset> frame_offsets(const std::vector<double>& frame_times, d
// Косинусные веса апертуры (порт sar_frame_weights); dt берётся относительно центра.
std::vector<double> frame_weights(const std::vector<FrameOffset>& offsets, double edge_weight);
// Сетка вокруг точки сцены (порт sar_build_grid_scene).
Grid build_grid_scene(const std::array<double, 3>& center_scene, const SARConfig&);
// Высота плоскости земли (радарные оси) на радарной глубине depth:
// y = (-ground_depth + depth*sin(pitch)) / cos(pitch) (порт bp_plane_y).
double ground_plane_y(double depth, double ground_depth_m, double pitch_rad);
// point_to_scene: оси радара → сцена (поворот на pitch, если ground).
// Сетка сцены по явным границам SARConfig (x/y/z_min..max, шаг voxel_m);
// оси — linspace c числом точек round(span/step)+1, как в ноутбуке.
Grid build_full_region_grid(const SARConfig&);
// Ноутбучные границы (sar_build_full_region_grid): X — размах антенн ±margin,
// Z — полоса, покрывающая радарные глубины [min_depth, max_depth] на всех
// высотах Y. Используются parity-инструментом для сверки с эталонными npz.
void notebook_grid_bounds(const GeometryConfig&, const BPConfig&, const SARConfig&,
double pitch_deg, double x_margin_m,
double& x_min, double& x_max, double& z_min, double& z_max);
// point_to_scene: оси радара → сцена (поворот на pitch + подъём на ground_depth,
// если ground; Y=0 сцены — земля).
std::array<double, 3> point_to_scene(double x, double y, double z, const MotionConfig&,
const SARConfig&);
@@ -19,13 +19,24 @@ struct ApertureFrame {
double weight;
};
// Сырой результат объёма (Python делает finish: abs+gate+gaussian+normalize).
// Сырой результат объёма (до finish_volume).
struct VolumeRaw {
std::vector<c64> acc; // [NY*NZ*NX]
std::vector<float> cnt; // [NY*NZ*NX]
int NX = 0, NY = 0, NZ = 0;
};
// Финальный объём (порт sar_finish_volume): |acc/cnt| → маска глубины →
// нормировка по Y-планам без краёв → гаусс по (Z,X) → нормировка → clip.
struct FinishedVolume {
std::vector<float> volume; // [NY*NZ*NX], нормирован на ~[0,1]
double level = 0.0; // масштаб, на который поделён объём
int NX = 0, NY = 0, NZ = 0;
};
FinishedVolume finish_volume(const VolumeRaw&, const Grid&, double pitch_deg,
const BPConfig&, const SARConfig&);
class SarEngine {
public:
SarEngine(const GeometryConfig&, const BPConfig&, const MotionConfig&, const SARConfig&);
+79 -27
View File
@@ -9,19 +9,69 @@
namespace gpr {
PairTimestamps compute_pair_timestamps(const MotionConfig& mot, const GeometryConfig& geom) {
// interleaved_tx_by_frequency: слот Tx = порядок в отсортированном output_idx.
std::vector<int> tx = geom.output_idx, rx = geom.input_idx;
std::sort(tx.begin(), tx.end());
std::sort(rx.begin(), rx.end());
int n_slots = (int)tx.size();
double frame_center = 0.5 * n_slots * mot.tx_sweep_time_s;
TxWindows uniform_tx_windows(const std::vector<int>& tx_order, double acq_s, int tx_per_segment) {
TxWindows windows;
if (tx_order.empty()) return windows;
int per_seg = std::max(1, tx_per_segment);
int n_seg = std::max(1, (int)((tx_order.size() + per_seg - 1) / per_seg));
double t_seg = acq_s / n_seg;
for (size_t k = 0; k < tx_order.size(); ++k) {
int seg = (int)k / per_seg;
windows[tx_order[k]] = {seg * t_seg, seg * t_seg + t_seg};
}
return windows;
}
PairTimestamps compute_pair_timestamps(const TxWindows& tx_windows_s, double acq_s,
const std::vector<int>& tx_order,
const GeometryConfig& geom, const MotionConfig& mot) {
// Аппаратный порядок: сперва tx_order (run.combos), затем прочие ключи окон.
std::vector<int> hw_order;
for (int tx : tx_order)
if (tx_windows_s.count(tx) && std::find(hw_order.begin(), hw_order.end(), tx) == hw_order.end())
hw_order.push_back(tx);
for (auto& [tx, _] : tx_windows_s)
if (std::find(hw_order.begin(), hw_order.end(), tx) == hw_order.end())
hw_order.push_back(tx);
// Tx с одинаковым окном — один сегмент; сегменты по возрастанию окна,
// слот — порядок появления в hw_order (порт группировки ноутбука).
std::vector<std::pair<std::pair<double, double>, std::vector<int>>> groups;
for (int tx : hw_order) {
auto win = tx_windows_s.at(tx);
auto found = std::find_if(groups.begin(), groups.end(),
[&](const auto& g) { return g.first == win; });
if (found == groups.end())
groups.push_back({win, {tx}});
else
found->second.push_back(tx);
}
std::sort(groups.begin(), groups.end(),
[](const auto& a, const auto& b) { return a.first < b.first; });
struct SegSlot { double t0, t1; int slot, n_slots; };
std::map<int, SegSlot> seg_of;
for (auto& [win, txs] : groups)
for (size_t slot = 0; slot < txs.size(); ++slot)
seg_of[txs[slot]] = {win.first, win.second, (int)slot, (int)txs.size()};
double cos_theta = std::cos(mot.look_angle_deg * M_PI / 180.0);
double t_ref = 0.5 * acq_s;
PairTimestamps ts;
for (int slot = 0; slot < (int)tx.size(); ++slot) {
int i_tx = slot; // индекс в отсортированном == логический
for (int j = 0; j < (int)rx.size(); ++j) {
ts[{i_tx, j}] = PairTiming{slot, n_slots, frame_center};
}
for (int i = 0; i < geom.n_tx; ++i) {
auto it = seg_of.find(geom.output_idx[i]);
if (it == seg_of.end())
throw std::runtime_error("pair timestamps: нет окна для Tx");
const SegSlot& s = it->second;
PairTiming timing;
timing.tx_interleave_slot = s.slot;
timing.tx_interleave_slots = s.n_slots;
timing.segment_start_s = s.t0;
timing.segment_time_s = s.t1 - s.t0;
timing.t_center_s = 0.5 * (s.t0 + s.t1);
double dz = mot.direction_sign * mot.speed_m_s * (timing.t_center_s - t_ref) * cos_theta;
timing.dtau_pair_s = 2.0 * dz / geom.v;
for (int j = 0; j < geom.n_rx; ++j) ts[{i, j}] = timing;
}
return ts;
}
@@ -37,7 +87,7 @@ static double median_of_diff(const std::vector<double>& f) {
std::vector<c128> apply_intra_sweep_phase_correction(
const std::vector<c128>& s21, const std::vector<double>& freq_full,
const PairTiming& info, const MotionConfig& cfg, double wave_speed,
double f_start, double f_stop, int n_slots_default) {
double f_start, double f_stop) {
std::vector<c128> out = s21;
int n_full = (int)freq_full.size();
// индексы в полосе
@@ -47,21 +97,23 @@ std::vector<c128> apply_intra_sweep_phase_correction(
if (idx.empty()) throw std::runtime_error("phase corr: пустая полоса");
if (n_full < 2) return out;
int n_slots = info.tx_interleave_slots ? info.tx_interleave_slots : n_slots_default;
// Сегмент проходит все n_full частот за своё окно; на каждой частоте по
// очереди работают его n_slots передатчиков, поэтому событие частоты i для
// слота slot имеет номер (n_slots*i + slot) из (n_slots*n_full).
int n_slots = std::max(1, info.tx_interleave_slots);
int slot = info.tx_interleave_slot;
double dt_base = cfg.tx_sweep_time_s / (n_full - 1);
// t_all по всему свипу; t_center — середина ПОЛНОГО свипа (как в preprocess.py)
double t_first = (double)(n_slots * 0 + slot) * dt_base;
double t_last = (double)(n_slots * (n_full - 1) + slot) * dt_base;
double t_center = 0.5 * (t_first + t_last);
double theta = cfg.look_angle_deg * M_PI / 180.0;
if (!cfg.apply_freq_phase_correction) return out;
double dt_event = info.segment_time_s / ((double)n_slots * n_full);
double cos_theta = std::cos(cfg.look_angle_deg * M_PI / 180.0);
double pair_dtau = cfg.apply_pair_motion_correction ? info.dtau_pair_s : 0.0;
if (!cfg.apply_freq_phase_correction && pair_dtau == 0.0) return out;
for (int i : idx) {
double t_abs = (double)(n_slots * i + slot) * dt_base;
double dt_intra = t_abs - t_center;
double delta_range = cfg.direction_sign * cfg.speed_m_s * std::cos(theta) * dt_intra;
double dtau = 2.0 * delta_range / wave_speed;
double phi = 2.0 * M_PI * freq_full[i] * dtau;
double dtau_intra = 0.0;
if (cfg.apply_freq_phase_correction) {
double t_abs = info.segment_start_s + (double)(n_slots * i + slot) * dt_event;
double dt_intra = t_abs - info.t_center_s;
dtau_intra = 2.0 * cfg.direction_sign * cfg.speed_m_s * cos_theta * dt_intra / wave_speed;
}
double phi = 2.0 * M_PI * freq_full[i] * (dtau_intra + pair_dtau);
double ang = cfg.motion_phase_sign * phi;
out[i] *= c128(std::cos(ang), std::sin(ang));
}
@@ -163,7 +215,7 @@ AScanSet preprocess_snapshot(
auto itt = ts.find(key);
if (itt != ts.end())
s = apply_intra_sweep_phase_correction(s, fr, itt->second, mot, geom.v,
bp.f_start, bp.f_stop, geom.n_tx);
bp.f_start, bp.f_stop);
std::vector<double> t;
std::vector<c128> h;
compute_ascan_bp(s, fr, bp.f_start, bp.f_stop, bp.window, bp.oversample, t, h);
@@ -1,5 +1,6 @@
#include "gpr/geometry.hpp"
#include <algorithm>
#include <cmath>
namespace gpr {
@@ -48,21 +49,56 @@ std::array<double, 3> point_to_scene(double x, double y, double z, const MotionC
const SARConfig& sar) {
if (!sar.coord_mode_ground) return {x, y, z};
double p = mot.pitch_deg * M_PI / 180.0;
double ys = y * std::cos(p) - z * std::sin(p);
double ys = y * std::cos(p) - z * std::sin(p) + sar.ground_depth_m;
double zs = y * std::sin(p) + z * std::cos(p);
return {x, ys, zs};
}
Grid build_grid_scene(const std::array<double, 3>& c, const SARConfig& sar) {
Grid g;
auto lin = [](double a, double b, int n, std::vector<double>& out) {
double ground_plane_y(double depth, double ground_depth_m, double pitch_rad) {
return (-ground_depth_m + depth * std::sin(pitch_rad)) / std::cos(pitch_rad);
}
Grid build_full_region_grid(const SARConfig& sar) {
// numpy: linspace(lo, hi, max(round(span/step)+1, 2)); nearbyint = банковское
// округление python round()
auto axis = [](double lo, double hi, double step, std::vector<double>& out) {
int n = std::max((int)std::nearbyint((hi - lo) / step) + 1, 2);
out.resize(n);
for (int i = 0; i < n; ++i) out[i] = a + (b - a) * i / (double)(n - 1);
for (int i = 0; i < n; ++i) out[i] = lo + (hi - lo) * i / (double)(n - 1);
};
lin(c[0] - sar.half_width_m[0], c[0] + sar.half_width_m[0], sar.grid[0], g.x);
lin(c[1] - sar.half_width_m[1], c[1] + sar.half_width_m[1], sar.grid[1], g.y);
lin(c[2] - sar.half_width_m[2], c[2] + sar.half_width_m[2], sar.grid[2], g.z);
Grid g;
axis(sar.x_min_m, sar.x_max_m, sar.voxel_m[0], g.x);
axis(sar.y_min_m, sar.y_max_m, sar.voxel_m[1], g.y);
axis(sar.z_min_m, sar.z_max_m, sar.voxel_m[2], g.z);
return g;
}
void notebook_grid_bounds(const GeometryConfig& geom, const BPConfig& bp, const SARConfig& sar,
double pitch_deg, double x_margin_m,
double& x_min, double& x_max, double& z_min, double& z_max) {
double x_lo = geom.x_tx[0], x_hi = geom.x_tx[0];
for (double x : geom.x_tx) { x_lo = std::min(x_lo, x); x_hi = std::max(x_hi, x); }
for (double x : geom.x_rx) { x_lo = std::min(x_lo, x); x_hi = std::max(x_hi, x); }
x_min = x_lo - x_margin_m;
x_max = x_hi + x_margin_m;
// Радарная z линейна по Z сцены при фиксированной высоте:
// z_rad(Y, Z=0) = -(Y - ground_depth)*sin(p), z = (depth - z_rad)/cos(p).
double p = pitch_deg * M_PI / 180.0;
double cos_p = std::cos(p), sin_p = std::sin(p);
double z_lo = 0.0, z_hi = 0.0;
bool first = true;
for (double y : {sar.y_min_m, sar.y_max_m}) {
for (double depth : {bp.min_depth, bp.max_depth}) {
double z_rad0 = -(y - sar.ground_depth_m) * sin_p;
double z = (depth - z_rad0) / cos_p;
if (first) { z_lo = z_hi = z; first = false; }
z_lo = std::min(z_lo, z);
z_hi = std::max(z_hi, z);
}
}
z_min = z_lo;
z_max = z_hi;
}
} // namespace gpr
@@ -1,70 +1,81 @@
// Исходник CUDA-ядра для NVRTC (порт bp_aperture_contribution из gpr_core/kernels.py).
// Self-contained complex<float> (без cupy-заголовков).
// Исходник CUDA-ядра для NVRTC (порт sar_frame_contribution из
// Multi_antenna_3D_clean_v1.ipynb). Self-contained complex<float>.
//
// Геометрия (координаты, дальности, tau, вес компенсации) считается в double:
// эталон ноутбука — float64, а в float32 позиция выборки A-скана (тысячи
// отсчётов) теряет ~3e-4 отсчёта и расходится с эталоном на ~1e-4. A-сканы
// остаются float32 — их значения O(1) и линейная интерполяция не теряет.
#pragma once
namespace gpr {
inline const char* APERTURE_KERNEL_SRC = R"CUDA(
struct cf { float x, y; };
__device__ __forceinline__ cf operator+(cf a, cf b){ return {a.x+b.x, a.y+b.y}; }
__device__ __forceinline__ cf operator-(cf a, cf b){ return {a.x-b.x, a.y-b.y}; }
__device__ __forceinline__ cf operator*(float s, cf a){ return {s*a.x, s*a.y}; }
// Весь взвешенный объём апертуры за один запуск: один поток на воксель,
// внутренний цикл по n_contrib = n_frames*n_pairs. acc/cnt пишутся один раз.
// Взвешенный объём апертуры: один поток на воксель, внутренний цикл по
// n_contrib = n_frames*n_pairs. Обрабатывается диапазон вокселей
// [gid_begin, gid_end) — хост запускает ядро чанками, переиспользуя небольшие
// acc/cnt-буферы (полный объём в VRAM не обязан помещаться: на десктопе с
// занятой другими процессами картой свободных может быть <200 МБ).
extern "C" __global__ void bp_aperture_contribution(
const cf* __restrict__ H,
const float t0, const float inv_dt, const int n_fft, const float inv_v,
const float* __restrict__ xt, const float* __restrict__ yt, const float* __restrict__ zt,
const float* __restrict__ xr, const float* __restrict__ yr, const float* __restrict__ zr,
const float* __restrict__ inv_geo_ref, const float* __restrict__ wgt,
const float* __restrict__ gx, const float* __restrict__ gy, const float* __restrict__ gz,
const double t0, const double inv_dt, const int n_fft, const double inv_v,
const double* __restrict__ xt, const double* __restrict__ yt, const double* __restrict__ zt,
const double* __restrict__ xr, const double* __restrict__ yr, const double* __restrict__ zr,
const double* __restrict__ inv_geo_ref, const float* __restrict__ wgt,
const double* __restrict__ gx, const double* __restrict__ gy, const double* __restrict__ gz,
const int NX, const int NY, const int NZ, const int n_contrib,
const float cos_p, const float sin_p,
const float min_depth, const float max_depth,
const float comp_pow, const float w_max, const int use_comp,
const long gid_begin, const long gid_end,
const double cos_p, const double sin_p, const double ground_depth,
const double min_depth, const double max_depth,
const double comp_pow, const double w_max, const int use_comp,
cf* __restrict__ acc, float* __restrict__ cnt)
{
const long total = (long)NX * NY * NZ;
for (long gid = blockIdx.x * blockDim.x + threadIdx.x;
gid < total; gid += (long)gridDim.x * blockDim.x) {
for (long gid = gid_begin + blockIdx.x * blockDim.x + threadIdx.x;
gid < gid_end; gid += (long)gridDim.x * blockDim.x) {
int ix = gid % NX;
int iz = (gid / NX) % NZ;
int iy = gid / ((long)NX * NZ);
float x = gx[ix];
float ys = gy[iy], zs = gz[iz];
float y_rad = ys * cos_p + zs * sin_p;
float z_rad = -ys * sin_p + zs * cos_p;
double x = gx[ix];
double zs = gz[iz];
// Сцена → оси радара: Y=0 сцены — земля, ground_depth ниже нижнего ряда
// антенн (порт sar_scene_to_radar_yz).
double y0 = gy[iy] - ground_depth;
double y_rad = y0 * cos_p + zs * sin_p;
double z_rad = -y0 * sin_p + zs * cos_p;
bool depth_ok = (z_rad >= min_depth) && (z_rad <= max_depth);
cf a = {0.f, 0.f};
double ax = 0.0, ay = 0.0;
float c = 0.f;
for (int i = 0; i < n_contrib; ++i) {
float w = wgt[i];
float dxt = x - xt[i], dyt = y_rad - yt[i], dzt = z_rad - zt[i];
float dxr = x - xr[i], dyr = y_rad - yr[i], dzr = z_rad - zr[i];
float Rt = sqrtf(dxt*dxt + dyt*dyt + dzt*dzt);
float Rr = sqrtf(dxr*dxr + dyr*dyr + dzr*dzr);
float pos = ((Rt + Rr) * inv_v - t0) * inv_dt;
bool inside = (pos >= 0.f) && (pos <= (float)(n_fft - 1));
if (!inside) continue;
int k = (int)floorf(pos);
double dxt = x - xt[i], dyt = y_rad - yt[i], dzt = z_rad - zt[i];
double dxr = x - xr[i], dyr = y_rad - yr[i], dzr = z_rad - zr[i];
double Rt = sqrt(dxt*dxt + dyt*dyt + dzt*dzt);
double Rr = sqrt(dxr*dxr + dyr*dyr + dzr*dzr);
double pos = ((Rt + Rr) * inv_v - t0) * inv_dt;
bool inside = (pos >= 0.0) && (pos <= (double)(n_fft - 1));
if (!inside || !depth_ok) continue;
int k = (int)floor(pos);
if (k < 0) k = 0;
if (k > n_fft - 2) k = n_fft - 2;
float f = pos - (float)k;
float f = (float)(pos - (double)k);
const cf* h = H + (long)i * n_fft;
cf hv = h[k] + f * (h[k+1] - h[k]);
float hx = h[k].x + f * (h[k+1].x - h[k].x);
float hy = h[k].y + f * (h[k+1].y - h[k].y);
double w = (double)wgt[i];
if (use_comp) {
float ratio = (1.f/(Rt*Rr + 1e-12f) + 1e-30f) * inv_geo_ref[i];
float denom = (comp_pow == 0.25f) ? sqrtf(sqrtf(ratio)) : powf(ratio, comp_pow);
float ww = 1.f/(denom + 1e-12f);
ww = fminf(fmaxf(ww, 0.f), w_max);
hv = ww * hv;
double ratio = (1.0/(Rt*Rr + 1e-12) + 1e-30) * inv_geo_ref[i];
double ww = 1.0/(pow(ratio, comp_pow) + 1e-12);
ww = fmin(fmax(ww, 0.0), w_max);
w *= ww;
}
if (depth_ok) { a = a + (w * hv); c += w; }
ax += w * (double)hx;
ay += w * (double)hy;
c += wgt[i];
}
acc[gid] = a;
cnt[gid] = c;
acc[gid - gid_begin].x = (float)ax;
acc[gid - gid_begin].y = (float)ay;
cnt[gid - gid_begin] = c;
}
}
)CUDA";
+162 -40
View File
@@ -31,11 +31,18 @@ VolumeRaw SarEngine::aperture_volume(const std::vector<ApertureFrame>& frames, c
int NX = grid.NX(), NY = grid.NY(), NZ = grid.NZ();
size_t nvox = (size_t)NX * NY * NZ;
double pitch = mot_.pitch_deg * M_PI / 180.0;
// Опорное расхождение пары — на плоскости земли на глубине comp_ref_depth.
double geo_ref_y = sar_.coord_mode_ground
? ground_plane_y(bp_.comp_ref_depth, sar_.ground_depth_m, pitch)
: 0.0;
// Плоские массивы по (кадр, пара)
std::vector<c64> H_stack((size_t)n_contrib * n_fft);
std::vector<float> xt(n_contrib), yt(n_contrib), zt(n_contrib);
std::vector<float> xr(n_contrib), yr(n_contrib), zr(n_contrib);
std::vector<float> igr(n_contrib), wgt(n_contrib);
std::vector<double> xt(n_contrib), yt(n_contrib), zt(n_contrib);
std::vector<double> xr(n_contrib), yr(n_contrib), zr(n_contrib);
std::vector<double> igr(n_contrib);
std::vector<float> wgt(n_contrib);
int idx = 0;
for (int f = 0; f < n_frames; ++f) {
const AScanSet& a = *frames[f].ascans;
@@ -49,62 +56,177 @@ VolumeRaw SarEngine::aperture_volume(const std::vector<ApertureFrame>& frames, c
const c128* src = a.pair_ptr(p);
c64* dst = H_stack.data() + (size_t)idx * n_fft;
for (int k = 0; k < n_fft; ++k) dst[k] = c64((float)src[k].real(), (float)src[k].imag());
xt[idx] = (float)(geom_.x_tx[i] + dx); yt[idx] = (float)(geom_.y_tx[i] + dy);
zt[idx] = (float)(geom_.z_tx[i] + dz);
xr[idx] = (float)(geom_.x_rx[j] + dx); yr[idx] = (float)(geom_.y_rx[j] + dy);
zr[idx] = (float)(geom_.z_rx[j] + dz);
igr[idx] = (float)(1.0 / pair_geo_ref(geom_, bp_, i, j, geom_.bp_plane_y));
xt[idx] = geom_.x_tx[i] + dx; yt[idx] = geom_.y_tx[i] + dy;
zt[idx] = geom_.z_tx[i] + dz;
xr[idx] = geom_.x_rx[j] + dx; yr[idx] = geom_.y_rx[j] + dy;
zr[idx] = geom_.z_rx[j] + dz;
igr[idx] = 1.0 / pair_geo_ref(geom_, bp_, i, j, geo_ref_y);
wgt[idx] = w;
++idx;
}
}
std::vector<float> gx(grid.x.begin(), grid.x.end());
std::vector<float> gy(grid.y.begin(), grid.y.end());
std::vector<float> gz(grid.z.begin(), grid.z.end());
float cos_p = 1.f, sin_p = 0.f;
double cos_p = 1.0, sin_p = 0.0, ground_depth = 0.0;
if (sar_.coord_mode_ground) {
double p = mot_.pitch_deg * M_PI / 180.0;
cos_p = (float)std::cos(p); sin_p = (float)std::sin(p);
cos_p = std::cos(pitch);
sin_p = std::sin(pitch);
ground_depth = sar_.ground_depth_m;
}
// Заливка на устройство
DeviceBuffer<c64> dH(H_stack.size()); dH.upload(H_stack.data(), H_stack.size());
DeviceBuffer<float> dxt(n_contrib), dyt(n_contrib), dzt(n_contrib);
DeviceBuffer<float> dxr(n_contrib), dyr(n_contrib), dzr(n_contrib);
DeviceBuffer<float> digr(n_contrib), dwgt(n_contrib);
DeviceBuffer<double> dxt(n_contrib), dyt(n_contrib), dzt(n_contrib);
DeviceBuffer<double> dxr(n_contrib), dyr(n_contrib), dzr(n_contrib);
DeviceBuffer<double> digr(n_contrib);
DeviceBuffer<float> dwgt(n_contrib);
dxt.upload(xt.data(), n_contrib); dyt.upload(yt.data(), n_contrib); dzt.upload(zt.data(), n_contrib);
dxr.upload(xr.data(), n_contrib); dyr.upload(yr.data(), n_contrib); dzr.upload(zr.data(), n_contrib);
digr.upload(igr.data(), n_contrib); dwgt.upload(wgt.data(), n_contrib);
DeviceBuffer<float> dgx(NX), dgy(NY), dgz(NZ);
dgx.upload(gx.data(), NX); dgy.upload(gy.data(), NY); dgz.upload(gz.data(), NZ);
DeviceBuffer<c64> dacc(nvox);
DeviceBuffer<float> dcnt(nvox);
DeviceBuffer<double> dgx(NX), dgy(NY), dgz(NZ);
dgx.upload(grid.x.data(), NX); dgy.upload(grid.y.data(), NY); dgz.upload(grid.z.data(), NZ);
float f_t0 = (float)t0, f_inv_dt = (float)(1.0 / dt), f_inv_v = (float)(1.0 / geom_.v);
float min_d = (float)bp_.min_depth, max_d = (float)bp_.max_depth;
float comp_pow = (float)bp_.comp_range_power, w_max = (float)bp_.comp_weight_max;
// Объём считается чанками с переиспользуемыми acc/cnt-буферами: полные
// буферы (~115 МБ на сетку 9.6 млн вокселей) могут не влезать в остаток
// VRAM, когда карту делят десктоп-процессы.
const size_t chunk_vox = std::min<size_t>(nvox, 2'000'000);
DeviceBuffer<c64> dacc(chunk_vox);
DeviceBuffer<float> dcnt(chunk_vox);
double inv_dt = 1.0 / dt, inv_v = 1.0 / geom_.v;
double min_d = bp_.min_depth, max_d = bp_.max_depth;
double comp_pow = bp_.comp_range_power, w_max = bp_.comp_weight_max;
int use_comp = sar_.use_compensation ? 1 : 0;
void* args[] = {
dH.addr(), &f_t0, &f_inv_dt, &n_fft, &f_inv_v,
dxt.addr(), dyt.addr(), dzt.addr(), dxr.addr(), dyr.addr(), dzr.addr(),
digr.addr(), dwgt.addr(), dgx.addr(), dgy.addr(), dgz.addr(),
&NX, &NY, &NZ, &n_contrib, &cos_p, &sin_p, &min_d, &max_d,
&comp_pow, &w_max, &use_comp, dacc.addr(), dcnt.addr()};
int block = 256;
int gridDim = (int)std::min<long>(65535, ((long)nvox + block - 1) / block);
cu_check(cuLaunchKernel(kernel_->func(), gridDim, 1, 1, block, 1, 1, 0, 0, args, nullptr),
"launch aperture");
cu_check(cuCtxSynchronize(), "sync");
VolumeRaw out;
out.NX = NX; out.NY = NY; out.NZ = NZ;
out.acc.resize(nvox); out.cnt.resize(nvox);
dacc.download(out.acc.data(), nvox);
dcnt.download(out.cnt.data(), nvox);
int block = 256;
for (size_t begin = 0; begin < nvox; begin += chunk_vox) {
const size_t end = std::min(nvox, begin + chunk_vox);
long gid_begin = (long)begin, gid_end = (long)end;
void* args[] = {
dH.addr(), &t0, &inv_dt, &n_fft, &inv_v,
dxt.addr(), dyt.addr(), dzt.addr(), dxr.addr(), dyr.addr(), dzr.addr(),
digr.addr(), dwgt.addr(), dgx.addr(), dgy.addr(), dgz.addr(),
&NX, &NY, &NZ, &n_contrib, &gid_begin, &gid_end,
&cos_p, &sin_p, &ground_depth, &min_d, &max_d,
&comp_pow, &w_max, &use_comp, dacc.addr(), dcnt.addr()};
int gridDim = (int)std::min<long>(65535, ((long)(end - begin) + block - 1) / block);
cu_check(cuLaunchKernel(kernel_->func(), gridDim, 1, 1, block, 1, 1, 0, 0, args, nullptr),
"launch aperture");
cu_check(cuCtxSynchronize(), "sync");
dacc.download(out.acc.data() + begin, end - begin);
dcnt.download(out.cnt.data() + begin, end - begin);
}
return out;
}
namespace {
// scipy.ndimage.gaussian_filter: ядро радиуса int(truncate*sigma+0.5), reflect.
std::vector<double> gaussian_kernel_1d(double sigma, int radius) {
std::vector<double> w(2 * radius + 1);
double sum = 0.0;
for (int i = -radius; i <= radius; ++i) {
w[i + radius] = std::exp(-0.5 * i * i / (sigma * sigma));
sum += w[i + radius];
}
for (auto& v : w) v /= sum;
return w;
}
// Корреляция вдоль оси с reflect-границами ('d c b a | a b c d').
inline int reflect_index(int i, int n) {
while (i < 0 || i >= n) {
if (i < 0) i = -i - 1;
if (i >= n) i = 2 * n - i - 1;
}
return i;
}
} // namespace
FinishedVolume finish_volume(const VolumeRaw& raw, const Grid& grid, double pitch_deg,
const BPConfig& bp, const SARConfig& sar) {
int NX = raw.NX, NY = raw.NY, NZ = raw.NZ;
size_t nvox = (size_t)NX * NY * NZ;
FinishedVolume out;
out.NX = NX; out.NY = NY; out.NZ = NZ;
std::vector<double> vol(nvox);
for (size_t i = 0; i < nvox; ++i)
vol[i] = std::abs(c128(raw.acc[i].real(), raw.acc[i].imag())) / ((double)raw.cnt[i] + 1e-12);
// Маска глубины в радарных осях (в ядре невалидные вкладов не получали,
// но |0/(0+eps)| уже 0; обнуление здесь повторяет ноутбук буквально).
if (sar.coord_mode_ground) {
double p = pitch_deg * M_PI / 180.0;
double cos_p = std::cos(p), sin_p = std::sin(p);
for (int iy = 0; iy < NY; ++iy) {
double y0 = grid.y[iy] - sar.ground_depth_m;
for (int iz = 0; iz < NZ; ++iz) {
double z_rad = -y0 * sin_p + grid.z[iz] * cos_p;
if (z_rad >= bp.min_depth && z_rad <= bp.max_depth) continue;
double* row = vol.data() + ((size_t)iy * NZ + iz) * NX;
std::fill(row, row + NX, 0.0);
}
}
}
// Плоскости Y для нормировки: без верхних/нижних norm_edge_margin_m.
std::vector<bool> ref_plane(NY, true);
double margin = sar.norm_edge_margin_m;
double y_span = grid.y.back() - grid.y.front();
if (margin > 0.0 && y_span > 2.5 * margin) {
for (int iy = 0; iy < NY; ++iy)
ref_plane[iy] = grid.y[iy] >= grid.y.front() + margin &&
grid.y[iy] <= grid.y.back() - margin;
}
auto ref_max = [&](const std::vector<double>& v) {
double m = 0.0;
for (int iy = 0; iy < NY; ++iy) {
if (!ref_plane[iy]) continue;
const double* plane = v.data() + (size_t)iy * NZ * NX;
for (size_t i = 0; i < (size_t)NZ * NX; ++i) m = std::max(m, plane[i]);
}
return m + 1e-12;
};
double level_raw = ref_max(vol);
for (auto& v : vol) v /= level_raw;
// Гаусс sigma=(0, s, s): сглаживание по Z и X внутри каждого Y-плана
// (порядок осей как у scipy: сначала ось 1, затем ось 2).
if (sar.smooth_sigma > 0.0) {
int radius = (int)(4.0 * sar.smooth_sigma + 0.5);
const auto kernel = gaussian_kernel_1d(sar.smooth_sigma, radius);
std::vector<double> tmp((size_t)NZ * NX);
for (int iy = 0; iy < NY; ++iy) {
double* plane = vol.data() + (size_t)iy * NZ * NX;
// вдоль Z (ось 1)
for (int iz = 0; iz < NZ; ++iz)
for (int ix = 0; ix < NX; ++ix) {
double s = 0.0;
for (int k = -radius; k <= radius; ++k)
s += kernel[k + radius] * plane[(size_t)reflect_index(iz + k, NZ) * NX + ix];
tmp[(size_t)iz * NX + ix] = s;
}
// вдоль X (ось 2)
for (int iz = 0; iz < NZ; ++iz)
for (int ix = 0; ix < NX; ++ix) {
double s = 0.0;
for (int k = -radius; k <= radius; ++k)
s += kernel[k + radius] * tmp[(size_t)iz * NX + reflect_index(ix + k, NX)];
plane[(size_t)iz * NX + ix] = s;
}
}
}
double level_smooth = ref_max(vol);
out.volume.resize(nvox);
for (size_t i = 0; i < nvox; ++i)
out.volume[i] = (float)std::max(vol[i] / level_smooth, 0.0);
out.level = level_raw * level_smooth;
return out;
}
@@ -90,21 +90,47 @@ struct ProcessingLiveConfig {
float gpr_visible_x_max_m = 2.0F;
float gpr_visible_z_min_m = 0.0F;
float gpr_visible_z_max_m = 14.0F;
// 3D SAR (sar3d processor) parameters; defaults mirror VNA_AlGo gpr_core.
// Frequency band, depth gate, antenna selection and motion (speed, look
// angle, direction, phase correction) are shared with the gpr_* fields.
std::uint32_t sar3d_frame_radius = 6U;
float sar3d_edge_frame_weight = 0.75F;
float sar3d_pitch_deg = 6.0F;
float sar3d_center_depth_m = 5.6F;
float sar3d_half_width_x_m = 1.8F;
float sar3d_half_width_y_m = 1.0F;
float sar3d_half_width_z_m = 1.8F;
std::uint32_t sar3d_grid_nx = 150U;
std::uint32_t sar3d_grid_ny = 120U;
std::uint32_t sar3d_grid_nz = 160U;
// 3D SAR (sar3d processor) parameters; defaults mirror the reference
// notebook Multi_antenna_3D_clean_v1.ipynb (cells 0/4). Frequency band,
// depth gate, antenna selection and motion (speed, look angle, direction,
// phase correction) are shared with the gpr_* fields.
std::uint32_t sar3d_frame_radius = 0U;
float sar3d_edge_frame_weight = 0.8F;
// Frame pitch for ground coordinates and inter-frame offsets. Override
// disabled -> the pitch follows gpr_look_angle_deg (notebook
// BP_GROUND_PITCH_DEG = None).
bool sar3d_pitch_override_enabled = false;
float sar3d_pitch_deg = 17.5F;
// Height of the LOWER antenna row above the ground; scene Y=0 is the ground.
float sar3d_ground_depth_m = 1.25F;
// Scene grid bounds (metres): X across, Y height above the ground, Z along
// the ground away from the radar. Voxels whose radar depth falls outside
// [gpr_min_depth, gpr_max_depth] are zeroed by the mask regardless.
// Doubles, not floats: the parity tool feeds the notebook-derived bounds
// through here and a float roundtrip shifts the grid by ~1e-7 m, flipping
// knife-edge voxels of the depth mask against the reference volumes.
double sar3d_x_min_m = -2.0;
double sar3d_x_max_m = 2.0;
double sar3d_y_min_m = -0.5;
double sar3d_y_max_m = 2.0;
double sar3d_z_min_m = 0.0;
double sar3d_z_max_m = 8.0;
float sar3d_voxel_x_m = 0.025F;
float sar3d_voxel_y_m = 0.02F;
float sar3d_voxel_z_m = 0.02F;
// Volume finalization: normalization ignores this much of the grid's Y
// edges (constant-range ridge), then Gaussian smoothing over (Z, X).
float sar3d_norm_edge_margin_m = 0.3F;
float sar3d_smooth_sigma = 1.5F;
// Geometric spreading compensation 1/(Rtx*Rrx)^power, capped.
float sar3d_comp_range_power = 0.3F;
// Fallback frame segmentation when per-trace capture stamps are absent.
std::uint32_t sar3d_tx_per_segment = 2U;
bool sar3d_use_compensation = true;
bool sar3d_normalize_pairs = true;
// Debug: when non-empty, the processor dumps every finished volume to this
// directory (raw float32 + json header). Used by the parity tool only.
std::string sar3d_volume_dump_dir{};
// Bottom-view projection angle around the vertical Y axis, degrees:
// 0 = front view (X horizontal), 90 = side view (Z horizontal),
// intermediate values project along the rotated horizontal direction.
@@ -392,10 +392,13 @@ void apply_legacy_gpr_algorithm_alias(ProcessingLiveConfig& config, const std::s
for (const auto& [key, target] : {
std::pair{"sar3d_edge_frame_weight", &config.sar3d_edge_frame_weight},
std::pair{"sar3d_pitch_deg", &config.sar3d_pitch_deg},
std::pair{"sar3d_center_depth_m", &config.sar3d_center_depth_m},
std::pair{"sar3d_half_width_x_m", &config.sar3d_half_width_x_m},
std::pair{"sar3d_half_width_y_m", &config.sar3d_half_width_y_m},
std::pair{"sar3d_half_width_z_m", &config.sar3d_half_width_z_m},
std::pair{"sar3d_ground_depth_m", &config.sar3d_ground_depth_m},
std::pair{"sar3d_voxel_x_m", &config.sar3d_voxel_x_m},
std::pair{"sar3d_voxel_y_m", &config.sar3d_voxel_y_m},
std::pair{"sar3d_voxel_z_m", &config.sar3d_voxel_z_m},
std::pair{"sar3d_norm_edge_margin_m", &config.sar3d_norm_edge_margin_m},
std::pair{"sar3d_smooth_sigma", &config.sar3d_smooth_sigma},
std::pair{"sar3d_comp_range_power", &config.sar3d_comp_range_power},
std::pair{"sar3d_bottom_view_deg", &config.sar3d_bottom_view_deg},
}) {
if (const auto found = root.find(key); found != root.end()) {
@@ -406,25 +409,37 @@ void apply_legacy_gpr_algorithm_alias(ProcessingLiveConfig& config, const std::s
}
}
for (const auto& [key, target] : {
std::pair{"sar3d_grid_nx", &config.sar3d_grid_nx},
std::pair{"sar3d_grid_ny", &config.sar3d_grid_ny},
std::pair{"sar3d_grid_nz", &config.sar3d_grid_nz},
std::pair{"sar3d_x_min_m", &config.sar3d_x_min_m},
std::pair{"sar3d_x_max_m", &config.sar3d_x_max_m},
std::pair{"sar3d_y_min_m", &config.sar3d_y_min_m},
std::pair{"sar3d_y_max_m", &config.sar3d_y_max_m},
std::pair{"sar3d_z_min_m", &config.sar3d_z_min_m},
std::pair{"sar3d_z_max_m", &config.sar3d_z_max_m},
}) {
if (const auto found = root.find(key); found != root.end()) {
*target = parse_u32_number(*found, std::string("processing.") + key);
if (!found->is_number()) {
throw std::runtime_error(std::string("processing.") + key + " must be number");
}
*target = found->get<double>();
}
}
if (const auto found = root.find("sar3d_use_compensation"); found != root.end()) {
if (!found->is_boolean()) {
throw std::runtime_error("processing.sar3d_use_compensation must be bool");
if (const auto found = root.find("sar3d_tx_per_segment"); found != root.end()) {
config.sar3d_tx_per_segment = parse_u32_number(*found, "processing.sar3d_tx_per_segment");
if (config.sar3d_tx_per_segment == 0U) {
throw std::runtime_error("processing.sar3d_tx_per_segment must be >= 1");
}
config.sar3d_use_compensation = found->get<bool>();
}
if (const auto found = root.find("sar3d_normalize_pairs"); found != root.end()) {
if (!found->is_boolean()) {
throw std::runtime_error("processing.sar3d_normalize_pairs must be bool");
for (const auto& [key, target] : {
std::pair{"sar3d_pitch_override_enabled", &config.sar3d_pitch_override_enabled},
std::pair{"sar3d_use_compensation", &config.sar3d_use_compensation},
std::pair{"sar3d_normalize_pairs", &config.sar3d_normalize_pairs},
}) {
if (const auto found = root.find(key); found != root.end()) {
if (!found->is_boolean()) {
throw std::runtime_error(std::string("processing.") + key + " must be bool");
}
*target = found->get<bool>();
}
config.sar3d_normalize_pairs = found->get<bool>();
}
if (const auto found = root.find("ignore_socket_speed"); found != root.end()) {
if (!found->is_boolean()) {
@@ -4,7 +4,10 @@
#include <bit>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <fstream>
#include <iostream>
#include <limits>
#include <map>
#include <optional>
#include <stdexcept>
@@ -19,7 +22,6 @@ namespace radar::processing {
namespace {
constexpr std::size_t kAscanCacheLimit = 64U;
constexpr double kVolumeEps = 1e-12;
void hash_mix(std::size_t& seed, std::uint64_t value) {
seed ^= value + 0x9e3779b97f4a7c15ULL + (seed << 6U) + (seed >> 2U);
@@ -130,6 +132,7 @@ struct Selection {
bp.pair_normalize = live_config.sar3d_normalize_pairs;
bp.pair_norm_depth_min = bp.min_depth;
bp.pair_norm_depth_max = bp.max_depth;
bp.comp_range_power = static_cast<double>(live_config.sar3d_comp_range_power);
return bp;
}
@@ -137,9 +140,13 @@ struct Selection {
gpr::MotionConfig motion{};
motion.speed_m_s = static_cast<double>(live_config.gpr_speed_m_s);
motion.look_angle_deg = static_cast<double>(live_config.gpr_look_angle_deg);
motion.pitch_deg = static_cast<double>(live_config.sar3d_pitch_deg);
// Notebook BP_GROUND_PITCH_DEG = None -> the frame pitch follows the look angle.
motion.pitch_deg = live_config.sar3d_pitch_override_enabled
? static_cast<double>(live_config.sar3d_pitch_deg)
: static_cast<double>(live_config.gpr_look_angle_deg);
motion.direction_sign = static_cast<double>(live_config.gpr_direction_sign);
motion.apply_freq_phase_correction = live_config.gpr_apply_freq_phase_correction;
motion.tx_per_segment = static_cast<int>(std::max<std::uint32_t>(1U, live_config.sar3d_tx_per_segment));
return motion;
}
@@ -147,20 +154,38 @@ struct Selection {
gpr::SARConfig sar{};
sar.frame_radius = static_cast<int>(std::max<std::uint32_t>(0U, live_config.sar3d_frame_radius));
sar.edge_frame_weight = static_cast<double>(live_config.sar3d_edge_frame_weight);
sar.half_width_m = {
static_cast<double>(live_config.sar3d_half_width_x_m),
static_cast<double>(live_config.sar3d_half_width_y_m),
static_cast<double>(live_config.sar3d_half_width_z_m),
};
sar.grid = {
static_cast<int>(std::max<std::uint32_t>(2U, live_config.sar3d_grid_nx)),
static_cast<int>(std::max<std::uint32_t>(2U, live_config.sar3d_grid_ny)),
static_cast<int>(std::max<std::uint32_t>(2U, live_config.sar3d_grid_nz)),
sar.ground_depth_m = static_cast<double>(live_config.sar3d_ground_depth_m);
sar.x_min_m = static_cast<double>(live_config.sar3d_x_min_m);
sar.x_max_m = static_cast<double>(live_config.sar3d_x_max_m);
sar.y_min_m = static_cast<double>(live_config.sar3d_y_min_m);
sar.y_max_m = static_cast<double>(live_config.sar3d_y_max_m);
sar.z_min_m = static_cast<double>(live_config.sar3d_z_min_m);
sar.z_max_m = static_cast<double>(live_config.sar3d_z_max_m);
sar.voxel_m = {
static_cast<double>(live_config.sar3d_voxel_x_m),
static_cast<double>(live_config.sar3d_voxel_y_m),
static_cast<double>(live_config.sar3d_voxel_z_m),
};
sar.norm_edge_margin_m = static_cast<double>(live_config.sar3d_norm_edge_margin_m);
sar.smooth_sigma = static_cast<double>(live_config.sar3d_smooth_sigma);
sar.use_compensation = live_config.sar3d_use_compensation;
return sar;
}
// Hardware Tx order: first occurrence of each output in run_combos. Deliberately
// ignores the live Tx selection so segment/slot assignment stays a hardware
// property (notebook tx_order_full).
[[nodiscard]] auto tx_order_full(const config::RunConfig& run_config) -> std::vector<int> {
std::vector<int> order{};
for (const auto& combo : run_config.run_combos) {
const int output = static_cast<int>(combo.output_pos);
if (std::find(order.begin(), order.end(), output) == order.end()) {
order.push_back(output);
}
}
return order;
}
[[nodiscard]] auto ascan_fingerprint(
const Selection& selection,
const gpr::GeometryConfig& geom,
@@ -184,6 +209,7 @@ struct Selection {
hash_mix(seed, motion.look_angle_deg);
hash_mix(seed, motion.direction_sign);
hash_mix(seed, static_cast<std::uint64_t>(motion.apply_freq_phase_correction ? 1U : 0U));
hash_mix(seed, static_cast<std::uint64_t>(motion.tx_per_segment));
return seed;
}
@@ -215,6 +241,7 @@ struct Selection {
hash_mix(seed, bp.comp_weight_max);
hash_mix(seed, bp.comp_ref_depth);
hash_mix(seed, motion.pitch_deg);
hash_mix(seed, sar.ground_depth_m);
hash_mix(seed, static_cast<std::uint64_t>(sar.use_compensation ? 1U : 0U));
return seed;
}
@@ -228,14 +255,80 @@ struct Selection {
return static_cast<double>(collection.monotonic_ns) * 1e-9;
}
// One Tx's interleaved sweep spans the whole capture window, so the per-Tx
// sweep time is capture_span / N_tx (same derivation as the coherent BP path).
[[nodiscard]] auto tx_sweep_time_s(const ipc::PreprocessedCollection& collection, int n_tx) -> double {
if (collection.capture_end_ns <= collection.capture_start_ns || n_tx <= 0) {
return 0.0;
// Real frame timing from the per-trace capture stamps: traces measured in one
// segment share one window on both producers, so the first window seen per Tx
// reconstructs the segments exactly (notebook snapshot_acquisition_timing).
// Invalid when any trace carries the zero "unmeasured" sentinel.
struct SnapshotTiming {
gpr::TxWindows windows{}; // seconds relative to the first trace start
double acq_s = 0.0;
double center_abs_s = 0.0; // absolute acquisition centre (frame time)
bool valid = false;
};
[[nodiscard]] auto snapshot_timing(const ipc::PreprocessedCollection& collection) -> SnapshotTiming {
SnapshotTiming timing{};
if (collection.traces.empty()) {
return timing;
}
const double span_s = static_cast<double>(collection.capture_end_ns - collection.capture_start_ns) * 1e-9;
return span_s / static_cast<double>(n_tx);
std::map<int, std::pair<std::uint64_t, std::uint64_t>> by_tx{};
std::uint64_t first = std::numeric_limits<std::uint64_t>::max();
std::uint64_t last = 0U;
for (const auto& trace : collection.traces) {
if (trace.capture_end_ns <= trace.capture_start_ns) {
return timing;
}
by_tx.try_emplace(
static_cast<int>(trace.combo.output_pos),
std::pair{trace.capture_start_ns, trace.capture_end_ns}
);
first = std::min(first, trace.capture_start_ns);
last = std::max(last, trace.capture_end_ns);
}
for (const auto& [tx, window] : by_tx) {
timing.windows[tx] = {
static_cast<double>(window.first - first) * 1e-9,
static_cast<double>(window.second - first) * 1e-9,
};
}
timing.acq_s = static_cast<double>(last - first) * 1e-9;
timing.center_abs_s = 0.5 * (static_cast<double>(first) + static_cast<double>(last)) * 1e-9;
timing.valid = true;
return timing;
}
// Debug dump for the parity tool: raw float32 volume [iy][iz][ix] plus a json
// header with the axes and level.
void dump_volume(
const std::string& dir,
std::uint64_t collection_id,
const gpr::FinishedVolume& finished,
const gpr::Grid& grid
) {
char name[64];
std::snprintf(name, sizeof(name), "/sar3d_volume_%llu",
static_cast<unsigned long long>(collection_id));
const std::string base = dir + name;
std::ofstream bin(base + ".f32", std::ios::binary | std::ios::trunc);
bin.write(
reinterpret_cast<const char*>(finished.volume.data()),
static_cast<std::streamsize>(finished.volume.size() * sizeof(float))
);
std::ofstream meta(base + ".json", std::ios::trunc);
meta.precision(17);
meta << "{\"nx\": " << finished.NX << ", \"ny\": " << finished.NY
<< ", \"nz\": " << finished.NZ << ", \"level\": " << finished.level;
const auto write_axis = [&meta](const char* key, const std::vector<double>& axis) {
meta << ", \"" << key << "\": [";
for (std::size_t index = 0U; index < axis.size(); ++index) {
meta << (index > 0U ? "," : "") << axis[index];
}
meta << "]";
};
write_axis("x_m", grid.x);
write_axis("y_m", grid.y);
write_axis("z_m", grid.z);
meta << "}\n";
}
[[nodiscard]] auto build_image_payload(
@@ -285,14 +378,26 @@ struct Sar3dProcessor::State {
const Selection& selection,
const gpr::GeometryConfig& geom,
const gpr::BPConfig& bp,
gpr::MotionConfig motion
const gpr::MotionConfig& motion,
const std::vector<int>& tx_order
) -> const CacheEntry& {
if (const auto found = ascan_cache.find(collection.collection_id); found != ascan_cache.end()) {
return found->second;
}
CacheEntry entry{};
entry.time_s = collection_time_s(collection);
auto timing = snapshot_timing(collection);
// The segmented model needs a window for every selected Tx; fall back to
// the uniform split otherwise (notebook median model).
if (timing.valid) {
for (const auto output : geom.output_idx) {
if (timing.windows.find(output) == timing.windows.end()) {
timing.valid = false;
break;
}
}
}
entry.time_s = timing.valid ? timing.center_abs_s : collection_time_s(collection);
std::map<gpr::PairKey, std::vector<gpr::c128>> s21{};
std::map<gpr::PairKey, std::vector<double>> freq{};
@@ -317,9 +422,25 @@ struct Sar3dProcessor::State {
}
if (s21.size() == selection.inputs.size() * selection.outputs.size()) {
motion.tx_sweep_time_s = tx_sweep_time_s(collection, geom.n_tx);
const auto timestamps = gpr::compute_pair_timestamps(motion, geom);
gpr::TxWindows windows{};
double acq_s = 0.0;
if (timing.valid) {
windows = std::move(timing.windows);
acq_s = timing.acq_s;
} else {
if (collection.capture_end_ns > collection.capture_start_ns) {
acq_s = static_cast<double>(collection.capture_end_ns - collection.capture_start_ns) * 1e-9;
}
auto order = tx_order;
for (const auto output : geom.output_idx) {
if (std::find(order.begin(), order.end(), output) == order.end()) {
order.push_back(output);
}
}
windows = gpr::uniform_tx_windows(order, acq_s, motion.tx_per_segment);
}
try {
const auto timestamps = gpr::compute_pair_timestamps(windows, acq_s, tx_order, geom, motion);
entry.ascans = std::make_shared<gpr::AScanSet>(gpr::preprocess_snapshot(
s21, freq, /*background=*/nullptr, motion, geom, bp, timestamps, bp.pair_normalize
));
@@ -399,12 +520,14 @@ auto Sar3dProcessor::process_collection(
}
window.push_back(&collection);
const auto tx_order = tx_order_full(run_config);
std::vector<std::shared_ptr<gpr::AScanSet>> ascan_sets{};
std::vector<double> frame_times{};
ascan_sets.reserve(window.size());
frame_times.reserve(window.size());
for (const auto* entry_collection : window) {
const auto& entry = state_->get_or_compute_ascans(*entry_collection, selection, geom, bp, motion);
const auto& entry =
state_->get_or_compute_ascans(*entry_collection, selection, geom, bp, motion, tx_order);
if (entry.ascans == nullptr) {
continue;
}
@@ -450,14 +573,17 @@ auto Sar3dProcessor::process_collection(
state_->engine_fp = engine_fp;
}
const auto center_scene =
gpr::point_to_scene(0.0, 0.0, static_cast<double>(live_config.sar3d_center_depth_m), motion, sar);
const gpr::Grid grid = gpr::build_grid_scene(center_scene, sar);
const gpr::VolumeRaw volume = state_->engine->aperture_volume(frames, grid);
const gpr::Grid grid = gpr::build_full_region_grid(sar);
const gpr::VolumeRaw raw = state_->engine->aperture_volume(frames, grid);
// Full notebook finalization: |acc/cnt|, depth gate, normalization
// without the Y-edge planes, Gaussian smoothing over (Z, X), clip.
const gpr::FinishedVolume finished =
gpr::finish_volume(raw, grid, motion.pitch_deg, bp, sar);
if (!live_config.sar3d_volume_dump_dir.empty()) {
dump_volume(live_config.sar3d_volume_dump_dir, collection.collection_id, finished, grid);
}
// Finalization (port of gpr_core sar_finish_volume, minus smoothing):
// |acc/(cnt+eps)|, ground-coordinate depth gate, normalization to peak.
//
// Two projections match the legacy render conventions (X across,
// Z range/"вдаль", Y up):
// - top view `sar3d_proj_xz`: max over Y, axes X (horizontal) / Z;
@@ -466,12 +592,9 @@ auto Sar3dProcessor::process_collection(
// `sar3d_bottom_view_deg` around Y: 0° = front view (u = X, max
// along Z), 90° = side view (u = Z, max along X), intermediate
// angles project along the rotated perpendicular.
const int nx = volume.NX;
const int ny = volume.NY;
const int nz = volume.NZ;
const double pitch_rad = motion.pitch_deg * M_PI / 180.0;
const double cos_p = std::cos(pitch_rad);
const double sin_p = std::sin(pitch_rad);
const int nx = finished.NX;
const int ny = finished.NY;
const int nz = finished.NZ;
const double theta_deg = std::clamp(static_cast<double>(live_config.sar3d_bottom_view_deg), 0.0, 90.0);
const double theta_rad = theta_deg * M_PI / 180.0;
@@ -501,24 +624,18 @@ auto Sar3dProcessor::process_collection(
const double u_start = u_axis.front();
const double inv_du = nu > 1 ? static_cast<double>(nu - 1) / (u_axis.back() - u_start) : 0.0;
// The volume is already normalized to ~[0, 1] by finish_volume, so the
// projections are published as-is; the scalar payload carries `level`
// (the absolute scale the volume was divided by).
std::vector<float> projection_xz(static_cast<std::size_t>(nz) * nx, 0.0F);
std::vector<float> projection_bottom(static_cast<std::size_t>(ny) * nu, 0.0F);
float peak = 0.0F;
for (int iy = 0; iy < ny; ++iy) {
for (int iz = 0; iz < nz; ++iz) {
const double z_rad = -grid.y[iy] * sin_p + grid.z[iz] * cos_p;
const bool depth_ok = z_rad >= bp.min_depth && z_rad <= bp.max_depth;
if (!depth_ok) {
continue;
}
const std::size_t row = (static_cast<std::size_t>(iy) * nz + iz) * nx;
const std::size_t bottom_row = static_cast<std::size_t>(iy) * nu;
const double z_sin_t = grid.z[iz] * sin_t;
for (int ix = 0; ix < nx; ++ix) {
const std::size_t voxel = row + static_cast<std::size_t>(ix);
const float value = std::abs(volume.acc[voxel]) /
(volume.cnt[voxel] + static_cast<float>(kVolumeEps));
peak = std::max(peak, value);
const float value = finished.volume[row + static_cast<std::size_t>(ix)];
auto& xz_cell = projection_xz[static_cast<std::size_t>(iz) * nx + ix];
xz_cell = std::max(xz_cell, value);
int bin = ix;
@@ -534,14 +651,6 @@ auto Sar3dProcessor::process_collection(
}
}
}
if (peak > 0.0F) {
for (auto& value : projection_xz) {
value /= peak;
}
for (auto& value : projection_bottom) {
value /= peak;
}
}
results.collection_payloads.push_back(
build_image_payload("sar3d_proj_xz", grid.x, grid.z, std::move(projection_xz))
@@ -552,7 +661,7 @@ auto Sar3dProcessor::process_collection(
ipc::ResultPayload peak_payload{};
peak_payload.processing_name = "sar3d_peak";
peak_payload.kind = ipc::ResultKind::ScalarF32;
peak_payload.scalar_value = peak;
peak_payload.scalar_value = static_cast<float>(finished.level);
results.collection_payloads.push_back(std::move(peak_payload));
} catch (const std::exception& exc) {
// A CUDA/NVRTC or data error must not kill the data_processor loop:
@@ -1,12 +1,18 @@
// Parity check for the sar3d processor against the VNA_AlGo reference path.
// Parity check of the sar3d processor against the reference notebook
// Multi_antenna_3D_clean_v1.ipynb.
//
// Loads the same .npy snapshot dataset the gpr_cpp `standalone_bp` example
// uses, feeds it through Sar3dProcessor (combo mapping, timestamps, A-scan
// cache, aperture, finalization) and prints the volume peak and the argmax of
// both projections. Compare against `standalone_bp <dir> <center>` output:
// the peak |acc/cnt| and the (X, Y, Z) of the maximum voxel must match.
// Feeds a preprocessed snapshot of an extracted replay dataset through
// Sar3dProcessor exactly the way the live pipeline would (geometry and combo
// order from the dataset's config_profile.json, per-trace capture windows from
// the snapshot's meta.json) and dumps the finished volume for comparison with
// the notebook's *_volume.npz export:
//
// Usage: sar3d_parity <preprocessed_dir> <center_snapshot_name>
// sar3d_parity <dataset_root> <center_snapshot> [--dump-dir <dir>]
// python_app/scripts/sar3d_volume_compare.py <dir>/sar3d_volume_<id> <ref.npz>
//
// <dataset_root> contains config_profile.json and preprocessed/<snapshot>/.
// Defaults mirror the parameters the reference volumes were built with
// (f 2..6 GHz, depth 2.5..8 m, speed 1.45 m/s, look angle 17.5 deg).
#include <algorithm>
#include <cctype>
#include <complex>
@@ -19,6 +25,10 @@
#include <string>
#include <vector>
#include <nlohmann/json.hpp>
#include "gpr/config.hpp"
#include "gpr/geometry.hpp"
#include "processing_live_config.hpp"
#include "run_config.hpp"
#include "sar3d_processor.hpp"
@@ -26,8 +36,7 @@
namespace {
// Minimal numpy .npy loader (v1.0, 1-D complex128/complex64/float64/float32),
// mirroring gpr_cpp/examples/standalone_bp.cpp.
// Minimal numpy .npy loader (v1.0, 1-D complex128/complex64/float64/float32).
struct NpyArray {
std::vector<double> f64{};
std::vector<std::complex<double>> c128{};
@@ -77,43 +86,53 @@ NpyArray load_npy(const std::string& path) {
return array;
}
// VNA_AlGo reference geometry (gpr_core/config.py) keyed by switch positions.
auto build_run_config() -> radar::config::RunConfig {
radar::config::RunConfig config{};
config.gpr.relative_permittivity = 1.0F;
config.gpr.tx_geometry = {
{.output_pos = 0U, .x_m = 0.785F, .y_m = 0.0F, .z_m = 0.0F},
{.output_pos = 1U, .x_m = -0.785F, .y_m = 0.0F, .z_m = 0.0F},
};
config.gpr.rx_geometry = {
{.input_pos = 0U, .x_m = -0.42F, .y_m = 0.0F, .z_m = 0.0F},
{.input_pos = 1U, .x_m = 0.48F, .y_m = 0.0F, .z_m = 0.0F},
{.input_pos = 2U, .x_m = 0.18F, .y_m = 0.0F, .z_m = 0.0F},
{.input_pos = 3U, .x_m = -0.12F, .y_m = 0.0F, .z_m = 0.0F},
};
return config;
}
// Live config matching the VNA_AlGo defaults (BPConfig/MotionConfig/SARConfig).
auto build_live_config() -> radar::processing::ProcessingLiveConfig {
// Live config matching the reference notebook run (npz meta parameters).
// The default program grid is a fixed window (X -2..2, Z 0..8), while the
// notebook derives its bounds from the antenna span and the depth gate — so
// the notebook bounds are recomputed here and passed explicitly.
auto build_live_config(const radar::config::RunConfig& run_config)
-> radar::processing::ProcessingLiveConfig {
radar::processing::ProcessingLiveConfig live{};
live.processor_mode = "sar3d";
live.gpr_start_freq_mhz = 2000.0F;
live.gpr_stop_freq_mhz = 6000.0F;
live.gpr_min_depth_m = 3.2F;
live.gpr_max_depth_m = 10.0F;
live.gpr_speed_m_s = 0.854F;
live.gpr_look_angle_deg = 6.0F;
live.gpr_min_depth_m = 2.5F;
live.gpr_max_depth_m = 8.0F;
live.gpr_speed_m_s = 1.45F;
live.gpr_look_angle_deg = 17.5F;
live.gpr_direction_sign = 1.0F;
live.gpr_apply_freq_phase_correction = true;
gpr::GeometryConfig geom{};
geom.input_idx.clear();
geom.output_idx.clear();
geom.rx_positions.clear();
geom.tx_positions.clear();
for (const auto& entry : run_config.gpr.rx_geometry) {
geom.input_idx.push_back(static_cast<int>(entry.input_pos));
geom.rx_positions.push_back({entry.x_m, entry.y_m, entry.z_m});
}
for (const auto& entry : run_config.gpr.tx_geometry) {
geom.output_idx.push_back(static_cast<int>(entry.output_pos));
geom.tx_positions.push_back({entry.x_m, entry.y_m, entry.z_m});
}
geom.finalize();
gpr::BPConfig bp{};
bp.min_depth = live.gpr_min_depth_m;
bp.max_depth = live.gpr_max_depth_m;
gpr::SARConfig sar{};
double x_min = 0.0, x_max = 0.0, z_min = 0.0, z_max = 0.0;
gpr::notebook_grid_bounds(geom, bp, sar, live.gpr_look_angle_deg, /*x_margin_m=*/2.0,
x_min, x_max, z_min, z_max);
live.sar3d_x_min_m = x_min;
live.sar3d_x_max_m = x_max;
live.sar3d_z_min_m = z_min;
live.sar3d_z_max_m = z_max;
return live;
}
constexpr double kTxSweepTimeS = 0.079;
constexpr int kTxCount = 2;
constexpr int kRxCount = 4;
auto load_collection(
const radar::config::RunConfig& run_config,
const std::string& root,
const std::string& snapshot,
std::uint64_t collection_id
@@ -121,46 +140,73 @@ auto load_collection(
radar::ipc::PreprocessedCollection collection{};
collection.collection_id = collection_id;
// Frame time from the `_ns<...>` folder suffix, identical to standalone_bp.
const auto ns_pos = snapshot.rfind("_ns");
if (ns_pos == std::string::npos) {
throw std::runtime_error("snapshot has no _ns suffix: " + snapshot);
}
const auto mid_ns = static_cast<std::uint64_t>(std::stoull(snapshot.substr(ns_pos + 3)));
// The processor derives tx_sweep_time from capture_span / n_tx and the frame
// time from the capture midpoint, so set the window symmetric around mid_ns.
const auto span_ns = static_cast<std::uint64_t>(kTxCount * kTxSweepTimeS * 1e9);
collection.capture_start_ns = mid_ns - span_ns / 2U;
collection.capture_end_ns = mid_ns + span_ns / 2U;
collection.monotonic_ns = mid_ns;
for (int rx = 0; rx < kRxCount; ++rx) {
for (int tx = 0; tx < kTxCount; ++tx) {
char path[512];
std::snprintf(path, sizeof(path), "%s/%s/i%d_o%d_s21.npy", root.c_str(), snapshot.c_str(), rx, tx);
const auto s21 = load_npy(path);
std::snprintf(path, sizeof(path), "%s/%s/i%d_o%d_freq.npy", root.c_str(), snapshot.c_str(), rx, tx);
const auto freq = load_npy(path);
radar::ipc::SweepTraceBlock trace{};
trace.combo = radar::ipc::ComboKey{
.input_pos = static_cast<std::uint32_t>(rx),
.output_pos = static_cast<std::uint32_t>(tx),
};
trace.frequency_hz.assign(freq.f64.begin(), freq.f64.end());
trace.s21.reserve(s21.c128.size());
for (const auto& sample : s21.c128) {
trace.s21.push_back(radar::ipc::Complex32{
.re = static_cast<float>(sample.real()),
.im = static_cast<float>(sample.imag()),
});
// Per-trace capture windows from the snapshot's meta.json; the segmented
// timing model is reconstructed from them exactly like in the notebook.
std::map<std::pair<std::uint32_t, std::uint32_t>, std::pair<std::uint64_t, std::uint64_t>> windows{};
{
std::ifstream in(root + "/preprocessed/" + snapshot + "/meta.json");
if (in) {
nlohmann::json meta{};
in >> meta;
for (const auto& trace : meta.value("traces", nlohmann::json::array())) {
if (!trace.contains("capture_start_ns") || !trace.contains("capture_end_ns")) {
continue;
}
windows[{trace.at("input").get<std::uint32_t>(), trace.at("output").get<std::uint32_t>()}] = {
trace.at("capture_start_ns").get<std::uint64_t>(),
trace.at("capture_end_ns").get<std::uint64_t>(),
};
}
trace.s11.resize(trace.s21.size());
trace.capture_start_ns = collection.capture_start_ns;
trace.capture_end_ns = collection.capture_end_ns;
collection.traces.push_back(std::move(trace));
}
}
std::uint64_t first = 0U;
std::uint64_t last = 0U;
for (const auto& [key, window] : windows) {
if (first == 0U || window.first < first) {
first = window.first;
}
last = std::max(last, window.second);
}
if (first == 0U) {
// No meta labels: fall back to the `_ns<...>` folder suffix.
const auto ns_pos = snapshot.rfind("_ns");
if (ns_pos == std::string::npos) {
throw std::runtime_error("snapshot has no meta.json and no _ns suffix: " + snapshot);
}
last = static_cast<std::uint64_t>(std::stoull(snapshot.substr(ns_pos + 3)));
first = last;
}
collection.capture_start_ns = first;
collection.capture_end_ns = last;
collection.monotonic_ns = (first + last) / 2U;
for (const auto& combo : run_config.run_combos) {
char path[512];
std::snprintf(path, sizeof(path), "%s/preprocessed/%s/i%u_o%u_s21.npy",
root.c_str(), snapshot.c_str(), combo.input_pos, combo.output_pos);
const auto s21 = load_npy(path);
std::snprintf(path, sizeof(path), "%s/preprocessed/%s/i%u_o%u_freq.npy",
root.c_str(), snapshot.c_str(), combo.input_pos, combo.output_pos);
const auto freq = load_npy(path);
radar::ipc::SweepTraceBlock trace{};
trace.combo = combo;
trace.frequency_hz.assign(freq.f64.begin(), freq.f64.end());
trace.s21.reserve(s21.c128.size());
for (const auto& sample : s21.c128) {
trace.s21.push_back(radar::ipc::Complex32{
.re = static_cast<float>(sample.real()),
.im = static_cast<float>(sample.imag()),
});
}
trace.s11.resize(trace.s21.size());
if (const auto found = windows.find({combo.input_pos, combo.output_pos}); found != windows.end()) {
trace.capture_start_ns = found->second.first;
trace.capture_end_ns = found->second.second;
}
collection.traces.push_back(std::move(trace));
}
return collection;
}
@@ -195,14 +241,21 @@ void print_projection_argmax(
int main(int argc, char** argv) {
if (argc < 3) {
std::printf("usage: %s <preprocessed_dir> <center_snapshot>\n", argv[0]);
std::printf("usage: %s <dataset_root> <center_snapshot> [--dump-dir <dir>]\n", argv[0]);
return 1;
}
const std::string root = argv[1];
const std::string center = argv[2];
std::string dump_dir{};
for (int index = 3; index + 1 < argc; ++index) {
if (std::string{argv[index]} == "--dump-dir") {
dump_dir = argv[index + 1];
}
}
const std::string preprocessed = root + "/preprocessed";
std::vector<std::string> snapshots{};
if (DIR* dir = opendir(root.c_str())) {
if (DIR* dir = opendir(preprocessed.c_str())) {
for (dirent* entry = nullptr; (entry = readdir(dir)) != nullptr;) {
const std::string name = entry->d_name;
if (!name.empty() && (std::isdigit(name[0]) != 0)) {
@@ -214,12 +267,19 @@ int main(int argc, char** argv) {
std::sort(snapshots.begin(), snapshots.end());
const auto center_it = std::find(snapshots.begin(), snapshots.end(), center);
if (center_it == snapshots.end()) {
std::printf("center %s not found\n", center.c_str());
std::printf("center %s not found in %s\n", center.c_str(), preprocessed.c_str());
return 1;
}
const auto center_index = static_cast<std::size_t>(center_it - snapshots.begin());
const auto live = build_live_config();
const auto run_config = radar::config::load_run_config(root + "/config_profile.json");
if (run_config.gpr.tx_geometry.empty() || run_config.gpr.rx_geometry.empty()) {
std::printf("FAIL: config_profile.json has no gpr geometry\n");
return 1;
}
auto live = build_live_config(run_config);
live.sar3d_volume_dump_dir = dump_dir;
const auto radius = static_cast<std::size_t>(live.sar3d_frame_radius);
if (center_index < radius || center_index + radius >= snapshots.size()) {
std::printf("not enough snapshots around center for radius %zu\n", radius);
@@ -228,14 +288,13 @@ int main(int argc, char** argv) {
// The processor centers the aperture `radius` frames before the newest
// collection, so feeding [center-r .. center+r] in order reproduces the
// symmetric aperture of standalone_bp around `center`.
// symmetric aperture around `center`.
std::vector<radar::ipc::PreprocessedCollection> history{};
for (std::size_t index = center_index - radius; index <= center_index + radius; ++index) {
history.push_back(load_collection(root, snapshots[index], index));
history.push_back(load_collection(run_config, root, snapshots[index], index));
}
std::printf("aperture: %zu frames around %s\n", history.size(), center.c_str());
const auto run_config = build_run_config();
radar::processing::Sar3dProcessor processor{};
radar::ipc::ResultCollection result{};
for (std::size_t index = 0U; index < history.size(); ++index) {
@@ -253,7 +312,7 @@ int main(int argc, char** argv) {
}
for (const auto& payload : result.collection_payloads) {
if (payload.kind == radar::ipc::ResultKind::ScalarF32) {
std::printf("%s: %.6e (compare with standalone_bp max |acc/cnt|)\n",
std::printf("%s (volume level): %.9e\n",
payload.processing_name.c_str(), static_cast<double>(payload.scalar_value));
}
}
@@ -269,6 +328,7 @@ int main(int argc, char** argv) {
std::span<const radar::ipc::PreprocessedCollection>(history.data(), history.size() - 1U);
for (const float angle_deg : {0.0F, 30.0F, 60.0F, 90.0F}) {
auto live_angle = live;
live_angle.sar3d_volume_dump_dir.clear();
live_angle.sar3d_bottom_view_deg = angle_deg;
const auto angle_result =
processor.process_collection(run_config, history.back(), last_previous, live_angle);
@@ -279,6 +339,11 @@ int main(int argc, char** argv) {
}
}
}
if (!dump_dir.empty()) {
std::printf("volume dumped to %s/sar3d_volume_%llu.{f32,json}\n",
dump_dir.c_str(),
static_cast<unsigned long long>(history.back().collection_id));
}
std::printf("OK: sar3d processor path completed\n");
return 0;
}
@@ -90,14 +90,22 @@ _WEB_LIVE_SCHEMA = [
("gpr_speed_m_s", "Motion", _GPR_LIKE_MODES, _dual("_gpr_speed_m_s", {"legacy_gpr": "_legacy_gpr_speed_m_s", "sar3d": "_sar3d_speed_m_s"})),
("sar3d_frame_radius", "3D SAR", ("sar3d",), _attr("_sar3d_frame_radius")),
("sar3d_edge_frame_weight", "3D SAR", ("sar3d",), _attr("_sar3d_edge_frame_weight")),
("sar3d_pitch_override_enabled", "3D SAR", ("sar3d",), _attr("_sar3d_pitch_override_enabled")),
("sar3d_pitch_deg", "3D SAR", ("sar3d",), _attr("_sar3d_pitch_deg")),
("sar3d_center_depth_m", "3D SAR", ("sar3d",), _attr("_sar3d_center_depth_m")),
("sar3d_half_width_x_m", "3D SAR", ("sar3d",), _attr("_sar3d_half_width_x_m")),
("sar3d_half_width_y_m", "3D SAR", ("sar3d",), _attr("_sar3d_half_width_y_m")),
("sar3d_half_width_z_m", "3D SAR", ("sar3d",), _attr("_sar3d_half_width_z_m")),
("sar3d_grid_nx", "3D SAR", ("sar3d",), _attr("_sar3d_grid_nx")),
("sar3d_grid_ny", "3D SAR", ("sar3d",), _attr("_sar3d_grid_ny")),
("sar3d_grid_nz", "3D SAR", ("sar3d",), _attr("_sar3d_grid_nz")),
("sar3d_ground_depth_m", "3D SAR", ("sar3d",), _attr("_sar3d_ground_depth_m")),
("sar3d_x_min_m", "3D SAR", ("sar3d",), _attr("_sar3d_x_min_m")),
("sar3d_x_max_m", "3D SAR", ("sar3d",), _attr("_sar3d_x_max_m")),
("sar3d_y_min_m", "3D SAR", ("sar3d",), _attr("_sar3d_y_min_m")),
("sar3d_y_max_m", "3D SAR", ("sar3d",), _attr("_sar3d_y_max_m")),
("sar3d_z_min_m", "3D SAR", ("sar3d",), _attr("_sar3d_z_min_m")),
("sar3d_z_max_m", "3D SAR", ("sar3d",), _attr("_sar3d_z_max_m")),
("sar3d_voxel_x_m", "3D SAR", ("sar3d",), _attr("_sar3d_voxel_x_m")),
("sar3d_voxel_y_m", "3D SAR", ("sar3d",), _attr("_sar3d_voxel_y_m")),
("sar3d_voxel_z_m", "3D SAR", ("sar3d",), _attr("_sar3d_voxel_z_m")),
("sar3d_norm_edge_margin_m", "3D SAR", ("sar3d",), _attr("_sar3d_norm_edge_margin_m")),
("sar3d_smooth_sigma", "3D SAR", ("sar3d",), _attr("_sar3d_smooth_sigma")),
("sar3d_comp_range_power", "3D SAR", ("sar3d",), _attr("_sar3d_comp_range_power")),
("sar3d_tx_per_segment", "3D SAR", ("sar3d",), _attr("_sar3d_tx_per_segment")),
("sar3d_use_compensation", "3D SAR", ("sar3d",), _attr("_sar3d_use_compensation")),
("sar3d_normalize_pairs", "3D SAR", ("sar3d",), _attr("_sar3d_normalize_pairs")),
("sar3d_bottom_view_deg", "3D SAR", ("sar3d",), _attr("_sar3d_bottom_view")),
@@ -553,9 +561,8 @@ class AppWindowLiveProcessingMixin:
f"depth={self._sar3d_min_depth_m.value():g}..{self._sar3d_max_depth_m.value():g} m, "
f"freq={self._sar3d_start_freq_mhz.value():g}..{self._sar3d_stop_freq_mhz.value():g} MHz, "
f"frame_radius={self._sar3d_frame_radius.value()}, "
f"grid={self._sar3d_grid_nx.value()}x{self._sar3d_grid_ny.value()}x{self._sar3d_grid_nz.value()}, "
f"center_depth={self._sar3d_center_depth_m.value():g} m, "
f"pitch={self._sar3d_pitch_deg.value():g} deg)"
f"voxel={self._sar3d_voxel_x_m.value():g}x{self._sar3d_voxel_y_m.value():g}x{self._sar3d_voxel_z_m.value():g} m, "
f"ground_depth={self._sar3d_ground_depth_m.value():g} m)"
)
def _clear_history_mode_caches(self) -> None:
@@ -350,14 +350,22 @@ class AppWindowConfigProfileIOMixin:
self._sar3d_apply_freq_phase_correction,
self._sar3d_frame_radius,
self._sar3d_edge_frame_weight,
self._sar3d_pitch_override_enabled,
self._sar3d_pitch_deg,
self._sar3d_center_depth_m,
self._sar3d_half_width_x_m,
self._sar3d_half_width_y_m,
self._sar3d_half_width_z_m,
self._sar3d_grid_nx,
self._sar3d_grid_ny,
self._sar3d_grid_nz,
self._sar3d_ground_depth_m,
self._sar3d_x_min_m,
self._sar3d_x_max_m,
self._sar3d_y_min_m,
self._sar3d_y_max_m,
self._sar3d_z_min_m,
self._sar3d_z_max_m,
self._sar3d_voxel_x_m,
self._sar3d_voxel_y_m,
self._sar3d_voxel_z_m,
self._sar3d_norm_edge_margin_m,
self._sar3d_smooth_sigma,
self._sar3d_comp_range_power,
self._sar3d_tx_per_segment,
self._sar3d_use_compensation,
self._sar3d_normalize_pairs,
self._save_count,
@@ -562,14 +570,22 @@ class AppWindowConfigProfileIOMixin:
self._sar3d_apply_freq_phase_correction.setChecked(bool(sar3d_state.apply_freq_phase_correction))
self._sar3d_frame_radius.setValue(int(sar3d_state.frame_radius))
self._sar3d_edge_frame_weight.setValue(float(sar3d_state.edge_frame_weight))
self._sar3d_pitch_override_enabled.setChecked(bool(sar3d_state.pitch_override_enabled))
self._sar3d_pitch_deg.setValue(float(sar3d_state.pitch_deg))
self._sar3d_center_depth_m.setValue(float(sar3d_state.center_depth_m))
self._sar3d_half_width_x_m.setValue(float(sar3d_state.half_width_x_m))
self._sar3d_half_width_y_m.setValue(float(sar3d_state.half_width_y_m))
self._sar3d_half_width_z_m.setValue(float(sar3d_state.half_width_z_m))
self._sar3d_grid_nx.setValue(int(sar3d_state.grid_nx))
self._sar3d_grid_ny.setValue(int(sar3d_state.grid_ny))
self._sar3d_grid_nz.setValue(int(sar3d_state.grid_nz))
self._sar3d_ground_depth_m.setValue(float(sar3d_state.ground_depth_m))
self._sar3d_x_min_m.setValue(float(sar3d_state.x_min_m))
self._sar3d_x_max_m.setValue(float(sar3d_state.x_max_m))
self._sar3d_y_min_m.setValue(float(sar3d_state.y_min_m))
self._sar3d_y_max_m.setValue(float(sar3d_state.y_max_m))
self._sar3d_z_min_m.setValue(float(sar3d_state.z_min_m))
self._sar3d_z_max_m.setValue(float(sar3d_state.z_max_m))
self._sar3d_voxel_x_m.setValue(float(sar3d_state.voxel_x_m))
self._sar3d_voxel_y_m.setValue(float(sar3d_state.voxel_y_m))
self._sar3d_voxel_z_m.setValue(float(sar3d_state.voxel_z_m))
self._sar3d_norm_edge_margin_m.setValue(float(sar3d_state.norm_edge_margin_m))
self._sar3d_smooth_sigma.setValue(float(sar3d_state.smooth_sigma))
self._sar3d_comp_range_power.setValue(float(sar3d_state.comp_range_power))
self._sar3d_tx_per_segment.setValue(int(sar3d_state.tx_per_segment))
self._sar3d_use_compensation.setChecked(bool(sar3d_state.use_compensation))
self._sar3d_normalize_pairs.setChecked(bool(sar3d_state.normalize_pairs))
self._set_combo_current_text(
@@ -438,14 +438,22 @@ class AppWindowConfigStateBuildersMixin:
apply_freq_phase_correction=bool(self._sar3d_apply_freq_phase_correction.isChecked()),
frame_radius=int(self._sar3d_frame_radius.value()),
edge_frame_weight=float(self._sar3d_edge_frame_weight.value()),
pitch_override_enabled=bool(self._sar3d_pitch_override_enabled.isChecked()),
pitch_deg=float(self._sar3d_pitch_deg.value()),
center_depth_m=float(self._sar3d_center_depth_m.value()),
half_width_x_m=float(self._sar3d_half_width_x_m.value()),
half_width_y_m=float(self._sar3d_half_width_y_m.value()),
half_width_z_m=float(self._sar3d_half_width_z_m.value()),
grid_nx=int(self._sar3d_grid_nx.value()),
grid_ny=int(self._sar3d_grid_ny.value()),
grid_nz=int(self._sar3d_grid_nz.value()),
ground_depth_m=float(self._sar3d_ground_depth_m.value()),
x_min_m=float(self._sar3d_x_min_m.value()),
x_max_m=float(self._sar3d_x_max_m.value()),
y_min_m=float(self._sar3d_y_min_m.value()),
y_max_m=float(self._sar3d_y_max_m.value()),
z_min_m=float(self._sar3d_z_min_m.value()),
z_max_m=float(self._sar3d_z_max_m.value()),
voxel_x_m=float(self._sar3d_voxel_x_m.value()),
voxel_y_m=float(self._sar3d_voxel_y_m.value()),
voxel_z_m=float(self._sar3d_voxel_z_m.value()),
norm_edge_margin_m=float(self._sar3d_norm_edge_margin_m.value()),
smooth_sigma=float(self._sar3d_smooth_sigma.value()),
comp_range_power=float(self._sar3d_comp_range_power.value()),
tx_per_segment=int(self._sar3d_tx_per_segment.value()),
use_compensation=bool(self._sar3d_use_compensation.isChecked()),
normalize_pairs=bool(self._sar3d_normalize_pairs.isChecked()),
bottom_view_deg=sar3d_bottom_view_deg_from_label(
@@ -660,49 +660,114 @@ def build_processing_group(owner) -> QGroupBox:
"Weight of the outermost aperture frames (1 = flat window, 0 = strong taper)."
)
owner._sar3d_pitch_override_enabled = QCheckBox("Override pitch")
owner._sar3d_pitch_override_enabled.setChecked(bool(sar3d_defaults.pitch_override_enabled))
owner._sar3d_pitch_override_enabled.setToolTip(
"When unchecked, the frame pitch follows the look angle (notebook "
"BP_GROUND_PITCH_DEG = None); check to set it separately below."
)
owner._sar3d_pitch_deg = QDoubleSpinBox()
owner._sar3d_pitch_deg.setDecimals(2)
owner._sar3d_pitch_deg.setRange(-90.0, 90.0)
owner._sar3d_pitch_deg.setSingleStep(0.5)
owner._sar3d_pitch_deg.setValue(float(sar3d_defaults.pitch_deg))
owner._sar3d_pitch_deg.setToolTip("Antenna pitch above the horizon (deg); tilts the imaged volume.")
owner._sar3d_pitch_deg.setToolTip(
"Frame pitch (deg) used for ground coordinates and aperture offsets when the override is on."
)
owner._sar3d_center_depth_m = QDoubleSpinBox()
owner._sar3d_center_depth_m.setDecimals(2)
owner._sar3d_center_depth_m.setRange(0.0, 50.0)
owner._sar3d_center_depth_m.setSingleStep(0.1)
owner._sar3d_center_depth_m.setValue(float(sar3d_defaults.center_depth_m))
owner._sar3d_center_depth_m.setToolTip("Range distance to the center of the imaged volume (m).")
owner._sar3d_ground_depth_m = QDoubleSpinBox()
owner._sar3d_ground_depth_m.setDecimals(2)
owner._sar3d_ground_depth_m.setRange(0.0, 10.0)
owner._sar3d_ground_depth_m.setSingleStep(0.05)
owner._sar3d_ground_depth_m.setValue(float(sar3d_defaults.ground_depth_m))
owner._sar3d_ground_depth_m.setToolTip(
"Height of the lower antenna row above the ground (m); scene Y=0 is the ground."
)
owner._sar3d_half_width_x_m = QDoubleSpinBox()
owner._sar3d_half_width_x_m.setDecimals(2)
owner._sar3d_half_width_x_m.setRange(0.01, 50.0)
owner._sar3d_half_width_x_m.setSingleStep(0.1)
owner._sar3d_half_width_x_m.setValue(float(sar3d_defaults.half_width_x_m))
owner._sar3d_x_min_m = QDoubleSpinBox()
owner._sar3d_x_min_m.setDecimals(2)
owner._sar3d_x_min_m.setRange(-50.0, 50.0)
owner._sar3d_x_min_m.setSingleStep(0.1)
owner._sar3d_x_min_m.setValue(float(sar3d_defaults.x_min_m))
owner._sar3d_half_width_y_m = QDoubleSpinBox()
owner._sar3d_half_width_y_m.setDecimals(2)
owner._sar3d_half_width_y_m.setRange(0.01, 50.0)
owner._sar3d_half_width_y_m.setSingleStep(0.1)
owner._sar3d_half_width_y_m.setValue(float(sar3d_defaults.half_width_y_m))
owner._sar3d_x_max_m = QDoubleSpinBox()
owner._sar3d_x_max_m.setDecimals(2)
owner._sar3d_x_max_m.setRange(-50.0, 50.0)
owner._sar3d_x_max_m.setSingleStep(0.1)
owner._sar3d_x_max_m.setValue(float(sar3d_defaults.x_max_m))
owner._sar3d_half_width_z_m = QDoubleSpinBox()
owner._sar3d_half_width_z_m.setDecimals(2)
owner._sar3d_half_width_z_m.setRange(0.01, 50.0)
owner._sar3d_half_width_z_m.setSingleStep(0.1)
owner._sar3d_half_width_z_m.setValue(float(sar3d_defaults.half_width_z_m))
owner._sar3d_z_min_m = QDoubleSpinBox()
owner._sar3d_z_min_m.setDecimals(2)
owner._sar3d_z_min_m.setRange(-50.0, 50.0)
owner._sar3d_z_min_m.setSingleStep(0.1)
owner._sar3d_z_min_m.setValue(float(sar3d_defaults.z_min_m))
owner._sar3d_grid_nx = QSpinBox()
owner._sar3d_grid_nx.setRange(2, 1024)
owner._sar3d_grid_nx.setValue(int(sar3d_defaults.grid_nx))
owner._sar3d_z_max_m = QDoubleSpinBox()
owner._sar3d_z_max_m.setDecimals(2)
owner._sar3d_z_max_m.setRange(-50.0, 50.0)
owner._sar3d_z_max_m.setSingleStep(0.1)
owner._sar3d_z_max_m.setValue(float(sar3d_defaults.z_max_m))
owner._sar3d_grid_ny = QSpinBox()
owner._sar3d_grid_ny.setRange(2, 1024)
owner._sar3d_grid_ny.setValue(int(sar3d_defaults.grid_ny))
owner._sar3d_y_min_m = QDoubleSpinBox()
owner._sar3d_y_min_m.setDecimals(2)
owner._sar3d_y_min_m.setRange(-10.0, 10.0)
owner._sar3d_y_min_m.setSingleStep(0.1)
owner._sar3d_y_min_m.setValue(float(sar3d_defaults.y_min_m))
owner._sar3d_grid_nz = QSpinBox()
owner._sar3d_grid_nz.setRange(2, 1024)
owner._sar3d_grid_nz.setValue(int(sar3d_defaults.grid_nz))
owner._sar3d_y_max_m = QDoubleSpinBox()
owner._sar3d_y_max_m.setDecimals(2)
owner._sar3d_y_max_m.setRange(-10.0, 10.0)
owner._sar3d_y_max_m.setSingleStep(0.1)
owner._sar3d_y_max_m.setValue(float(sar3d_defaults.y_max_m))
owner._sar3d_voxel_x_m = QDoubleSpinBox()
owner._sar3d_voxel_x_m.setDecimals(3)
owner._sar3d_voxel_x_m.setRange(0.005, 1.0)
owner._sar3d_voxel_x_m.setSingleStep(0.005)
owner._sar3d_voxel_x_m.setValue(float(sar3d_defaults.voxel_x_m))
owner._sar3d_voxel_y_m = QDoubleSpinBox()
owner._sar3d_voxel_y_m.setDecimals(3)
owner._sar3d_voxel_y_m.setRange(0.005, 1.0)
owner._sar3d_voxel_y_m.setSingleStep(0.005)
owner._sar3d_voxel_y_m.setValue(float(sar3d_defaults.voxel_y_m))
owner._sar3d_voxel_z_m = QDoubleSpinBox()
owner._sar3d_voxel_z_m.setDecimals(3)
owner._sar3d_voxel_z_m.setRange(0.005, 1.0)
owner._sar3d_voxel_z_m.setSingleStep(0.005)
owner._sar3d_voxel_z_m.setValue(float(sar3d_defaults.voxel_z_m))
owner._sar3d_norm_edge_margin_m = QDoubleSpinBox()
owner._sar3d_norm_edge_margin_m.setDecimals(2)
owner._sar3d_norm_edge_margin_m.setRange(0.0, 5.0)
owner._sar3d_norm_edge_margin_m.setSingleStep(0.05)
owner._sar3d_norm_edge_margin_m.setValue(float(sar3d_defaults.norm_edge_margin_m))
owner._sar3d_norm_edge_margin_m.setToolTip(
"Grid Y edges excluded from volume normalization (constant-range ridge)."
)
owner._sar3d_smooth_sigma = QDoubleSpinBox()
owner._sar3d_smooth_sigma.setDecimals(2)
owner._sar3d_smooth_sigma.setRange(0.0, 10.0)
owner._sar3d_smooth_sigma.setSingleStep(0.25)
owner._sar3d_smooth_sigma.setValue(float(sar3d_defaults.smooth_sigma))
owner._sar3d_smooth_sigma.setToolTip("Gaussian smoothing sigma over (Z, X) planes, voxels.")
owner._sar3d_comp_range_power = QDoubleSpinBox()
owner._sar3d_comp_range_power.setDecimals(2)
owner._sar3d_comp_range_power.setRange(0.0, 2.0)
owner._sar3d_comp_range_power.setSingleStep(0.05)
owner._sar3d_comp_range_power.setValue(float(sar3d_defaults.comp_range_power))
owner._sar3d_comp_range_power.setToolTip("Geometric spreading compensation power.")
owner._sar3d_tx_per_segment = QSpinBox()
owner._sar3d_tx_per_segment.setRange(1, 64)
owner._sar3d_tx_per_segment.setValue(int(sar3d_defaults.tx_per_segment))
owner._sar3d_tx_per_segment.setToolTip(
"Fallback Tx-per-segment split of a frame when per-trace capture stamps are absent."
)
owner._sar3d_use_compensation = QCheckBox("Use per-voxel compensation")
owner._sar3d_use_compensation.setChecked(bool(sar3d_defaults.use_compensation))
@@ -738,19 +803,27 @@ def build_processing_group(owner) -> QGroupBox:
owner._sar3d_apply_freq_phase_correction,
("Frame radius", owner._sar3d_frame_radius),
("Edge frame weight", owner._sar3d_edge_frame_weight),
owner._sar3d_pitch_override_enabled,
("Pitch deg", owner._sar3d_pitch_deg),
("Center depth m", owner._sar3d_center_depth_m),
("Half width X m", owner._sar3d_half_width_x_m),
("Half width Y m", owner._sar3d_half_width_y_m),
("Half width Z m", owner._sar3d_half_width_z_m),
("Grid Nx", owner._sar3d_grid_nx),
("Grid Ny", owner._sar3d_grid_ny),
("Grid Nz", owner._sar3d_grid_nz),
("Ground depth m", owner._sar3d_ground_depth_m),
("X min m", owner._sar3d_x_min_m),
("X max m", owner._sar3d_x_max_m),
("Y min m", owner._sar3d_y_min_m),
("Y max m", owner._sar3d_y_max_m),
("Z min m", owner._sar3d_z_min_m),
("Z max m", owner._sar3d_z_max_m),
("Voxel X m", owner._sar3d_voxel_x_m),
("Voxel Y m", owner._sar3d_voxel_y_m),
("Voxel Z m", owner._sar3d_voxel_z_m),
("Norm edge margin m", owner._sar3d_norm_edge_margin_m),
("Smooth sigma", owner._sar3d_smooth_sigma),
("Comp range power", owner._sar3d_comp_range_power),
("Tx per segment", owner._sar3d_tx_per_segment),
("Bottom view", owner._sar3d_bottom_view),
owner._sar3d_use_compensation,
owner._sar3d_normalize_pairs,
],
split_index=11,
split_index=13,
)
owner._processing_mode_pages.addWidget(sar3d_page)
@@ -839,14 +912,22 @@ def build_processing_group(owner) -> QGroupBox:
owner._sar3d_apply_freq_phase_correction.toggled.connect(owner._on_processing_live_settings_changed)
owner._sar3d_frame_radius.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_edge_frame_weight.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_pitch_override_enabled.toggled.connect(owner._on_processing_live_settings_changed)
owner._sar3d_pitch_deg.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_center_depth_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_half_width_x_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_half_width_y_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_half_width_z_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_grid_nx.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_grid_ny.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_grid_nz.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_ground_depth_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_x_min_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_x_max_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_z_min_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_z_max_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_y_min_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_y_max_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_voxel_x_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_voxel_y_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_voxel_z_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_norm_edge_margin_m.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_smooth_sigma.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_comp_range_power.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_tx_per_segment.valueChanged.connect(owner._on_processing_live_settings_changed)
owner._sar3d_use_compensation.toggled.connect(owner._on_processing_live_settings_changed)
owner._sar3d_normalize_pairs.toggled.connect(owner._on_processing_live_settings_changed)
owner._sar3d_bottom_view.currentTextChanged.connect(owner._on_processing_live_settings_changed)
+105 -24
View File
@@ -533,39 +533,102 @@ def gui_profile_from_dict(payload: dict[str, Any]) -> GuiProfileModel:
gui.processing.sar3d.edge_frame_weight,
"gui.processing.sar3d",
),
pitch_override_enabled=_optional_bool(
sar3d_object,
"pitch_override_enabled",
gui.processing.sar3d.pitch_override_enabled,
"gui.processing.sar3d",
),
pitch_deg=_optional_float(
sar3d_object,
"pitch_deg",
gui.processing.sar3d.pitch_deg,
"gui.processing.sar3d",
),
center_depth_m=_optional_float(
ground_depth_m=_optional_float(
sar3d_object,
"center_depth_m",
gui.processing.sar3d.center_depth_m,
"ground_depth_m",
gui.processing.sar3d.ground_depth_m,
"gui.processing.sar3d",
),
half_width_x_m=_optional_float(
x_min_m=_optional_float(
sar3d_object,
"half_width_x_m",
gui.processing.sar3d.half_width_x_m,
"x_min_m",
gui.processing.sar3d.x_min_m,
"gui.processing.sar3d",
),
half_width_y_m=_optional_float(
x_max_m=_optional_float(
sar3d_object,
"half_width_y_m",
gui.processing.sar3d.half_width_y_m,
"x_max_m",
gui.processing.sar3d.x_max_m,
"gui.processing.sar3d",
),
half_width_z_m=_optional_float(
y_min_m=_optional_float(
sar3d_object,
"half_width_z_m",
gui.processing.sar3d.half_width_z_m,
"y_min_m",
gui.processing.sar3d.y_min_m,
"gui.processing.sar3d",
),
y_max_m=_optional_float(
sar3d_object,
"y_max_m",
gui.processing.sar3d.y_max_m,
"gui.processing.sar3d",
),
z_min_m=_optional_float(
sar3d_object,
"z_min_m",
gui.processing.sar3d.z_min_m,
"gui.processing.sar3d",
),
z_max_m=_optional_float(
sar3d_object,
"z_max_m",
gui.processing.sar3d.z_max_m,
"gui.processing.sar3d",
),
voxel_x_m=_optional_float(
sar3d_object,
"voxel_x_m",
gui.processing.sar3d.voxel_x_m,
"gui.processing.sar3d",
),
voxel_y_m=_optional_float(
sar3d_object,
"voxel_y_m",
gui.processing.sar3d.voxel_y_m,
"gui.processing.sar3d",
),
voxel_z_m=_optional_float(
sar3d_object,
"voxel_z_m",
gui.processing.sar3d.voxel_z_m,
"gui.processing.sar3d",
),
norm_edge_margin_m=_optional_float(
sar3d_object,
"norm_edge_margin_m",
gui.processing.sar3d.norm_edge_margin_m,
"gui.processing.sar3d",
),
smooth_sigma=_optional_float(
sar3d_object,
"smooth_sigma",
gui.processing.sar3d.smooth_sigma,
"gui.processing.sar3d",
),
comp_range_power=_optional_float(
sar3d_object,
"comp_range_power",
gui.processing.sar3d.comp_range_power,
"gui.processing.sar3d",
),
tx_per_segment=_optional_int(
sar3d_object,
"tx_per_segment",
gui.processing.sar3d.tx_per_segment,
"gui.processing.sar3d",
),
grid_nx=_optional_int(sar3d_object, "grid_nx", gui.processing.sar3d.grid_nx, "gui.processing.sar3d"),
grid_ny=_optional_int(sar3d_object, "grid_ny", gui.processing.sar3d.grid_ny, "gui.processing.sar3d"),
grid_nz=_optional_int(sar3d_object, "grid_nz", gui.processing.sar3d.grid_nz, "gui.processing.sar3d"),
use_compensation=_optional_bool(
sar3d_object,
"use_compensation",
@@ -628,8 +691,18 @@ def gui_profile_from_dict(payload: dict[str, Any]) -> GuiProfileModel:
raise ValueError("gui.processing.sar3d.frame_radius must be >= 0")
if not 0.0 <= gui.processing.sar3d.edge_frame_weight <= 1.0:
raise ValueError("gui.processing.sar3d.edge_frame_weight must be within [0, 1]")
if min(gui.processing.sar3d.grid_nx, gui.processing.sar3d.grid_ny, gui.processing.sar3d.grid_nz) < 2:
raise ValueError("gui.processing.sar3d.grid_nx/ny/nz must each be >= 2")
if min(gui.processing.sar3d.voxel_x_m, gui.processing.sar3d.voxel_y_m, gui.processing.sar3d.voxel_z_m) <= 0.0:
raise ValueError("gui.processing.sar3d.voxel_x/y/z_m must each be > 0")
if gui.processing.sar3d.x_max_m <= gui.processing.sar3d.x_min_m:
raise ValueError("gui.processing.sar3d.x_max_m must be > x_min_m")
if gui.processing.sar3d.y_max_m <= gui.processing.sar3d.y_min_m:
raise ValueError("gui.processing.sar3d.y_max_m must be > y_min_m")
if gui.processing.sar3d.z_max_m <= gui.processing.sar3d.z_min_m:
raise ValueError("gui.processing.sar3d.z_max_m must be > z_min_m")
if gui.processing.sar3d.ground_depth_m < 0.0:
raise ValueError("gui.processing.sar3d.ground_depth_m must be >= 0")
if gui.processing.sar3d.tx_per_segment < 1:
raise ValueError("gui.processing.sar3d.tx_per_segment must be >= 1")
data_actions_object = _as_dict(gui_object.get("data_actions"), "gui.data_actions")
gui.data_actions = GuiDataActionsStateModel(
@@ -795,14 +868,22 @@ def gui_profile_to_dict(model: GuiProfileModel) -> dict[str, Any]:
"apply_freq_phase_correction": gui.processing.sar3d.apply_freq_phase_correction,
"frame_radius": gui.processing.sar3d.frame_radius,
"edge_frame_weight": gui.processing.sar3d.edge_frame_weight,
"pitch_override_enabled": gui.processing.sar3d.pitch_override_enabled,
"pitch_deg": gui.processing.sar3d.pitch_deg,
"center_depth_m": gui.processing.sar3d.center_depth_m,
"half_width_x_m": gui.processing.sar3d.half_width_x_m,
"half_width_y_m": gui.processing.sar3d.half_width_y_m,
"half_width_z_m": gui.processing.sar3d.half_width_z_m,
"grid_nx": gui.processing.sar3d.grid_nx,
"grid_ny": gui.processing.sar3d.grid_ny,
"grid_nz": gui.processing.sar3d.grid_nz,
"ground_depth_m": gui.processing.sar3d.ground_depth_m,
"x_min_m": gui.processing.sar3d.x_min_m,
"x_max_m": gui.processing.sar3d.x_max_m,
"y_min_m": gui.processing.sar3d.y_min_m,
"y_max_m": gui.processing.sar3d.y_max_m,
"z_min_m": gui.processing.sar3d.z_min_m,
"z_max_m": gui.processing.sar3d.z_max_m,
"voxel_x_m": gui.processing.sar3d.voxel_x_m,
"voxel_y_m": gui.processing.sar3d.voxel_y_m,
"voxel_z_m": gui.processing.sar3d.voxel_z_m,
"norm_edge_margin_m": gui.processing.sar3d.norm_edge_margin_m,
"smooth_sigma": gui.processing.sar3d.smooth_sigma,
"comp_range_power": gui.processing.sar3d.comp_range_power,
"tx_per_segment": gui.processing.sar3d.tx_per_segment,
"use_compensation": gui.processing.sar3d.use_compensation,
"normalize_pairs": gui.processing.sar3d.normalize_pairs,
"bottom_view_deg": gui.processing.sar3d.bottom_view_deg,
+24 -16
View File
@@ -128,31 +128,39 @@ class GuiSar3dStateModel:
"""UI-only defaults for the 3D SAR (sar3d) live settings.
The band, depth gate, antenna selection and motion knobs mirror the shared
gpr_* live-config fields; the sar3d_* knobs mirror the VNA_AlGo gpr_core
defaults used by the C++ sar3d processor.
gpr_* live-config fields; the sar3d_* knobs mirror the reference notebook
(Multi_antenna_3D_clean_v1.ipynb) defaults used by the C++ sar3d processor.
"""
input_positions: str = ""
output_positions: str = ""
min_depth_m: float = 2.0
max_depth_m: float = 14.0
start_freq_mhz: float = 3000.0
min_depth_m: float = 2.5
max_depth_m: float = 8.0
start_freq_mhz: float = 2000.0
stop_freq_mhz: float = 6000.0
speed_m_s: float = 0.0
ignore_socket_speed_enabled: bool = False
look_angle_deg: float = 0.0
look_angle_deg: float = 17.5
direction_sign: float = 1.0
apply_freq_phase_correction: bool = True
frame_radius: int = 6
edge_frame_weight: float = 0.75
pitch_deg: float = 6.0
center_depth_m: float = 5.6
half_width_x_m: float = 1.8
half_width_y_m: float = 1.0
half_width_z_m: float = 1.8
grid_nx: int = 150
grid_ny: int = 120
grid_nz: int = 160
frame_radius: int = 0
edge_frame_weight: float = 0.8
pitch_override_enabled: bool = False
pitch_deg: float = 17.5
ground_depth_m: float = 1.25
x_min_m: float = -2.0
x_max_m: float = 2.0
y_min_m: float = -0.5
y_max_m: float = 2.0
z_min_m: float = 0.0
z_max_m: float = 8.0
voxel_x_m: float = 0.025
voxel_y_m: float = 0.02
voxel_z_m: float = 0.02
norm_edge_margin_m: float = 0.3
smooth_sigma: float = 1.5
comp_range_power: float = 0.3
tx_per_segment: int = 2
use_compensation: bool = True
normalize_pairs: bool = True
# Bottom-view projection angle around the vertical Y axis, degrees:
@@ -72,19 +72,35 @@ class ProcessingLiveConfig:
gpr_visible_x_max_m: float = 2.0
gpr_visible_z_min_m: float = 0.0
gpr_visible_z_max_m: float = 14.0
# 3D SAR (sar3d processor) parameters; defaults mirror VNA_AlGo gpr_core.
# The frequency band, depth gate, antenna selection and motion parameters
# are shared with the gpr_* fields above.
sar3d_frame_radius: int = 6
sar3d_edge_frame_weight: float = 0.75
sar3d_pitch_deg: float = 6.0
sar3d_center_depth_m: float = 5.6
sar3d_half_width_x_m: float = 1.8
sar3d_half_width_y_m: float = 1.0
sar3d_half_width_z_m: float = 1.8
sar3d_grid_nx: int = 150
sar3d_grid_ny: int = 120
sar3d_grid_nz: int = 160
# 3D SAR (sar3d processor) parameters; defaults mirror the reference
# notebook Multi_antenna_3D_clean_v1.ipynb. The frequency band, depth gate,
# antenna selection and motion parameters are shared with the gpr_* fields.
sar3d_frame_radius: int = 0
sar3d_edge_frame_weight: float = 0.8
# Frame pitch override; disabled -> the pitch follows gpr_look_angle_deg.
sar3d_pitch_override_enabled: bool = False
sar3d_pitch_deg: float = 17.5
# Height of the lower antenna row above the ground (scene Y=0 = ground).
sar3d_ground_depth_m: float = 1.25
# Scene grid bounds (metres): X across, Y height above the ground, Z along
# the ground away from the radar. Voxels outside the radar depth gate are
# zeroed by the mask regardless of the Z window.
sar3d_x_min_m: float = -2.0
sar3d_x_max_m: float = 2.0
sar3d_y_min_m: float = -0.5
sar3d_y_max_m: float = 2.0
sar3d_z_min_m: float = 0.0
sar3d_z_max_m: float = 8.0
sar3d_voxel_x_m: float = 0.025
sar3d_voxel_y_m: float = 0.02
sar3d_voxel_z_m: float = 0.02
# Volume finalization: Y-edge margin excluded from normalization, Gaussian
# smoothing sigma over (Z, X) planes.
sar3d_norm_edge_margin_m: float = 0.3
sar3d_smooth_sigma: float = 1.5
sar3d_comp_range_power: float = 0.3
# Fallback frame segmentation when per-trace capture stamps are absent.
sar3d_tx_per_segment: int = 2
sar3d_use_compensation: bool = True
sar3d_normalize_pairs: bool = True
# Bottom-view projection angle around the vertical Y axis, degrees:
@@ -156,14 +172,22 @@ class ProcessingLiveConfig:
"gpr_visible_z_max_m": float(self.gpr_visible_z_max_m),
"sar3d_frame_radius": int(self.sar3d_frame_radius),
"sar3d_edge_frame_weight": float(self.sar3d_edge_frame_weight),
"sar3d_pitch_override_enabled": bool(self.sar3d_pitch_override_enabled),
"sar3d_pitch_deg": float(self.sar3d_pitch_deg),
"sar3d_center_depth_m": float(self.sar3d_center_depth_m),
"sar3d_half_width_x_m": float(self.sar3d_half_width_x_m),
"sar3d_half_width_y_m": float(self.sar3d_half_width_y_m),
"sar3d_half_width_z_m": float(self.sar3d_half_width_z_m),
"sar3d_grid_nx": int(self.sar3d_grid_nx),
"sar3d_grid_ny": int(self.sar3d_grid_ny),
"sar3d_grid_nz": int(self.sar3d_grid_nz),
"sar3d_ground_depth_m": float(self.sar3d_ground_depth_m),
"sar3d_x_min_m": float(self.sar3d_x_min_m),
"sar3d_x_max_m": float(self.sar3d_x_max_m),
"sar3d_y_min_m": float(self.sar3d_y_min_m),
"sar3d_y_max_m": float(self.sar3d_y_max_m),
"sar3d_z_min_m": float(self.sar3d_z_min_m),
"sar3d_z_max_m": float(self.sar3d_z_max_m),
"sar3d_voxel_x_m": float(self.sar3d_voxel_x_m),
"sar3d_voxel_y_m": float(self.sar3d_voxel_y_m),
"sar3d_voxel_z_m": float(self.sar3d_voxel_z_m),
"sar3d_norm_edge_margin_m": float(self.sar3d_norm_edge_margin_m),
"sar3d_smooth_sigma": float(self.sar3d_smooth_sigma),
"sar3d_comp_range_power": float(self.sar3d_comp_range_power),
"sar3d_tx_per_segment": int(self.sar3d_tx_per_segment),
"sar3d_use_compensation": bool(self.sar3d_use_compensation),
"sar3d_normalize_pairs": bool(self.sar3d_normalize_pairs),
"sar3d_bottom_view_deg": float(self.sar3d_bottom_view_deg),
+35 -5
View File
@@ -20,6 +20,7 @@ never see time jump backwards. ``collection_id`` increases forever.
from __future__ import annotations
import argparse
import json
import logging
from pathlib import Path
import re
@@ -75,6 +76,22 @@ def _load_dataset(dataset_dir: Path, combos: list[ComboKey], expected_points: in
timestamps_ns: list[int] = []
for name, timestamp_ns in named:
folder = dataset_dir / name
# Recorded per-trace capture windows from the snapshot's meta.json: the
# sar3d segmented timing model reconstructs the acquisition segments
# from them. Missing/partial meta leaves the zero "unmeasured" sentinel
# and the consumer falls back to its uniform-segment model.
windows: dict[tuple[int, int], tuple[int, int]] = {}
meta_path = folder / "meta.json"
if meta_path.is_file():
try:
meta = json.loads(meta_path.read_text())
for entry in meta.get("traces", []):
windows[(int(entry["input"]), int(entry["output"]))] = (
int(entry["capture_start_ns"]),
int(entry["capture_end_ns"]),
)
except (OSError, ValueError, KeyError, TypeError):
windows = {}
traces: list[TraceData] = []
for combo in combos:
s21_path = folder / f"i{combo.input}_o{combo.output}_s21.npy"
@@ -94,12 +111,15 @@ def _load_dataset(dataset_dir: Path, combos: list[ComboKey], expected_points: in
f"{s21_path}: {s21.size} points, but radar.sweep.points = "
f"{expected_points}. Align the profile's sweep with the dataset."
)
window = windows.get((combo.input, combo.output), (0, 0))
traces.append(
TraceData(
combo=combo,
frequency_hz=freq,
s11=np.zeros_like(s21),
s21=s21,
capture_start_ns=window[0],
capture_end_ns=window[1],
)
)
snapshots.append(traces)
@@ -114,28 +134,38 @@ def _build_collection(
collection_id: int,
frame_ns: int,
) -> SweepCollection:
"""Wrap one preloaded snapshot in a SweepCollection at the given frame time."""
"""Wrap one preloaded snapshot in a SweepCollection at the given frame time.
Recorded per-trace capture windows are re-anchored by the same shift as the
frame time, preserving the real intra-frame segment layout. Traces without
recorded windows get a synthetic window spanning the median frame gap.
"""
span_ns = dataset.median_gap_ns
start_ns = frame_ns - span_ns // 2
end_ns = frame_ns + span_ns // 2
shift_ns = frame_ns - dataset.timestamps_ns[index]
traces = []
for source in dataset.snapshots[index]:
if source.capture_end_ns > source.capture_start_ns:
capture = (source.capture_start_ns + shift_ns, source.capture_end_ns + shift_ns)
else:
capture = (start_ns, end_ns)
traces.append(
TraceData(
combo=source.combo,
frequency_hz=source.frequency_hz,
s11=source.s11,
s21=source.s21,
capture_start_ns=start_ns,
capture_end_ns=end_ns,
capture_start_ns=capture[0],
capture_end_ns=capture[1],
)
)
return SweepCollection(
collection_id=collection_id,
monotonic_ns=frame_ns,
traces=traces,
capture_start_ns=start_ns,
capture_end_ns=end_ns,
capture_start_ns=min(t.capture_start_ns for t in traces),
capture_end_ns=max(t.capture_end_ns for t in traces),
)
+21 -7
View File
@@ -84,14 +84,28 @@ def _sar3d_live_config(profile_path: Path) -> ProcessingLiveConfig:
live.sar3d_edge_frame_weight = float(
sar3d.get("edge_frame_weight", live.sar3d_edge_frame_weight)
)
live.sar3d_pitch_override_enabled = bool(
sar3d.get("pitch_override_enabled", live.sar3d_pitch_override_enabled)
)
live.sar3d_pitch_deg = float(sar3d.get("pitch_deg", live.sar3d_pitch_deg))
live.sar3d_center_depth_m = float(sar3d.get("center_depth_m", live.sar3d_center_depth_m))
live.sar3d_half_width_x_m = float(sar3d.get("half_width_x_m", live.sar3d_half_width_x_m))
live.sar3d_half_width_y_m = float(sar3d.get("half_width_y_m", live.sar3d_half_width_y_m))
live.sar3d_half_width_z_m = float(sar3d.get("half_width_z_m", live.sar3d_half_width_z_m))
live.sar3d_grid_nx = int(sar3d.get("grid_nx", live.sar3d_grid_nx))
live.sar3d_grid_ny = int(sar3d.get("grid_ny", live.sar3d_grid_ny))
live.sar3d_grid_nz = int(sar3d.get("grid_nz", live.sar3d_grid_nz))
live.sar3d_ground_depth_m = float(sar3d.get("ground_depth_m", live.sar3d_ground_depth_m))
live.sar3d_x_min_m = float(sar3d.get("x_min_m", live.sar3d_x_min_m))
live.sar3d_x_max_m = float(sar3d.get("x_max_m", live.sar3d_x_max_m))
live.sar3d_y_min_m = float(sar3d.get("y_min_m", live.sar3d_y_min_m))
live.sar3d_y_max_m = float(sar3d.get("y_max_m", live.sar3d_y_max_m))
live.sar3d_z_min_m = float(sar3d.get("z_min_m", live.sar3d_z_min_m))
live.sar3d_z_max_m = float(sar3d.get("z_max_m", live.sar3d_z_max_m))
live.sar3d_voxel_x_m = float(sar3d.get("voxel_x_m", live.sar3d_voxel_x_m))
live.sar3d_voxel_y_m = float(sar3d.get("voxel_y_m", live.sar3d_voxel_y_m))
live.sar3d_voxel_z_m = float(sar3d.get("voxel_z_m", live.sar3d_voxel_z_m))
live.sar3d_norm_edge_margin_m = float(
sar3d.get("norm_edge_margin_m", live.sar3d_norm_edge_margin_m)
)
live.sar3d_smooth_sigma = float(sar3d.get("smooth_sigma", live.sar3d_smooth_sigma))
live.sar3d_comp_range_power = float(
sar3d.get("comp_range_power", live.sar3d_comp_range_power)
)
live.sar3d_tx_per_segment = int(sar3d.get("tx_per_segment", live.sar3d_tx_per_segment))
live.sar3d_use_compensation = bool(sar3d.get("use_compensation", live.sar3d_use_compensation))
live.sar3d_normalize_pairs = bool(sar3d.get("normalize_pairs", live.sar3d_normalize_pairs))
live.sar3d_bottom_view_deg = float(sar3d.get("bottom_view_deg", live.sar3d_bottom_view_deg))
+389
View File
@@ -0,0 +1,389 @@
"""Reference implementation of the sar3d volume from Multi_antenna_3D_clean_v1.ipynb.
Executable specification for the C++ port: recomputes a single-frame
(aperture radius 0) 3D SAR volume from a preprocessed snapshot exactly the
way the notebook does, then compares it against a reference *_volume.npz
exported by the notebook.
Usage:
python sar3d_reference_check.py <dataset_root> <reference_volume.npz>
<dataset_root> is an extracted replay dataset (contains config_profile.json
and preprocessed/<snapshot>/i*_o*_{s21,freq}.npy + meta.json). The snapshot
name and all physics parameters are taken from the npz meta.
Exits non-zero when the recomputed volume differs from the reference by more
than the tolerance. The only expected difference is the Gaussian smoothing
(manual separable filter vs scipy) at the 1e-6 level.
"""
from __future__ import annotations
import argparse
import json
import sys
from pathlib import Path
import numpy as np
TOLERANCE = 2e-6
C_LIGHT = 3e8
# Notebook constants not recorded in the npz meta (cell 0 / cell 4).
BP_OVERSAMPLE = 8
BP_WINDOW = True
PAIR_NORM_PERCENTILE = 50.0
PAIR_NORM_EPS = 1e-15
COMP_RANGE_POWER = 0.3
COMP_RANGE_WEIGHT_MAX = 5.0
COMP_REF_DEPTH = 5.0
BP_SMOOTH_SIGMA = 1.5
SAR_Y_RANGE_M = (-0.5, 2.0)
SAR_NORM_EDGE_MARGIN_M = 0.3
X_MARGIN_M = 2.0
MOTION_PHASE_SIGN = -1.0
DIRECTION_SIGN = +1.0
def gaussian_kernel1d(sigma: float, radius: int) -> np.ndarray:
x = np.arange(-radius, radius + 1, dtype=float)
phi = np.exp(-0.5 / (sigma * sigma) * x * x)
return phi / phi.sum()
def correlate1d_reflect(arr: np.ndarray, weights: np.ndarray, axis: int) -> np.ndarray:
r = (len(weights) - 1) // 2
a = np.moveaxis(arr, axis, -1)
pad = np.concatenate([a[..., r - 1::-1], a, a[..., :-r - 1:-1]], axis=-1)
out = np.zeros_like(a)
for k, w in enumerate(weights):
out += w * pad[..., k:k + a.shape[-1]]
return np.moveaxis(out, -1, axis)
def gaussian_filter(arr: np.ndarray, sigma: tuple[float, ...], truncate: float = 4.0) -> np.ndarray:
"""scipy.ndimage.gaussian_filter equivalent (mode='reflect')."""
out = np.asarray(arr, dtype=float).copy()
for axis, s in enumerate(sigma):
if s <= 0:
continue
radius = int(truncate * float(s) + 0.5)
out = correlate1d_reflect(out, gaussian_kernel1d(float(s), radius), axis)
return out
def load_geometry(profile_path: Path) -> dict:
profile = json.loads(profile_path.read_text())
gpr = profile["gpr"]
tx_geom = {int(e["output_pos"]): (e["x_m"], e.get("y_m", 0.0), e.get("z_m", 0.0))
for e in gpr["tx_geometry"]}
rx_geom = {int(e["input_pos"]): (e["x_m"], e.get("y_m", 0.0), e.get("z_m", 0.0))
for e in gpr["rx_geometry"]}
combos = profile["run"]["combos"]
pair_order = [(int(c["output"]), int(c["input"])) for c in combos]
tx_order_full: list[int] = []
for tx, _ in pair_order:
if tx not in tx_order_full:
tx_order_full.append(tx)
output_idx = sorted(tx_geom)
input_idx = sorted(rx_geom)
tx_xyz = np.array([tx_geom[i] for i in output_idx], dtype=float)
rx_xyz = np.array([rx_geom[i] for i in input_idx], dtype=float)
return {
"output_idx": output_idx, "input_idx": input_idx,
"tx_xyz": tx_xyz, "rx_xyz": rx_xyz,
"pair_order": pair_order, "tx_order_full": tx_order_full,
}
def snapshot_timing(snapshot_dir: Path) -> dict | None:
"""Per-Tx capture windows from the snapshot's meta.json trace labels."""
try:
meta = json.loads((snapshot_dir / "meta.json").read_text())
except OSError:
return None
traces = [t for t in meta.get("traces", [])
if "capture_start_ns" in t and "capture_end_ns" in t and "output" in t]
if not traces:
return None
first = min(int(t["capture_start_ns"]) for t in traces)
last = max(int(t["capture_end_ns"]) for t in traces)
by_tx: dict[int, tuple[int, int]] = {}
for t in traces:
by_tx.setdefault(int(t["output"]), (int(t["capture_start_ns"]), int(t["capture_end_ns"])))
return {
"acq_s": (last - first) * 1e-9,
"center_abs_s": 0.5 * (first + last) * 1e-9,
"tx_window_s": {tx: ((a - first) * 1e-9, (b - first) * 1e-9) for tx, (a, b) in by_tx.items()},
}
def pair_timestamps(geometry: dict, timing: dict, speed_m_s: float, look_angle_deg: float,
wave_speed: float) -> dict:
"""Segment/slot assignment and coarse pair dtau (notebook compute_pair_timestamps_for_snapshot)."""
windows = timing["tx_window_s"]
hw_order = [tx for tx in geometry["tx_order_full"] if tx in windows]
hw_order += [tx for tx in sorted(windows) if tx not in hw_order]
groups: dict[tuple[float, float], list[int]] = {}
for tx in hw_order:
groups.setdefault(windows[tx], []).append(tx)
seg_of = {}
for seg_idx, win in enumerate(sorted(groups)):
for slot, tx in enumerate(groups[win]):
seg_of[tx] = (seg_idx, slot, len(groups[win]))
cos_theta = np.cos(np.radians(look_angle_deg))
t_ref = 0.5 * timing["acq_s"]
tx_map = {p: k for k, p in enumerate(geometry["output_idx"])}
rx_map = {p: k for k, p in enumerate(geometry["input_idx"])}
ts = {}
for tx_phys, rx_phys in geometry["pair_order"]:
_, slot, n_slots = seg_of[tx_phys]
t0, t1 = windows[tx_phys]
t_center = 0.5 * (t0 + t1)
dz_motion = DIRECTION_SIGN * speed_m_s * (t_center - t_ref) * cos_theta
ts[(tx_map[tx_phys], rx_map[rx_phys])] = {
"slot": slot, "n_slots": n_slots,
"segment_start_s": t0, "segment_time_s": t1 - t0,
"t_center_s": t_center,
"dtau_motion_s": 2.0 * dz_motion / wave_speed,
}
return ts
def timing_abs_center(snapshot_dir: Path) -> float:
timing = snapshot_timing(snapshot_dir)
if timing is None:
raise RuntimeError(f"meta.json trace labels required in {snapshot_dir}")
return timing["center_abs_s"]
def load_mimo(snapshot_dir: Path, geometry: dict) -> tuple[dict, dict]:
s21_data, freq_data = {}, {}
for f in snapshot_dir.glob("i*_o*_s21.npy"):
parts = f.stem.split("_")
rx_phys = int(parts[0][1:])
tx_phys = int(parts[1][1:])
if tx_phys not in geometry["output_idx"] or rx_phys not in geometry["input_idx"]:
continue
key = (geometry["output_idx"].index(tx_phys), geometry["input_idx"].index(rx_phys))
s21_data[key] = np.load(f)
freq_data[key] = np.asarray(np.load(f.parent / f"i{rx_phys}_o{tx_phys}_freq.npy"), dtype=float)
return s21_data, freq_data
def apply_phase_correction(s21: np.ndarray, freq_full: np.ndarray, info: dict,
f_start: float, f_stop: float, speed_m_s: float,
look_angle_deg: float, wave_speed: float) -> np.ndarray:
"""Coarse pair + intra-sweep motion correction as one phase ramp before the IFFT."""
idx = np.flatnonzero((freq_full >= f_start) & (freq_full <= f_stop))
n_full = len(freq_full)
dt_event = info["segment_time_s"] / (info["n_slots"] * n_full)
t_abs = info["segment_start_s"] + (info["n_slots"] * idx + info["slot"]) * dt_event
dt_intra = t_abs - info["t_center_s"]
dtau_intra = 2.0 * DIRECTION_SIGN * speed_m_s * np.cos(np.radians(look_angle_deg)) * dt_intra / wave_speed
total_dtau = dtau_intra + info["dtau_motion_s"]
phi = 2.0 * np.pi * freq_full[idx] * total_dtau
out = np.array(s21, dtype=np.complex128, copy=True)
out[idx] *= np.exp(1j * MOTION_PHASE_SIGN * phi)
return out
def compute_ascan(s21: np.ndarray, freq: np.ndarray, f_start: float, f_stop: float):
"""Oversampled complex A-scan with the notebook's frequency-shifted IFFT."""
mask = (freq >= f_start) & (freq <= f_stop)
freq_cut = freq[mask]
s21_cut = s21[mask]
df = float(np.median(np.diff(freq_cut)))
n = len(freq_cut)
k0 = int(round(freq_cut[0] / df))
n_fft = (1 << int(np.ceil(np.log2(2 * (k0 + n - 1))))) * BP_OVERSAMPLE
dt = 1.0 / (n_fft * df)
s = s21_cut * np.hanning(n) if BP_WINDOW else s21_cut.copy()
spectrum = np.zeros(n_fft, dtype=np.complex128)
spectrum[k0:k0 + n] = s
return np.arange(n_fft, dtype=float) * dt, np.fft.ifft(spectrum)
def interp_complex(tau: np.ndarray, t_axis: np.ndarray, h_axis: np.ndarray) -> np.ndarray:
n = t_axis.size
dt = (float(t_axis[-1]) - float(t_axis[0])) / (n - 1)
pos = (tau - float(t_axis[0])) / dt
inside = (pos >= 0.0) & (pos <= n - 1)
k = np.clip(np.floor(pos), 0, n - 2).astype(np.intp)
frac = pos - k
out = h_axis[k] + frac * (h_axis[k + 1] - h_axis[k])
return np.where(inside, out, 0.0 + 0.0j)
def axis_linspace(lo: float, hi: float, step: float) -> np.ndarray:
return np.linspace(lo, hi, max(int(round((hi - lo) / step)) + 1, 2))
def compute_volume(dataset_root: Path, meta: dict) -> dict:
"""Single-frame volume via the notebook's ANIMATION path.
The reference npz files were exported from animation frames: the base grid
is shifted by the radar displacement between the frame and the scene-origin
snapshot (meta reference_snapshot), antennas get the same offset, and the
depth window follows the radar (depth_shift). Physically it all cancels for
an aperture of radius 0, but replicating the exact float arithmetic is what
makes the comparison bit-tight at the depth-mask knife edge.
"""
geometry = load_geometry(dataset_root / "config_profile.json")
snapshot_dir = dataset_root / "preprocessed" / meta["snapshot"]
eps_r = 1.0
wave_speed = C_LIGHT / np.sqrt(eps_r)
f_start = float(meta["f_start_hz"])
f_stop = float(meta["f_stop_hz"])
min_depth = float(meta["min_depth_m"])
max_depth = float(meta["max_depth_m"])
speed = float(meta["speed_m_s"])
look_angle = float(meta["look_angle_deg"])
ground_depth = float(meta["ground_depth_m"])
voxel = [float(a) for a in meta["voxel_m"]]
x_tx, y_tx, z_tx = geometry["tx_xyz"].T
x_rx, y_rx, z_rx = geometry["rx_xyz"].T
n_tx, n_rx = len(x_tx), len(x_rx)
timing = snapshot_timing(snapshot_dir)
if timing is None:
raise RuntimeError(f"meta.json trace labels required in {snapshot_dir}")
ts = pair_timestamps(geometry, timing, speed, look_angle, wave_speed)
s21_data, freq_data = load_mimo(snapshot_dir, geometry)
ascan_h, ascan_t = {}, {}
for key, s21 in s21_data.items():
corrected = apply_phase_correction(
s21, freq_data[key], ts[key], f_start, f_stop, speed, look_angle, wave_speed)
t_axis, h = compute_ascan(corrected, freq_data[key], f_start, f_stop)
z_axis = t_axis * wave_speed / 2
m = (z_axis >= min_depth) & (z_axis <= max_depth)
scale = float(np.percentile(np.abs(h)[m], PAIR_NORM_PERCENTILE))
if not np.isfinite(scale) or scale <= PAIR_NORM_EPS:
scale = 1.0
ascan_h[key] = h / (scale + PAIR_NORM_EPS)
ascan_t[key] = t_axis
pitch = np.radians(look_angle)
cos_p, sin_p = np.cos(pitch), np.sin(pitch)
def scene_to_radar_yz(y_scene, z_scene):
y0 = y_scene - ground_depth
return y0 * cos_p + z_scene * sin_p, -y0 * sin_p + z_scene * cos_p
x_ant = np.concatenate([x_tx, x_rx])
y_lo, y_hi = SAR_Y_RANGE_M
z_span = [(d - float(scene_to_radar_yz(y, 0.0)[1])) / cos_p
for y in (y_lo, y_hi) for d in (min_depth, max_depth)]
gx = axis_linspace(float(x_ant.min() - X_MARGIN_M), float(x_ant.max() + X_MARGIN_M), voxel[0])
gy = axis_linspace(y_lo, y_hi, voxel[1])
gz = axis_linspace(min(z_span), max(z_span), voxel[2])
# Displacement of this frame relative to the scene-origin snapshot
# (notebook sar_frame_offsets + anim_frame_displacement).
off_dy = off_dz = 0.0
disp_y = disp_z = 0.0
reference_name = meta.get("reference_snapshot", meta["snapshot"])
if reference_name != meta["snapshot"]:
ref_timing = snapshot_timing(dataset_root / "preprocessed" / reference_name)
if ref_timing is None:
raise RuntimeError(f"meta.json trace labels required in {reference_name}")
dt = timing_abs_center(snapshot_dir) - ref_timing["center_abs_s"]
horizontal = DIRECTION_SIGN * speed * dt
off_dy = float(horizontal * np.sin(pitch))
off_dz = float(horizontal * np.cos(pitch))
# Scene displacement = difference of two transformed points; the ground
# lift cancels, leaving essentially (0, ~0, horizontal).
disp_y = (off_dy * cos_p - off_dz * sin_p + ground_depth) - ground_depth
disp_z = off_dy * sin_p + off_dz * cos_p
gy_s = gy + disp_y
gz_s = gz + disp_z
depth_shift = off_dz
grid_xx, grid_zz = np.meshgrid(gx, gz_s)
ground_plane_y = (-ground_depth + COMP_REF_DEPTH * sin_p) / cos_p
def pair_geo_ref(i, j):
xc = (x_tx[i] + x_rx[j]) / 2.0
r_tx = np.sqrt((xc - x_tx[i]) ** 2 + (ground_plane_y - y_tx[i]) ** 2 + (COMP_REF_DEPTH - z_tx[i]) ** 2)
r_rx = np.sqrt((xc - x_rx[j]) ** 2 + (ground_plane_y - y_rx[j]) ** 2 + (COMP_REF_DEPTH - z_rx[j]) ** 2)
return 1.0 / (r_tx * r_rx + 1e-12) + 1e-30
geo_ref = {(i, j): pair_geo_ref(i, j) for i in range(n_tx) for j in range(n_rx)}
dx2_tx = [(grid_xx - x_tx[i]) ** 2 for i in range(n_tx)]
dx2_rx = [(grid_xx - x_rx[j]) ** 2 for j in range(n_rx)]
acc_vol = np.zeros((gy.size, gz.size, gx.size), dtype=np.complex128)
cnt_vol = np.zeros(acc_vol.shape, dtype=float)
for iy, y_plane in enumerate(gy_s):
y_rad, z_rad = scene_to_radar_yz(y_plane, grid_zz)
depth_ok = (z_rad - depth_shift >= min_depth) & (z_rad - depth_shift <= max_depth)
r_tx_all = [np.sqrt(dx2_tx[i] + (y_rad - (y_tx[i] + off_dy)) ** 2 + (z_rad - (z_tx[i] + off_dz)) ** 2)
for i in range(n_tx)]
r_rx_all = [np.sqrt(dx2_rx[j] + (y_rad - (y_rx[j] + off_dy)) ** 2 + (z_rad - (z_rx[j] + off_dz)) ** 2)
for j in range(n_rx)]
acc = np.zeros_like(grid_xx, dtype=np.complex128)
cnt = np.zeros_like(grid_xx, dtype=float)
for i in range(n_tx):
for j in range(n_rx):
key = (i, j)
if key not in ascan_h:
continue
r_tx, r_rx = r_tx_all[i], r_rx_all[j]
tau = (r_tx + r_rx) / wave_speed
t_axis = ascan_t[key]
valid = depth_ok & (tau >= t_axis[0]) & (tau <= t_axis[-1])
h = interp_complex(tau, t_axis, ascan_h[key])
ratio = (1.0 / (r_tx * r_rx + 1e-12) + 1e-30) / geo_ref[key]
weight = np.clip(1.0 / (ratio ** COMP_RANGE_POWER + 1e-12), 0.0, COMP_RANGE_WEIGHT_MAX)
acc += np.where(valid, h * weight, 0.0 + 0.0j)
cnt += valid
acc_vol[iy] = acc
cnt_vol[iy] = cnt
volume = np.abs(acc_vol / (cnt_vol + 1e-12))
z_rad_planes = np.asarray([scene_to_radar_yz(yp, grid_zz)[1] - depth_shift for yp in gy_s])
volume = np.where((z_rad_planes >= min_depth) & (z_rad_planes <= max_depth), volume, 0.0)
ref_planes = (gy_s >= gy_s[0] + SAR_NORM_EDGE_MARGIN_M) & (gy_s <= gy_s[-1] - SAR_NORM_EDGE_MARGIN_M)
level_raw = volume[ref_planes].max() + 1e-12
volume = gaussian_filter(volume / level_raw, sigma=(0.0, BP_SMOOTH_SIGMA, BP_SMOOTH_SIGMA))
level_smooth = volume[ref_planes].max() + 1e-12
volume = np.clip(volume / level_smooth, 0.0, None).astype(np.float32)
return {"volume": volume, "level": float(level_raw * level_smooth), "x": gx, "y": gy, "z": gz}
def main() -> int:
parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
parser.add_argument("dataset_root", type=Path,
help="extracted dataset dir (config_profile.json + preprocessed/)")
parser.add_argument("reference_npz", type=Path, help="reference *_volume.npz from the notebook")
args = parser.parse_args()
ref = np.load(args.reference_npz)
meta = json.loads(str(ref["meta"]))
print(f"snapshot: {meta['snapshot']}")
result = compute_volume(args.dataset_root, meta)
volume, ref_volume = result["volume"], ref["volume"]
axes_ok = all(np.allclose(result[k], ref[f"{k}_m"]) for k in ("x", "y", "z"))
level_rel = abs(result["level"] - meta["level"]) / abs(meta["level"])
print(f"axes match: {axes_ok}; level rel diff: {level_rel:.3e} "
f"(mine {result['level']:.9f}, ref {meta['level']:.9f})")
if volume.shape != ref_volume.shape:
print(f"FAIL: shape mismatch {volume.shape} vs {ref_volume.shape}")
return 1
diff = np.abs(volume.astype(np.float64) - ref_volume.astype(np.float64))
argmax = np.unravel_index(int(np.argmax(volume)), volume.shape)
ref_argmax = np.unravel_index(int(np.argmax(ref_volume)), ref_volume.shape)
print(f"max |diff| = {diff.max():.3e}, mean = {diff.mean():.3e}")
print(f"argmax: mine {argmax}, ref {ref_argmax}")
ok = axes_ok and diff.max() <= TOLERANCE and argmax == ref_argmax and level_rel <= TOLERANCE
print("OK" if ok else "FAIL")
return 0 if ok else 1
if __name__ == "__main__":
sys.exit(main())
@@ -0,0 +1,89 @@
"""Compare a sar3d volume dumped by the C++ processor with a notebook npz export.
Usage:
python sar3d_volume_compare.py <dump_base> <reference_volume.npz>
<dump_base> is the path prefix written by `sar3d_parity --dump-dir` (or by the
processor with `sar3d_volume_dump_dir` set): <dump_base>.f32 holds the raw
float32 volume [iy][iz][ix], <dump_base>.json its axes and level.
The C++ path accumulates the backprojection in float32 (CUDA kernel), so the
comparison is tolerance-based rather than bitwise; geometry is computed in
double, keeping the expected difference at the 1e-5 level. Voxels within the
smoothing radius of the depth-mask boundary are reported separately: there the
float32 knife-edge can flip whole boundary voxels (the same effect exists
between two float64 runs with shifted grids).
"""
from __future__ import annotations
import argparse
import json
import sys
from pathlib import Path
import numpy as np
MAX_DIFF_TOLERANCE = 1e-4
LEVEL_REL_TOLERANCE = 1e-4
# Gaussian smoothing radius (sigma 1.5 -> radius 6) around the depth-mask edge.
EDGE_EXCLUSION_VOXELS = 7
def main() -> int:
parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
parser.add_argument("dump_base", type=Path, help="path prefix of the .f32/.json dump")
parser.add_argument("reference_npz", type=Path, help="reference *_volume.npz from the notebook")
args = parser.parse_args()
header = json.loads(Path(str(args.dump_base) + ".json").read_text())
shape = (header["ny"], header["nz"], header["nx"])
volume = np.fromfile(str(args.dump_base) + ".f32", dtype=np.float32).reshape(shape)
ref = np.load(args.reference_npz)
ref_volume = ref["volume"]
meta = json.loads(str(ref["meta"]))
print(f"snapshot: {meta['snapshot']}")
if volume.shape != ref_volume.shape:
print(f"FAIL: shape mismatch {volume.shape} vs {ref_volume.shape}")
return 1
axes_ok = all(
np.allclose(np.asarray(header[key]), ref[key])
for key in ("x_m", "y_m", "z_m")
)
level_rel = abs(header["level"] - meta["level"]) / abs(meta["level"])
print(f"axes match: {axes_ok}; level rel diff: {level_rel:.3e} "
f"(dump {header['level']:.9f}, ref {meta['level']:.9f})")
diff = np.abs(volume.astype(np.float64) - ref_volume.astype(np.float64))
# Depth-mask boundary band: voxels whose (iy, iz) sit within the smoothing
# radius of a masked/unmasked transition. Reconstructed from the reference
# parameters exactly like the processor computes the mask.
pitch = np.radians(float(meta.get("ground_pitch_deg") or meta["look_angle_deg"]))
y0 = np.asarray(ref["y_m"]) - float(meta["ground_depth_m"])
z_rad = -y0[:, None] * np.sin(pitch) + np.asarray(ref["z_m"])[None, :] * np.cos(pitch)
mask = (z_rad >= float(meta["min_depth_m"])) & (z_rad <= float(meta["max_depth_m"]))
edge = np.zeros_like(mask)
r = EDGE_EXCLUSION_VOXELS
for shift in range(-r, r + 1):
edge |= np.roll(mask, shift, axis=1) != mask
interior = ~edge
argmax = np.unravel_index(int(np.argmax(volume)), volume.shape)
ref_argmax = np.unravel_index(int(np.argmax(ref_volume)), ref_volume.shape)
interior_diff = diff[interior[:, :, None] & np.ones(volume.shape, dtype=bool)]
print(f"max |diff| (all voxels) = {diff.max():.3e}")
print(f"max |diff| (off depth-mask edge) = {interior_diff.max():.3e}, "
f"mean = {interior_diff.mean():.3e}")
print(f"argmax: dump {argmax}, ref {ref_argmax}")
ok = (axes_ok and level_rel <= LEVEL_REL_TOLERANCE and argmax == ref_argmax
and interior_diff.max() <= MAX_DIFF_TOLERANCE)
print("OK" if ok else "FAIL")
return 0 if ok else 1
if __name__ == "__main__":
sys.exit(main())