awe/radar_system
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

sar3d: port the Multi_antenna_3D_clean_v1 notebook algorithm

Browse Source 
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>
This commit is contained in:
awe
2026-09-25 16:42:01 +03:00
co-authored by Claude Fable 5
parent e9a0b43ed3
commit b8cb572bc4
17 changed files with 1264 additions and 353 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Makefile
+6 -3
View File
@@ -145,11 +145,14 @@ all: $(TARGETS)
kamil_adc_collector: $(BIN_DIR)/kamil_adc_collector kamil_adc_collector: $(BIN_DIR)/kamil_adc_collector
# On-demand parity check of the sar3d processor against the gpr_cpp # On-demand parity check of the sar3d processor against the reference notebook
# standalone_bp example (same .npy dataset). Not part of `all`. # 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)/, \ SAR3D_PARITY_OBJS := $(addprefix $(BUILD_DIR)/, \
data_acq_and_processing/processing/tools/sar3d_parity_main.o \ 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 sar3d_parity: $(BIN_DIR)/sar3d_parity
data_acq_and_processing/gpr_cpp/README.md
+28 -30
View File
@@ -1,23 +1,24 @@
# gpr_cpp — C++/CUDA счёт-ядро MIMO GPR 3D SAR # gpr_cpp — C++/CUDA счёт-ядро MIMO GPR 3D SAR
> Вендорено из репозитория VNA_AlGo (`gpr_cpp/`, коммит `544dc59`); python-биндинги > Инфраструктура (сборка, NVRTC, PocketFFT, npy) вендорена из репозитория
> (`bindings/`, pybind11) при переносе не взяты — сборка здесь только C++ > VNA_AlGo (`gpr_cpp/`, коммит `544dc59`) без python-биндингов
> (`GPR_BUILD_PYTHON=OFF`). Собирается из корневого Makefile radar_system > (`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`). > (цель `$(GPR_CPP_LIB)` вызывает `build_cmake.sh`).
Порт вычислительного ядра из Python/CuPy (`gpr_core`) в самостоятельную ## Что внутри
C++/CUDA библиотеку для интеграции в C++ проект.
Алгоритм — см. `ALGORITHM.md` в репозитории VNA_AlGo.
## Что внутри (граница C++ ↔ Python) **В C++ (`libgpr_core_cpp`):** сегментированный тайминг кадра по per-trace
меткам захвата (`compute_pair_timestamps`), препроцессинг A-сканов (bg-subtract,
**В C++ (`libgpr_core_cpp`):** геометрия, препроцессинг A-сканов (bg-subtract, полная фазовая коррекция движения intra+pair до IFFT, k0-вставка + Ханн + IFFT
фазовая коррекция движения, k0-вставка + Ханн + IFFT через PocketFFT, нормировка), через PocketFFT, нормировка пар), полнообластная ground-сетка
3D SAR апертурный объём (одно CUDA-ядро на всю апертуру). Возвращает **сырые** (`build_full_region_grid`), 3D SAR апертурный объём (одно CUDA-ядро на всю
`acc`(complex64)/`cnt`(float32). апертуру; геометрия в double, A-сканы float32) и финализация
(`finish_volume`: |acc/cnt| → маска глубины → нормировка без Y-краёв → гаусс
**Остаётся в Python:** загрузка `.npy`, финализация объёма (abs+gate+gaussian+ по (Z,X) → нормировка; возвращает объём [0,1] + level).
нормировка), CLEAN-детекция, трек, оценка скорости, plotly.
## Ключевое решение: CUDA без nvcc (NVRTC + Driver API) ## Ключевое решение: CUDA без nvcc (NVRTC + Driver API)
@@ -42,14 +43,14 @@ libnvrtc прописан — бинарники/модуль грузятся
## Публичный C++ API ([include/gpr/](include/gpr/)) ## Публичный C++ API ([include/gpr/](include/gpr/))
```cpp ```cpp
#include "gpr/config.hpp" // GeometryConfig/BPConfig/MotionConfig/SARConfig (дефолты = config.py) #include "gpr/config.hpp" // GeometryConfig/BPConfig/MotionConfig/SARConfig (дефолты = ноутбук)
#include "gpr/ascan.hpp" // preprocess_snapshot(...) -> AScanSet #include "gpr/ascan.hpp" // compute_pair_timestamps, preprocess_snapshot(...) -> AScanSet
#include "gpr/sar3d.hpp" // SarEngine::aperture_volume(frames, grid) -> VolumeRaw{acc,cnt} #include "gpr/sar3d.hpp" // SarEngine::aperture_volume -> VolumeRaw; finish_volume -> FinishedVolume
#include "gpr/geometry.hpp" // frame_offsets, frame_weights, build_grid_scene, point_to_scene #include "gpr/geometry.hpp" // frame_offsets, frame_weights, build_full_region_grid
gpr::SarEngine engine(geom, bp, mot, sar); // компилит ядро под устройство gpr::SarEngine engine(geom, bp, mot, sar); // компилит ядро под устройство
auto vol = engine.aperture_volume(frames, grid); // один запуск ядра на апертуру auto raw = engine.aperture_volume(frames, grid); // один запуск ядра на апертуру
// vol.acc / vol.cnt — сырые; финализацию (gaussian/normalize) делает вызывающий auto fin = gpr::finish_volume(raw, grid, mot.pitch_deg, bp, sar); // объём [0,1] + level
``` ```
Пример без Python: [examples/standalone_bp.cpp](examples/standalone_bp.cpp) — грузит Пример без Python: [examples/standalone_bp.cpp](examples/standalone_bp.cpp) — грузит
S21 из `.npy`, строит объём, печатает максимум. 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 acc, cnt = eng.aperture_volume([(a, dx, dy, dz, w), ...], gx, gy, gz) # numpy
``` ```
## Паритет (проверено на x86 T600) ## Паритет с ноутбуком (проверено на x86 T600, датасет 20260922_mine_fast_true)
| Этап | Порог | Факт | Кадры 0034/0035/0036 против референсных `*_volume.npz` ноутбука
|---|---|---| (`make sar3d_parity` + `python_app/scripts/sar3d_volume_compare.py`):
| config дефолты vs dataclass | точно | OK | max|Δ| ≤ 1.6e-6 вне ножевой границы маски глубины, argmax и level совпадают
| A-сканы raw H vs numpy (float64) | <1e-9 | 1.7e-12 | (level rel ≤ 7e-7). Python-эталон той же цепочки:
| A-сканы норм. H (FFT-ULP через медиану) | <1e-4 | 2.5e-6 | `python_app/scripts/sar3d_reference_check.py`.
| **3D объём C++ ядро vs cupy** (float32) | <5e-3 | **3e-7, argmax идентичен** |
Тесты: [../tests/test_parity_cpp.py](../tests/test_parity_cpp.py) (C++ vs cupy).
## Статус ## Статус
data_acq_and_processing/gpr_cpp/examples/standalone_bp.cpp
+21 -10
View File
@@ -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); SarEngine engine(geom, bp, mot, sar);
// держим AScanSet живыми // держим AScanSet живыми
@@ -122,22 +135,20 @@ int main(int argc, char** argv) {
frames.push_back(ApertureFrame{a.get(), offs[f], wts[f]}); frames.push_back(ApertureFrame{a.get(), offs[f], wts[f]});
} }
// центр сетки: под центральным снапшотом на глубине 5.6 м (сцена) // полнообластная ground-сетка (как одиночный 3D в ноутбуке)
auto center_scene = point_to_scene(0.0, 0.0, 5.6, mot, sar); Grid grid = build_full_region_grid(geom, bp, sar, mot.pitch_deg);
Grid grid = build_grid_scene(center_scene, sar);
printf("сетка %dx%dx%d, счёт объёма...\n", grid.NX(), grid.NY(), grid.NZ()); printf("сетка %dx%dx%d, счёт объёма...\n", grid.NX(), grid.NY(), grid.NZ());
VolumeRaw vol = engine.aperture_volume(frames, grid); 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; double vmax = 0; size_t amax = 0;
for (size_t i = 0; i < vol.acc.size(); ++i) { for (size_t i = 0; i < fin.volume.size(); ++i)
double v = std::abs(vol.acc[i]) / (vol.cnt[i] + 1e-12); if (fin.volume[i] > vmax) { vmax = fin.volume[i]; amax = i; }
if (v > vmax) { vmax = v; amax = i; }
}
int NX = vol.NX, NZ = vol.NZ; int NX = vol.NX, NZ = vol.NZ;
int ix = amax % NX, iz = (amax / NX) % NZ, iy = amax / ((size_t)NX * 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", printf(" X=%.1f см, Y=%.1f см, Z=%.1f см\n",
grid.x[ix] * 100, grid.y[iy] * 100, grid.z[iz] * 100); grid.x[ix] * 100, grid.y[iy] * 100, grid.z[iz] * 100);
printf("OK: автономный C++ путь отработал без Python.\n"); printf("OK: автономный C++ путь отработал без Python.\n");
data_acq_and_processing/gpr_cpp/include/gpr/ascan.hpp
+29 -6
View File
@@ -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 #pragma once
#include <map> #include <map>
#include <string> #include <string>
@@ -12,15 +14,34 @@ namespace gpr {
using PairKey = std::pair<int, int>; 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 { struct PairTiming {
int tx_interleave_slot = 0; int tx_interleave_slot = 0;
int tx_interleave_slots = 1; 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>; 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). // A-сканы одного снапшота, общая ось времени у всех пар (assert).
struct AScanSet { struct AScanSet {
@@ -32,7 +53,7 @@ struct AScanSet {
const c128* pair_ptr(int p) const { return H.data() + (size_t)p * n_fft; } 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. // s21/freq/background — по ключам (i_tx,i_rx); background может быть nullptr.
AScanSet preprocess_snapshot( AScanSet preprocess_snapshot(
const std::map<PairKey, std::vector<c128>>& s21, const std::map<PairKey, std::vector<c128>>& s21,
@@ -42,10 +63,12 @@ AScanSet preprocess_snapshot(
const PairTimestamps&, bool normalize = true); 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( std::vector<c128> apply_intra_sweep_phase_correction(
const std::vector<c128>& s21, const std::vector<double>& freq_full, const std::vector<c128>& s21, const std::vector<double>& freq_full,
const PairTiming&, const MotionConfig&, double wave_speed, 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-параметр. // Возвращает t (ось), h (комплексный A-скан). n_fft через out-параметр.
void compute_ascan_bp(const std::vector<c128>& s21, const std::vector<double>& freq, void compute_ascan_bp(const std::vector<c128>& s21, const std::vector<double>& freq,
data_acq_and_processing/gpr_cpp/include/gpr/config.hpp
+44 -35
View File
@@ -1,4 +1,5 @@
// Конфигурация счёт-ядра. Дефолты 1:1 повторяют gpr_core/config.py. // Конфигурация счёт-ядра. Дефолты 1:1 повторяют Multi_antenna_3D_clean_v1.ipynb
// (ячейки 0 и 4) — эталонный ноутбук 3D SAR.
#pragma once #pragma once
#include <algorithm> #include <algorithm>
#include <array> #include <array>
@@ -11,13 +12,13 @@ namespace gpr {
struct GeometryConfig { struct GeometryConfig {
std::vector<int> input_idx{0, 1, 2, 3}; std::vector<int> input_idx{0, 1, 2, 3};
std::vector<int> output_idx{0, 1}; std::vector<int> output_idx{0, 1, 2, 3};
// (x,y,z) метры; порядок по физическим индексам как в config.py // (x,y,z) метры; порядок соответствует input_idx/output_idx
std::vector<std::array<double, 3>> tx_positions{{0.785, 0.0, 0.0}, {-0.785, 0.0, 0.0}}; 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{ 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 eps_r = 1.0;
double bp_plane_y = 0.0;
// --- производные (заполняются finalize) --- // --- производные (заполняются finalize) ---
std::vector<double> x_tx, y_tx, z_tx; // отсортированы по phys-индексу std::vector<double> x_tx, y_tx, z_tx; // отсортированы по phys-индексу
@@ -50,53 +51,61 @@ struct GeometryConfig {
}; };
struct BPConfig { struct BPConfig {
double f_start = 20 * 1e8; // 2 ГГц double f_start = 2.0e9;
double f_stop = 6.0e9; double f_stop = 6.0e9;
double min_depth = 3.2; double min_depth = 2.5;
double max_depth = 10.0; double max_depth = 8.0;
double grid_x_pad = 2.0;
double z_min = 0.1;
int nx = 300;
int nz = 300;
int oversample = 8; int oversample = 8;
bool window = true; bool window = true; // окно Ханна на полосе
bool pair_normalize = true; bool pair_normalize = true;
double pair_norm_percentile = 50.0; double pair_norm_percentile = 50.0;
double pair_norm_depth_min = 3.2; double pair_norm_depth_min = 2.5;
double pair_norm_depth_max = 10.0; double pair_norm_depth_max = 8.0;
double pair_norm_eps = 1e-15; 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_weight_max = 5.0;
double comp_ref_depth = 5.0; double comp_ref_depth = 5.0;
double smooth_sigma = 1.5;
bool bg_subtract = true;
double cf_eps = 1e-12;
}; };
struct MotionConfig { struct MotionConfig {
double speed_m_s = 0.854; double speed_m_s = 1.45;
double look_angle_deg = 6.0; double look_angle_deg = 17.5;
double pitch_deg = 6.0; // Наклон рамы для ground-координат и смещений кадров; в ноутбуке
double tx_sweep_time_s = 0.079; // BP_GROUND_PITCH_DEG=None означает «равен look_angle».
double tx_switch_time_s = 1e-5; double pitch_deg = 17.5;
// frequency_timing_mode всегда interleaved_tx_by_frequency
double direction_sign = +1.0; double direction_sign = +1.0;
// Внутрисвиповая (по частоте) фазовая коррекция движения.
bool apply_freq_phase_correction = true; bool apply_freq_phase_correction = true;
// Грубая межсегментная коррекция (dtau сегмента относительно центра кадра),
// вносится той же фазовой моделью до IFFT (режим 'frequency_domain').
bool apply_pair_motion_correction = true;
double motion_phase_sign = -1.0; double motion_phase_sign = -1.0;
// Fallback-разбиение кадра на сегменты, когда per-trace меток захвата нет.
double frame_time_s(int n_tx) const { int tx_per_segment = 2;
return n_tx * tx_sweep_time_s + std::max(0, n_tx - 1) * tx_switch_time_s;
}
}; };
struct SARConfig { struct SARConfig {
int frame_radius = 6; int frame_radius = 0; // 0 = один снимок (физическая вертикальная база)
double edge_frame_weight = 0.75; double edge_frame_weight = 0.8; // вес крайнего кадра апертуры
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' bool coord_mode_ground = true; // 'ground' vs 'radar_local'
bool use_compensation = true; bool use_compensation = true;
bool normalize_each_frame = true; // Высота НИЖНЕГО ряда антенн над землёй; Y=0 сцены — поверхность земли.
double ground_depth_m = 1.25;
// Полнообластная сетка: X — размах антенн ±x_margin_m, Y — высота над
// землёй, Z — полоса, покрывающая радарные глубины [min,max]_depth на всех
// высотах Y (см. build_full_region_grid).
double x_margin_m = 2.0;
double y_min_m = -0.5;
double y_max_m = 2.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 } // namespace gpr
data_acq_and_processing/gpr_cpp/include/gpr/geometry.hpp
+11 -3
View File
@@ -18,10 +18,18 @@ std::vector<FrameOffset> frame_offsets(const std::vector<double>& frame_times, d
// Косинусные веса апертуры (порт sar_frame_weights); dt берётся относительно центра. // Косинусные веса апертуры (порт sar_frame_weights); dt берётся относительно центра.
std::vector<double> frame_weights(const std::vector<FrameOffset>& offsets, double edge_weight); std::vector<double> frame_weights(const std::vector<FrameOffset>& offsets, double edge_weight);
// Сетка вокруг точки сцены (порт sar_build_grid_scene). // Высота плоскости земли (радарные оси) на радарной глубине depth:
Grid build_grid_scene(const std::array<double, 3>& center_scene, const SARConfig&); // 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). // Полнообластная сетка сцены (порт sar_build_full_region_grid):
// X — размах антенн ±x_margin, Y — [y_min, y_max] над землёй, Z — полоса,
// покрывающая радарные глубины [min_depth, max_depth] на всех высотах Y.
Grid build_full_region_grid(const GeometryConfig&, const BPConfig&, const SARConfig&,
double pitch_deg);
// point_to_scene: оси радара → сцена (поворот на pitch + подъём на ground_depth,
// если ground; Y=0 сцены — земля).
std::array<double, 3> point_to_scene(double x, double y, double z, const MotionConfig&, std::array<double, 3> point_to_scene(double x, double y, double z, const MotionConfig&,
const SARConfig&); const SARConfig&);
data_acq_and_processing/gpr_cpp/include/gpr/sar3d.hpp
+12 -1
View File
@@ -19,13 +19,24 @@ struct ApertureFrame {
double weight; double weight;
}; };
// Сырой результат объёма (Python делает finish: abs+gate+gaussian+normalize). // Сырой результат объёма (до finish_volume).
struct VolumeRaw { struct VolumeRaw {
std::vector<c64> acc; // [NY*NZ*NX] std::vector<c64> acc; // [NY*NZ*NX]
std::vector<float> cnt; // [NY*NZ*NX] std::vector<float> cnt; // [NY*NZ*NX]
int NX = 0, NY = 0, NZ = 0; 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 { class SarEngine {
public: public:
SarEngine(const GeometryConfig&, const BPConfig&, const MotionConfig&, const SARConfig&); SarEngine(const GeometryConfig&, const BPConfig&, const MotionConfig&, const SARConfig&);
data_acq_and_processing/gpr_cpp/src/ascan.cpp
+79 -27
View File
@@ -9,19 +9,69 @@
namespace gpr { namespace gpr {
PairTimestamps compute_pair_timestamps(const MotionConfig& mot, const GeometryConfig& geom) { TxWindows uniform_tx_windows(const std::vector<int>& tx_order, double acq_s, int tx_per_segment) {
// interleaved_tx_by_frequency: слот Tx = порядок в отсортированном output_idx. TxWindows windows;
std::vector<int> tx = geom.output_idx, rx = geom.input_idx; if (tx_order.empty()) return windows;
std::sort(tx.begin(), tx.end()); int per_seg = std::max(1, tx_per_segment);
std::sort(rx.begin(), rx.end()); int n_seg = std::max(1, (int)((tx_order.size() + per_seg - 1) / per_seg));
int n_slots = (int)tx.size(); double t_seg = acq_s / n_seg;
double frame_center = 0.5 * n_slots * mot.tx_sweep_time_s; for (size_t k = 0; k < tx_order.size(); ++k) {
PairTimestamps ts; int seg = (int)k / per_seg;
for (int slot = 0; slot < (int)tx.size(); ++slot) { windows[tx_order[k]] = {seg * t_seg, seg * t_seg + t_seg};
int i_tx = slot; // индекс в отсортированном == логический
for (int j = 0; j < (int)rx.size(); ++j) {
ts[{i_tx, j}] = PairTiming{slot, n_slots, frame_center};
} }
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 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; return ts;
} }
@@ -37,7 +87,7 @@ static double median_of_diff(const std::vector<double>& f) {
std::vector<c128> 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 std::vector<c128>& s21, const std::vector<double>& freq_full,
const PairTiming& info, const MotionConfig& cfg, double wave_speed, 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; std::vector<c128> out = s21;
int n_full = (int)freq_full.size(); 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 (idx.empty()) throw std::runtime_error("phase corr: пустая полоса");
if (n_full < 2) return out; 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; int slot = info.tx_interleave_slot;
double dt_base = cfg.tx_sweep_time_s / (n_full - 1); double dt_event = info.segment_time_s / ((double)n_slots * n_full);
// t_all по всему свипу; t_center — середина ПОЛНОГО свипа (как в preprocess.py) double cos_theta = std::cos(cfg.look_angle_deg * M_PI / 180.0);
double t_first = (double)(n_slots * 0 + slot) * dt_base; double pair_dtau = cfg.apply_pair_motion_correction ? info.dtau_pair_s : 0.0;
double t_last = (double)(n_slots * (n_full - 1) + slot) * dt_base; if (!cfg.apply_freq_phase_correction && pair_dtau == 0.0) return out;
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;
for (int i : idx) { for (int i : idx) {
double t_abs = (double)(n_slots * i + slot) * dt_base; double dtau_intra = 0.0;
double dt_intra = t_abs - t_center; if (cfg.apply_freq_phase_correction) {
double delta_range = cfg.direction_sign * cfg.speed_m_s * std::cos(theta) * dt_intra; double t_abs = info.segment_start_s + (double)(n_slots * i + slot) * dt_event;
double dtau = 2.0 * delta_range / wave_speed; double dt_intra = t_abs - info.t_center_s;
double phi = 2.0 * M_PI * freq_full[i] * dtau; 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; double ang = cfg.motion_phase_sign * phi;
out[i] *= c128(std::cos(ang), std::sin(ang)); out[i] *= c128(std::cos(ang), std::sin(ang));
} }
@@ -163,7 +215,7 @@ AScanSet preprocess_snapshot(
auto itt = ts.find(key); auto itt = ts.find(key);
if (itt != ts.end()) if (itt != ts.end())
s = apply_intra_sweep_phase_correction(s, fr, itt->second, mot, geom.v, 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<double> t;
std::vector<c128> h; std::vector<c128> h;
compute_ascan_bp(s, fr, bp.f_start, bp.f_stop, bp.window, bp.oversample, t, h); compute_ascan_bp(s, fr, bp.f_start, bp.f_stop, bp.window, bp.oversample, t, h);
data_acq_and_processing/gpr_cpp/src/geometry.cpp
+42 -8
View File
@@ -1,5 +1,6 @@
#include "gpr/geometry.hpp" #include "gpr/geometry.hpp"
#include <algorithm>
#include <cmath> #include <cmath>
namespace gpr { namespace gpr {
@@ -48,20 +49,53 @@ std::array<double, 3> point_to_scene(double x, double y, double z, const MotionC
const SARConfig& sar) { const SARConfig& sar) {
if (!sar.coord_mode_ground) return {x, y, z}; if (!sar.coord_mode_ground) return {x, y, z};
double p = mot.pitch_deg * M_PI / 180.0; 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); double zs = y * std::sin(p) + z * std::cos(p);
return {x, ys, zs}; return {x, ys, zs};
} }
Grid build_grid_scene(const std::array<double, 3>& c, const SARConfig& sar) { double ground_plane_y(double depth, double ground_depth_m, double pitch_rad) {
Grid g; return (-ground_depth_m + depth * std::sin(pitch_rad)) / std::cos(pitch_rad);
auto lin = [](double a, double b, int n, std::vector<double>& out) { }
Grid build_full_region_grid(const GeometryConfig& geom, const BPConfig& bp, const SARConfig& sar,
double pitch_deg) {
// 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); 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); double x_lo = geom.x_tx[0], x_hi = geom.x_tx[0];
lin(c[1] - sar.half_width_m[1], c[1] + sar.half_width_m[1], sar.grid[1], g.y); for (double x : geom.x_tx) { x_lo = std::min(x_lo, x); x_hi = std::max(x_hi, x); }
lin(c[2] - sar.half_width_m[2], c[2] + sar.half_width_m[2], sar.grid[2], g.z); for (double x : geom.x_rx) { x_lo = std::min(x_lo, x); x_hi = std::max(x_hi, x); }
x_lo -= sar.x_margin_m;
x_hi += sar.x_margin_m;
Grid g;
axis(x_lo, x_hi, sar.voxel_m[0], g.x);
axis(sar.y_min_m, sar.y_max_m, sar.voxel_m[1], g.y);
if (sar.coord_mode_ground) {
// Радарная 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);
}
}
axis(z_lo, z_hi, sar.voxel_m[2], g.z);
} else {
axis(bp.min_depth, bp.max_depth, sar.voxel_m[2], g.z);
}
return g; return g;
} }
data_acq_and_processing/gpr_cpp/src/kernel_src.hpp
+44 -36
View File
@@ -1,28 +1,30 @@
// Исходник CUDA-ядра для NVRTC (порт bp_aperture_contribution из gpr_core/kernels.py). // Исходник CUDA-ядра для NVRTC (порт sar_frame_contribution из
// Self-contained complex<float> (без cupy-заголовков). // 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 #pragma once
namespace gpr { namespace gpr {
inline const char* APERTURE_KERNEL_SRC = R"CUDA( inline const char* APERTURE_KERNEL_SRC = R"CUDA(
struct cf { float x, y; }; 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. acc/cnt пишутся один раз.
extern "C" __global__ void bp_aperture_contribution( extern "C" __global__ void bp_aperture_contribution(
const cf* __restrict__ H, const cf* __restrict__ H,
const float t0, const float inv_dt, const int n_fft, const float inv_v, const double t0, const double inv_dt, const int n_fft, const double inv_v,
const float* __restrict__ xt, const float* __restrict__ yt, const float* __restrict__ zt, const double* __restrict__ xt, const double* __restrict__ yt, const double* __restrict__ zt,
const float* __restrict__ xr, const float* __restrict__ yr, const float* __restrict__ zr, const double* __restrict__ xr, const double* __restrict__ yr, const double* __restrict__ zr,
const float* __restrict__ inv_geo_ref, const float* __restrict__ wgt, const double* __restrict__ inv_geo_ref, const float* __restrict__ wgt,
const float* __restrict__ gx, const float* __restrict__ gy, const float* __restrict__ gz, 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 int NX, const int NY, const int NZ, const int n_contrib,
const float cos_p, const float sin_p, const double cos_p, const double sin_p, const double ground_depth,
const float min_depth, const float max_depth, const double min_depth, const double max_depth,
const float comp_pow, const float w_max, const int use_comp, const double comp_pow, const double w_max, const int use_comp,
cf* __restrict__ acc, float* __restrict__ cnt) cf* __restrict__ acc, float* __restrict__ cnt)
{ {
const long total = (long)NX * NY * NZ; const long total = (long)NX * NY * NZ;
@@ -31,39 +33,45 @@ extern "C" __global__ void bp_aperture_contribution(
int ix = gid % NX; int ix = gid % NX;
int iz = (gid / NX) % NZ; int iz = (gid / NX) % NZ;
int iy = gid / ((long)NX * NZ); int iy = gid / ((long)NX * NZ);
float x = gx[ix]; double x = gx[ix];
float ys = gy[iy], zs = gz[iz]; double zs = gz[iz];
float y_rad = ys * cos_p + zs * sin_p; // Сцена → оси радара: Y=0 сцены — земля, ground_depth ниже нижнего ряда
float z_rad = -ys * sin_p + zs * cos_p; // антенн (порт 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); 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; float c = 0.f;
for (int i = 0; i < n_contrib; ++i) { for (int i = 0; i < n_contrib; ++i) {
float w = wgt[i]; double dxt = x - xt[i], dyt = y_rad - yt[i], dzt = z_rad - zt[i];
float 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];
float dxr = x - xr[i], dyr = y_rad - yr[i], dzr = z_rad - zr[i]; double Rt = sqrt(dxt*dxt + dyt*dyt + dzt*dzt);
float Rt = sqrtf(dxt*dxt + dyt*dyt + dzt*dzt); double Rr = sqrt(dxr*dxr + dyr*dyr + dzr*dzr);
float Rr = sqrtf(dxr*dxr + dyr*dyr + dzr*dzr); double pos = ((Rt + Rr) * inv_v - t0) * inv_dt;
float pos = ((Rt + Rr) * inv_v - t0) * inv_dt; bool inside = (pos >= 0.0) && (pos <= (double)(n_fft - 1));
bool inside = (pos >= 0.f) && (pos <= (float)(n_fft - 1)); if (!inside || !depth_ok) continue;
if (!inside) continue; int k = (int)floor(pos);
int k = (int)floorf(pos);
if (k < 0) k = 0; if (k < 0) k = 0;
if (k > n_fft - 2) k = n_fft - 2; 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; 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) { if (use_comp) {
float ratio = (1.f/(Rt*Rr + 1e-12f) + 1e-30f) * inv_geo_ref[i]; double ratio = (1.0/(Rt*Rr + 1e-12) + 1e-30) * inv_geo_ref[i];
float denom = (comp_pow == 0.25f) ? sqrtf(sqrtf(ratio)) : powf(ratio, comp_pow); double ww = 1.0/(pow(ratio, comp_pow) + 1e-12);
float ww = 1.f/(denom + 1e-12f); ww = fmin(fmax(ww, 0.0), w_max);
ww = fminf(fmaxf(ww, 0.f), w_max); w *= ww;
hv = ww * hv;
} }
if (depth_ok) { a = a + (w * hv); c += w; } ax += w * (double)hx;
ay += w * (double)hy;
c += wgt[i];
} }
acc[gid] = a; acc[gid].x = (float)ax;
acc[gid].y = (float)ay;
cnt[gid] = c; cnt[gid] = c;
} }
} }
data_acq_and_processing/gpr_cpp/src/sar3d.cpp
+139 -25
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(); int NX = grid.NX(), NY = grid.NY(), NZ = grid.NZ();
size_t nvox = (size_t)NX * NY * 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<c64> H_stack((size_t)n_contrib * n_fft);
std::vector<float> xt(n_contrib), yt(n_contrib), zt(n_contrib); std::vector<double> xt(n_contrib), yt(n_contrib), zt(n_contrib);
std::vector<float> xr(n_contrib), yr(n_contrib), zr(n_contrib); std::vector<double> xr(n_contrib), yr(n_contrib), zr(n_contrib);
std::vector<float> igr(n_contrib), wgt(n_contrib); std::vector<double> igr(n_contrib);
std::vector<float> wgt(n_contrib);
int idx = 0; int idx = 0;
for (int f = 0; f < n_frames; ++f) { for (int f = 0; f < n_frames; ++f) {
const AScanSet& a = *frames[f].ascans; const AScanSet& a = *frames[f].ascans;
@@ -49,49 +56,47 @@ VolumeRaw SarEngine::aperture_volume(const std::vector<ApertureFrame>& frames, c
const c128* src = a.pair_ptr(p); const c128* src = a.pair_ptr(p);
c64* dst = H_stack.data() + (size_t)idx * n_fft; 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()); 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); xt[idx] = geom_.x_tx[i] + dx; yt[idx] = geom_.y_tx[i] + dy;
zt[idx] = (float)(geom_.z_tx[i] + dz); zt[idx] = geom_.z_tx[i] + dz;
xr[idx] = (float)(geom_.x_rx[j] + dx); yr[idx] = (float)(geom_.y_rx[j] + dy); xr[idx] = geom_.x_rx[j] + dx; yr[idx] = geom_.y_rx[j] + dy;
zr[idx] = (float)(geom_.z_rx[j] + dz); zr[idx] = geom_.z_rx[j] + dz;
igr[idx] = (float)(1.0 / pair_geo_ref(geom_, bp_, i, j, geom_.bp_plane_y)); igr[idx] = 1.0 / pair_geo_ref(geom_, bp_, i, j, geo_ref_y);
wgt[idx] = w; wgt[idx] = w;
++idx; ++idx;
} }
} }
std::vector<float> gx(grid.x.begin(), grid.x.end()); double cos_p = 1.0, sin_p = 0.0, ground_depth = 0.0;
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;
if (sar_.coord_mode_ground) { if (sar_.coord_mode_ground) {
double p = mot_.pitch_deg * M_PI / 180.0; cos_p = std::cos(pitch);
cos_p = (float)std::cos(p); sin_p = (float)std::sin(p); 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<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<double> dxt(n_contrib), dyt(n_contrib), dzt(n_contrib);
DeviceBuffer<float> dxr(n_contrib), dyr(n_contrib), dzr(n_contrib); DeviceBuffer<double> dxr(n_contrib), dyr(n_contrib), dzr(n_contrib);
DeviceBuffer<float> digr(n_contrib), dwgt(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); 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); 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); digr.upload(igr.data(), n_contrib); dwgt.upload(wgt.data(), n_contrib);
DeviceBuffer<float> dgx(NX), dgy(NY), dgz(NZ); DeviceBuffer<double> dgx(NX), dgy(NY), dgz(NZ);
dgx.upload(gx.data(), NX); dgy.upload(gy.data(), NY); dgz.upload(gz.data(), NZ); dgx.upload(grid.x.data(), NX); dgy.upload(grid.y.data(), NY); dgz.upload(grid.z.data(), NZ);
DeviceBuffer<c64> dacc(nvox); DeviceBuffer<c64> dacc(nvox);
DeviceBuffer<float> dcnt(nvox); DeviceBuffer<float> dcnt(nvox);
float f_t0 = (float)t0, f_inv_dt = (float)(1.0 / dt), f_inv_v = (float)(1.0 / geom_.v); double inv_dt = 1.0 / dt, inv_v = 1.0 / geom_.v;
float min_d = (float)bp_.min_depth, max_d = (float)bp_.max_depth; double min_d = bp_.min_depth, max_d = bp_.max_depth;
float comp_pow = (float)bp_.comp_range_power, w_max = (float)bp_.comp_weight_max; double comp_pow = bp_.comp_range_power, w_max = bp_.comp_weight_max;
int use_comp = sar_.use_compensation ? 1 : 0; int use_comp = sar_.use_compensation ? 1 : 0;
void* args[] = { void* args[] = {
dH.addr(), &f_t0, &f_inv_dt, &n_fft, &f_inv_v, dH.addr(), &t0, &inv_dt, &n_fft, &inv_v,
dxt.addr(), dyt.addr(), dzt.addr(), dxr.addr(), dyr.addr(), dzr.addr(), dxt.addr(), dyt.addr(), dzt.addr(), dxr.addr(), dyr.addr(), dzr.addr(),
digr.addr(), dwgt.addr(), dgx.addr(), dgy.addr(), dgz.addr(), digr.addr(), dwgt.addr(), dgx.addr(), dgy.addr(), dgz.addr(),
&NX, &NY, &NZ, &n_contrib, &cos_p, &sin_p, &min_d, &max_d, &NX, &NY, &NZ, &n_contrib, &cos_p, &sin_p, &ground_depth, &min_d, &max_d,
&comp_pow, &w_max, &use_comp, dacc.addr(), dcnt.addr()}; &comp_pow, &w_max, &use_comp, dacc.addr(), dcnt.addr()};
int block = 256; int block = 256;
@@ -108,4 +113,113 @@ VolumeRaw SarEngine::aperture_volume(const std::vector<ApertureFrame>& frames, c
return out; 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;
}
} // namespace gpr } // namespace gpr
data_acq_and_processing/processing/data_processor/include/processing_live_config.hpp
+33 -13
View File
@@ -90,21 +90,41 @@ struct ProcessingLiveConfig {
float gpr_visible_x_max_m = 2.0F; float gpr_visible_x_max_m = 2.0F;
float gpr_visible_z_min_m = 0.0F; float gpr_visible_z_min_m = 0.0F;
float gpr_visible_z_max_m = 14.0F; float gpr_visible_z_max_m = 14.0F;
// 3D SAR (sar3d processor) parameters; defaults mirror VNA_AlGo gpr_core. // 3D SAR (sar3d processor) parameters; defaults mirror the reference
// Frequency band, depth gate, antenna selection and motion (speed, look // notebook Multi_antenna_3D_clean_v1.ipynb (cells 0/4). Frequency band,
// angle, direction, phase correction) are shared with the gpr_* fields. // depth gate, antenna selection and motion (speed, look angle, direction,
std::uint32_t sar3d_frame_radius = 6U; // phase correction) are shared with the gpr_* fields.
float sar3d_edge_frame_weight = 0.75F; std::uint32_t sar3d_frame_radius = 0U;
float sar3d_pitch_deg = 6.0F; float sar3d_edge_frame_weight = 0.8F;
float sar3d_center_depth_m = 5.6F; // Frame pitch for ground coordinates and inter-frame offsets. Override
float sar3d_half_width_x_m = 1.8F; // disabled -> the pitch follows gpr_look_angle_deg (notebook
float sar3d_half_width_y_m = 1.0F; // BP_GROUND_PITCH_DEG = None).
float sar3d_half_width_z_m = 1.8F; bool sar3d_pitch_override_enabled = false;
std::uint32_t sar3d_grid_nx = 150U; float sar3d_pitch_deg = 17.5F;
std::uint32_t sar3d_grid_ny = 120U; // Height of the LOWER antenna row above the ground; scene Y=0 is the ground.
std::uint32_t sar3d_grid_nz = 160U; float sar3d_ground_depth_m = 1.25F;
// Full-region grid: X spans the antennas +- x_margin, Y is height above the
// ground, Z covers radar depths [gpr_min_depth, gpr_max_depth] over the Y
// range. Voxel pitch per axis.
float sar3d_x_margin_m = 2.0F;
float sar3d_y_min_m = -0.5F;
float sar3d_y_max_m = 2.0F;
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_use_compensation = true;
bool sar3d_normalize_pairs = 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: // Bottom-view projection angle around the vertical Y axis, degrees:
// 0 = front view (X horizontal), 90 = side view (Z horizontal), // 0 = front view (X horizontal), 90 = side view (Z horizontal),
// intermediate values project along the rotated horizontal direction. // intermediate values project along the rotated horizontal direction.
data_acq_and_processing/processing/data_processor/src/processing_live_config.cpp
+21 -18
View File
@@ -392,10 +392,16 @@ void apply_legacy_gpr_algorithm_alias(ProcessingLiveConfig& config, const std::s
for (const auto& [key, target] : { for (const auto& [key, target] : {
std::pair{"sar3d_edge_frame_weight", &config.sar3d_edge_frame_weight}, std::pair{"sar3d_edge_frame_weight", &config.sar3d_edge_frame_weight},
std::pair{"sar3d_pitch_deg", &config.sar3d_pitch_deg}, std::pair{"sar3d_pitch_deg", &config.sar3d_pitch_deg},
std::pair{"sar3d_center_depth_m", &config.sar3d_center_depth_m}, std::pair{"sar3d_ground_depth_m", &config.sar3d_ground_depth_m},
std::pair{"sar3d_half_width_x_m", &config.sar3d_half_width_x_m}, std::pair{"sar3d_x_margin_m", &config.sar3d_x_margin_m},
std::pair{"sar3d_half_width_y_m", &config.sar3d_half_width_y_m}, std::pair{"sar3d_y_min_m", &config.sar3d_y_min_m},
std::pair{"sar3d_half_width_z_m", &config.sar3d_half_width_z_m}, std::pair{"sar3d_y_max_m", &config.sar3d_y_max_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}, std::pair{"sar3d_bottom_view_deg", &config.sar3d_bottom_view_deg},
}) { }) {
if (const auto found = root.find(key); found != root.end()) { if (const auto found = root.find(key); found != root.end()) {
@@ -405,26 +411,23 @@ void apply_legacy_gpr_algorithm_alias(ProcessingLiveConfig& config, const std::s
*target = static_cast<float>(found->get<double>()); *target = static_cast<float>(found->get<double>());
} }
} }
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");
}
}
for (const auto& [key, target] : { for (const auto& [key, target] : {
std::pair{"sar3d_grid_nx", &config.sar3d_grid_nx}, std::pair{"sar3d_pitch_override_enabled", &config.sar3d_pitch_override_enabled},
std::pair{"sar3d_grid_ny", &config.sar3d_grid_ny}, std::pair{"sar3d_use_compensation", &config.sar3d_use_compensation},
std::pair{"sar3d_grid_nz", &config.sar3d_grid_nz}, std::pair{"sar3d_normalize_pairs", &config.sar3d_normalize_pairs},
}) { }) {
if (const auto found = root.find(key); found != root.end()) { if (const auto found = root.find(key); found != root.end()) {
*target = parse_u32_number(*found, std::string("processing.") + key);
}
}
if (const auto found = root.find("sar3d_use_compensation"); found != root.end()) {
if (!found->is_boolean()) { if (!found->is_boolean()) {
throw std::runtime_error("processing.sar3d_use_compensation must be bool"); throw std::runtime_error(std::string("processing.") + key + " must be bool");
} }
config.sar3d_use_compensation = found->get<bool>(); *target = 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");
}
config.sar3d_normalize_pairs = found->get<bool>();
} }
if (const auto found = root.find("ignore_socket_speed"); found != root.end()) { if (const auto found = root.find("ignore_socket_speed"); found != root.end()) {
if (!found->is_boolean()) { if (!found->is_boolean()) {
data_acq_and_processing/processing/processors/src/sar3d_processor.cpp
+161 -55
View File
@@ -4,7 +4,10 @@
#include <bit> #include <bit>
#include <cmath> #include <cmath>
#include <cstdint> #include <cstdint>
#include <cstdio>
#include <fstream>
#include <iostream> #include <iostream>
#include <limits>
#include <map> #include <map>
#include <optional> #include <optional>
#include <stdexcept> #include <stdexcept>
@@ -19,7 +22,6 @@ namespace radar::processing {
namespace { namespace {
constexpr std::size_t kAscanCacheLimit = 64U; constexpr std::size_t kAscanCacheLimit = 64U;
constexpr double kVolumeEps = 1e-12;
void hash_mix(std::size_t& seed, std::uint64_t value) { void hash_mix(std::size_t& seed, std::uint64_t value) {
seed ^= value + 0x9e3779b97f4a7c15ULL + (seed << 6U) + (seed >> 2U); seed ^= value + 0x9e3779b97f4a7c15ULL + (seed << 6U) + (seed >> 2U);
@@ -130,6 +132,7 @@ struct Selection {
bp.pair_normalize = live_config.sar3d_normalize_pairs; bp.pair_normalize = live_config.sar3d_normalize_pairs;
bp.pair_norm_depth_min = bp.min_depth; bp.pair_norm_depth_min = bp.min_depth;
bp.pair_norm_depth_max = bp.max_depth; bp.pair_norm_depth_max = bp.max_depth;
bp.comp_range_power = static_cast<double>(live_config.sar3d_comp_range_power);
return bp; return bp;
} }
@@ -137,9 +140,13 @@ struct Selection {
gpr::MotionConfig motion{}; gpr::MotionConfig motion{};
motion.speed_m_s = static_cast<double>(live_config.gpr_speed_m_s); 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.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.direction_sign = static_cast<double>(live_config.gpr_direction_sign);
motion.apply_freq_phase_correction = live_config.gpr_apply_freq_phase_correction; 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; return motion;
} }
@@ -147,20 +154,35 @@ struct Selection {
gpr::SARConfig sar{}; gpr::SARConfig sar{};
sar.frame_radius = static_cast<int>(std::max<std::uint32_t>(0U, live_config.sar3d_frame_radius)); 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.edge_frame_weight = static_cast<double>(live_config.sar3d_edge_frame_weight);
sar.half_width_m = { sar.ground_depth_m = static_cast<double>(live_config.sar3d_ground_depth_m);
static_cast<double>(live_config.sar3d_half_width_x_m), sar.x_margin_m = static_cast<double>(live_config.sar3d_x_margin_m);
static_cast<double>(live_config.sar3d_half_width_y_m), sar.y_min_m = static_cast<double>(live_config.sar3d_y_min_m);
static_cast<double>(live_config.sar3d_half_width_z_m), sar.y_max_m = static_cast<double>(live_config.sar3d_y_max_m);
}; sar.voxel_m = {
sar.grid = { static_cast<double>(live_config.sar3d_voxel_x_m),
static_cast<int>(std::max<std::uint32_t>(2U, live_config.sar3d_grid_nx)), static_cast<double>(live_config.sar3d_voxel_y_m),
static_cast<int>(std::max<std::uint32_t>(2U, live_config.sar3d_grid_ny)), static_cast<double>(live_config.sar3d_voxel_z_m),
static_cast<int>(std::max<std::uint32_t>(2U, live_config.sar3d_grid_nz)),
}; };
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; sar.use_compensation = live_config.sar3d_use_compensation;
return sar; 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( [[nodiscard]] auto ascan_fingerprint(
const Selection& selection, const Selection& selection,
const gpr::GeometryConfig& geom, const gpr::GeometryConfig& geom,
@@ -184,6 +206,7 @@ struct Selection {
hash_mix(seed, motion.look_angle_deg); hash_mix(seed, motion.look_angle_deg);
hash_mix(seed, motion.direction_sign); 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.apply_freq_phase_correction ? 1U : 0U));
hash_mix(seed, static_cast<std::uint64_t>(motion.tx_per_segment));
return seed; return seed;
} }
@@ -215,6 +238,7 @@ struct Selection {
hash_mix(seed, bp.comp_weight_max); hash_mix(seed, bp.comp_weight_max);
hash_mix(seed, bp.comp_ref_depth); hash_mix(seed, bp.comp_ref_depth);
hash_mix(seed, motion.pitch_deg); 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)); hash_mix(seed, static_cast<std::uint64_t>(sar.use_compensation ? 1U : 0U));
return seed; return seed;
} }
@@ -228,14 +252,80 @@ struct Selection {
return static_cast<double>(collection.monotonic_ns) * 1e-9; return static_cast<double>(collection.monotonic_ns) * 1e-9;
} }
// One Tx's interleaved sweep spans the whole capture window, so the per-Tx // Real frame timing from the per-trace capture stamps: traces measured in one
// sweep time is capture_span / N_tx (same derivation as the coherent BP path). // segment share one window on both producers, so the first window seen per Tx
[[nodiscard]] auto tx_sweep_time_s(const ipc::PreprocessedCollection& collection, int n_tx) -> double { // reconstructs the segments exactly (notebook snapshot_acquisition_timing).
if (collection.capture_end_ns <= collection.capture_start_ns || n_tx <= 0) { // Invalid when any trace carries the zero "unmeasured" sentinel.
return 0.0; 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; std::map<int, std::pair<std::uint64_t, std::uint64_t>> by_tx{};
return span_s / static_cast<double>(n_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( [[nodiscard]] auto build_image_payload(
@@ -285,14 +375,26 @@ struct Sar3dProcessor::State {
const Selection& selection, const Selection& selection,
const gpr::GeometryConfig& geom, const gpr::GeometryConfig& geom,
const gpr::BPConfig& bp, const gpr::BPConfig& bp,
gpr::MotionConfig motion const gpr::MotionConfig& motion,
const std::vector<int>& tx_order
) -> const CacheEntry& { ) -> const CacheEntry& {
if (const auto found = ascan_cache.find(collection.collection_id); found != ascan_cache.end()) { if (const auto found = ascan_cache.find(collection.collection_id); found != ascan_cache.end()) {
return found->second; return found->second;
} }
CacheEntry entry{}; 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<gpr::c128>> s21{};
std::map<gpr::PairKey, std::vector<double>> freq{}; std::map<gpr::PairKey, std::vector<double>> freq{};
@@ -317,9 +419,25 @@ struct Sar3dProcessor::State {
} }
if (s21.size() == selection.inputs.size() * selection.outputs.size()) { if (s21.size() == selection.inputs.size() * selection.outputs.size()) {
motion.tx_sweep_time_s = tx_sweep_time_s(collection, geom.n_tx); gpr::TxWindows windows{};
const auto timestamps = gpr::compute_pair_timestamps(motion, geom); 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 { 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( entry.ascans = std::make_shared<gpr::AScanSet>(gpr::preprocess_snapshot(
s21, freq, /*background=*/nullptr, motion, geom, bp, timestamps, bp.pair_normalize s21, freq, /*background=*/nullptr, motion, geom, bp, timestamps, bp.pair_normalize
)); ));
@@ -399,12 +517,14 @@ auto Sar3dProcessor::process_collection(
} }
window.push_back(&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<std::shared_ptr<gpr::AScanSet>> ascan_sets{};
std::vector<double> frame_times{}; std::vector<double> frame_times{};
ascan_sets.reserve(window.size()); ascan_sets.reserve(window.size());
frame_times.reserve(window.size()); frame_times.reserve(window.size());
for (const auto* entry_collection : window) { 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) { if (entry.ascans == nullptr) {
continue; continue;
} }
@@ -450,14 +570,17 @@ auto Sar3dProcessor::process_collection(
state_->engine_fp = engine_fp; state_->engine_fp = engine_fp;
} }
const auto center_scene = const gpr::Grid grid = gpr::build_full_region_grid(geom, bp, sar, motion.pitch_deg);
gpr::point_to_scene(0.0, 0.0, static_cast<double>(live_config.sar3d_center_depth_m), motion, sar); const gpr::VolumeRaw raw = state_->engine->aperture_volume(frames, grid);
const gpr::Grid grid = gpr::build_grid_scene(center_scene, sar); // Full notebook finalization: |acc/cnt|, depth gate, normalization
const gpr::VolumeRaw volume = state_->engine->aperture_volume(frames, grid); // 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, // Two projections match the legacy render conventions (X across,
// Z range/"вдаль", Y up): // Z range/"вдаль", Y up):
// - top view `sar3d_proj_xz`: max over Y, axes X (horizontal) / Z; // - top view `sar3d_proj_xz`: max over Y, axes X (horizontal) / Z;
@@ -466,12 +589,9 @@ auto Sar3dProcessor::process_collection(
// `sar3d_bottom_view_deg` around Y: 0° = front view (u = X, max // `sar3d_bottom_view_deg` around Y: 0° = front view (u = X, max
// along Z), 90° = side view (u = Z, max along X), intermediate // along Z), 90° = side view (u = Z, max along X), intermediate
// angles project along the rotated perpendicular. // angles project along the rotated perpendicular.
const int nx = volume.NX; const int nx = finished.NX;
const int ny = volume.NY; const int ny = finished.NY;
const int nz = volume.NZ; const int nz = finished.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 double theta_deg = std::clamp(static_cast<double>(live_config.sar3d_bottom_view_deg), 0.0, 90.0); 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; const double theta_rad = theta_deg * M_PI / 180.0;
@@ -501,24 +621,18 @@ auto Sar3dProcessor::process_collection(
const double u_start = u_axis.front(); 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; 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_xz(static_cast<std::size_t>(nz) * nx, 0.0F);
std::vector<float> projection_bottom(static_cast<std::size_t>(ny) * nu, 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 iy = 0; iy < ny; ++iy) {
for (int iz = 0; iz < nz; ++iz) { 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 row = (static_cast<std::size_t>(iy) * nz + iz) * nx;
const std::size_t bottom_row = static_cast<std::size_t>(iy) * nu; const std::size_t bottom_row = static_cast<std::size_t>(iy) * nu;
const double z_sin_t = grid.z[iz] * sin_t; const double z_sin_t = grid.z[iz] * sin_t;
for (int ix = 0; ix < nx; ++ix) { for (int ix = 0; ix < nx; ++ix) {
const std::size_t voxel = row + static_cast<std::size_t>(ix); const float value = finished.volume[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);
auto& xz_cell = projection_xz[static_cast<std::size_t>(iz) * nx + ix]; auto& xz_cell = projection_xz[static_cast<std::size_t>(iz) * nx + ix];
xz_cell = std::max(xz_cell, value); xz_cell = std::max(xz_cell, value);
int bin = ix; int bin = ix;
@@ -534,14 +648,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( results.collection_payloads.push_back(
build_image_payload("sar3d_proj_xz", grid.x, grid.z, std::move(projection_xz)) build_image_payload("sar3d_proj_xz", grid.x, grid.z, std::move(projection_xz))
@@ -552,7 +658,7 @@ auto Sar3dProcessor::process_collection(
ipc::ResultPayload peak_payload{}; ipc::ResultPayload peak_payload{};
peak_payload.processing_name = "sar3d_peak"; peak_payload.processing_name = "sar3d_peak";
peak_payload.kind = ipc::ResultKind::ScalarF32; 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)); results.collection_payloads.push_back(std::move(peak_payload));
} catch (const std::exception& exc) { } catch (const std::exception& exc) {
// A CUDA/NVRTC or data error must not kill the data_processor loop: // A CUDA/NVRTC or data error must not kill the data_processor loop:
data_acq_and_processing/processing/tools/sar3d_parity_main.cpp
+96 -63
View File
@@ -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 // Feeds a preprocessed snapshot of an extracted replay dataset through
// uses, feeds it through Sar3dProcessor (combo mapping, timestamps, A-scan // Sar3dProcessor exactly the way the live pipeline would (geometry and combo
// cache, aperture, finalization) and prints the volume peak and the argmax of // order from the dataset's config_profile.json, per-trace capture windows from
// both projections. Compare against `standalone_bp <dir> <center>` output: // the snapshot's meta.json) and dumps the finished volume for comparison with
// the peak |acc/cnt| and the (X, Y, Z) of the maximum voxel must match. // 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 <algorithm>
#include <cctype> #include <cctype>
#include <complex> #include <complex>
@@ -19,6 +25,8 @@
#include <string> #include <string>
#include <vector> #include <vector>
#include <nlohmann/json.hpp>
#include "processing_live_config.hpp" #include "processing_live_config.hpp"
#include "run_config.hpp" #include "run_config.hpp"
#include "sar3d_processor.hpp" #include "sar3d_processor.hpp"
@@ -26,8 +34,7 @@
namespace { namespace {
// Minimal numpy .npy loader (v1.0, 1-D complex128/complex64/float64/float32), // Minimal numpy .npy loader (v1.0, 1-D complex128/complex64/float64/float32).
// mirroring gpr_cpp/examples/standalone_bp.cpp.
struct NpyArray { struct NpyArray {
std::vector<double> f64{}; std::vector<double> f64{};
std::vector<std::complex<double>> c128{}; std::vector<std::complex<double>> c128{};
@@ -77,43 +84,23 @@ NpyArray load_npy(const std::string& path) {
return array; return array;
} }
// VNA_AlGo reference geometry (gpr_core/config.py) keyed by switch positions. // Live config matching the reference notebook run (npz meta parameters).
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 { auto build_live_config() -> radar::processing::ProcessingLiveConfig {
radar::processing::ProcessingLiveConfig live{}; radar::processing::ProcessingLiveConfig live{};
live.processor_mode = "sar3d"; live.processor_mode = "sar3d";
live.gpr_start_freq_mhz = 2000.0F; live.gpr_start_freq_mhz = 2000.0F;
live.gpr_stop_freq_mhz = 6000.0F; live.gpr_stop_freq_mhz = 6000.0F;
live.gpr_min_depth_m = 3.2F; live.gpr_min_depth_m = 2.5F;
live.gpr_max_depth_m = 10.0F; live.gpr_max_depth_m = 8.0F;
live.gpr_speed_m_s = 0.854F; live.gpr_speed_m_s = 1.45F;
live.gpr_look_angle_deg = 6.0F; live.gpr_look_angle_deg = 17.5F;
live.gpr_direction_sign = 1.0F; live.gpr_direction_sign = 1.0F;
live.gpr_apply_freq_phase_correction = true; live.gpr_apply_freq_phase_correction = true;
return live; return live;
} }
constexpr double kTxSweepTimeS = 0.079;
constexpr int kTxCount = 2;
constexpr int kRxCount = 4;
auto load_collection( auto load_collection(
const radar::config::RunConfig& run_config,
const std::string& root, const std::string& root,
const std::string& snapshot, const std::string& snapshot,
std::uint64_t collection_id std::uint64_t collection_id
@@ -121,32 +108,58 @@ auto load_collection(
radar::ipc::PreprocessedCollection collection{}; radar::ipc::PreprocessedCollection collection{};
collection.collection_id = collection_id; collection.collection_id = collection_id;
// Frame time from the `_ns<...>` folder suffix, identical to standalone_bp. // 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>(),
};
}
}
}
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"); const auto ns_pos = snapshot.rfind("_ns");
if (ns_pos == std::string::npos) { if (ns_pos == std::string::npos) {
throw std::runtime_error("snapshot has no _ns suffix: " + snapshot); throw std::runtime_error("snapshot has no meta.json and no _ns suffix: " + snapshot);
} }
const auto mid_ns = static_cast<std::uint64_t>(std::stoull(snapshot.substr(ns_pos + 3))); last = 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 first = last;
// 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 = first;
collection.capture_start_ns = mid_ns - span_ns / 2U; collection.capture_end_ns = last;
collection.capture_end_ns = mid_ns + span_ns / 2U; collection.monotonic_ns = (first + last) / 2U;
collection.monotonic_ns = mid_ns;
for (int rx = 0; rx < kRxCount; ++rx) { for (const auto& combo : run_config.run_combos) {
for (int tx = 0; tx < kTxCount; ++tx) {
char path[512]; char path[512];
std::snprintf(path, sizeof(path), "%s/%s/i%d_o%d_s21.npy", root.c_str(), snapshot.c_str(), rx, tx); 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); 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); 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); const auto freq = load_npy(path);
radar::ipc::SweepTraceBlock trace{}; radar::ipc::SweepTraceBlock trace{};
trace.combo = radar::ipc::ComboKey{ trace.combo = combo;
.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.frequency_hz.assign(freq.f64.begin(), freq.f64.end());
trace.s21.reserve(s21.c128.size()); trace.s21.reserve(s21.c128.size());
for (const auto& sample : s21.c128) { for (const auto& sample : s21.c128) {
@@ -156,10 +169,11 @@ auto load_collection(
}); });
} }
trace.s11.resize(trace.s21.size()); trace.s11.resize(trace.s21.size());
trace.capture_start_ns = collection.capture_start_ns; if (const auto found = windows.find({combo.input_pos, combo.output_pos}); found != windows.end()) {
trace.capture_end_ns = collection.capture_end_ns; trace.capture_start_ns = found->second.first;
collection.traces.push_back(std::move(trace)); trace.capture_end_ns = found->second.second;
} }
collection.traces.push_back(std::move(trace));
} }
return collection; return collection;
} }
@@ -195,14 +209,21 @@ void print_projection_argmax(
int main(int argc, char** argv) { int main(int argc, char** argv) {
if (argc < 3) { 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; return 1;
} }
const std::string root = argv[1]; const std::string root = argv[1];
const std::string center = argv[2]; 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{}; 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;) { for (dirent* entry = nullptr; (entry = readdir(dir)) != nullptr;) {
const std::string name = entry->d_name; const std::string name = entry->d_name;
if (!name.empty() && (std::isdigit(name[0]) != 0)) { if (!name.empty() && (std::isdigit(name[0]) != 0)) {
@@ -214,28 +235,34 @@ int main(int argc, char** argv) {
std::sort(snapshots.begin(), snapshots.end()); std::sort(snapshots.begin(), snapshots.end());
const auto center_it = std::find(snapshots.begin(), snapshots.end(), center); const auto center_it = std::find(snapshots.begin(), snapshots.end(), center);
if (center_it == snapshots.end()) { 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; return 1;
} }
const auto center_index = static_cast<std::size_t>(center_it - snapshots.begin()); const auto center_index = static_cast<std::size_t>(center_it - snapshots.begin());
const auto live = build_live_config(); auto live = build_live_config();
live.sar3d_volume_dump_dir = dump_dir;
const auto radius = static_cast<std::size_t>(live.sar3d_frame_radius); const auto radius = static_cast<std::size_t>(live.sar3d_frame_radius);
if (center_index < radius || center_index + radius >= snapshots.size()) { if (center_index < radius || center_index + radius >= snapshots.size()) {
std::printf("not enough snapshots around center for radius %zu\n", radius); std::printf("not enough snapshots around center for radius %zu\n", radius);
return 1; return 1;
} }
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;
}
// The processor centers the aperture `radius` frames before the newest // The processor centers the aperture `radius` frames before the newest
// collection, so feeding [center-r .. center+r] in order reproduces the // 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{}; std::vector<radar::ipc::PreprocessedCollection> history{};
for (std::size_t index = center_index - radius; index <= center_index + radius; ++index) { 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()); 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::processing::Sar3dProcessor processor{};
radar::ipc::ResultCollection result{}; radar::ipc::ResultCollection result{};
for (std::size_t index = 0U; index < history.size(); ++index) { for (std::size_t index = 0U; index < history.size(); ++index) {
@@ -253,7 +280,7 @@ int main(int argc, char** argv) {
} }
for (const auto& payload : result.collection_payloads) { for (const auto& payload : result.collection_payloads) {
if (payload.kind == radar::ipc::ResultKind::ScalarF32) { 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)); payload.processing_name.c_str(), static_cast<double>(payload.scalar_value));
} }
} }
@@ -269,6 +296,7 @@ int main(int argc, char** argv) {
std::span<const radar::ipc::PreprocessedCollection>(history.data(), history.size() - 1U); std::span<const radar::ipc::PreprocessedCollection>(history.data(), history.size() - 1U);
for (const float angle_deg : {0.0F, 30.0F, 60.0F, 90.0F}) { for (const float angle_deg : {0.0F, 30.0F, 60.0F, 90.0F}) {
auto live_angle = live; auto live_angle = live;
live_angle.sar3d_volume_dump_dir.clear();
live_angle.sar3d_bottom_view_deg = angle_deg; live_angle.sar3d_bottom_view_deg = angle_deg;
const auto angle_result = const auto angle_result =
processor.process_collection(run_config, history.back(), last_previous, live_angle); processor.process_collection(run_config, history.back(), last_previous, live_angle);
@@ -279,6 +307,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"); std::printf("OK: sar3d processor path completed\n");
return 0; return 0;
} }
python_app/scripts/sar3d_reference_check.py
+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())
python_app/scripts/sar3d_volume_compare.py
+89
View File
@@ -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())