diff --git a/Makefile b/Makefile index b25b138..0578cf9 100644 --- a/Makefile +++ b/Makefile @@ -109,7 +109,7 @@ TARGETS := \ $(BIN_DIR)/data_preprocessor \ $(BIN_DIR)/data_processor -.PHONY: all clean rph_radar_collector check-native check-sar3d +.PHONY: all clean rph_radar_collector check-native check-sar3d check-gpr all: $(TARGETS) @@ -119,19 +119,23 @@ check-native: $(NATIVE_BUILD_DIR)/native_regressions $(BIN_DIR)/collector_regres LD_LIBRARY_PATH=$(abspath $(NATIVE_BUILD_DIR)/rph_sdk) $(BIN_DIR)/collector_regressions $(NATIVE_BUILD_DIR)/native_regressions run_config_examples/run_config_simulator.example.json -SAR3D_TEST_SOURCES := data_acq_and_processing/tests/sar3d_fixture.cpp \ +BACKPROJECTION_TEST_SOURCES := data_acq_and_processing/tests/backprojection_fixture.cpp \ data_acq_and_processing/common_cpp/ipc/src/shared_types.cpp \ data_acq_and_processing/processing/processors/src/sar3d_processor.cpp \ - data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp -SAR3D_TEST_OBJS := $(addprefix $(NATIVE_BUILD_DIR)/,$(SAR3D_TEST_SOURCES:.cpp=.o)) -DEPFILES += $(SAR3D_TEST_OBJS:.o=.d) + data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp \ + data_acq_and_processing/processing/processors/src/gpr_processor.cpp +BACKPROJECTION_TEST_OBJS := $(addprefix $(NATIVE_BUILD_DIR)/,$(BACKPROJECTION_TEST_SOURCES:.cpp=.o)) +DEPFILES += $(BACKPROJECTION_TEST_OBJS:.o=.d) # The sar3d processor on the real CUDA backend against the notebook reference; # without a usable GPU the test checks the warning and skips the rest. -check-sar3d: $(NATIVE_BUILD_DIR)/sar3d_fixture - SAR3D_FIXTURE="$(abspath $<)" $(PYTHON) data_acq_and_processing/tests/test_sar3d.py +check-sar3d: $(NATIVE_BUILD_DIR)/backprojection_fixture + BACKPROJECTION_FIXTURE="$(abspath $<)" $(PYTHON) data_acq_and_processing/tests/test_backprojection.py Sar3dTest -$(NATIVE_BUILD_DIR)/sar3d_fixture: $(SAR3D_TEST_OBJS) +check-gpr: $(NATIVE_BUILD_DIR)/backprojection_fixture + BACKPROJECTION_FIXTURE="$(abspath $<)" $(PYTHON) data_acq_and_processing/tests/test_backprojection.py GprTest + +$(NATIVE_BUILD_DIR)/backprojection_fixture: $(BACKPROJECTION_TEST_OBJS) $(CXX) $(NATIVE_CXXFLAGS) $^ -o $@ $(LDFLAGS) $(SAR3D_LDFLAGS) $(NATIVE_BUILD_DIR)/data_acq_and_processing/processing/processors/src/sar3d_gpu.o: NATIVE_CXXFLAGS += $(SAR3D_CXXFLAGS) diff --git a/README.md b/README.md index 3c2d6f3..ef46c90 100644 --- a/README.md +++ b/README.md @@ -66,9 +66,10 @@ controls, processing settings and a copy of the live plot through an embedded we - **pass_through** — magnitude and phase of S21 for each combo. - **bscan** — depth-versus-sweep B-scan heatmap. -- **gpr** — coherent back-projection GPR imaging with object detection. -- **legacy_gpr** — ellipse-intersection MIMO localization (point and extended/region modes). -- **sar3d** — 3D SAR of the reference notebook `Multi_antenna_3D_clean_v1.ipynb` +- **gpr** — coherent back-projection on the ground plane, as in `Horns_motion_4x4_ground.py`, + with object detection. Radar height and look angle define the slice in radar X/Z axes; + the wave speed is 3×10⁸ m/s (air). +- **sar3d** — 3D SAR of the reference notebook `Multi_antenna_demonstrate_depthnorm.ipynb` on an NVIDIA GPU through CUDA: front (X-Y), side (Z-Y), top (X-Z) and oblique maximum projections of the imaged volume. Its settings carry the notebook's parameter names; the scene window is @@ -77,6 +78,28 @@ controls, processing settings and a copy of the live plot through an embedded we (default 2 cm). The radar motion is corrected from each trace's capture timestamps, so the sweep timing needs no parameters. +Both coherent GPR and SAR apply `PAIR_DEPTH_EQUALIZE` (on by default) before the +per-pair median normalization. The notebook's 1.5 m sliding median estimates each +pair's depth-dependent level; positive gains bring the pairs to their geometric +mean while preserving phase and the common depth envelope. **Subtract mean of previous +collections** subtracts up to **Mean count** preceding spectra for each pair before +motion correction and the IFFT; zero count disables it. SAR uses a causal aperture: +`SAR_FRAME_RADIUS` previous frames plus the current frame, with the current scene origin. +The processor retains enough history for the aperture and each frame's preceding +mean; transport ring capacity does not shorten this background window. +Scene width `SAR_X_SIZE_M` and length `SAR_Z_SIZE_M` use 0 for the automatic notebook +bounds; positive values resize the axes about their automatic centres. `SAR_Y_RANGE_M` +sets the height limits. The depth gate remains independent of these scene bounds. + +SAR **Zoom (%)** shows the full scene at 100%; increasing it shrinks all four +views uniformly towards the brightest response. **SAR_CLOUD_THRESHOLD** hides +values below the selected fraction of the normalization maximum. Both controls +redraw the existing result without recalculating it. Tx/Rx markers appear on the +front, side and top projections in ground coordinates. A hollow square marks the +peak; its size follows the contiguous half-maximum width along the volume axes +(FWHM, at least one voxel). This measures the imaged response rather than the +reflector's physical dimensions. + The desktop uses a light theme with dark plot canvases. **Mode** stays below the **Actions** buttons while the settings sections scroll; the settings panel and runtime log can be collapsed with their arrow buttons. SAR3D status, peak @@ -88,8 +111,8 @@ first, then input, and wraps to the first recorded pair; it is enabled with at l two pairs. Desktop and browser selectors redraw the cached B-scan immediately. Tx/Rx antenna positions are edited in tables with switch position and X/Y/Z in -metres, on both the desktop and web interface. Coordinate and legacy-GPR medium -changes apply to the running processor and recompute its retained result. Profiles +metres, on both the desktop and web interface. Coordinate changes apply to the +running processor and recompute its retained result. Profiles preserve all three coordinates. A new table row becomes active after X is entered; desktop, web and profile geometry use one Python validator; C++ validates its own input. GPR/SAR views wait for a result marked with the current settings revision, @@ -431,6 +454,7 @@ backup to `radar_system`. Restart the service if it is normally used. QT_QPA_PLATFORM=offscreen .venv/bin/python -m unittest discover -s python_app/tests make check-native make check-sar3d PYTHON=.venv/bin/python +make check-gpr ASAN_OPTIONS=detect_leaks=0 make check-native NATIVE_SANITIZERS=address,undefined ``` @@ -441,24 +465,33 @@ also exercise mutated frames, ring concurrency, locator input and ADC whole-swee Firmware host tests check cart timer rollover and SBUS packing; they do not validate physical timing. Real devices, GPIO and the CUDA backend need their own hardware run. -`check-sar3d` runs the real CUDA processor (`data_acq_and_processing/tests/test_sar3d.py`). +`check-sar3d` runs the real CUDA processor (`data_acq_and_processing/tests/test_backprojection.py`). It compares it with `sar3d_reference.npz`: compact float32/complex64 inputs, the three -maximum projections of the notebook's volume, and its joint peak/level for four -compensation/normalization and timing cases. The expected values come from unchanged -functions in cells 0 and 4 of the original notebook. The same test checks the published +maximum projections of the notebook's volume, its joint peak/level and response +widths for six compensation/normalization, equalization, background, scene and +timing cases. The expected values come from unchanged functions in cells 1 and 14 +of the original notebook, with each frame's previous-collection mean subtracted +from its input spectra by NumPy. The same test checks the published status for invalid input, a result too large for the ring slot, an empty and a very weak signal. Without a usable GPU it checks only the CUDA warning and skips the rest; the notebook itself is not needed to run it. To regenerate the reference, pass the original notebook (it is kept outside the repository): ```bash -.venv/bin/python data_acq_and_processing/tests/generate_sar3d_reference.py \ - --notebook /path/to/Multi_antenna_3D_clean_v1.ipynb +.venv/bin/python data_acq_and_processing/tests/generate_backprojection_reference.py \ + --notebook /path/to/Multi_antenna_demonstrate_depthnorm.ipynb \ + --gpr-script /path/to/Horns_motion_4x4_ground.py ``` -The archive records cell-source SHA256 hashes and NumPy/SciPy versions. Another build +The generator also writes `gpr_reference.npz` from the unchanged Horns coherent BP +functions, with the notebook equalizer inserted before scalar normalization. +`make check-gpr` compares the CPU processor with those images, including motion and +previous-collection subtraction. Neither source file is needed to run the tests. + +The archives record source SHA256 hashes and NumPy/SciPy versions. Another build can be checked with `make check-sar3d BUILD_DIR=/path/to/build` or by running -`test_sar3d.py` with `SAR3D_FIXTURE=/path/to/build/native_tests/sar3d_fixture`. +`test_backprojection.py Sar3dTest` with +`BACKPROJECTION_FIXTURE=/path/to/build/native_tests/backprojection_fixture`. For older GPUs, pass the same `BUILD_DIR` and NVRTC 12 `CUDA_ROOTS` used to build. For development, install `pyright` into the venv, then run diff --git a/data_acq_and_processing/common_cpp/config/include/run_config.hpp b/data_acq_and_processing/common_cpp/config/include/run_config.hpp index 0d79507..9bd5e4b 100644 --- a/data_acq_and_processing/common_cpp/config/include/run_config.hpp +++ b/data_acq_and_processing/common_cpp/config/include/run_config.hpp @@ -92,7 +92,6 @@ struct GprRxGeometry { struct GprConfig { // Initial GPR settings; the data_processor can override them through live config. - float relative_permittivity = 1.0F; std::vector tx_geometry{}; std::vector rx_geometry{}; }; diff --git a/data_acq_and_processing/common_cpp/config/src/run_config.cpp b/data_acq_and_processing/common_cpp/config/src/run_config.cpp index 674295d..fb44aa0 100644 --- a/data_acq_and_processing/common_cpp/config/src/run_config.cpp +++ b/data_acq_and_processing/common_cpp/config/src/run_config.cpp @@ -277,7 +277,6 @@ void validate_combo_key(const ipc::ComboKey& key, const RunConfig& config) { [[nodiscard]] auto parse_gpr_config(const Json& object, const RunConfig& run_config) -> GprConfig { const auto* gpr_obj = as_object(object, "gpr"); GprConfig config{}; - config.relative_permittivity = optional_f32(*gpr_obj, "relative_permittivity", 1.0F); if (const auto* tx_value = optional_field(*gpr_obj, "tx_geometry"); tx_value != nullptr) { const auto* tx_array = as_array(*tx_value, "gpr.tx_geometry"); @@ -375,10 +374,6 @@ void validate_combo_key(const ipc::ComboKey& key, const RunConfig& config) { } // namespace void validate_gpr_config(const GprConfig& config, const RunConfig& run_config) { - if (!std::isfinite(config.relative_permittivity) || !(config.relative_permittivity > 0.0F)) { - throw std::runtime_error("gpr.relative_permittivity must be > 0"); - } - std::unordered_set seen_tx_positions{}; for (const auto& entry : config.tx_geometry) { if (!std::isfinite(entry.x_m) || !std::isfinite(entry.y_m) || !std::isfinite(entry.z_m)) { diff --git a/data_acq_and_processing/processing/data_processor/include/data_processor.hpp b/data_acq_and_processing/processing/data_processor/include/data_processor.hpp index e209e57..381dcff 100644 --- a/data_acq_and_processing/processing/data_processor/include/data_processor.hpp +++ b/data_acq_and_processing/processing/data_processor/include/data_processor.hpp @@ -46,12 +46,7 @@ class DataProcessor { // value exists, so the speed from the live config stays in force. [[nodiscard]] auto socket_speed(const ProcessingLiveConfig& live_config) const -> std::optional; - // Convert one processed collection into a locator filter spec, derived - // from the live config and current processor mode. - [[nodiscard]] auto build_locator_filter(const ProcessingLiveConfig& live_config) const - -> radar::locator::FilterParams; - - void publish_locator(const ipc::ResultCollection& collection, const ProcessingLiveConfig& live_config); + void publish_locator(const ipc::ResultCollection& collection); const config::RunConfig& config_; config::RunConfig processing_config_; diff --git a/data_acq_and_processing/processing/data_processor/include/processing_live_config.hpp b/data_acq_and_processing/processing/data_processor/include/processing_live_config.hpp index 446baaa..e3d2def 100644 --- a/data_acq_and_processing/processing/data_processor/include/processing_live_config.hpp +++ b/data_acq_and_processing/processing/data_processor/include/processing_live_config.hpp @@ -10,19 +10,6 @@ namespace radar::processing { -enum class LegacyGprMode { - Point, - Extended, -}; - -enum class LegacyGprReferenceMode { - // t_ref = midpoint between the first and last event centers. - FrameCenter, - // t_ref = center of the first event. Useful when motion offsets should be - // accumulated from frame start, e.g. for tagging frames by their head time. - FirstTxEvent, -}; - struct HistoryTarget { std::uint64_t collection_id = 0; std::uint64_t capture_start_ns = 0; @@ -36,15 +23,13 @@ struct ProcessingLiveConfig { float bscan_gain = 1.0F; float bscan_start_freq_mhz = 100.0F; float bscan_stop_freq_mhz = 8800.0F; - LegacyGprMode legacy_gpr_mode = LegacyGprMode::Point; - // Absent in legacy live files: use the initial run_config geometry and medium. + // Absent geometry overrides use the initial run_config geometry. std::optional gpr{}; std::vector gpr_input_positions{}; std::vector gpr_output_positions{}; float gpr_min_depth_m = 2.0F; float gpr_max_depth_m = 14.0F; float gpr_range_comp_power = 0.1F; - float gpr_comp_power = 0.2F; // BP object-detection stop level, as a fraction of the global peak; peaks // below it are not objects. float gpr_object_min_frac = 0.7F; @@ -54,29 +39,15 @@ struct ProcessingLiveConfig { std::string gpr_motion_mode = "int_minus"; float gpr_speed_m_s = 0.0F; float gpr_look_angle_deg = 0.0F; - // Motion model parameters for the legacy GPR pipeline. The direction sign - // selects which way later Tx-events appear deeper (+1) or shallower (-1). - // Intra-sweep phase correction compensates the motion that happens *inside* - // one Tx-sweep before the IFFT — it is independent of the per-pair coarse - // tau shift and can be disabled without affecting the rest of the pipeline. - // The reference mode picks the anchor used to compute dt_ref per event. float gpr_direction_sign = 1.0F; - bool gpr_apply_freq_phase_correction = true; - LegacyGprReferenceMode gpr_reference_mode = LegacyGprReferenceMode::FrameCenter; - float gpr_snr_thresh = 4.5F; - float gpr_snr_comp_max = 25.0F; float gpr_start_freq_mhz = 3000.0F; float gpr_stop_freq_mhz = 6000.0F; bool gpr_background_subtract_enabled = true; std::uint32_t gpr_background_mean_count = 10U; bool gpr_remove_sidelobe_objects_enabled = true; - // BP image is computed in the y=imaging_plane_y_m slice of the 3D grid. - // Default 0 keeps legacy 1D antenna layouts imaging in the antenna plane. - float gpr_imaging_plane_y_m = 0.0F; - // Coherent BP object visibility (window + the draw limits below) is applied in - // the processor itself, with no score threshold. Only legacy GPR thresholds, - // on a pair count. - float legacy_gpr_min_visible_pair_count = 0.0F; + // Ground plane in radar coordinates, at this height and the configured pitch. + float gpr_bp_ground_depth_m = 1.25F; + bool pair_depth_equalize = true; std::uint32_t gpr_max_detected_objects_to_draw = 0; std::uint32_t gpr_draw_top_m_objects = 0; // Cross-frame approach filter: an object is shown only once it persists as a @@ -89,7 +60,7 @@ struct ProcessingLiveConfig { float gpr_visible_z_min_m = 0.0F; float gpr_visible_z_max_m = 14.0F; // 3D SAR (sar3d processor) parameters, named after the reference notebook - // Multi_antenna_3D_clean_v1.ipynb (SAR_FRAME_RADIUS, BP_GROUND_DEPTH_M, ...). + // Multi_antenna_demonstrate_depthnorm.ipynb (SAR_FRAME_RADIUS, BP_GROUND_DEPTH_M, ...). // The band, depth gate, antenna selection, speed and look angle are the // shared gpr_* fields above. std::uint32_t sar3d_frame_radius = 0U; @@ -99,6 +70,10 @@ struct ProcessingLiveConfig { float sar3d_comp_range_power = 0.3F; bool sar3d_use_compensation = true; bool sar3d_normalize_each_frame = true; + float sar3d_x_size_m = 0.0F; + float sar3d_y_min_m = -0.5F; + float sar3d_y_max_m = 2.0F; + float sar3d_z_size_m = 0.0F; // When true, the data_processor ignores socket-supplied `vlc` updates and // keeps using `gpr_speed_m_s` from this file. Mirrored from the GUI's // "ignore socket speed" checkbox. diff --git a/data_acq_and_processing/processing/data_processor/src/data_processor.cpp b/data_acq_and_processing/processing/data_processor/src/data_processor.cpp index 14679cc..f612cb9 100644 --- a/data_acq_and_processing/processing/data_processor/src/data_processor.cpp +++ b/data_acq_and_processing/processing/data_processor/src/data_processor.cpp @@ -22,10 +22,15 @@ namespace { constexpr const char* kDefaultProcessorMode = "pass_through"; constexpr std::size_t kBscanReplayWindow = 50U; -[[nodiscard]] auto replay_history_limit(const config::RunConfig& config) -> std::size_t { +[[nodiscard]] auto processing_history_limit(const config::RunConfig& config, const ProcessingLiveConfig& live) -> std::size_t { const auto preprocessed_capacity = static_cast(std::max(1U, config.rings.preprocessed.capacity)); const auto results_capacity = static_cast(std::max(1U, config.rings.results.capacity)); - return std::min({preprocessed_capacity, results_capacity, kBscanReplayWindow}); + std::size_t required = 1U; + if (live.processor_mode == "gpr" || live.processor_mode == "sar3d") { + if (live.gpr_background_subtract_enabled) { required += live.gpr_background_mean_count; } + if (live.processor_mode == "sar3d") { required += live.sar3d_frame_radius; } + } + return std::max(required, std::min({preprocessed_capacity, results_capacity, kBscanReplayWindow})); } void publish_result_collection( @@ -68,7 +73,6 @@ DataProcessor::DataProcessor( void DataProcessor::run(const std::atomic& stop_requested) { std::vector bytes{}; std::vector preprocessed_history{}; - const std::size_t history_limit = replay_history_limit(config_); std::uint64_t last_replayed_revision = live_config_loader_.revision(); runtime::CollectionErrors errors; std::uint64_t oversized_results = 0; @@ -92,6 +96,7 @@ void DataProcessor::run(const std::atomic& stop_requested) { processing_config_.gpr = live_config.gpr.value_or(config_.gpr); } auto& processor = resolve_processor(live_config); + const auto history_limit = processing_history_limit(config_, live_config); const bool socket_speed_dirty = current_socket_speed != last_reprocessed_socket_speed; const auto was_removed = [&](const ipc::PreprocessedCollection& collection) { @@ -107,6 +112,9 @@ void DataProcessor::run(const std::atomic& stop_requested) { }; if (live_revision != last_replayed_revision || socket_speed_dirty) { std::erase_if(preprocessed_history, was_removed); + if (preprocessed_history.size() > history_limit) { + preprocessed_history.erase(preprocessed_history.begin(), preprocessed_history.end() - history_limit); + } if (!live_config.reprocess_current_result) { // Socket-fed speed updates should affect only future preprocessed collections, @@ -120,7 +128,7 @@ void DataProcessor::run(const std::atomic& stop_requested) { live_config ); publish_result_collection(replay_result, results_ring_, oversized_results); - publish_locator(replay_result, live_config); + publish_locator(replay_result); } } else if (!preprocessed_history.empty()) { const auto replay_result = process_collection( @@ -130,7 +138,7 @@ void DataProcessor::run(const std::atomic& stop_requested) { live_config ); publish_result_collection(replay_result, results_ring_, oversized_results); - publish_locator(replay_result, live_config); + publish_locator(replay_result); } last_replayed_revision = live_revision; //record the speed we just reprocessed with so an @@ -156,8 +164,8 @@ void DataProcessor::run(const std::atomic& stop_requested) { preprocessed.capture_start_ns <= preprocessed_history.back().capture_start_ns; if (!fed_back) { preprocessed_history.push_back(std::move(preprocessed)); - while (preprocessed_history.size() > history_limit) { - preprocessed_history.erase(preprocessed_history.begin()); + if (preprocessed_history.size() > history_limit) { + preprocessed_history.erase(preprocessed_history.begin(), preprocessed_history.end() - history_limit); } } const auto& current = fed_back ? preprocessed : preprocessed_history.back(); @@ -175,7 +183,7 @@ void DataProcessor::run(const std::atomic& stop_requested) { ); publish_result_collection(result_collection, results_ring_, oversized_results); if (!fed_back) { - publish_locator(result_collection, live_config); + publish_locator(result_collection); } errors.progress(); continue; @@ -237,40 +245,13 @@ auto DataProcessor::socket_speed(const ProcessingLiveConfig& live_config) const return locator_server_->latest_socket_speed(); } -auto DataProcessor::build_locator_filter(const ProcessingLiveConfig& live_config) const - -> radar::locator::FilterParams { - const std::string requested_mode = - live_config.processor_mode.empty() ? default_processor_mode_ : live_config.processor_mode; - - radar::locator::FilterParams filter{}; - if (requested_mode == "legacy_gpr") { - // Legacy GPR emits its objects unfiltered, so the socket applies the legacy - // rule here: clip to the visible window and threshold on the pair count. The - // GUI disables "draw top N" for legacy, so we skip it on the wire to match. - filter.visible_bounds = radar::locator::VisibleBounds{ - .x_min = live_config.gpr_visible_x_min_m, - .x_max = live_config.gpr_visible_x_max_m, - .z_min = live_config.gpr_visible_z_min_m, - .z_max = live_config.gpr_visible_z_max_m, - }; - filter.min_score = live_config.legacy_gpr_min_visible_pair_count; - filter.draw_limits.reset(); - } - // Coherent BP already emits the FINAL visible object set from the processor - // (window + N/M, no score threshold), so the - // socket forwards it verbatim. The default FilterParams passes everything through - // (it only drops non-finite rows), keeping the filtering logic in one place. - return filter; -} - void DataProcessor::publish_locator( - const ipc::ResultCollection& collection, - const ProcessingLiveConfig& live_config + const ipc::ResultCollection& collection ) { - if (locator_server_ == nullptr || !locator_server_->is_running()) { + if (locator_server_ == nullptr) { return; } - locator_server_->publish(collection, build_locator_filter(live_config)); + locator_server_->publish(collection); } auto create_default_processors() -> ProcessorRegistry { @@ -287,10 +268,6 @@ auto create_default_processors() -> ProcessorRegistry { auto processor = std::make_unique(); processors.emplace(processor->name(), std::move(processor)); } - { - auto processor = std::make_unique(); - processors.emplace(processor->name(), std::move(processor)); - } { auto processor = std::make_unique(); processors.emplace(processor->name(), std::move(processor)); diff --git a/data_acq_and_processing/processing/data_processor/src/processing_live_config.cpp b/data_acq_and_processing/processing/data_processor/src/processing_live_config.cpp index c6465a4..c783972 100644 --- a/data_acq_and_processing/processing/data_processor/src/processing_live_config.cpp +++ b/data_acq_and_processing/processing/data_processor/src/processing_live_config.cpp @@ -19,28 +19,6 @@ namespace { using Json = nlohmann::json; -[[nodiscard]] auto parse_legacy_gpr_mode(const std::string& value, const std::string& field_name) -> LegacyGprMode { - if (value == "point") { - return LegacyGprMode::Point; - } - if (value == "extended") { - return LegacyGprMode::Extended; - } - throw std::runtime_error(field_name + " must be one of: point, extended"); -} - -[[nodiscard]] auto parse_legacy_gpr_reference_mode( - const std::string& value, const std::string& field_name -) -> LegacyGprReferenceMode { - if (value == "frame_center") { - return LegacyGprReferenceMode::FrameCenter; - } - if (value == "first_tx_event") { - return LegacyGprReferenceMode::FirstTxEvent; - } - throw std::runtime_error(field_name + " must be one of: frame_center, first_tx_event"); -} - [[nodiscard]] auto parse_gpr_score_mode(const std::string& value, const std::string& field_name) -> std::string { if (value == "peak" || value == "combined") { return value; @@ -55,24 +33,6 @@ using Json = nlohmann::json; throw std::runtime_error(field_name + " must be one of: int_minus, int_focus"); } -void apply_legacy_gpr_algorithm_alias(ProcessingLiveConfig& config, const std::string& value) { - if (value == "backprojection") { - return; - } - if (value == "legacy_point") { - config.legacy_gpr_mode = LegacyGprMode::Point; - } else if (value == "legacy_extended") { - config.legacy_gpr_mode = LegacyGprMode::Extended; - } else { - throw std::runtime_error( - "processing.gpr_algorithm must be one of: backprojection, legacy_point, legacy_extended" - ); - } - if (config.processor_mode.empty() || config.processor_mode == "gpr") { - config.processor_mode = "legacy_gpr"; - } -} - [[nodiscard]] auto parse_f32(const Json& value, const std::string& field_name) -> float { if (!value.is_number()) { throw std::runtime_error(field_name + " must be number"); @@ -192,11 +152,8 @@ auto parse_geometry(const Json& value, const char* position_field, std::uint32_t if (const auto found = root.find("settings_revision"); found != root.end()) { config.settings_revision = parse_u64_number(*found, "processing.settings_revision"); } - if (root.contains("relative_permittivity") || root.contains("tx_geometry") || root.contains("rx_geometry")) { + if (root.contains("tx_geometry") || root.contains("rx_geometry")) { config.gpr = run_config.gpr; - if (const auto found = root.find("relative_permittivity"); found != root.end()) { - config.gpr->relative_permittivity = parse_f32(*found, "processing.relative_permittivity"); - } if (const auto found = root.find("tx_geometry"); found != root.end()) { config.gpr->tx_geometry = parse_geometry(*found, "output_pos", &config::GprTxGeometry::output_pos); } @@ -206,25 +163,16 @@ auto parse_geometry(const Json& value, const char* position_field, std::uint32_t config::validate_gpr_config(*config.gpr, run_config); } assign_optional(root, "processor_mode", config.processor_mode); + if (!config.processor_mode.empty() && config.processor_mode != "pass_through" + && config.processor_mode != "bscan" && config.processor_mode != "gpr" && config.processor_mode != "sar3d") { + throw std::runtime_error("Unsupported processing.processor_mode: " + config.processor_mode); + } assign_optional(root, "bscan_axis", config.bscan_axis); assign_optional(root, "bscan_cut_m", config.bscan_cut_m); assign_optional(root, "bscan_max_depth_m", config.bscan_max_depth_m); assign_optional(root, "bscan_gain", config.bscan_gain); assign_optional(root, "bscan_start_freq_mhz", config.bscan_start_freq_mhz); assign_optional(root, "bscan_stop_freq_mhz", config.bscan_stop_freq_mhz); - if (const auto found = root.find("legacy_gpr_mode"); found != root.end()) { - if (!found->is_string()) { - throw std::runtime_error("processing.legacy_gpr_mode must be string"); - } - config.legacy_gpr_mode = - parse_legacy_gpr_mode(found->get(), "processing.legacy_gpr_mode"); - } - if (const auto found = root.find("gpr_algorithm"); found != root.end()) { - if (!found->is_string()) { - throw std::runtime_error("processing.gpr_algorithm must be string"); - } - apply_legacy_gpr_algorithm_alias(config, found->get()); - } if (const auto found = root.find("gpr_input_positions"); found != root.end()) { config.gpr_input_positions = parse_u32_array(*found, "processing.gpr_input_positions"); } @@ -235,7 +183,6 @@ auto parse_geometry(const Json& value, const char* position_field, std::uint32_t assign_optional(root, "gpr_max_depth_m", config.gpr_max_depth_m); assign_optional(root, "gpr_range_comp_power", config.gpr_range_comp_power); assign_optional(root, "gpr_object_min_frac", config.gpr_object_min_frac); - assign_optional(root, "gpr_comp_power", config.gpr_comp_power); if (const auto found = root.find("gpr_score_mode"); found != root.end()) { if (!found->is_string()) { throw std::runtime_error("processing.gpr_score_mode must be string"); @@ -251,16 +198,6 @@ auto parse_geometry(const Json& value, const char* position_field, std::uint32_t assign_optional(root, "gpr_speed_m_s", config.gpr_speed_m_s); assign_optional(root, "gpr_look_angle_deg", config.gpr_look_angle_deg); assign_optional(root, "gpr_direction_sign", config.gpr_direction_sign); - assign_optional(root, "gpr_apply_freq_phase_correction", config.gpr_apply_freq_phase_correction); - if (const auto found = root.find("gpr_reference_mode"); found != root.end()) { - if (!found->is_string()) { - throw std::runtime_error("processing.gpr_reference_mode must be string"); - } - config.gpr_reference_mode = - parse_legacy_gpr_reference_mode(found->get(), "processing.gpr_reference_mode"); - } - assign_optional(root, "gpr_snr_thresh", config.gpr_snr_thresh); - assign_optional(root, "gpr_snr_comp_max", config.gpr_snr_comp_max); assign_optional(root, "gpr_start_freq_mhz", config.gpr_start_freq_mhz); assign_optional(root, "gpr_stop_freq_mhz", config.gpr_stop_freq_mhz); assign_optional(root, "gpr_background_subtract_enabled", config.gpr_background_subtract_enabled); @@ -268,7 +205,8 @@ auto parse_geometry(const Json& value, const char* position_field, std::uint32_t config.gpr_background_mean_count = parse_u32_number(*found, "processing.gpr_background_mean_count"); } assign_optional(root, "gpr_remove_sidelobe_objects_enabled", config.gpr_remove_sidelobe_objects_enabled); - assign_optional(root, "gpr_imaging_plane_y_m", config.gpr_imaging_plane_y_m); + assign_optional(root, "gpr_bp_ground_depth_m", config.gpr_bp_ground_depth_m); + assign_optional(root, "pair_depth_equalize", config.pair_depth_equalize); for (const auto& [key, target] : { std::pair{"gpr_visible_x_min_m", &config.gpr_visible_x_min_m}, std::pair{"gpr_visible_x_max_m", &config.gpr_visible_x_max_m}, @@ -277,7 +215,6 @@ auto parse_geometry(const Json& value, const char* position_field, std::uint32_t }) { assign_optional(root, key, *target); } - assign_optional(root, "legacy_gpr_min_visible_pair_count", config.legacy_gpr_min_visible_pair_count); if (const auto found = root.find("gpr_max_detected_objects_to_draw"); found != root.end()) { config.gpr_max_detected_objects_to_draw = parse_u32_number(*found, "processing.gpr_max_detected_objects_to_draw"); @@ -297,6 +234,10 @@ auto parse_geometry(const Json& value, const char* position_field, std::uint32_t std::pair{"sar3d_voxel_m", &config.sar3d_voxel_m}, std::pair{"sar3d_bp_ground_depth_m", &config.sar3d_bp_ground_depth_m}, std::pair{"sar3d_comp_range_power", &config.sar3d_comp_range_power}, + std::pair{"sar3d_x_size_m", &config.sar3d_x_size_m}, + std::pair{"sar3d_y_min_m", &config.sar3d_y_min_m}, + std::pair{"sar3d_y_max_m", &config.sar3d_y_max_m}, + std::pair{"sar3d_z_size_m", &config.sar3d_z_size_m}, }) { assign_optional(root, key, *target); } @@ -306,6 +247,10 @@ auto parse_geometry(const Json& value, const char* position_field, std::uint32_t }) { assign_optional(root, key, *target); } + if (config.sar3d_x_size_m < 0.0F || config.sar3d_z_size_m < 0.0F || + config.sar3d_y_max_m <= config.sar3d_y_min_m) { + throw std::runtime_error("SAR scene sizes must be non-negative and SAR_Y_RANGE_M must increase"); + } assign_optional(root, "ignore_socket_speed", config.ignore_socket_speed); assign_optional(root, "reprocess_current_result", config.reprocess_current_result); if (const auto found = root.find("history_clear_before_ns"); found != root.end()) { diff --git a/data_acq_and_processing/processing/locator/include/locator/payload_builder.hpp b/data_acq_and_processing/processing/locator/include/locator/payload_builder.hpp index ed13015..a51ccd5 100644 --- a/data_acq_and_processing/processing/locator/include/locator/payload_builder.hpp +++ b/data_acq_and_processing/processing/locator/include/locator/payload_builder.hpp @@ -1,7 +1,6 @@ #pragma once #include -#include #include #include @@ -10,47 +9,9 @@ namespace radar::locator { -// Inclusive bounds used to clip locator observations to a visible window. -struct VisibleBounds { - float x_min = 0.0F; - float x_max = 0.0F; - float z_min = 0.0F; - float z_max = 0.0F; -}; - -// Limits mirroring the GUI semantics: if the number of detected objects exceeds -// `max_detected_objects`, the result is intentionally empty (matches Python -// reference). Otherwise, at most `draw_top_objects` rows are emitted. -struct DrawLimits { - std::uint32_t max_detected_objects = 0; - std::uint32_t draw_top_objects = 0; -}; - -// Filter parameters for the payload builder. None of these fields couple to GUI -// state; they are resolved by the data_processor from the live processing config -// and passed in explicitly so the builder remains a pure function. -struct FilterParams { - float min_score = 0.0F; - std::optional visible_bounds{}; - std::optional draw_limits{}; -}; - -// Extract observations from one result collection. -// -// Mirrors `python_app/orchestration/gpr_locator.py::locator_observations_from_collection`: -// * Looks for a `gpr_points` TableF32 payload first, then falls back to -// `gpr_region_centers`. Both layouts encode `[x_m, z_m, score, ...]` rows. -// * Drops rows with non-finite coordinates and those below `min_score`. -// * If `visible_bounds` are provided, drops rows outside the inclusive window. -// * If `draw_limits` are provided and the surviving count exceeds -// `max_detected_objects`, returns an empty vector (the GUI's "too many to -// trust" heuristic). Otherwise, keeps the first `draw_top_objects` rows. -// -// The function never throws; malformed payloads degrade to an empty result. -[[nodiscard]] auto observations_from_collection( - const ipc::ResultCollection& collection, - const FilterParams& filter -) -> std::vector; +// Extract the final coherent-GPR objects, dropping malformed or non-finite rows. +[[nodiscard]] auto observations_from_collection(const ipc::ResultCollection& collection) + -> std::vector; // Serialize a list of observations into the JSON payload format expected by // locator clients: `{"ver": , "tim": "HH:MM:SS.mmm", "sts": 1, "obs": [...]}`. diff --git a/data_acq_and_processing/processing/locator/include/locator/tcp_server.hpp b/data_acq_and_processing/processing/locator/include/locator/tcp_server.hpp index 1a8f1f0..80af34c 100644 --- a/data_acq_and_processing/processing/locator/include/locator/tcp_server.hpp +++ b/data_acq_and_processing/processing/locator/include/locator/tcp_server.hpp @@ -138,7 +138,7 @@ class TcpServer { // Build and broadcast one locator packet derived from a GPR result // collection. Never blocks for I/O. Filter parameters are supplied by the // caller so the builder remains state-free. - void publish(const ipc::ResultCollection& collection, const FilterParams& filter); + void publish(const ipc::ResultCollection& collection); // Most recent `vlc` value received from any connected client, or nullopt // if no client has ever submitted one since startup. Lock-free. diff --git a/data_acq_and_processing/processing/locator/src/payload_builder.cpp b/data_acq_and_processing/processing/locator/src/payload_builder.cpp index 4611a99..fe4969c 100644 --- a/data_acq_and_processing/processing/locator/src/payload_builder.cpp +++ b/data_acq_and_processing/processing/locator/src/payload_builder.cpp @@ -19,14 +19,12 @@ namespace { using Json = nlohmann::json; constexpr const char* kGprPointsName = "gpr_points"; -constexpr const char* kGprRegionCentersName = "gpr_region_centers"; constexpr std::uint32_t kMinTableColumns = 3U; // [x_m, z_m, score, ...] // Locate the first TableF32 payload matching one of the GPR object names. // Returns nullptr if no usable payload is present in the collection. [[nodiscard]] auto find_object_table(const ipc::ResultCollection& collection) -> const ipc::ResultPayload* { - const ipc::ResultPayload* fallback = nullptr; for (const auto& payload : collection.collection_payloads) { if (payload.kind != ipc::ResultKind::TableF32) { continue; @@ -37,11 +35,8 @@ constexpr std::uint32_t kMinTableColumns = 3U; // [x_m, z_m, score, ...] if (payload.processing_name == kGprPointsName) { return &payload; } - if (fallback == nullptr && payload.processing_name == kGprRegionCentersName) { - fallback = &payload; - } } - return fallback; + return nullptr; } // Quantise to one decimal place (0.1 m = decimetres). Equivalent to Python's @@ -51,29 +46,6 @@ constexpr std::uint32_t kMinTableColumns = 3U; // [x_m, z_m, score, ...] return std::round(value * 10.0) / 10.0; } -[[nodiscard]] auto passes_basic_filter( - float x_m, - float z_m, - float score, - float min_score -) -> bool { - if (!std::isfinite(x_m) || !std::isfinite(z_m) || !std::isfinite(score)) { - return false; - } - return score >= min_score; -} - -[[nodiscard]] auto passes_visible_bounds( - float x_m, - float z_m, - const VisibleBounds& bounds -) -> bool { - return x_m >= bounds.x_min - && x_m <= bounds.x_max - && z_m >= bounds.z_min - && z_m <= bounds.z_max; -} - // Render the current wall-clock time as "HH:MM:SS.mmm". Uses localtime to match // the Python reference behaviour (which calls datetime.now() with no tz info). [[nodiscard]] auto format_timestamp_now() -> std::string { @@ -107,66 +79,20 @@ void append_u32_little_endian(std::vector& buffer, std::uint32_t v } // namespace -auto observations_from_collection( - const ipc::ResultCollection& collection, - const FilterParams& filter -) -> std::vector { +auto observations_from_collection(const ipc::ResultCollection& collection) -> std::vector { const auto* payload = find_object_table(collection); if (payload == nullptr) { return {}; } - - const std::uint32_t columns = payload->table_columns; - if (columns == 0U) { - return {}; - } - const std::size_t row_count = payload->table_values.size() / columns; - if (row_count == 0U) { - return {}; - } - - std::vector visible; - visible.reserve(row_count); - for (std::size_t row = 0; row < row_count; ++row) { - const std::size_t base = row * static_cast(columns); - const float x_m = payload->table_values[base + 0U]; - const float z_m = payload->table_values[base + 1U]; - const float score = payload->table_values[base + 2U]; - - if (!passes_basic_filter(x_m, z_m, score, filter.min_score)) { - continue; - } - if (filter.visible_bounds.has_value() - && !passes_visible_bounds(x_m, z_m, *filter.visible_bounds)) { - continue; - } - visible.push_back({.crs = x_m, .dst = z_m}); - } - - if (!filter.draw_limits.has_value()) { - std::vector quantised; - quantised.reserve(visible.size()); - for (const auto& observation : visible) { - quantised.push_back({ - .crs = quantise_to_decimetres(observation.crs), - .dst = quantise_to_decimetres(observation.dst), - }); - } - return quantised; - } - - const auto& limits = *filter.draw_limits; - if (visible.size() > limits.max_detected_objects) { - return {}; - } - const std::size_t kept = std::min(visible.size(), limits.draw_top_objects); + const std::size_t columns = payload->table_columns; std::vector result; - result.reserve(kept); - for (std::size_t index = 0; index < kept; ++index) { - result.push_back({ - .crs = quantise_to_decimetres(visible[index].crs), - .dst = quantise_to_decimetres(visible[index].dst), - }); + for (std::size_t base = 0; base + columns <= payload->table_values.size(); base += columns) { + const auto x = payload->table_values[base]; + const auto z = payload->table_values[base + 1]; + const auto score = payload->table_values[base + 2]; + if (std::isfinite(x) && std::isfinite(z) && std::isfinite(score)) { + result.push_back({.crs = quantise_to_decimetres(x), .dst = quantise_to_decimetres(z)}); + } } return result; } diff --git a/data_acq_and_processing/processing/locator/src/tcp_server.cpp b/data_acq_and_processing/processing/locator/src/tcp_server.cpp index a65b070..e9c94bd 100644 --- a/data_acq_and_processing/processing/locator/src/tcp_server.cpp +++ b/data_acq_and_processing/processing/locator/src/tcp_server.cpp @@ -461,12 +461,12 @@ auto TcpServer::is_running() const -> bool { return running_.load(std::memory_order_acquire); } -void TcpServer::publish(const ipc::ResultCollection& collection, const FilterParams& filter) { +void TcpServer::publish(const ipc::ResultCollection& collection) { if (!running_.load(std::memory_order_acquire)) { return; } - const auto observations = observations_from_collection(collection, filter); + const auto observations = observations_from_collection(collection); const auto payload_json = build_payload_json(observations, config_.protocol_version); auto packet = encode_packet(payload_json, config_.device_id); cache_latest_packet(packet); diff --git a/data_acq_and_processing/processing/processors/include/gpr_processor.hpp b/data_acq_and_processing/processing/processors/include/gpr_processor.hpp index 18567b6..e5299ba 100644 --- a/data_acq_and_processing/processing/processors/include/gpr_processor.hpp +++ b/data_acq_and_processing/processing/processors/include/gpr_processor.hpp @@ -21,15 +21,4 @@ class GprProcessor final : public ProcessorInterface { ObjectApproachFilter approach_filter_{}; }; -class LegacyGprProcessor final : public ProcessorInterface { - public: - [[nodiscard]] auto name() const -> std::string override; - [[nodiscard]] auto process_collection( - const config::RunConfig& run_config, - const ipc::PreprocessedCollection& collection, - std::span previous_collections, - const ProcessingLiveConfig& live_config - ) -> ipc::ResultCollection override; -}; - } // namespace radar::processing diff --git a/data_acq_and_processing/processing/processors/include/sar3d_processor.hpp b/data_acq_and_processing/processing/processors/include/sar3d_processor.hpp index a203d31..48c3b6d 100644 --- a/data_acq_and_processing/processing/processors/include/sar3d_processor.hpp +++ b/data_acq_and_processing/processing/processors/include/sar3d_processor.hpp @@ -8,7 +8,7 @@ namespace radar::processing { -// 3D SAR of the reference notebook Multi_antenna_3D_clean_v1.ipynb on an NVIDIA +// 3D SAR of the reference notebook Multi_antenna_demonstrate_depthnorm.ipynb on an NVIDIA // GPU. Publishes max-projections of the normalized volume as ImageF32 payloads: // sar3d_front (X-Y), sar3d_side (Z-Y), sar3d_top (X-Z) and sar3d_oblique (the // notebook's 3D scene view, axes in the plane of sar3d_oblique_basis), plus diff --git a/data_acq_and_processing/processing/processors/src/backprojection.hpp b/data_acq_and_processing/processing/processors/src/backprojection.hpp new file mode 100644 index 0000000..77882fe --- /dev/null +++ b/data_acq_and_processing/processing/processors/src/backprojection.hpp @@ -0,0 +1,231 @@ +#pragma once + +#include "shared_types.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace radar::processing::bp { + +constexpr double kWaveSpeed = 3e8; // v in both reference algorithms +using Background = std::unordered_map>, ipc::ComboKeyHash>; + +// Mean of up to count previous collections, separately for each physical pair. +inline auto background_mean(std::span previous, bool enabled, + std::size_t count) -> Background { + Background result; + if (!enabled || count == 0U) { return result; } + std::unordered_map counts; + for (const auto& collection : previous.last(std::min(count, previous.size()))) { + for (const auto& trace : collection.traces) { + const auto key = trace.combo; + auto& sum = result[key]; + if (sum.empty()) { sum.resize(trace.s21.size()); } + if (sum.size() != trace.s21.size()) { continue; } + for (std::size_t i = 0; i < sum.size(); ++i) { + sum[i] += std::complex(trace.s21[i].re, trace.s21[i].im); + } + ++counts[key]; + } + } + for (auto& [key, samples] : result) { + for (auto& sample : samples) { sample /= static_cast(counts.at(key)); } + } + return result; +} + +inline auto subtract_background(const ipc::SweepTraceBlock& trace, const Background& background) + -> std::vector> { + std::vector> samples; + samples.reserve(trace.s21.size()); + const auto found = background.find(trace.combo); + const bool subtract = found != background.end() && found->second.size() == trace.s21.size(); + for (std::size_t i = 0; i < trace.s21.size(); ++i) { + samples.emplace_back(std::complex(trace.s21[i].re, trace.s21[i].im) - + (subtract ? found->second[i] : std::complex{})); + } + return samples; +} + +// When the samples of a trace were measured, relative to the acquisition centre +// of its collection (compute_pair_timestamps_for_snapshot and +// compute_frequency_sample_times). Traces sharing a capture window come from +// one sweep: at every frequency point its stimulus stages - the outputs, which +// are numbered in device stage order - take turns while all inputs receive. +struct TraceClock { + double start_s = 0.0; // sweep start + double event_s = 0.0; // sweep duration / (stages * points) + std::size_t stage = 0U; + std::size_t stages = 1U; + + [[nodiscard]] auto sample_time_s(std::size_t point) const -> double { + return start_s + (static_cast(stages) * point + stage) * event_s; + } +}; + +struct CollectionTiming { + std::uint64_t center_ns = 0U; // centre of the collection's capture window + std::vector clocks; // per collection trace +}; + +// Timing from the per-trace capture windows (snapshot_acquisition_timing); every +// trace takes part, so the stages do not depend on the INPUT_IDX/OUTPUT_IDX choice. +[[nodiscard]] inline auto collection_timing(const ipc::PreprocessedCollection& collection) -> CollectionTiming { + if (collection.traces.empty()) { + throw std::runtime_error("collection has no traces"); + } + std::map, std::vector> sweep_outputs{}; + std::uint64_t first = std::numeric_limits::max(); + std::uint64_t last = 0U; + for (const auto& trace : collection.traces) { + if (trace.capture_end_ns <= trace.capture_start_ns) { + throw std::runtime_error("collection has traces without a capture window"); + } + if (trace.frequency_hz.empty()) { + throw std::runtime_error("collection has a trace without frequency points"); + } + sweep_outputs[{trace.capture_start_ns, trace.capture_end_ns}].push_back(trace.combo.output_pos); + first = std::min(first, trace.capture_start_ns); + last = std::max(last, trace.capture_end_ns); + } + for (auto& [window, outputs] : sweep_outputs) { + std::sort(outputs.begin(), outputs.end()); + outputs.erase(std::unique(outputs.begin(), outputs.end()), outputs.end()); + } + + CollectionTiming timing{}; + const auto duration_ns = last - first; + timing.center_ns = first + duration_ns / 2U; + for (const auto& trace : collection.traces) { + const auto& stages = sweep_outputs.at({trace.capture_start_ns, trace.capture_end_ns}); + const auto stage = std::find(stages.begin(), stages.end(), trace.combo.output_pos) - stages.begin(); + TraceClock clock{}; + clock.start_s = 1e-9 * (static_cast(trace.capture_start_ns - first) - + 0.5 * static_cast(duration_ns)); + clock.stage = stage; + clock.stages = stages.size(); + clock.event_s = 1e-9 * static_cast(trace.capture_end_ns - trace.capture_start_ns) / + static_cast(clock.stages) / trace.frequency_hz.size(); + timing.clocks.push_back(clock); + } + return timing; +} + +// apply_intra_sweep_phase_correction: one phase model includes the pair offset +// and the intra-sweep motion. Focus mode removes its constant and linear terms. +inline void apply_motion_phase(std::vector>& samples, const std::vector& frequency, + const TraceClock& clock, double speed, double pitch, double direction, bool focus, + double f_start, double f_stop) { + if (speed == 0.0) { return; } + std::vector points; + std::vector phase; + const double factor = 4.0 * std::numbers::pi * direction * speed * std::cos(pitch) / kWaveSpeed; + double sum_f = 0.0, sum_phase = 0.0; + for (std::size_t i = 0; i < frequency.size(); ++i) { + if (frequency[i] < f_start || frequency[i] > f_stop) { continue; } + points.push_back(i); + phase.push_back(factor * frequency[i] * clock.sample_time_s(i)); + sum_f += frequency[i]; + sum_phase += phase.back(); + } + if (points.empty()) { return; } + if (focus) { + const double mean_f = sum_f / points.size(), mean_phase = sum_phase / points.size(); + double cross = 0.0, square = 0.0; + for (std::size_t i = 0; i < points.size(); ++i) { + const double f = frequency[points[i]] - mean_f; + cross += f * phase[i]; square += f * f; + } + const double slope = square > 0.0 ? cross / square : 0.0; + for (std::size_t i = 0; i < points.size(); ++i) { + phase[i] -= mean_phase + slope * (frequency[points[i]] - mean_f); + } + } + for (std::size_t i = 0; i < points.size(); ++i) { samples[points[i]] *= std::polar(1.0, -phase[i]); } +} + +// pair_depth_eq_z_max: cover the far scene corner, aperture and half a median window. +inline auto equalize_depth_limit(double max_depth, double x_min, double x_max, double antenna_x_max) -> double { + const double x_far = std::max(std::abs(x_min), std::abs(x_max)) + antenna_x_max; + return std::sqrt(max_depth * max_depth + x_far * x_far + 2.5 * 2.5) + 1.5 + 0.75; +} + +// bp_plane_y, Horns_motion_4x4_ground.py: the ground plane in radar axes. +inline auto ground_y(double z, double height, double pitch) -> double { + if (std::abs(std::cos(pitch)) < 1e-6) { throw std::runtime_error("Ground BP is undefined at a look angle of 90 degrees"); } + return (-height + z * std::sin(pitch)) / std::cos(pitch); +} + +// equalize_pairs_by_depth, Multi_antenna_demonstrate_depthnorm.ipynb, cell 1. +// All time axes start at zero; nearest median filtering uses the upper median +// for even windows, as scipy.ndimage.median_filter does. +template +void equalize_pairs_by_depth(const std::vector>>& pairs, double dz, double z_max) { + if (pairs.size() < 2U) { return; } + std::size_t length = pairs.front().size(); + for (const auto pair : pairs) { length = std::min(length, pair.size()); } + const double step = std::max(1.0, std::nearbyint(0.10 / dz)); + if (step > length) { return; } + const auto d = static_cast(step); + const auto n = static_cast(std::max(3.0, std::nearbyint(1.5 / (dz * d)))); + const auto nc = static_cast(std::min(static_cast(length), std::ceil(z_max / dz))) / d; + if (nc < 2U) { return; } + std::vector> levels(pairs.size(), std::vector(nc)); + std::vector common(nc, 0.0), envelope(nc), window(n); + for (std::size_t a = 0; a < pairs.size(); ++a) { + for (std::size_t i = 0; i < nc; ++i) { + double sum = 0.0; + for (std::size_t k = 0; k < d; ++k) { sum += std::abs(std::complex(pairs[a][i * d + k])); } + envelope[i] = sum / static_cast(d); + } + for (std::size_t i = 0; i < nc; ++i) { + for (std::size_t k = 0; k < n; ++k) { + const auto index = std::clamp(static_cast(i) + + static_cast(k) - static_cast(n / 2), 0, nc - 1); + window[k] = envelope[index]; + } + std::nth_element(window.begin(), window.begin() + n / 2, window.end()); + levels[a][i] = std::log(std::max(window[n / 2], 1e-30)); + common[i] += levels[a][i] / static_cast(pairs.size()); + } + } + for (std::size_t a = 0; a < pairs.size(); ++a) { + for (std::size_t i = 0; i < nc; ++i) { levels[a][i] = std::exp(common[i] - levels[a][i]); } + for (std::size_t i = 0; i < nc * d; ++i) { + const double pos = (static_cast(i) - (static_cast(d) - 1.0) / 2.0) / d; + const auto j = static_cast(std::clamp(std::floor(pos), 0.0, static_cast(nc - 2))); + const double w = std::clamp(pos - j, 0.0, 1.0); + pairs[a][i] *= static_cast(levels[a][j] * (1.0 - w) + levels[a][j + 1] * w); + } + } +} + +// Scalar median normalization follows the depth equalizer in both references. +template +void normalize_pair(std::span> samples, double dz, double min_depth, double max_depth) { + std::vector magnitudes; + for (std::size_t i = 0; i < samples.size(); ++i) { + if (i * dz >= min_depth && i * dz <= max_depth) { magnitudes.push_back(std::abs(std::complex(samples[i]))); } + } + if (magnitudes.empty()) { + for (const auto sample : samples) { magnitudes.push_back(std::abs(std::complex(sample))); } + } + if (magnitudes.empty()) { return; } + const auto middle = magnitudes.begin() + magnitudes.size() / 2; + std::nth_element(magnitudes.begin(), middle, magnitudes.end()); + double scale = *middle; + if (magnitudes.size() % 2 == 0) { scale = 0.5 * (scale + *std::max_element(magnitudes.begin(), middle)); } + if (!std::isfinite(scale) || scale <= 1e-15) { scale = 1.0; } + for (auto& sample : samples) { sample /= static_cast(scale + 1e-15); } +} + +} // namespace radar::processing::bp diff --git a/data_acq_and_processing/processing/processors/src/gpr_backprojection_processor.ipp b/data_acq_and_processing/processing/processors/src/gpr_backprojection_processor.ipp index 034a537..7501820 100644 --- a/data_acq_and_processing/processing/processors/src/gpr_backprojection_processor.ipp +++ b/data_acq_and_processing/processing/processors/src/gpr_backprojection_processor.ipp @@ -1,5 +1,4 @@ constexpr double kPi = 3.14159265358979323846; -constexpr double kSpeedOfLightMetersPerSec = 299'792'458.0; constexpr double kXMarginM = 2.0; constexpr double kGridZMinM = 0.1; constexpr std::size_t kGridWidth = 300U; @@ -7,13 +6,8 @@ constexpr std::size_t kGridHeight = 300U; constexpr std::size_t kAscanOversample = 8U; constexpr double kRangeWeightMax = 5.0; -constexpr double kAngleWeightMax = 2.0; -constexpr double kTotalWeightMax = 8.0; constexpr double kCompensationReferenceDepthM = 5.0; -constexpr bool kPairNormalize = true; -constexpr double kPairNormPercentile = 50.0; -constexpr double kPairNormEps = 1e-15; constexpr double kSmoothSigma = 1.5; // Gaussian-kernel half-width in sigmas, matching scipy.ndimage.gaussian_filter's @@ -55,15 +49,8 @@ using PairKey = std::uint64_t; // IntMinus — full intra-frequency correction to the interleaved-sweep centre. // IntFocus — remove the constant + linear-in-frequency part of phi(f), leaving // the focusing residual without a net Z shift. -enum class MotionMode { IntMinus, IntFocus }; - -[[nodiscard]] auto parse_motion_mode(const std::string& value) -> MotionMode { - return value == "int_focus" ? MotionMode::IntFocus : MotionMode::IntMinus; -} - struct GeometrySelection { - // Per-local-index Tx/Rx antenna coordinates in metres. y_/z_ default to 0 - // for legacy configs so the imaging plane coincides with the antennas. + // Per-local-index Tx/Rx antenna coordinates in metres. std::vector input_positions{}; std::vector output_positions{}; std::vector x_tx{}; @@ -76,11 +63,6 @@ struct GeometrySelection { std::unordered_map output_local_by_pos{}; }; -struct BackgroundAccumulator { - std::vector> sum{}; - std::size_t count = 0U; -}; - struct SelectedTrace { ipc::ComboKey combo{}; std::uint32_t tx_local_index = 0U; @@ -268,20 +250,6 @@ void fft_inplace(std::vector>& values, bool inverse) { return axis; } -[[nodiscard]] auto clamp_index(std::ptrdiff_t value, std::size_t limit) -> std::size_t { - if (limit == 0U) { - return 0U; - } - if (value < 0) { - return 0U; - } - const auto max_index = static_cast(limit - 1U); - if (value > max_index) { - return limit - 1U; - } - return static_cast(value); -} - // scipy 'reflect' boundary (half-sample symmetric): the signal is mirrored about // the outer edge of the first/last sample, so index -1 maps to 0, -2 to 1, n to // n-1, and so on. Matches scipy.ndimage.gaussian_filter's default mode. @@ -314,7 +282,7 @@ void normalize_in_place(std::vector& values) { return; } for (auto& value : values) { - value /= maximum; + value /= maximum + 1e-12; } } @@ -545,67 +513,10 @@ void validate_collection_trace_order( } } -[[nodiscard]] auto build_background_mean( - std::span previous_collections, - const GeometrySelection& selection, - const ProcessingLiveConfig& live_config -) -> std::unordered_map>> { - std::unordered_map>> result{}; - if (!live_config.gpr_background_subtract_enabled || live_config.gpr_background_mean_count == 0U) { - return result; - } - - const std::size_t mean_count = static_cast(live_config.gpr_background_mean_count); - const std::size_t start_index = - previous_collections.size() > mean_count ? previous_collections.size() - mean_count : 0U; - - std::unordered_map accumulators{}; - for (std::size_t collection_index = start_index; collection_index < previous_collections.size(); ++collection_index) { - const auto& collection = previous_collections[collection_index]; - for (const auto& trace : collection.traces) { - const auto output_it = selection.output_local_by_pos.find(trace.combo.output_pos); - const auto input_it = selection.input_local_by_pos.find(trace.combo.input_pos); - if (output_it == selection.output_local_by_pos.end() || input_it == selection.input_local_by_pos.end()) { - continue; - } - - const auto key = make_pair_key(output_it->second, input_it->second); - auto& accumulator = accumulators[key]; - if (accumulator.sum.empty()) { - accumulator.sum.assign(trace.s21.size(), std::complex(0.0, 0.0)); - } - if (accumulator.sum.size() != trace.s21.size()) { - continue; - } - - for (std::size_t sample_index = 0U; sample_index < trace.s21.size(); ++sample_index) { - const auto& sample = trace.s21[sample_index]; - accumulator.sum[sample_index] += std::complex(sample.re, sample.im); - } - accumulator.count += 1U; - } - } - - for (auto& [key, accumulator] : accumulators) { - if (accumulator.count == 0U) { - continue; - } - - auto& mean_trace = result[key]; - mean_trace = std::move(accumulator.sum); - const double inverse_count = 1.0 / static_cast(accumulator.count); - for (auto& sample : mean_trace) { - sample *= inverse_count; - } - } - - return result; -} - [[nodiscard]] auto collect_selected_traces( const ipc::PreprocessedCollection& collection, const GeometrySelection& selection, - const std::unordered_map>>& background_mean + const bp::Background& background_mean ) -> std::vector { std::vector traces{}; traces.reserve(collection.traces.size()); @@ -637,15 +548,7 @@ void validate_collection_trace_order( throw std::runtime_error("GPR requires unique selected combos; duplicate combo: " + combo_to_string(trace.combo)); } - const auto background_it = background_mean.find(key); - selected.s21.reserve(trace.s21.size()); - for (std::size_t sample_index = 0U; sample_index < trace.s21.size(); ++sample_index) { - std::complex sample(trace.s21[sample_index].re, trace.s21[sample_index].im); - if (background_it != background_mean.end() && background_it->second.size() == trace.s21.size()) { - sample -= background_it->second[sample_index]; - } - selected.s21.push_back(sample); - } + selected.s21 = bp::subtract_background(trace, background_mean); traces.push_back(std::move(selected)); } @@ -653,112 +556,6 @@ void validate_collection_trace_order( return traces; } -// Frequency-domain speed correction applied to S21 *before* the IFFT. The measurement is interleaved Tx-by-frequency -// — Tx0(f0), Tx1(f0), Tx0(f1), ... — so each Tx occupies one interleave slot and -// its frequency k is sampled at global step (slot_count*k + slot). Each in-band -// frequency's phase is rotated back to the centre of the (full) interleaved sweep. -// -// Unlike the offline Python script, every physical input is taken from the live -// system: `speed_mps` is the socket-fed velocity, `look_angle_deg`/`direction_sign` -// are config, and `tx_sweep_time_s` is derived per-collection from the acquisition -// timestamps (capture_span / Tx count). A zero speed or non-positive sweep time -// makes this a no-op. -void apply_intra_sweep_motion_correction( - SelectedTrace& trace, - std::size_t slot, - std::size_t slot_count, - double start_hz, - double stop_hz, - double velocity_mps, - double speed_mps, - double look_angle_deg, - double direction_sign, - double tx_sweep_time_s, - MotionMode motion_mode -) { - const std::size_t n_full = trace.frequency_hz.size(); - if (n_full < 2U || trace.s21.size() != n_full || slot_count == 0U - || !(velocity_mps > 0.0) || !(tx_sweep_time_s > 0.0)) { - return; - } - - const double low_hz = std::min(start_hz, stop_hz); - const double high_hz = std::max(start_hz, stop_hz); - - // Absolute sample time of full-sweep frequency index k for this Tx slot. - const double dt_base = tx_sweep_time_s / static_cast(n_full - 1U); - const auto sample_time_s = [&](std::size_t index) { - return ((static_cast(slot_count) * static_cast(index)) + static_cast(slot)) * dt_base; - }; - // t_center is the midpoint of the *full* interleaved sweep for this slot, - // i.e. between the first and last full-array frequencies (not the cut band). - const double t_center_s = 0.5 * (sample_time_s(0U) + sample_time_s(n_full - 1U)); - - const double theta_rad = (look_angle_deg * kPi) / 180.0; - const double motion_factor = direction_sign * speed_mps * std::cos(theta_rad); - - std::vector band_indices{}; - std::vector band_frequency_hz{}; - std::vector phi{}; - band_indices.reserve(n_full); - band_frequency_hz.reserve(n_full); - phi.reserve(n_full); - - bool any_nonzero = false; - for (std::size_t index = 0U; index < n_full; ++index) { - const double frequency = trace.frequency_hz[index]; - if (frequency < low_hz || frequency > high_hz) { - continue; - } - const double dt_intra_s = sample_time_s(index) - t_center_s; - const double delta_range_m = motion_factor * dt_intra_s; - const double delta_path_m = 2.0 * delta_range_m; - const double dtau_s = delta_path_m / velocity_mps; - if (dtau_s != 0.0) { - any_nonzero = true; - } - band_indices.push_back(index); - band_frequency_hz.push_back(frequency); - phi.push_back(2.0 * kPi * frequency * dtau_s); - } - - if (!any_nonzero || band_indices.empty()) { - return; - } - - // int_focus removes the constant + linear-in-frequency component of phi(f). - // With f_rel = f - mean(f) the design columns [1, f_rel] are orthogonal, so - // the least-squares fit is intercept = mean(phi), slope = /. - if (motion_mode == MotionMode::IntFocus) { - const auto count = static_cast(band_frequency_hz.size()); - double mean_frequency = 0.0; - double mean_phi = 0.0; - for (std::size_t i = 0U; i < band_frequency_hz.size(); ++i) { - mean_frequency += band_frequency_hz[i]; - mean_phi += phi[i]; - } - mean_frequency /= count; - mean_phi /= count; - - double cross = 0.0; - double f_rel_sq = 0.0; - for (std::size_t i = 0U; i < band_frequency_hz.size(); ++i) { - const double f_rel = band_frequency_hz[i] - mean_frequency; - cross += f_rel * phi[i]; - f_rel_sq += f_rel * f_rel; - } - const double slope = (f_rel_sq > 0.0) ? (cross / f_rel_sq) : 0.0; - for (std::size_t i = 0U; i < phi.size(); ++i) { - const double f_rel = band_frequency_hz[i] - mean_frequency; - phi[i] -= mean_phi + (slope * f_rel); - } - } - - for (std::size_t i = 0U; i < band_indices.size(); ++i) { - trace.s21[band_indices[i]] *= std::polar(1.0, -phi[i]); - } -} - [[nodiscard]] auto compute_ascan( const SelectedTrace& trace, double start_hz, @@ -804,7 +601,7 @@ void apply_intra_sweep_motion_correction( return result; } - const auto start_bin = static_cast(std::llround(frequency_hz.front() / df_hz)); + const auto start_bin = static_cast(std::nearbyint(frequency_hz.front() / df_hz)); if (start_bin < 0) { return result; } @@ -842,50 +639,6 @@ void apply_intra_sweep_motion_correction( return result; } -void normalize_pair_ascans( - std::unordered_map& ascans_by_pair, - double velocity_mps, - double min_depth_m, - double max_depth_m -) { - if (!kPairNormalize) { - return; - } - - for (auto& item : ascans_by_pair) { - auto& ascan = item.second; - if (ascan.samples.empty() || ascan.time_s.size() != ascan.samples.size()) { - continue; - } - - std::vector amplitudes{}; - amplitudes.reserve(ascan.samples.size()); - for (std::size_t index = 0U; index < ascan.samples.size(); ++index) { - const double depth_m = 0.5 * velocity_mps * ascan.time_s[index]; - if (depth_m < min_depth_m || depth_m > max_depth_m) { - continue; - } - amplitudes.push_back(std::abs(ascan.samples[index])); - } - if (amplitudes.empty()) { - amplitudes.reserve(ascan.samples.size()); - for (const auto& sample : ascan.samples) { - amplitudes.push_back(std::abs(sample)); - } - } - - double scale = percentile_copy(std::move(amplitudes), kPairNormPercentile); - if (!std::isfinite(scale) || !(scale > kPairNormEps)) { - scale = 1.0; - } - - const double denominator = scale + kPairNormEps; - for (auto& sample : ascan.samples) { - sample /= denominator; - } - } -} - [[nodiscard]] auto distance_3d(double dx, double dy, double dz) -> double { return std::sqrt((dx * dx) + (dy * dy) + (dz * dz)); } @@ -894,7 +647,7 @@ void normalize_pair_ascans( const GeometrySelection& selection, double max_depth_m, double min_z_m, - double imaging_plane_y_m + double ground_depth_m, double pitch_rad ) -> GridDefinition { GridDefinition grid{}; if (selection.x_tx.empty() || selection.x_rx.empty() || !(max_depth_m > min_z_m)) { @@ -922,14 +675,14 @@ void normalize_pair_ascans( for (std::size_t tx_index = 0U; tx_index < selection.x_tx.size(); ++tx_index) { grid.tx_distance_grids[tx_index][cell_index] = distance_3d( x_value - selection.x_tx[tx_index], - imaging_plane_y_m - selection.y_tx[tx_index], + bp::ground_y(z_value, ground_depth_m, pitch_rad) - selection.y_tx[tx_index], z_value - selection.z_tx[tx_index] ); } for (std::size_t rx_index = 0U; rx_index < selection.x_rx.size(); ++rx_index) { grid.rx_distance_grids[rx_index][cell_index] = distance_3d( x_value - selection.x_rx[rx_index], - imaging_plane_y_m - selection.y_rx[rx_index], + bp::ground_y(z_value, ground_depth_m, pitch_rad) - selection.y_rx[rx_index], z_value - selection.z_rx[rx_index] ); } @@ -946,10 +699,10 @@ void normalize_pair_ascans( const GeometrySelection& selection, double max_depth_m, double min_z_m, - double imaging_plane_y_m + double ground_depth_m, double pitch_rad ) -> std::vector { - std::vector signature; - signature.reserve((selection.x_tx.size() + selection.x_rx.size()) * 3U + 3U); + std::vector signature{static_cast(selection.x_tx.size()), static_cast(selection.x_rx.size())}; + signature.reserve((selection.x_tx.size() + selection.x_rx.size()) * 3U + 6U); const auto append = [&](const std::vector& axis) { signature.insert(signature.end(), axis.begin(), axis.end()); }; @@ -961,7 +714,8 @@ void normalize_pair_ascans( append(selection.z_rx); signature.push_back(max_depth_m); signature.push_back(min_z_m); - signature.push_back(imaging_plane_y_m); + signature.push_back(ground_depth_m); + signature.push_back(pitch_rad); return signature; } @@ -983,67 +737,18 @@ void normalize_pair_ascans( return (ascan.samples[index] * (1.0 - frac)) + (ascan.samples[index + 1U] * frac); } -[[nodiscard]] auto attenuation_components( - double r_tx, - double r_rx, - double dz_tx, - double dz_rx -) -> std::pair { - // Antennas boresight along +Z, so cos(theta) = (z_pixel - z_antenna) / R. - const double geo = 1.0 / ((r_tx * r_rx) + 1e-12); - const double angle = - std::pow(dz_tx / (r_tx + 1e-12), 2.0) * - std::pow(dz_rx / (r_rx + 1e-12), 2.0); - return {geo + 1e-30, angle + 1e-30}; -} - -[[nodiscard]] auto attenuation_components_at_ref_depth( - std::uint32_t tx_index, - std::uint32_t rx_index, - const GeometrySelection& selection, - double imaging_plane_y_m -) -> std::pair { - const double x_center = 0.5 * (selection.x_tx[tx_index] + selection.x_rx[rx_index]); - const double dz_tx = kCompensationReferenceDepthM - selection.z_tx[tx_index]; - const double dz_rx = kCompensationReferenceDepthM - selection.z_rx[rx_index]; - const double r_tx = distance_3d( - x_center - selection.x_tx[tx_index], - imaging_plane_y_m - selection.y_tx[tx_index], - dz_tx - ); - const double r_rx = distance_3d( - x_center - selection.x_rx[rx_index], - imaging_plane_y_m - selection.y_rx[rx_index], - dz_rx - ); - return attenuation_components(r_tx, r_rx, dz_tx, dz_rx); -} - -[[nodiscard]] auto compensation_weight( - double r_tx, - double r_rx, - double dz_tx, - double dz_rx, - double geo_ref, - double angle_ref, - double range_power, - double angle_power +// attenuation_geo_at_ref_depth, Horns_motion_4x4_ground.py (no angle compensation). +[[nodiscard]] auto attenuation_geo_at_ref_depth( + std::uint32_t tx, std::uint32_t rx, const GeometrySelection& selection, + double ground_depth_m, double pitch_rad ) -> double { - const auto [geo, angle] = attenuation_components(r_tx, r_rx, dz_tx, dz_rx); - const double geo_norm = geo / geo_ref; - const double angle_norm = angle / angle_ref; - - const double range_weight = std::clamp( - 1.0 / (std::pow(geo_norm, range_power) + 1e-12), - 0.0, - kRangeWeightMax - ); - const double angle_weight = std::clamp( - 1.0 / (std::pow(angle_norm, angle_power) + 1e-12), - 0.0, - kAngleWeightMax - ); - return std::clamp(range_weight * angle_weight, 0.0, kTotalWeightMax); + const double x_center = 0.5 * (selection.x_tx[tx] + selection.x_rx[rx]); + const double y = bp::ground_y(kCompensationReferenceDepthM, ground_depth_m, pitch_rad); + const double r_tx = distance_3d(x_center - selection.x_tx[tx], y - selection.y_tx[tx], + kCompensationReferenceDepthM - selection.z_tx[tx]); + const double r_rx = distance_3d(x_center - selection.x_rx[rx], y - selection.y_rx[rx], + kCompensationReferenceDepthM - selection.z_rx[rx]); + return 1.0 / (r_tx * r_rx + 1e-12) + 1e-30; } [[nodiscard]] auto backproject_coherent( @@ -1051,12 +756,11 @@ void normalize_pair_ascans( const std::unordered_map& ascans_by_pair, const GridDefinition& grid, const GeometrySelection& selection, - double imaging_plane_y_m, + double ground_depth_m, double pitch_rad, double velocity_mps, double min_depth_m, double max_depth_m, - double range_power, - double angle_power + double range_power ) -> BpMap { const std::size_t width = grid.x_grid.size(); const std::size_t height = grid.z_grid.size(); @@ -1076,9 +780,6 @@ void normalize_pair_ascans( const std::vector* rx_distances; const AscanResult* ascan; double geo_ref; - double angle_ref; - double z_tx_ant; - double z_rx_ant; }; std::vector contributions; @@ -1089,20 +790,17 @@ void normalize_pair_ascans( if (ascan_it == ascans_by_pair.end() || ascan_it->second.time_s.empty()) { continue; } - const auto [geo_ref, angle_ref] = attenuation_components_at_ref_depth( + const double geo_ref = attenuation_geo_at_ref_depth( trace.tx_local_index, trace.rx_local_index, selection, - imaging_plane_y_m + ground_depth_m, pitch_rad ); contributions.push_back(PairContribution{ &grid.tx_distance_grids[trace.tx_local_index], &grid.rx_distance_grids[trace.rx_local_index], &ascan_it->second, geo_ref, - angle_ref, - selection.z_tx[trace.tx_local_index], - selection.z_rx[trace.rx_local_index], }); } @@ -1133,16 +831,9 @@ void normalize_pair_ascans( } const auto sample = interpolate_complex(ascan, tau_s); - const double weight = compensation_weight( - r_tx, - r_rx, - z_m - contribution.z_tx_ant, - z_m - contribution.z_rx_ant, - contribution.geo_ref, - contribution.angle_ref, - range_power, - angle_power - ); + const double geo = 1.0 / (r_tx * r_rx + 1e-12) + 1e-30; + const double weight = std::clamp(1.0 / (std::pow(geo / contribution.geo_ref, range_power) + 1e-12), + 0.0, kRangeWeightMax); coherent_sum += sample * weight; incoherent_sum += std::abs(sample) * weight; contribution_count += 1.0; @@ -1151,8 +842,8 @@ void normalize_pair_ascans( if (!(contribution_count > 0.0)) { continue; } - coherent_sum /= contribution_count; - incoherent_sum /= contribution_count; + coherent_sum /= contribution_count + 1e-12; + incoherent_sum /= contribution_count + 1e-12; result.coherent[cell_index] = coherent_sum; result.incoherent[cell_index] = incoherent_sum; result.image[cell_index] = std::abs(coherent_sum); @@ -1503,7 +1194,7 @@ void apply_depth_gate( const ObjectRecord& object, const std::vector& selected_traces, const GeometrySelection& selection, - double imaging_plane_y_m + double ground_depth_m, double pitch_rad ) -> std::vector { std::vector signature{}; signature.reserve(selected_traces.size()); @@ -1512,12 +1203,12 @@ void apply_depth_gate( const auto rx = trace.rx_local_index; const double r_tx = distance_3d( object.x_m - selection.x_tx[tx], - imaging_plane_y_m - selection.y_tx[tx], + bp::ground_y(object.z_m, ground_depth_m, pitch_rad) - selection.y_tx[tx], object.z_m - selection.z_tx[tx] ); const double r_rx = distance_3d( object.x_m - selection.x_rx[rx], - imaging_plane_y_m - selection.y_rx[rx], + bp::ground_y(object.z_m, ground_depth_m, pitch_rad) - selection.y_rx[rx], object.z_m - selection.z_rx[rx] ); signature.push_back(0.5 * (r_tx + r_rx)); @@ -1529,12 +1220,12 @@ void mark_sidelobe_candidates( std::vector& objects, const std::vector& selected_traces, const GeometrySelection& selection, - double imaging_plane_y_m + double ground_depth_m, double pitch_rad ) { std::vector> signatures{}; signatures.reserve(objects.size()); for (const auto& object : objects) { - signatures.push_back(bistatic_depth_signature(object, selected_traces, selection, imaging_plane_y_m)); + signatures.push_back(bistatic_depth_signature(object, selected_traces, selection, ground_depth_m, pitch_rad)); } for (std::size_t object_index = 0U; object_index < objects.size(); ++object_index) { @@ -1849,15 +1540,15 @@ void add_bp_score_metrics( } validate_collection_trace_order(run_config, collection); - const auto background_mean = build_background_mean(previous_collections, selection, live_config); + const auto background_mean = bp::background_mean(previous_collections, live_config.gpr_background_subtract_enabled, + live_config.gpr_background_mean_count); auto selected_traces = collect_selected_traces(collection, selection, background_mean); if (selected_traces.empty()) { return results; } - // Coherent BP fixes the medium to vacuum/air (eps_r = 1); the configurable - // relative_permittivity is a legacy-GPR-only knob. - const double velocity_mps = kSpeedOfLightMetersPerSec; + // Coherent BP propagates in air, as in Horns_motion_4x4_ground.py. + const double velocity_mps = bp::kWaveSpeed; const double start_hz = static_cast(live_config.gpr_start_freq_mhz) * 1'000'000.0; const double stop_hz = static_cast(live_config.gpr_stop_freq_mhz) * 1'000'000.0; const double min_depth_m = static_cast(live_config.gpr_min_depth_m); @@ -1866,36 +1557,12 @@ void add_bp_score_metrics( return results; } - // Intra-sweep speed correction before the IFFT. The Tx interleave slot is the - // Tx local index (outputs are sorted), with slot_count = number of selected Tx. - // tx_sweep_time is derived per-collection from the acquisition timestamps: - // one Tx's interleaved sweep spans the whole frame, so capture_span / N_tx. - // Speed comes from the socket-fed live config; look-angle/direction from config. - const auto motion_mode = parse_motion_mode(live_config.gpr_motion_mode); - const std::size_t motion_slot_count = selection.output_positions.size(); - const double capture_span_s = - collection.capture_end_ns > collection.capture_start_ns - ? static_cast(collection.capture_end_ns - collection.capture_start_ns) * 1e-9 - : 0.0; - const double tx_sweep_time_s = - motion_slot_count > 0U ? capture_span_s / static_cast(motion_slot_count) : 0.0; - const double motion_speed_mps = static_cast(live_config.gpr_speed_m_s); - const double motion_look_angle_deg = static_cast(live_config.gpr_look_angle_deg); - const double motion_direction_sign = static_cast(live_config.gpr_direction_sign); + const auto timing = bp::collection_timing(collection); + const double motion_look_angle_deg = live_config.gpr_look_angle_deg; for (auto& trace : selected_traces) { - apply_intra_sweep_motion_correction( - trace, - trace.tx_local_index, - motion_slot_count, - start_hz, - stop_hz, - velocity_mps, - motion_speed_mps, - motion_look_angle_deg, - motion_direction_sign, - tx_sweep_time_s, - motion_mode - ); + bp::apply_motion_phase(trace.s21, collection.traces[trace.run_order].frequency_hz, timing.clocks[trace.run_order], + live_config.gpr_speed_m_s, motion_look_angle_deg * kPi / 180.0, + live_config.gpr_direction_sign, live_config.gpr_motion_mode == "int_focus", start_hz, stop_hz); } std::unordered_map ascans_by_pair{}; @@ -1912,9 +1579,30 @@ void add_bp_score_metrics( if (ascans_by_pair.empty()) { return results; } - normalize_pair_ascans(ascans_by_pair, velocity_mps, min_depth_m, max_depth_m); - - const double imaging_plane_y_m = static_cast(live_config.gpr_imaging_plane_y_m); + std::vector>> pairs; + for (const auto& trace : selected_traces) { + const auto found = ascans_by_pair.find(make_pair_key(trace.tx_local_index, trace.rx_local_index)); + if (found != ascans_by_pair.end()) { pairs.emplace_back(found->second.samples); } + } + const double dt = ascans_by_pair.begin()->second.dt_s; + double x_min = selection.x_tx.front(), x_max = x_min, max_abs_x = 0.0; + for (const auto* xs : {&selection.x_tx, &selection.x_rx}) { + for (const double x : *xs) { + x_min = std::min(x_min, x); x_max = std::max(x_max, x); max_abs_x = std::max(max_abs_x, std::abs(x)); + } + } + if (live_config.pair_depth_equalize) { + for (const auto& [key, ascan] : ascans_by_pair) { + if (ascan.dt_s != dt) { throw std::runtime_error("PAIR_DEPTH_EQUALIZE requires a common frequency grid"); } + } + bp::equalize_pairs_by_depth(pairs, dt * velocity_mps / 2.0, + bp::equalize_depth_limit(max_depth_m, x_min - kXMarginM, x_max + kXMarginM, max_abs_x)); + } + for (auto& [key, ascan] : ascans_by_pair) { + bp::normalize_pair(std::span(ascan.samples), ascan.dt_s * velocity_mps / 2.0, min_depth_m, max_depth_m); + } + const double ground_depth_m = live_config.gpr_bp_ground_depth_m; + const double pitch_rad = motion_look_angle_deg * kPi / 180.0; // The imaging grid (90k cells x 6 distance fields) only depends on antenna // geometry and the depth window, which are constant across a run, so it is @@ -1923,9 +1611,9 @@ void add_bp_score_metrics( // backprojection workers join before it returns), so a local static is safe. static std::vector cached_grid_signature; static GridDefinition cached_grid; - const auto signature = grid_signature(selection, max_depth_m, kGridZMinM, imaging_plane_y_m); + const auto signature = grid_signature(selection, max_depth_m, kGridZMinM, ground_depth_m, pitch_rad); if (cached_grid.x_grid.empty() || signature != cached_grid_signature) { - cached_grid = build_grid(selection, max_depth_m, kGridZMinM, imaging_plane_y_m); + cached_grid = build_grid(selection, max_depth_m, kGridZMinM, ground_depth_m, pitch_rad); cached_grid_signature = signature; } const GridDefinition& grid = cached_grid; @@ -1938,12 +1626,11 @@ void add_bp_score_metrics( ascans_by_pair, grid, selection, - imaging_plane_y_m, + ground_depth_m, pitch_rad, velocity_mps, min_depth_m, max_depth_m, - std::max(0.0, static_cast(live_config.gpr_range_comp_power)), - 0.0 // angle compensation is fixed off for coherent BP + std::max(0.0, static_cast(live_config.gpr_range_comp_power)) ); if (bp.image.empty()) { return results; @@ -1956,7 +1643,7 @@ void add_bp_score_metrics( auto objects = find_bp_objects(display_map, grid, static_cast(live_config.gpr_object_min_frac)); add_local_prominence_metrics(objects, display_map, grid, min_depth_m, max_depth_m); add_incoherent_support_metrics(objects, incoherent_display_map, bp.coherence_factor); - mark_sidelobe_candidates(objects, selected_traces, selection, imaging_plane_y_m); + mark_sidelobe_candidates(objects, selected_traces, selection, ground_depth_m, pitch_rad); add_bp_score_metrics(objects, live_config); if (live_config.gpr_remove_sidelobe_objects_enabled) { diff --git a/data_acq_and_processing/processing/processors/src/gpr_legacy_processor.ipp b/data_acq_and_processing/processing/processors/src/gpr_legacy_processor.ipp deleted file mode 100644 index 46072f4..0000000 --- a/data_acq_and_processing/processing/processors/src/gpr_legacy_processor.ipp +++ /dev/null @@ -1,1066 +0,0 @@ -// Legacy MIMO GPR — ellipse-intersection localizer with motion compensation. -// -// This translation unit is included from `gpr_processor.cpp` *after* -// `gpr_backprojection_processor.ipp`, which defines the shared building blocks -// (kPi, SelectedTrace, GeometrySelection, distance_3d, fft_inplace, ...). Do -// not include this file directly. -// -// Pipeline overview (mirrors the Python reference Ellips_motion_remake_2.py): -// 1. Pre-process each pair: optional background subtraction (already done in -// collect_selected_traces); optional intra-sweep phase correction of S21 -// before the IFFT — compensates radar displacement that happens *inside* -// one sweep where the frequencies are stepped linearly in time. -// 2. IFFT each pair → A-scan; find peaks above an SNR threshold within the -// depth gate; record both raw and attenuation-compensated SNR. -// 3. Apply coarse per-event motion correction so every peak carries a -// motion-corrected (tau_corr, z_corr) on top of the apparent values. -// 4. CLEAN-style iterative ellipse intersection: build a soft Gaussian-shell -// accumulator, take the strongest pixel, count agreeing pairs, suppress -// its depth band, repeat. -// 5. Optional extended-mode region detection for diffuse reflectors. -// -// Sweep-event model — radar topology matters: -// * Matrix radars (`librevna_multi`) fire one Tx at a time and -// receive on all Rx channels in parallel. A run of 8 traces is just 2 -// Tx-events; all (tx_k, *) pairs share one timestamp. -// * Sequential radars (single librevna with switches, rph_radar) -// measure each pair separately. A run of 8 traces is 8 separate sweeps; -// every pair has its own timestamp. -// Both cases reduce to: there are N events in the frame, each event lasts -// `event_duration_s = (capture_end_ns - capture_start_ns) / N`. What changes -// is how event indices are assigned to pairs — by Tx for matrix radars, by -// trace run-order for sequential ones. `event_duration_s` is derived from -// the collection metadata, never from a live-config knob. - -constexpr double kLegacyGridZMinM = 0.20; -constexpr double kLegacySmoothSigma = 3.0; -constexpr double kLegacyCleanSuppressRadiusM = 0.07; -constexpr double kLegacyCleanThresholdFrac = 0.05; -constexpr std::size_t kLegacyMaxObjects = 15U; -constexpr double kLegacyExtendedThresholdFrac = 0.75; -constexpr double kLegacyExtendedMinAreaCm2 = 2.0; -constexpr double kLegacyAttenuationReferenceDepthM = 3.0; - -struct LegacyAscanResult { - std::vector time_s{}; - std::vector depth_m{}; - std::vector amplitude{}; - double bandwidth_hz = 0.0; -}; - -struct LegacyPeakRecord { - double z_app = 0.0; - double tau = 0.0; - double tau_corr = 0.0; - double z_corr = 0.0; - double snr_raw = 0.0; - double snr_comp = 0.0; -}; - -struct LegacyPointRecord { - double x_m = 0.0; - double z_m = 0.0; - double score = 0.0; -}; - -struct LegacyRegionRecord { - double x_m = 0.0; - double z_m = 0.0; - double score = 0.0; - double pixel_count = 0.0; - std::vector mask{}; -}; - -// Per-pair timing snapshot. For matrix radars, all (tx_k, *) pairs share the -// same row; for sequential radars every pair has a distinct row. Indexing by -// pair keeps the rest of the pipeline ignorant of the radar topology. -struct LegacyEventTiming { - std::size_t event_index = 0U; // 0-based order of the event in the frame - double t_start_s = 0.0; // start of this event sweep relative to frame - double t_center_s = 0.0; // center of this event sweep - double dt_ref_s = 0.0; // t_center_s - t_frame_ref_s - double dz_motion_m = 0.0; // direction_sign * speed * dt_ref * cos(theta) - double dtau_motion_s = 0.0; // 2 * dz_motion_m / velocity -}; - -struct LegacyMotionTiming { - double event_duration_s = 0.0; // (capture_end - capture_start) / num_events - double cos_look_angle = 1.0; - double direction_sign = 1.0; - double speed_m_s = 0.0; - bool apply_intra_sweep_phase = false; - bool parallel_rx_per_tx_event = false; - std::unordered_map by_pair{}; -}; - -enum class LegacyPeakDomain { - Apparent, - Corrected, -}; - -[[nodiscard]] auto find_legacy_peak_indices( - const std::vector& values, - std::size_t start_index, - std::size_t stop_index, - double threshold, - std::size_t min_distance -) -> std::vector { - if (stop_index <= start_index + 2U) { - return {}; - } - - std::vector candidates{}; - for (std::size_t index = start_index + 1U; index + 1U < stop_index; ++index) { - const bool below_threshold = values[index] < threshold; - const bool not_rising = values[index] <= values[index - 1U]; - const bool still_rising = values[index] < values[index + 1U]; - if (below_threshold || not_rising || still_rising) { - continue; - } - candidates.push_back(index); - } - - std::sort(candidates.begin(), candidates.end(), [&](std::size_t left, std::size_t right) { - return values[left] > values[right]; - }); - - std::vector selected{}; - for (const auto index : candidates) { - bool keep = true; - for (const auto accepted : selected) { - const auto distance = accepted > index ? accepted - index : index - accepted; - if (distance < min_distance) { - keep = false; - break; - } - } - if (keep) { - selected.push_back(index); - } - } - - std::sort(selected.begin(), selected.end()); - return selected; -} - -// Expected geo*pattern attenuation for a target directly under the virtual -// pair center at depth `z_app`. Uses full 3D antenna coordinates so this -// generalizes to non-coplanar antenna layouts; boresight is taken along +Z. -[[nodiscard]] auto legacy_attenuation_at_depth( - std::size_t tx_index, - std::size_t rx_index, - double z_app, - const GeometrySelection& selection, - double imaging_plane_y_m -) -> double { - const double x_center = 0.5 * (selection.x_tx[tx_index] + selection.x_rx[rx_index]); - const double r_tx = distance_3d( - x_center - selection.x_tx[tx_index], - imaging_plane_y_m - selection.y_tx[tx_index], - z_app - selection.z_tx[tx_index] - ); - const double r_rx = distance_3d( - x_center - selection.x_rx[rx_index], - imaging_plane_y_m - selection.y_rx[rx_index], - z_app - selection.z_rx[rx_index] - ); - const double cos_tx = (z_app - selection.z_tx[tx_index]) / (r_tx + 1e-12); - const double cos_rx = (z_app - selection.z_rx[rx_index]) / (r_rx + 1e-12); - const double geo = 1.0 / ((r_tx * r_rx) + 1e-12); - const double pattern = (cos_tx * cos_tx) * (cos_rx * cos_rx); - return (geo * pattern) + 1e-30; -} - -[[nodiscard]] auto lower_bound_index(const std::vector& axis, double value) -> std::size_t { - const auto found = std::lower_bound(axis.begin(), axis.end(), value); - return static_cast(std::distance(axis.begin(), found)); -} - -[[nodiscard]] auto compute_legacy_ascan( - const SelectedTrace& trace, - double start_hz, - double stop_hz, - double velocity_mps -) -> LegacyAscanResult { - LegacyAscanResult result{}; - if (trace.frequency_hz.size() != trace.s21.size()) { - return result; - } - - const double low_hz = std::min(start_hz, stop_hz); - const double high_hz = std::max(start_hz, stop_hz); - std::vector frequency_hz{}; - std::vector> s21{}; - frequency_hz.reserve(trace.frequency_hz.size()); - s21.reserve(trace.s21.size()); - for (std::size_t index = 0U; index < trace.frequency_hz.size(); ++index) { - const double frequency_value = trace.frequency_hz[index]; - if (frequency_value < low_hz || frequency_value > high_hz) { - continue; - } - frequency_hz.push_back(frequency_value); - s21.push_back(trace.s21[index]); - } - - if (frequency_hz.size() < 2U) { - return result; - } - - const std::size_t point_count = frequency_hz.size(); - const double df_hz = (frequency_hz.back() - frequency_hz.front()) / static_cast(point_count - 1U); - if (!(df_hz > 0.0)) { - return result; - } - - const auto start_bin = static_cast(std::llround(frequency_hz.front() / df_hz)); - if (start_bin < 0) { - return result; - } - - const auto start_index = static_cast(start_bin); - const std::size_t min_fft_len = 2U * (start_index + point_count - 1U); - const std::size_t fft_len = next_power_of_two(min_fft_len); - if (fft_len < min_fft_len || start_index > fft_len || point_count > (fft_len - start_index)) { - return result; - } - - std::vector> spectrum(fft_len, std::complex(0.0, 0.0)); - for (std::size_t index = 0U; index < point_count; ++index) { - const double window = point_count > 1U - ? 0.5 - (0.5 * std::cos((2.0 * kPi * static_cast(index)) / static_cast(point_count - 1U))) - : 1.0; - spectrum[start_index + index] = s21[index] * window; - } - - fft_inplace(spectrum, true); - - result.bandwidth_hz = frequency_hz.back() - frequency_hz.front(); - const double dt_s = 1.0 / (static_cast(fft_len) * df_hz); - result.time_s.resize(fft_len, 0.0); - result.depth_m.resize(fft_len, 0.0); - result.amplitude.resize(fft_len, 0.0); - for (std::size_t index = 0U; index < fft_len; ++index) { - result.time_s[index] = static_cast(index) * dt_s; - result.depth_m[index] = result.time_s[index] * velocity_mps * 0.5; - result.amplitude[index] = std::abs(spectrum[index]); - } - return result; -} - -[[nodiscard]] auto legacy_peak_depth_for_domain(const LegacyPeakRecord& peak, LegacyPeakDomain domain) -> double { - return domain == LegacyPeakDomain::Corrected ? peak.z_corr : peak.z_app; -} - -[[nodiscard]] auto legacy_peak_tau_for_domain(const LegacyPeakRecord& peak, LegacyPeakDomain domain) -> double { - return domain == LegacyPeakDomain::Corrected ? peak.tau_corr : peak.tau; -} - -[[nodiscard]] auto is_legacy_depth_excluded( - double z_value, - const std::vector>& ranges -) -> bool { - for (const auto& [low, high] : ranges) { - if (z_value >= low && z_value <= high) { - return true; - } - } - return false; -} - -// True when `z_value` should participate in the ellipse vote — both the static -// depth gate ([min, max]) and the per-step CLEAN-suppression bands must allow it. -// Mirrors Python's `_peak_in_work_depth` + the per-step `excl_z` filter. -[[nodiscard]] auto is_legacy_depth_active( - double z_value, - double min_depth_m, - double max_depth_m, - const std::vector>& excluded_ranges -) -> bool { - if (z_value < min_depth_m || z_value > max_depth_m) { - return false; - } - return !is_legacy_depth_excluded(z_value, excluded_ranges); -} - -[[nodiscard]] auto is_matrix_radar_model(const std::string& model) -> bool { - return model == "librevna_multi"; -} - -// Assign an `event_index` to every selected pair. The mapping depends on the -// radar topology: -// * Matrix radar — all (tx_k, *) pairs share one event, ordered by the Tx's -// first appearance in run order. So 8 traces with 2 Tx's give 2 events. -// * Sequential radar — every pair is its own event, ordered by run order. -// So 8 traces give 8 events. -[[nodiscard]] auto assign_event_indices( - const std::vector& selected_traces, - bool matrix_radar -) -> std::pair, std::size_t> { - std::unordered_map event_index_by_pair{}; - event_index_by_pair.reserve(selected_traces.size()); - - if (matrix_radar) { - std::unordered_map event_by_tx{}; - std::size_t next_event = 0U; - for (const auto& trace : selected_traces) { - const auto [event_it, inserted] = event_by_tx.try_emplace(trace.tx_local_index, next_event); - if (inserted) { - ++next_event; - } - event_index_by_pair[make_pair_key(trace.tx_local_index, trace.rx_local_index)] = event_it->second; - } - return {std::move(event_index_by_pair), next_event}; - } - - // Sequential mode: rank traces by their run_order so the event index is a - // dense 0..N-1 sequence regardless of any holes in run_order. - std::vector> ordered{}; - ordered.reserve(selected_traces.size()); - for (const auto& trace : selected_traces) { - ordered.emplace_back(trace.run_order, make_pair_key(trace.tx_local_index, trace.rx_local_index)); - } - std::sort(ordered.begin(), ordered.end(), - [](const auto& lhs, const auto& rhs) { return lhs.first < rhs.first; }); - for (std::size_t event_index = 0U; event_index < ordered.size(); ++event_index) { - event_index_by_pair[ordered[event_index].second] = event_index; - } - return {std::move(event_index_by_pair), ordered.size()}; -} - -// Build per-pair motion timing. `event_duration_s` is derived from collection -// metadata as `(capture_end_ns - capture_start_ns) / num_events` — it is the -// duration of one sweep event in the frame, never a live-config knob. If the -// motion model is disabled (speed = 0 and phase correction off), the function -// still returns one row per pair so downstream code can index uniformly. -[[nodiscard]] auto compute_legacy_motion_timing( - const std::vector& selected_traces, - bool matrix_radar, - std::uint64_t capture_start_ns, - std::uint64_t capture_end_ns, - const ProcessingLiveConfig& live_config, - double velocity_mps -) -> LegacyMotionTiming { - LegacyMotionTiming timing{}; - timing.speed_m_s = static_cast(live_config.gpr_speed_m_s); - timing.direction_sign = static_cast(live_config.gpr_direction_sign); - timing.cos_look_angle = std::cos((static_cast(live_config.gpr_look_angle_deg) * kPi) / 180.0); - timing.apply_intra_sweep_phase = live_config.gpr_apply_freq_phase_correction; - timing.parallel_rx_per_tx_event = matrix_radar; - - if (selected_traces.empty()) { - return timing; - } - - auto [event_index_by_pair, num_events] = assign_event_indices(selected_traces, matrix_radar); - if (num_events == 0U) { - return timing; - } - - const bool speed_meaningful = std::abs(timing.speed_m_s) > 1e-12; - const bool model_active = speed_meaningful || timing.apply_intra_sweep_phase; - - // No motion and no phase correction — populate with zeroed rows and bail. - if (!model_active) { - for (const auto& [pair_key, event_index] : event_index_by_pair) { - timing.by_pair.emplace(pair_key, LegacyEventTiming{.event_index = event_index}); - } - return timing; - } - - if (capture_end_ns <= capture_start_ns) { - throw std::runtime_error( - "Legacy GPR motion model requires valid capture_start_ns/capture_end_ns metadata" - ); - } - - const double total_span_s = static_cast(capture_end_ns - capture_start_ns) * 1e-9; - const double event_duration_s = total_span_s / static_cast(num_events); - if (!(event_duration_s > 0.0)) { - throw std::runtime_error("Legacy GPR motion model requires positive per-event duration"); - } - timing.event_duration_s = event_duration_s; - - // Reference anchor for dt_ref: either the midpoint between the first and - // last event centers (frame_center) or just the first event center - // (first_tx_event). Matches Python's `MOTION_CONFIG.reference_mode`. - const double first_center_s = 0.5 * event_duration_s; - const double last_center_s = (static_cast(num_events) - 0.5) * event_duration_s; - const double t_ref_s = live_config.gpr_reference_mode == LegacyGprReferenceMode::FirstTxEvent - ? first_center_s - : 0.5 * (first_center_s + last_center_s); - - const double motion_factor = speed_meaningful - ? timing.direction_sign * timing.speed_m_s * timing.cos_look_angle - : 0.0; - - for (const auto& [pair_key, event_index] : event_index_by_pair) { - LegacyEventTiming row{}; - row.event_index = event_index; - row.t_start_s = static_cast(event_index) * event_duration_s; - row.t_center_s = row.t_start_s + (0.5 * event_duration_s); - row.dt_ref_s = row.t_center_s - t_ref_s; - row.dz_motion_m = motion_factor * row.dt_ref_s; - row.dtau_motion_s = (2.0 * row.dz_motion_m) / velocity_mps; - timing.by_pair.emplace(pair_key, row); - } - - return timing; -} - -// Compensate for the radar moving while a single sweep is being recorded. -// Frequencies inside one sweep are stepped linearly in time, so each frequency -// is sampled from a slightly different antenna position. The correction shifts -// each frequency's phase back to the event center; after that the IFFT -// produces an A-scan as if the whole sweep were captured at one position. -void apply_intra_sweep_phase_correction( - SelectedTrace& trace, - const LegacyEventTiming& timing, - const LegacyMotionTiming& motion, - double velocity_mps -) { - if (!motion.apply_intra_sweep_phase || !(motion.event_duration_s > 0.0)) { - return; - } - if (!(std::abs(motion.speed_m_s) > 1e-12)) { - return; // No motion → zero phase shift, no-op. - } - const std::size_t point_count = trace.frequency_hz.size(); - if (point_count < 2U || trace.s21.size() != point_count) { - return; - } - - const double dt_freq_s = motion.event_duration_s / static_cast(point_count - 1U); - const double motion_factor = motion.direction_sign * motion.speed_m_s * motion.cos_look_angle; - - for (std::size_t index = 0U; index < point_count; ++index) { - const double t_abs_s = timing.t_start_s + (static_cast(index) * dt_freq_s); - const double dt_intra_s = t_abs_s - timing.t_center_s; - const double delta_path_m = 2.0 * motion_factor * dt_intra_s; - const double phi = (2.0 * kPi * trace.frequency_hz[index] * delta_path_m) / velocity_mps; - trace.s21[index] *= std::polar(1.0, phi); - } -} - -void apply_legacy_motion_correction( - std::unordered_map>& peaks_by_pair, - const LegacyMotionTiming& motion, - double velocity_mps -) { - for (auto& [key, peaks] : peaks_by_pair) { - const auto timing_it = motion.by_pair.find(key); - if (timing_it == motion.by_pair.end()) { - throw std::runtime_error("Missing motion timing for selected legacy GPR pair"); - } - const double dtau_motion_s = timing_it->second.dtau_motion_s; - for (auto& peak : peaks) { - peak.tau_corr = peak.tau + dtau_motion_s; - peak.z_corr = 0.5 * velocity_mps * peak.tau_corr; - } - } -} - -[[nodiscard]] auto build_legacy_accumulator( - const GridDefinition& grid, - const std::unordered_map>& peaks_by_pair, - const std::vector>& exclude_ranges, - double velocity_mps, - double shell_sigma_m, - double min_depth_m, - double max_depth_m, - const std::vector& x_tx, - const std::vector& x_rx, - LegacyPeakDomain domain -) -> std::vector { - const std::size_t width = grid.x_grid.size(); - const std::size_t height = grid.z_grid.size(); - std::vector accumulator(width * height, 0.0); - if (!(shell_sigma_m > 0.0)) { - return accumulator; - } - - for (std::size_t tx_index = 0U; tx_index < x_tx.size(); ++tx_index) { - for (std::size_t rx_index = 0U; rx_index < x_rx.size(); ++rx_index) { - const auto peak_it = peaks_by_pair.find( - make_pair_key(static_cast(tx_index), static_cast(rx_index)) - ); - if (peak_it == peaks_by_pair.end()) { - continue; - } - - const auto& tx_grid = grid.tx_distance_grids[tx_index]; - const auto& rx_grid = grid.rx_distance_grids[rx_index]; - for (const auto& peak : peak_it->second) { - const double depth = legacy_peak_depth_for_domain(peak, domain); - if (!is_legacy_depth_active(depth, min_depth_m, max_depth_m, exclude_ranges)) { - continue; - } - - const double range_total = velocity_mps * legacy_peak_tau_for_domain(peak, domain); - for (std::size_t cell_index = 0U; cell_index < accumulator.size(); ++cell_index) { - const double residual = tx_grid[cell_index] + rx_grid[cell_index] - range_total; - const double shell = std::exp(-0.5 * std::pow(residual / shell_sigma_m, 2.0)); - accumulator[cell_index] += shell * peak.snr_comp; - } - } - } - } - - return accumulator; -} - -[[nodiscard]] auto count_legacy_agreeing_ellipses( - double x_est, - double z_est, - const std::unordered_map>& peaks_by_pair, - const std::vector>& exclude_ranges, - const std::vector& x_tx, - const std::vector& x_rx, - double velocity_mps, - double shell_sigma_m, - double min_depth_m, - double max_depth_m, - LegacyPeakDomain domain -) -> double { - std::size_t count = 0U; - for (std::size_t tx_index = 0U; tx_index < x_tx.size(); ++tx_index) { - for (std::size_t rx_index = 0U; rx_index < x_rx.size(); ++rx_index) { - const auto peak_it = peaks_by_pair.find( - make_pair_key(static_cast(tx_index), static_cast(rx_index)) - ); - if (peak_it == peaks_by_pair.end()) { - continue; - } - - for (const auto& peak : peak_it->second) { - const double depth = legacy_peak_depth_for_domain(peak, domain); - if (!is_legacy_depth_active(depth, min_depth_m, max_depth_m, exclude_ranges)) { - continue; - } - - const double r_tx = std::hypot(x_est - x_tx[tx_index], z_est); - const double r_rx = std::hypot(x_est - x_rx[rx_index], z_est); - if (std::abs((r_tx + r_rx) - (velocity_mps * legacy_peak_tau_for_domain(peak, domain))) < - shell_sigma_m * 6.0) { - count += 1U; - break; - } - } - } - } - - return static_cast(count); -} - -[[nodiscard]] auto find_legacy_centroid( - const std::vector& values, - std::size_t width, - std::size_t height, - std::size_t row, - std::size_t col, - std::size_t radius_z, - std::size_t radius_x, - const std::vector& x_grid, - const std::vector& z_grid -) -> std::pair { - if (values.empty()) { - return {x_grid[col], z_grid[row]}; - } - - const auto row_start = row > radius_z ? row - radius_z : 0U; - const auto row_stop = std::min(height, row + radius_z + 1U); - const auto col_start = col > radius_x ? col - radius_x : 0U; - const auto col_stop = std::min(width, col + radius_x + 1U); - - double weight_sum = 0.0; - double row_weighted_sum = 0.0; - double col_weighted_sum = 0.0; - for (std::size_t sample_row = row_start; sample_row < row_stop; ++sample_row) { - for (std::size_t sample_col = col_start; sample_col < col_stop; ++sample_col) { - const double weight = values[(sample_row * width) + sample_col]; - weight_sum += weight; - row_weighted_sum += static_cast(sample_row) * weight; - col_weighted_sum += static_cast(sample_col) * weight; - } - } - - if (!(weight_sum > 0.0)) { - return {x_grid[col], z_grid[row]}; - } - - const auto centroid_row = - clamp_index(static_cast(std::llround(row_weighted_sum / weight_sum)), height); - const auto centroid_col = - clamp_index(static_cast(std::llround(col_weighted_sum / weight_sum)), width); - return {x_grid[centroid_col], z_grid[centroid_row]}; -} - -[[nodiscard]] auto clean_legacy_find_points( - const GridDefinition& grid, - const std::unordered_map>& peaks_by_pair, - const std::vector& x_tx, - const std::vector& x_rx, - double velocity_mps, - double shell_sigma_m, - double min_depth_m, - double max_depth_m, - LegacyPeakDomain domain -) -> std::pair, std::vector> { - std::vector found{}; - const auto accumulator = build_legacy_accumulator( - grid, peaks_by_pair, {}, velocity_mps, shell_sigma_m, min_depth_m, max_depth_m, x_tx, x_rx, domain - ); - const double initial_max = max_value(accumulator); - if (!(initial_max > 0.0) || grid.x_grid.size() < 2U || grid.z_grid.size() < 2U) { - return {found, gaussian_filter_2d(accumulator, grid.x_grid.size(), grid.z_grid.size(), kLegacySmoothSigma)}; - } - - // Suppression radius rounds down to mirror Python's `int(0.07 / dx)`. - const double dx = grid.x_grid[1] - grid.x_grid[0]; - const double dz = grid.z_grid[1] - grid.z_grid[0]; - const auto radius_x = - static_cast(std::max(1.0, std::floor(kLegacyCleanSuppressRadiusM / std::max(dx, 1e-6)))); - const auto radius_z = - static_cast(std::max(1.0, std::floor(kLegacyCleanSuppressRadiusM / std::max(dz, 1e-6)))); - - std::vector> excluded_ranges{}; - for (std::size_t step = 0U; step < kLegacyMaxObjects; ++step) { - const auto current = build_legacy_accumulator( - grid, - peaks_by_pair, - excluded_ranges, - velocity_mps, - shell_sigma_m, - min_depth_m, - max_depth_m, - x_tx, - x_rx, - domain - ); - const auto smoothed = gaussian_filter_2d(current, grid.x_grid.size(), grid.z_grid.size(), kLegacySmoothSigma); - const double smoothed_max = max_value(smoothed); - if (!(smoothed_max > (kLegacyCleanThresholdFrac * initial_max))) { - break; - } - - const auto max_it = std::max_element(smoothed.begin(), smoothed.end()); - const auto max_index = static_cast(std::distance(smoothed.begin(), max_it)); - const auto peak_row = max_index / grid.x_grid.size(); - const auto peak_col = max_index % grid.x_grid.size(); - const auto [x_est, z_est] = find_legacy_centroid( - smoothed, - grid.x_grid.size(), - grid.z_grid.size(), - peak_row, - peak_col, - radius_z, - radius_x, - grid.x_grid, - grid.z_grid - ); - - found.push_back( - LegacyPointRecord{ - .x_m = x_est, - .z_m = z_est, - .score = count_legacy_agreeing_ellipses( - x_est, - z_est, - peaks_by_pair, - excluded_ranges, - x_tx, - x_rx, - velocity_mps, - shell_sigma_m, - min_depth_m, - max_depth_m, - domain - ), - } - ); - - // Collect depths of all peaks consistent with the just-detected point; - // they form the next exclusion band so subsequent CLEAN steps cannot - // re-pick the same target. - std::vector matched_depths{}; - for (std::size_t tx_index = 0U; tx_index < x_tx.size(); ++tx_index) { - for (std::size_t rx_index = 0U; rx_index < x_rx.size(); ++rx_index) { - const auto peak_it = peaks_by_pair.find( - make_pair_key(static_cast(tx_index), static_cast(rx_index)) - ); - if (peak_it == peaks_by_pair.end()) { - continue; - } - for (const auto& peak : peak_it->second) { - const double depth = legacy_peak_depth_for_domain(peak, domain); - if (!is_legacy_depth_active(depth, min_depth_m, max_depth_m, excluded_ranges)) { - continue; - } - const double r_tx = std::hypot(x_est - x_tx[tx_index], z_est); - const double r_rx = std::hypot(x_est - x_rx[rx_index], z_est); - if (std::abs((r_tx + r_rx) - (velocity_mps * legacy_peak_tau_for_domain(peak, domain))) < - shell_sigma_m * 3.0) { - matched_depths.push_back(depth); - } - } - } - } - - if (!matched_depths.empty()) { - const auto [min_it, max_it_depth] = std::minmax_element(matched_depths.begin(), matched_depths.end()); - excluded_ranges.emplace_back(*min_it - shell_sigma_m, *max_it_depth + shell_sigma_m); - } - } - - return {found, gaussian_filter_2d(accumulator, grid.x_grid.size(), grid.z_grid.size(), kLegacySmoothSigma)}; -} - -[[nodiscard]] auto extended_legacy_find_regions( - const GridDefinition& grid, - const std::unordered_map>& peaks_by_pair, - const std::vector& x_tx, - const std::vector& x_rx, - double velocity_mps, - double shell_sigma_m, - double min_depth_m, - double max_depth_m -) -> std::pair, std::vector> { - std::vector regions{}; - const auto accumulator = build_legacy_accumulator( - grid, - peaks_by_pair, - {}, - velocity_mps, - shell_sigma_m, - min_depth_m, - max_depth_m, - x_tx, - x_rx, - LegacyPeakDomain::Apparent - ); - const auto smoothed = gaussian_filter_2d(accumulator, grid.x_grid.size(), grid.z_grid.size(), kLegacySmoothSigma); - const double smoothed_max = max_value(smoothed); - if (!(smoothed_max > 0.0) || grid.x_grid.size() < 2U || grid.z_grid.size() < 2U) { - return {regions, smoothed}; - } - - const double dx_cm = std::abs(grid.x_grid[1] - grid.x_grid[0]) * 100.0; - const double dz_cm = std::abs(grid.z_grid[1] - grid.z_grid[0]) * 100.0; - const double pixel_area_cm2 = std::max(dx_cm * dz_cm, 1e-6); - const auto min_pixels = - static_cast(std::max(1.0, std::floor(kLegacyExtendedMinAreaCm2 / pixel_area_cm2))); - - const std::size_t width = grid.x_grid.size(); - const std::size_t height = grid.z_grid.size(); - const double threshold = kLegacyExtendedThresholdFrac * smoothed_max; - std::vector visited(width * height, 0U); - - for (std::size_t row = 0U; row < height; ++row) { - for (std::size_t col = 0U; col < width; ++col) { - const auto start_index = (row * width) + col; - if (visited[start_index] != 0U || smoothed[start_index] <= threshold) { - continue; - } - - std::vector stack{start_index}; - std::vector component{}; - visited[start_index] = 1U; - - while (!stack.empty()) { - const auto cell_index = stack.back(); - stack.pop_back(); - component.push_back(cell_index); - - const auto cell_row = cell_index / width; - const auto cell_col = cell_index % width; - const std::pair offsets[] = { - {-1, 0}, - {1, 0}, - {0, -1}, - {0, 1}, - }; - - for (const auto& [row_offset, col_offset] : offsets) { - const auto next_row = static_cast(cell_row) + row_offset; - const auto next_col = static_cast(cell_col) + col_offset; - if (next_row < 0 || next_col < 0) { - continue; - } - const bool row_out_of_bounds = next_row >= static_cast(height); - const bool col_out_of_bounds = next_col >= static_cast(width); - if (row_out_of_bounds || col_out_of_bounds) { - continue; - } - const auto next_index = - (static_cast(next_row) * width) + static_cast(next_col); - if (visited[next_index] != 0U || smoothed[next_index] <= threshold) { - continue; - } - visited[next_index] = 1U; - stack.push_back(next_index); - } - } - - if (component.size() < min_pixels) { - continue; - } - - LegacyRegionRecord region{}; - region.mask.assign(width * height, 0.0F); - double weight_sum = 0.0; - double x_weight_sum = 0.0; - double z_weight_sum = 0.0; - for (const auto cell_index : component) { - const auto cell_row = cell_index / width; - const auto cell_col = cell_index % width; - const double weight = smoothed[cell_index]; - weight_sum += weight; - x_weight_sum += grid.x_grid[cell_col] * weight; - z_weight_sum += grid.z_grid[cell_row] * weight; - region.mask[cell_index] = 1.0F; - } - - if (!(weight_sum > 0.0)) { - continue; - } - - region.x_m = x_weight_sum / weight_sum; - region.z_m = z_weight_sum / weight_sum; - region.score = count_legacy_agreeing_ellipses( - region.x_m, - region.z_m, - peaks_by_pair, - {}, - x_tx, - x_rx, - velocity_mps, - shell_sigma_m, - min_depth_m, - max_depth_m, - LegacyPeakDomain::Apparent - ); - region.pixel_count = static_cast(component.size()); - regions.push_back(std::move(region)); - } - } - - return {regions, smoothed}; -} - -[[nodiscard]] auto process_legacy_gpr( - const config::RunConfig& run_config, - const ipc::PreprocessedCollection& collection, - std::span previous_collections, - const ProcessingLiveConfig& live_config -) -> ipc::ResultCollection { - ipc::ResultCollection results{}; - - const auto selection = build_geometry_selection(run_config, live_config); - if (selection.input_positions.empty() || selection.output_positions.empty()) { - return results; - } - - validate_collection_trace_order(run_config, collection); - const auto background_mean = build_background_mean(previous_collections, selection, live_config); - auto selected_traces = collect_selected_traces(collection, selection, background_mean); - if (selected_traces.empty()) { - return results; - } - - const double velocity_mps = - kSpeedOfLightMetersPerSec / std::sqrt(std::max(1e-6, static_cast(run_config.gpr.relative_permittivity))); - const double start_hz = static_cast(live_config.gpr_start_freq_mhz) * 1'000'000.0; - const double stop_hz = static_cast(live_config.gpr_stop_freq_mhz) * 1'000'000.0; - const double min_depth_m = static_cast(live_config.gpr_min_depth_m); - const double max_depth_m = static_cast(live_config.gpr_max_depth_m); - if (!(max_depth_m > min_depth_m)) { - return results; - } - const double imaging_plane_y_m = static_cast(live_config.gpr_imaging_plane_y_m); - - // Motion timing is computed once per collection. Matrix radars get one - // event per Tx (parallel Rx); sequential radars get one event per pair. - const bool matrix_radar = is_matrix_radar_model(run_config.radar_model); - const auto motion_timing = compute_legacy_motion_timing( - selected_traces, - matrix_radar, - collection.capture_start_ns, - collection.capture_end_ns, - live_config, - velocity_mps - ); - - // Intra-sweep phase correction (frequency-domain) — happens BEFORE the IFFT - // because it modifies the S21 spectrum that compute_legacy_ascan consumes. - if (motion_timing.apply_intra_sweep_phase) { - for (auto& trace : selected_traces) { - const auto pair_key = make_pair_key(trace.tx_local_index, trace.rx_local_index); - const auto timing_it = motion_timing.by_pair.find(pair_key); - if (timing_it == motion_timing.by_pair.end()) { - continue; - } - apply_intra_sweep_phase_correction(trace, timing_it->second, motion_timing, velocity_mps); - } - } - - std::unordered_map ascans_by_pair{}; - double bandwidth_hz = 0.0; - for (const auto& trace : selected_traces) { - auto ascan = compute_legacy_ascan(trace, start_hz, stop_hz, velocity_mps); - if (ascan.amplitude.empty() || !(ascan.bandwidth_hz > 0.0)) { - continue; - } - bandwidth_hz = std::max(bandwidth_hz, ascan.bandwidth_hz); - ascans_by_pair.emplace(make_pair_key(trace.tx_local_index, trace.rx_local_index), std::move(ascan)); - } - if (ascans_by_pair.empty() || !(bandwidth_hz > 0.0)) { - return results; - } - - const auto grid = build_grid(selection, max_depth_m, kLegacyGridZMinM, imaging_plane_y_m); - if (grid.x_grid.empty() || grid.z_grid.empty()) { - return results; - } - - const double shell_sigma_m = velocity_mps / bandwidth_hz * 0.5; - const double snr_thresh = std::max(0.0, static_cast(live_config.gpr_snr_thresh)); - const double snr_comp_max = std::max(0.0, static_cast(live_config.gpr_snr_comp_max)); - const double comp_power = std::max(0.0, static_cast(live_config.gpr_comp_power)); - - std::unordered_map> peaks_by_pair{}; - for (std::size_t tx_index = 0U; tx_index < selection.x_tx.size(); ++tx_index) { - for (std::size_t rx_index = 0U; rx_index < selection.x_rx.size(); ++rx_index) { - const auto key = make_pair_key(static_cast(tx_index), static_cast(rx_index)); - const auto ascan_it = ascans_by_pair.find(key); - if (ascan_it == ascans_by_pair.end()) { - continue; - } - - const auto& ascan = ascan_it->second; - if (ascan.depth_m.size() < 3U || ascan.amplitude.size() < 3U) { - continue; - } - - const auto min_index = lower_bound_index(ascan.depth_m, min_depth_m); - const auto max_index = lower_bound_index(ascan.depth_m, max_depth_m); - if (max_index <= min_index + 2U || max_index > ascan.amplitude.size()) { - continue; - } - - const auto noise_begin = ascan.amplitude.begin() + static_cast(min_index); - const auto noise_end = ascan.amplitude.begin() + static_cast(max_index); - const double noise = median_copy(std::vector(noise_begin, noise_end)); - const double z_step = std::max(ascan.depth_m[1] - ascan.depth_m[0], 1e-6); - const std::size_t min_distance = static_cast(std::max( - 4.0, - std::floor(((velocity_mps / (2.0 * bandwidth_hz)) / z_step) * 0.7) - )); - const auto peak_indices = - find_legacy_peak_indices(ascan.amplitude, min_index, max_index, noise * snr_thresh, min_distance); - - const double attenuation_ref = legacy_attenuation_at_depth( - tx_index, rx_index, kLegacyAttenuationReferenceDepthM, selection, imaging_plane_y_m - ); - - auto& peaks = peaks_by_pair[key]; - peaks.reserve(peak_indices.size()); - for (const auto peak_index : peak_indices) { - const double z_app = ascan.depth_m[peak_index]; - const double snr_raw = ascan.amplitude[peak_index] / std::max(noise, 1e-12); - const double attenuation = - legacy_attenuation_at_depth(tx_index, rx_index, z_app, selection, imaging_plane_y_m); - const double snr_comp = std::min( - snr_raw / (std::pow(attenuation / attenuation_ref, comp_power) + 1e-12), - snr_comp_max - ); - peaks.push_back( - LegacyPeakRecord{ - .z_app = z_app, - .tau = ascan.time_s[peak_index], - .tau_corr = ascan.time_s[peak_index], - .z_corr = z_app, - .snr_raw = snr_raw, - .snr_comp = snr_comp, - } - ); - } - } - } - - if (peaks_by_pair.empty()) { - return results; - } - - if (live_config.legacy_gpr_mode == LegacyGprMode::Extended) { - const auto [regions, smoothed_accumulator] = extended_legacy_find_regions( - grid, - peaks_by_pair, - selection.x_tx, - selection.x_rx, - velocity_mps, - shell_sigma_m, - min_depth_m, - max_depth_m - ); - results.collection_payloads.push_back( - build_image_payload("gpr_accumulator", grid.x_grid, grid.z_grid, smoothed_accumulator) - ); - - std::vector> region_rows{}; - region_rows.reserve(regions.size()); - for (std::size_t index = 0U; index < regions.size(); ++index) { - const auto& region = regions[index]; - region_rows.push_back( - { - static_cast(region.x_m), - static_cast(region.z_m), - static_cast(region.score), - static_cast(region.pixel_count), - } - ); - results.collection_payloads.push_back( - build_image_payload( - "gpr_region_mask_" + std::to_string(index), - grid.x_grid, - grid.z_grid, - std::vector(region.mask.begin(), region.mask.end()) - ) - ); - } - results.collection_payloads.push_back(build_table_payload("gpr_region_centers", region_rows, 4U)); - return results; - } - - // Apply coarse per-Tx-event motion correction to the peak set produced - // above. After this step every peak carries both apparent and motion- - // corrected (tau, depth) values; the CLEAN search uses the corrected domain. - apply_legacy_motion_correction(peaks_by_pair, motion_timing, velocity_mps); - - const auto [points, smoothed_accumulator] = clean_legacy_find_points( - grid, - peaks_by_pair, - selection.x_tx, - selection.x_rx, - velocity_mps, - shell_sigma_m, - min_depth_m, - max_depth_m, - LegacyPeakDomain::Corrected - ); - results.collection_payloads.push_back( - build_image_payload("gpr_accumulator", grid.x_grid, grid.z_grid, smoothed_accumulator) - ); - - std::vector> point_rows{}; - point_rows.reserve(points.size()); - for (const auto& point : points) { - point_rows.push_back( - { - static_cast(point.x_m), - static_cast(point.z_m), - static_cast(point.score), - } - ); - } - results.collection_payloads.push_back(build_table_payload("gpr_points", point_rows, 3U)); - return results; -} diff --git a/data_acq_and_processing/processing/processors/src/gpr_processor.cpp b/data_acq_and_processing/processing/processors/src/gpr_processor.cpp index 43172ce..6882c4c 100644 --- a/data_acq_and_processing/processing/processors/src/gpr_processor.cpp +++ b/data_acq_and_processing/processing/processors/src/gpr_processor.cpp @@ -1,4 +1,5 @@ #include "gpr_processor.hpp" +#include "backprojection.hpp" #include #include @@ -22,7 +23,6 @@ namespace { // Algorithm bodies stay private to this translation unit; the public processor remains a small dispatcher. #include "gpr_backprojection_processor.ipp" -#include "gpr_legacy_processor.ipp" } // namespace @@ -39,17 +39,4 @@ auto GprProcessor::process_collection( return process_backprojection_gpr(run_config, collection, previous_collections, live_config, approach_filter_); } -auto LegacyGprProcessor::name() const -> std::string { - return "legacy_gpr"; -} - -auto LegacyGprProcessor::process_collection( - const config::RunConfig& run_config, - const ipc::PreprocessedCollection& collection, - std::span previous_collections, - const ProcessingLiveConfig& live_config -) -> ipc::ResultCollection { - return process_legacy_gpr(run_config, collection, previous_collections, live_config); -} - } // namespace radar::processing diff --git a/data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp b/data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp index 113ea6c..dec6acb 100644 --- a/data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp +++ b/data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp @@ -226,6 +226,34 @@ extern "C" __global__ void peak_planes( indices[y] = index; } +// Contiguous half-maximum interval through the joint peak on each scene axis. +// Interpolate the two crossings; use one voxel as the resolution floor. +extern "C" __global__ void peak_widths( + const float* volume, const double* gx, const double* gy, const double* gz, + const int nx, const int ny, const int nz, const int px, const int py, const int pz, + float* widths) +{ + const int a = (int)(blockIdx.x * blockDim.x + threadIdx.x); + if (a >= 3) return; + const int count[3] = {nx, ny, nz}, stride[3] = {1, nx*nz, nx}, peak[3] = {px, py, pz}; + const double* axis = a == 0 ? gx : a == 1 ? gy : gz; + const int origin = (py*nz + pz)*nx + px, k = peak[a], step = stride[a]; + const float half = .5f * volume[origin]; + int low = k, high = k; + while (low > 0 && volume[origin + (low-k-1)*step] >= half) --low; + while (high+1 < count[a] && volume[origin + (high-k+1)*step] >= half) ++high; + double left = axis[low], right = axis[high]; + if (low > 0) { + const float inside = volume[origin + (low-k)*step], outside = volume[origin + (low-k-1)*step]; + left -= (axis[low]-axis[low-1]) * (inside-half)/(inside-outside); + } + if (high+1 < count[a]) { + const float inside = volume[origin + (high-k)*step], outside = volume[origin + (high-k+1)*step]; + right += (axis[high+1]-axis[high]) * (inside-half)/(inside-outside); + } + widths[a] = (float)fmax(right-left, axis[1]-axis[0]); +} + // Orthographic max projection. The volume is non-negative, so float bit // patterns order like the values and an integer atomicMax does the job. The // pixel is found in double, so it does not hinge on float rounding. @@ -417,6 +445,7 @@ struct Gpu::Impl { CUfunction smooth = module.function("smooth"); CUfunction max_project = module.function("max_project"); CUfunction peak_planes = module.function("peak_planes"); + CUfunction peak_widths = module.function("peak_widths"); CUfunction project_oblique = module.function("project_oblique"); int smooth_radius = 0; DeviceBuffer taps; @@ -510,6 +539,13 @@ auto Gpu::Impl::project(const VolumeRequest& request) -> Projections { const auto plane_levels = download(levels, static_cast(ny)); const auto plane_indices = download(indices, static_cast(ny)); choose_peak(result, request, plane_levels, plane_indices); + if (result.level > 0.0F) { + const auto widths = gpu.peak_stats.reserve(3U * sizeof(float)); + launch(gpu.peak_widths, 3, volume, gx, gy, gz, nx, ny, nz, + static_cast(result.peak[2]), static_cast(result.peak[0]), static_cast(result.peak[1]), widths); + const auto measured = download(widths, 3U); + std::copy(measured.begin(), measured.end(), result.target_width.begin()); + } launch(gpu.max_project, front_size, volume, front, front_size, nx, plane, 1, nz, nx); launch(gpu.max_project, side_size, volume, side, side_size, nz, plane, nx, nx, 1); launch(gpu.max_project, top_size, volume, top, top_size, nx, nx, 1, ny, plane); diff --git a/data_acq_and_processing/processing/processors/src/sar3d_gpu.hpp b/data_acq_and_processing/processing/processors/src/sar3d_gpu.hpp index 7236a34..623c15c 100644 --- a/data_acq_and_processing/processing/processors/src/sar3d_gpu.hpp +++ b/data_acq_and_processing/processing/processors/src/sar3d_gpu.hpp @@ -10,7 +10,7 @@ // CUDA backend of the sar3d processor: A-scans of the band spectra, and the // backprojection of an aperture into the 3D volume, smoothed and projected. -// Only the projections and peak leave the backend. +// Only the projections, peak and its axial half-maximum widths leave the backend. namespace radar::processing::sar3d { class CudaUnavailable : public std::runtime_error { @@ -70,6 +70,7 @@ struct Projections { float level = 0.0F; float raw_level = 0.0F; std::array peak{}; // joint voxel maximum, [y, z, x] + std::array target_width{}; // contiguous FWHM through the peak, [X, Y, Z], metres }; // CUDA context and runtime-compiled kernels. A missing device raises diff --git a/data_acq_and_processing/processing/processors/src/sar3d_processor.cpp b/data_acq_and_processing/processing/processors/src/sar3d_processor.cpp index c26a6c3..09a7457 100644 --- a/data_acq_and_processing/processing/processors/src/sar3d_processor.cpp +++ b/data_acq_and_processing/processing/processors/src/sar3d_processor.cpp @@ -1,4 +1,5 @@ #include "sar3d_processor.hpp" +#include "backprojection.hpp" #include #include @@ -23,8 +24,8 @@ #include "sar3d_gpu.hpp" -// Port of the 3D SAR part of the reference notebook Multi_antenna_3D_clean_v1.ipynb -// (cells 0 and 4). Comments name the notebook functions and constants each step +// Port of the 3D SAR part of the reference notebook Multi_antenna_demonstrate_depthnorm.ipynb +// (cells 1 and 14). Comments name the notebook functions and constants each step // follows; the parameters keep the notebook names in the GUI. namespace radar::processing { namespace { @@ -32,10 +33,8 @@ namespace { using std::numbers::pi; // Notebook constants that are not live parameters. -constexpr double kWaveSpeed = 3e8; // v, eps_r = 1 +constexpr double kWaveSpeed = bp::kWaveSpeed; // v, eps_r = 1 constexpr std::size_t kOversample = 8U; // BP_OVERSAMPLE -constexpr double kPairNormPercentile = 50.0; // PAIR_NORM_PERCENTILE -constexpr double kPairNormEps = 1e-15; // PAIR_NORM_EPS constexpr double kCompRefDepth = 5.0; // COMP_REF_DEPTH constexpr float kCompRangeWeightMax = 5.0F; // COMP_RANGE_WEIGHT_MAX constexpr double kEdgeFrameWeight = 0.8; // SAR_EDGE_FRAME_WEIGHT @@ -45,9 +44,6 @@ constexpr double kNormEdgeMargin = 0.3; // SAR_NORM_EDGE_MARGIN_M // (x_ant.min() - 2 ... x_ant.max() + 2), Y is the height above the ground // (SAR_Y_RANGE_M), Z along the ground follows the depth gate. constexpr double kSceneXMargin = 2.0; -constexpr double kSceneYMin = -0.5; -constexpr double kSceneYMax = 2.0; -static_assert(kSceneYMax - kSceneYMin > 2.5 * kNormEdgeMargin, "else sar_norm_reference_mask keeps every plane"); // The 3D view looks along the camera of the notebook's 3D scene // (sar_scene_layout: Plotly eye (1.7, -2.2, 0.85) in its axes X, Z, Y). constexpr std::array kCameraEyeXzy{1.7, -2.2, 0.85}; @@ -71,10 +67,17 @@ struct Settings { bool use_compensation = true; // SAR_USE_COMPENSATION bool normalize_each_frame = true; // SAR_NORMALIZE_EACH_FRAME + bool pair_depth_equalize = true; // PAIR_DEPTH_EQUALIZE + bool subtract_mean = true; + std::size_t mean_count = 10U; + double x_size_m = 0.0, z_size_m = 0.0; // 0: notebook's automatic extent + double y_min_m = -0.5, y_max_m = 2.0; // SAR_Y_RANGE_M + // Everything the cached A-scans depend on. [[nodiscard]] auto ascan_key() const { return std::make_tuple(inputs, outputs, f_start_hz, f_stop_hz, min_depth_m, max_depth_m, speed_m_s, - look_angle_rad, normalize_each_frame); + look_angle_rad, normalize_each_frame, pair_depth_equalize, subtract_mean, mean_count, + tx, rx, x_size_m); } }; using AscanKey = decltype(std::declval().ascan_key()); @@ -127,6 +130,16 @@ void select_antennas( settings.comp_range_power = live.sar3d_comp_range_power; settings.use_compensation = live.sar3d_use_compensation; settings.normalize_each_frame = live.sar3d_normalize_each_frame; + settings.pair_depth_equalize = live.pair_depth_equalize; + settings.subtract_mean = live.gpr_background_subtract_enabled; + settings.mean_count = live.gpr_background_mean_count; + settings.x_size_m = live.sar3d_x_size_m; + settings.z_size_m = live.sar3d_z_size_m; + settings.y_min_m = live.sar3d_y_min_m; + settings.y_max_m = live.sar3d_y_max_m; + if (settings.x_size_m < 0.0 || settings.z_size_m < 0.0 || !(settings.y_max_m > settings.y_min_m)) { + throw std::runtime_error("SAR scene sizes must be non-negative and SAR_Y_RANGE_M must increase"); + } if (!(settings.f_stop_hz > settings.f_start_hz)) { throw std::runtime_error("F_STOP must be above F_START"); } @@ -152,73 +165,6 @@ void select_antennas( return *lower_it + (position - static_cast(lower)) * (upper - *lower_it); } -// When the samples of a trace were measured, relative to the acquisition centre -// of its collection (compute_pair_timestamps_for_snapshot and -// compute_frequency_sample_times). Traces sharing a capture window come from -// one sweep: at every frequency point its stimulus stages - the outputs, which -// are numbered in device stage order - take turns while all inputs receive. -struct TraceClock { - double start_s = 0.0; // sweep start - double event_s = 0.0; // sweep duration / (stages * points) - std::size_t stage = 0U; - std::size_t stages = 1U; - - [[nodiscard]] auto sample_time_s(std::size_t point) const -> double { - return start_s + (static_cast(stages) * point + stage) * event_s; - } -}; - -struct CollectionTiming { - std::uint64_t center_ns = 0U; // centre of the collection's capture window - std::vector clocks; // per collection trace -}; - -// Timing from the per-trace capture windows (snapshot_acquisition_timing); every -// trace takes part, so the stages do not depend on the INPUT_IDX/OUTPUT_IDX choice. -[[nodiscard]] auto collection_timing(const ipc::PreprocessedCollection& collection) -> CollectionTiming { - if (collection.traces.empty()) { - throw std::runtime_error("collection has no traces"); - } - std::map, std::vector> sweep_outputs{}; - std::uint64_t first = std::numeric_limits::max(); - std::uint64_t last = 0U; - for (const auto& trace : collection.traces) { - if (trace.capture_end_ns <= trace.capture_start_ns) { - throw std::runtime_error("collection has traces without a capture window"); - } - if (trace.frequency_hz.empty()) { - throw std::runtime_error("collection has a trace without frequency points"); - } - sweep_outputs[{trace.capture_start_ns, trace.capture_end_ns}].push_back(trace.combo.output_pos); - first = std::min(first, trace.capture_start_ns); - last = std::max(last, trace.capture_end_ns); - } - for (auto& [window, outputs] : sweep_outputs) { - std::sort(outputs.begin(), outputs.end()); - outputs.erase(std::unique(outputs.begin(), outputs.end()), outputs.end()); - } - - CollectionTiming timing{}; - const auto duration_ns = last - first; - timing.center_ns = first + duration_ns / 2U; - for (const auto& trace : collection.traces) { - const auto& stages = sweep_outputs.at({trace.capture_start_ns, trace.capture_end_ns}); - const auto stage = std::find(stages.begin(), stages.end(), trace.combo.output_pos) - stages.begin(); - TraceClock clock{}; - clock.start_s = 1e-9 * (static_cast(trace.capture_start_ns - first) - - 0.5 * static_cast(duration_ns)); - clock.stage = stage; - clock.stages = stages.size(); - clock.event_s = 1e-9 * static_cast(trace.capture_end_ns - trace.capture_start_ns) / - static_cast(clock.stages) / trace.frequency_hz.size(); - timing.clocks.push_back(clock); - } - return timing; -} - -// Band of one trace, ready for the IFFT (apply_intra_sweep_phase_correction and -// compute_ascan_bp up to the IFFT): the radar motion removed as a phase, then -// the Hann window over the band. struct Band { std::size_t first_bin = 0U; // k0, the IFFT bin of the first band point std::size_t size = 0U; // IFFT length: 2(k0 + n - 1) rounded up to a power of two, x BP_OVERSAMPLE @@ -226,7 +172,8 @@ struct Band { std::vector> spectrum{}; }; -[[nodiscard]] auto trace_band(const ipc::SweepTraceBlock& trace, const TraceClock& clock, const Settings& settings) +[[nodiscard]] auto trace_band(const ipc::SweepTraceBlock& trace, const bp::TraceClock& clock, const Settings& settings, + const bp::Background& background) -> Band { const auto& frequency = trace.frequency_hz; if (trace.s21.size() != frequency.size()) { @@ -269,47 +216,24 @@ struct Band { // from the collection centre at the sample time becomes the phase // exp(-j 4 pi f offset / v): the pair and the intra-sweep corrections in one // (MOTION_PHASE_SIGN = -1, direction +1). - const double phase_per_hz_s = -4.0 * pi * settings.speed_m_s * std::cos(settings.look_angle_rad) / kWaveSpeed; + auto samples = bp::subtract_background(trace, background); + bp::apply_motion_phase(samples, frequency, clock, settings.speed_m_s, settings.look_angle_rad, + 1.0, false, settings.f_start_hz, settings.f_stop_hz); const auto last = static_cast(points.size() - 1U); for (std::size_t index = 0U; index < points.size(); ++index) { const std::size_t point = points[index]; const double hann = 0.5 - 0.5 * std::cos(2.0 * pi * static_cast(index) / last); - const std::complex sample(trace.s21[point].re, trace.s21[point].im); + const auto sample = samples[point]; if (!std::isfinite(sample.real()) || !std::isfinite(sample.imag())) { throw std::runtime_error("S21 contains a non-finite sample"); } band.spectrum.emplace_back( - sample * std::polar(hann, phase_per_hz_s * frequency[point] * clock.sample_time_s(point)) + sample * hann ); } return band; } -// SAR_NORMALIZE_EACH_FRAME: divide a channel by the PAIR_NORM_PERCENTILE of |h| -// inside the depth gate (sar_pair_ascans). -void normalize_ascan(std::span> ascan, double dt_s, const Settings& settings) { - std::vector magnitudes{}; - for (std::size_t index = 0U; index < ascan.size(); ++index) { - const double depth = static_cast(index) * dt_s * kWaveSpeed / 2.0; - if (depth >= settings.min_depth_m && depth <= settings.max_depth_m) { - magnitudes.push_back(std::abs(std::complex(ascan[index]))); - } - } - if (magnitudes.empty()) { - for (const auto& sample : ascan) { - magnitudes.push_back(std::abs(std::complex(sample))); - } - } - double scale = percentile(std::move(magnitudes), kPairNormPercentile); - if (!std::isfinite(scale) || scale <= kPairNormEps) { - scale = 1.0; - } - const auto divisor = static_cast(scale + kPairNormEps); - for (auto& sample : ascan) { - sample /= divisor; - } -} - // A-scans of every selected Tx/Rx pair of one collection, (Tx, Rx) order. struct Frame { std::uint64_t center_ns = 0U; @@ -321,9 +245,11 @@ struct Frame { [[nodiscard]] auto compute_frame( const ipc::PreprocessedCollection& collection, const Settings& settings, + std::span previous, sar3d::Gpu& backend ) -> Frame { - const auto timing = collection_timing(collection); + const auto timing = bp::collection_timing(collection); + const auto background = bp::background_mean(previous, settings.subtract_mean, settings.mean_count); std::vector bands{}; for (const auto output : settings.outputs) { for (const auto input : settings.inputs) { @@ -337,7 +263,7 @@ struct Frame { ); } const auto& clock = timing.clocks[static_cast(trace - traces.begin())]; - bands.push_back(trace_band(*trace, clock, settings)); + bands.push_back(trace_band(*trace, clock, settings, background)); } } const Band& first = bands.front(); @@ -354,10 +280,27 @@ struct Frame { frame.samples = first.size; frame.dt_s = 1.0 / (static_cast(first.size) * first.df_hz); frame.ascans = backend.inverse_dft(spectra, first.spectrum.size(), first.first_bin, first.size); - if (settings.normalize_each_frame) { - for (std::size_t offset = 0U; offset < frame.ascans.size(); offset += frame.samples) { - normalize_ascan(std::span(frame.ascans).subspan(offset, frame.samples), frame.dt_s, settings); + std::vector>> pairs; + for (std::size_t offset = 0; offset < frame.ascans.size(); offset += frame.samples) { + pairs.push_back(std::span(frame.ascans).subspan(offset, frame.samples)); + } + const double dz = frame.dt_s * kWaveSpeed / 2.0; + if (settings.pair_depth_equalize) { + double x_min = settings.tx.front()[0], x_max = x_min, max_abs_x = 0.0; + for (const auto* antennas : {&settings.tx, &settings.rx}) { + for (const auto& antenna : *antennas) { + x_min = std::min(x_min, antenna[0]); + x_max = std::max(x_max, antenna[0]); + max_abs_x = std::max(max_abs_x, std::abs(antenna[0])); + } } + const double half_size = settings.x_size_m > 0.0 ? settings.x_size_m / 2.0 : (x_max - x_min) / 2.0 + kSceneXMargin; + const double center = (x_min + x_max) / 2.0; + bp::equalize_pairs_by_depth(pairs, dz, + bp::equalize_depth_limit(settings.max_depth_m, center - half_size, center + half_size, max_abs_x)); + } + if (settings.normalize_each_frame) { + for (auto pair : pairs) { bp::normalize_pair(pair, dz, settings.min_depth_m, settings.max_depth_m); } } return frame; } @@ -406,7 +349,7 @@ struct Grid { } double z_min = std::numeric_limits::infinity(); double z_max = -std::numeric_limits::infinity(); - for (const double height : {kSceneYMin, kSceneYMax}) { + for (const double height : {settings.y_min_m, settings.y_max_m}) { const double depth_at_zero = -(height - settings.ground_depth_m) * sin_p; for (const double depth : {settings.min_depth_m, settings.max_depth_m}) { const double z = (depth - depth_at_zero) / cos_p; @@ -414,18 +357,23 @@ struct Grid { z_max = std::max(z_max, z); } } + const auto resize = [](double& low, double& high, double size) { + if (size > 0.0) { const double center = (low + high) / 2.0; low = center - size / 2.0; high = center + size / 2.0; } + }; + resize(x_min, x_max, settings.x_size_m); + resize(z_min, z_max, settings.z_size_m); // Match sar_build_full_region_grid and reject oversized products before // allocating the axes. Every downstream kernel uses signed 32-bit indices. const auto axis_size = [&](double span) { return checked_axis_size(std::max(std::nearbyint(span / settings.voxel_m) + 1.0, 2.0)); }; const auto nx = axis_size(x_max - x_min); - const auto ny = axis_size(kSceneYMax - kSceneYMin); + const auto ny = axis_size(settings.y_max_m - settings.y_min_m); const auto nz = axis_size(z_max - z_min); if (nx > static_cast(std::numeric_limits::max()) / ny / nz) { throw std::runtime_error("SAR grid exceeds 32-bit kernel indexing; increase SAR_VOXEL_M"); } - return {linspace(x_min, x_max, nx), linspace(kSceneYMin, kSceneYMax, ny), + return {linspace(x_min, x_max, nx), linspace(settings.y_min_m, settings.y_max_m, ny), linspace(z_min, z_max, nz)}; } @@ -584,6 +532,7 @@ struct Sar3dProcessor::State { } struct CachedFrame { + std::vector background; std::shared_ptr value; std::string error; }; @@ -594,21 +543,28 @@ struct Sar3dProcessor::State { std::string last_error{}; std::uint64_t repeated_errors = 0U; - auto frame(const ipc::PreprocessedCollection& collection, const Settings& settings, sar3d::Gpu& backend) + auto frame(const ipc::PreprocessedCollection& collection, const Settings& settings, + std::span previous, sar3d::Gpu& backend) -> std::shared_ptr { - if (const auto found = frames.find(sweep_key(collection)); found != frames.end()) { + std::vector background; + if (settings.subtract_mean) { + for (const auto& member : previous.last(std::min(previous.size(), settings.mean_count))) { + background.push_back(sweep_key(member)); + } + } + if (const auto found = frames.find(sweep_key(collection)); found != frames.end() && found->second.background == background) { if (found->second.value == nullptr) { report_error(found->second.error); } return found->second.value; } std::shared_ptr frame{}; std::string failure; try { - frame = std::make_shared(compute_frame(collection, settings, backend)); + frame = std::make_shared(compute_frame(collection, settings, previous, backend)); } catch (const std::exception& error) { failure = error.what(); report_error(failure); } - frames.emplace(sweep_key(collection), CachedFrame{frame, std::move(failure)}); + frames.insert_or_assign(sweep_key(collection), CachedFrame{std::move(background), frame, std::move(failure)}); return frame; } @@ -631,22 +587,20 @@ struct Sar3dProcessor::State { frames.clear(); frames_key = settings.ascan_key(); } - const auto recent = previous.last(std::min(previous.size(), 2U * settings.frame_radius)); + const auto recent = previous.last(std::min(previous.size(), settings.frame_radius)); std::set retained{sweep_key(collection)}; for (const auto& member : recent) { retained.insert(sweep_key(member)); } // Prune even when the newest frame is invalid, so a bad input stream // cannot grow the failed-frame cache forever. std::erase_if(frames, [&](const auto& entry) { return !retained.contains(entry.first); }); - const auto newest = frame(collection, settings, backend); + const auto newest = frame(collection, settings, previous, backend); if (newest == nullptr) { return {status_payload(2.0F, last_error)}; } - // Aperture: the newest 2r+1 collections, centred r collections before - // the newest one, so the image trails the input by r frames. Collections - // that left the aperture are dropped from the cache. + // SAR_CAUSAL_ONLY: current collection and r previous ones, current scene origin. std::vector> aperture{}; - for (const auto& member : recent) { - auto member_frame = frame(member, settings, backend); + for (std::size_t i = previous.size() - recent.size(); i < previous.size(); ++i) { + auto member_frame = frame(previous[i], settings, previous.first(i), backend); const bool same_axis = member_frame != nullptr && member_frame->samples == newest->samples && member_frame->dt_s == newest->dt_s; if (same_axis) { @@ -654,9 +608,7 @@ struct Sar3dProcessor::State { } } aperture.push_back(newest); - const std::size_t centre_index = - aperture.size() > settings.frame_radius ? aperture.size() - 1U - settings.frame_radius : 0U; - const Frame& centre = *aperture[centre_index]; + const Frame& centre = *newest; std::vector offsets_s{}; for (const auto& member : aperture) { @@ -761,6 +713,8 @@ struct Sar3dProcessor::State { basis.insert(basis.end(), oblique.view.v.begin(), oblique.view.v.end()); payloads.push_back(table_payload("sar3d_oblique_basis", 3U, std::move(basis))); if (level > 0.0F) { + payloads.push_back(table_payload("sar3d_target_width", 3U, + {images.target_width[0], images.target_width[1], images.target_width[2]})); payloads.push_back(table_payload("sar3d_peak", 4U, {static_cast(normalization), static_cast(grid.x[peak_x]), static_cast(grid.y[peak_y]), static_cast(grid.z[peak_z])})); diff --git a/data_acq_and_processing/tests/backprojection_fixture.cpp b/data_acq_and_processing/tests/backprojection_fixture.cpp new file mode 100644 index 0000000..2d5d8fa --- /dev/null +++ b/data_acq_and_processing/tests/backprojection_fixture.cpp @@ -0,0 +1,113 @@ +#include "sar3d_processor.hpp" +#include "gpr_processor.hpp" + +#include +#include +#include +#include +#include + +using Json = nlohmann::json; +using namespace radar; + +namespace { +auto complex_samples(const Json& samples) -> std::vector> { + std::vector> values; + for (const auto& sample : samples) { + values.emplace_back(sample.at(0).get(), sample.at(1).get()); + } + return values; +} + +auto collection(const Json& source) -> ipc::PreprocessedCollection { + ipc::PreprocessedCollection result; + result.collection_id = source.at("id"); + result.monotonic_ns = source.at("timestamp"); + for (const auto& entry : source.at("traces")) { + ipc::SweepTraceBlock trace; + trace.combo = {entry.at("input"), entry.at("output")}; + trace.capture_start_ns = entry.at("start"); + trace.capture_end_ns = entry.at("end"); + trace.frequency_hz = entry.at("frequency").get>(); + for (const auto& sample : complex_samples(entry.at("s21"))) { + trace.s21.push_back({sample.real(), sample.imag()}); + } + if (result.traces.empty() || trace.capture_start_ns < result.capture_start_ns) { + result.capture_start_ns = trace.capture_start_ns; + } + result.capture_end_ns = std::max(result.capture_end_ns, trace.capture_end_ns); + result.traces.push_back(std::move(trace)); + } + return result; +} +} + +// Runs GPR or sar3d on one request or a batch, retaining processor caches. Each request holds the antenna +// geometry, the live settings and the latest and previous collections; the +// published payloads go to stdout as JSON. +int main() { + try { + const Json source = Json::parse(std::cin); + const bool batch = source.is_array(); + const Json requests = batch ? source : Json::array({source}); + Json results = Json::array(); + processing::Sar3dProcessor sar; + processing::GprProcessor gpr; + for (const auto& input : requests) { + config::RunConfig config; + config.rings.results.slot_size_bytes = input.value("slot", config.rings.results.slot_size_bytes); + for (const auto& row : input.at("tx")) { + config.gpr.tx_geometry.push_back({row.at(0), row.at(1), row.at(2), row.at(3)}); + } + for (const auto& row : input.at("rx")) { + config.gpr.rx_geometry.push_back({row.at(0), row.at(1), row.at(2), row.at(3)}); + } + processing::ProcessingLiveConfig live; + live.gpr_start_freq_mhz = input.at("f_start").get() / 1e6; + live.gpr_stop_freq_mhz = input.at("f_stop").get() / 1e6; + live.gpr_min_depth_m = input.at("min_depth"); + live.gpr_max_depth_m = input.at("max_depth"); + live.gpr_speed_m_s = input.at("speed"); + live.gpr_look_angle_deg = input.at("pitch"); + live.sar3d_voxel_m = input.at("voxel"); + live.sar3d_bp_ground_depth_m = input.at("ground"); + live.sar3d_frame_radius = input.at("radius"); + live.sar3d_comp_range_power = input.at("power"); + live.sar3d_use_compensation = input.at("compensation"); + live.sar3d_normalize_each_frame = input.at("normalize"); + live.pair_depth_equalize = input.value("equalize", true); + live.gpr_background_subtract_enabled = input.value("subtract_mean", false); + live.gpr_background_mean_count = input.value("mean_count", 10U); + live.sar3d_x_size_m = input.value("x_size", 0.0F); + live.sar3d_y_min_m = input.value("y_min", -0.5F); + live.sar3d_y_max_m = input.value("y_max", 2.0F); + live.sar3d_z_size_m = input.value("z_size", 0.0F); + live.gpr_bp_ground_depth_m = input.at("ground"); + live.gpr_range_comp_power = input.at("power"); + live.gpr_motion_mode = input.value("motion", "int_minus"); + live.gpr_object_approach_min_frames = 1; + live.gpr_remove_sidelobe_objects_enabled = false; + std::vector previous; + for (const auto& member : input.at("previous")) { + previous.push_back(collection(member)); + } + const auto latest = collection(input.at("latest")); + for (const auto& trace : latest.traces) { config.run_combos.push_back(trace.combo); } + processing::ProcessorInterface& processor = input.value("mode", "sar3d") == "gpr" + ? static_cast(gpr) : static_cast(sar); + const auto result = processor.process_collection(config, latest, previous, live); + Json payloads = Json::array(); + for (const auto& payload : result.collection_payloads) { + payloads.push_back({{"name", payload.processing_name}, + {"x", payload.image_x_axis}, {"y", payload.image_y_axis}, + {"image", payload.image_values}, {"table", payload.table_values}, + {"columns", payload.table_columns}}); + } + results.push_back(Json{{"payloads", payloads}}); + } + std::cout << (batch ? results : results.at(0)) << '\n'; + } catch (const std::exception& error) { + std::cout << Json{{"error", error.what()}} << '\n'; + return 1; + } +} diff --git a/data_acq_and_processing/tests/generate_sar3d_reference.py b/data_acq_and_processing/tests/generate_backprojection_reference.py similarity index 52% rename from data_acq_and_processing/tests/generate_sar3d_reference.py rename to data_acq_and_processing/tests/generate_backprojection_reference.py index 76ee957..db90345 100644 --- a/data_acq_and_processing/tests/generate_sar3d_reference.py +++ b/data_acq_and_processing/tests/generate_backprojection_reference.py @@ -1,12 +1,12 @@ -"""Generate small SAR3D references by executing notebook cells 0 and 4. +"""Generate SAR3D/GPR references from notebook cells 1/14 and Horns coherent BP. -Requires the original, untracked Multi_antenna_3D_clean_v1.ipynb; a missing notebook +Requires the original, untracked Multi_antenna_demonstrate_depthnorm.ipynb; a missing notebook is an error. NumPy/SciPy execute its function/class ASTs unchanged. Only parameter values are replaced; plotting imports and the interactive top-level runs are excluded. No production SAR implementation or hand-transcribed oracle is used. Run from the repository root: - .venv/bin/python data_acq_and_processing/tests/generate_sar3d_reference.py + .venv/bin/python data_acq_and_processing/tests/generate_backprojection_reference.py """ from __future__ import annotations @@ -24,14 +24,14 @@ ROOT = Path(__file__).resolve().parents[2] REFERENCE = Path(__file__).with_name("sar3d_reference.npz") -def notebook_functions(notebook: Path, parameters: dict) -> tuple[dict, dict]: - cells = json.loads(notebook.read_text(encoding="utf-8"))["cells"] +def reference_functions(source_path: Path, parameters: dict, cells: tuple[int, ...] = (1, 14)) -> tuple[dict, dict]: + contents = source_path.read_text(encoding="utf-8") + sources = [(str(i), "".join(json.loads(contents)["cells"][i]["source"])) for i in cells] if cells else [("script", contents)] namespace = {"__name__": __name__} hashes = {} - for index in (0, 4): - source = "".join(cells[index]["source"]) + for index, source in sources: hashes[str(index)] = hashlib.sha256(source.encode()).hexdigest() - tree = ast.parse(source, filename=f"{notebook.name}:cell{index}") + tree = ast.parse(source, filename=f"{source_path.name}:cell{index}") definitions = [] for node in tree.body: if isinstance(node, (ast.FunctionDef, ast.ClassDef, ast.ImportFrom)): @@ -45,23 +45,24 @@ def notebook_functions(notebook: Path, parameters: dict) -> tuple[dict, dict]: except (ValueError, TypeError): pass # Dataset loading and the interactive run are supplied below. namespace.update(parameters) - exec(compile(ast.Module(definitions, type_ignores=[]), f"{notebook.name}:cell{index}", + exec(compile(ast.Module(definitions, type_ignores=[]), f"{source_path.name}:cell{index}", "exec", dont_inherit=True), namespace) return namespace, hashes -def scene() -> dict: +def scene(f_start=1e9, f_stop=4e9, points=201) -> dict: """Deterministic two-echo input with motion, vertical baselines and four pairs.""" request = dict( tx=[[0, -.4, 0, 0], [1, .45, .2, .05]], - rx=[[0, -.2, .1, 0], [1, .3, .35, .02]], f_start=1e9, f_stop=4e9, + rx=[[0, -.2, .1, 0], [1, .3, .35, .02]], f_start=f_start, f_stop=f_stop, min_depth=1., max_depth=3., speed=.4, pitch=13., ground=.3, - voxel=.4, radius=1, power=1.5, compensation=True, normalize=True, + voxel=.4, radius=2, power=1.5, compensation=True, normalize=True, + equalize=True, subtract_mean=False, mean_count=1, motion="int_minus", ) # Match the float32 settings/geometry accepted by the native fixture. for key in ("tx", "rx"): request[key] = np.asarray(request[key], dtype=np.float32).astype(float).tolist() - frequency = np.linspace(1e9, 4e9, 201, dtype=np.float32) + frequency = np.linspace(f_start, f_stop, points, dtype=np.float32) members = [] for frame, offset in enumerate((0, 180_000_000, 510_000_000)): traces = [] @@ -91,14 +92,17 @@ def write_dataset(directory: Path, request: dict) -> list[Path]: ] (directory / "config_profile.json").write_text(json.dumps(profile), encoding="utf-8") snapshots = [] - for index, member in enumerate([*request["previous"], request["latest"]]): + members = [*request["previous"], request["latest"]] + for index, member in enumerate(members): + background = previous_mean({**request, "previous": members[:index]}) snapshot = directory / f'{index:04d}_id{member["id"]}_ns{member["timestamp"]}' snapshot.mkdir() metadata = {"traces": []} for trace in member["traces"]: stem = f'i{trace["input"]}_o{trace["output"]}' - values = np.asarray(trace["s21"], dtype=np.float32) - np.save(snapshot / f"{stem}_s21.npy", values[:, 0] + 1j * values[:, 1]) + values = np.asarray(trace["s21"]) + samples = values[:, 0] + 1j * values[:, 1] - background.get((trace["output"], trace["input"]), 0.) + np.save(snapshot / f"{stem}_s21.npy", samples) np.save(snapshot / f"{stem}_freq.npy", np.asarray(trace["frequency"], dtype=np.float32)) metadata["traces"].append(dict(input=trace["input"], output=trace["output"], capture_start_ns=trace["start"], capture_end_ns=trace["end"])) @@ -109,7 +113,7 @@ def write_dataset(directory: Path, request: dict) -> list[Path]: def reference(notebook: Path, request: dict, directory: Path) -> tuple[dict, dict]: snapshots = write_dataset(directory, request) - center = snapshots[max(len(snapshots) - 1 - request["radius"], 0)] + center = snapshots[-1] parameters = dict( INPUT_IDX=[int(row[0]) for row in request["rx"]], OUTPUT_IDX=[int(row[0]) for row in request["tx"]], @@ -120,9 +124,10 @@ def reference(notebook: Path, request: dict, directory: Path) -> tuple[dict, dic BP_GROUND_DEPTH_M=request["ground"], COMP_RANGE_POWER=request["power"], SAR_VOXEL_M=(request["voxel"],) * 3, SAR_FRAME_RADIUS=request["radius"], SAR_USE_COMPENSATION=request["compensation"], SAR_NORMALIZE_EACH_FRAME=request["normalize"], + PAIR_DEPTH_EQUALIZE=request["equalize"], SAR_SCAN_ROOT=directory, DATA_PATH=center, BG_PATH=directory, v=3e8, ) - ns, hashes = notebook_functions(notebook, parameters) + ns, hashes = reference_functions(notebook, parameters) tx, rx, order, tx_order, pair_count = ns["load_antenna_config_profile"]( directory / "config_profile.json", ns["INPUT_IDX"], ns["OUTPUT_IDX"]) ns["N_tx"], ns["N_rx"] = len(tx), len(rx) @@ -130,7 +135,9 @@ def reference(notebook: Path, request: dict, directory: Path) -> tuple[dict, dic for axis, values in zip("xyz", ns["positions_to_xyz"](positions).T): ns[f"{axis}_{kind}"] = values antenna_x = np.concatenate((ns["x_tx"], ns["x_rx"])) - ns["x_grid_bp"] = np.linspace(antenna_x.min() - 2, antenna_x.max() + 2, ns["NX_BP"]) + ns["x_ant"] = antenna_x + ns["x_min"], ns["x_max"] = antenna_x.min() - 2, antenna_x.max() + 2 + ns["x_grid_bp"] = np.linspace(ns["x_min"], ns["x_max"], ns["NX_BP"]) sweep_time, _, _ = ns["estimate_sweep_time_from_meta"](directory, n_pairs=pair_count) config = ns["MotionConfig"]( speed_m_s=ns["SPEED_M_S"], look_angle_deg=ns["LOOK_ANGLE_DEG"], @@ -143,23 +150,118 @@ def reference(notebook: Path, request: dict, directory: Path) -> tuple[dict, dic motion_phase_sign=ns["MOTION_PHASE_SIGN"], ) timestamps, _ = ns["compute_pair_timestamps_for_snapshot"](center, config) - result = ns["sar_compute_volume"](center, None, config, timestamps, verbose=False) + grid = ns["sar_build_full_region_grid"](config, y_range_m=(request.get("y_min", -.5), request.get("y_max", 2.))) + for axis in ("x", "z"): + size = request.get(f"{axis}_size", 0.) + if size: + center_axis = .5 * (grid[axis][0] + grid[axis][-1]) + grid[axis] = np.linspace(center_axis - size/2, center_axis + size/2, max(round(size/request["voxel"])+1, 2)) + grid["XX"], grid["ZZ"] = np.meshgrid(grid["x"], grid["z"]) + ns["x_min"], ns["x_max"] = grid["x"][0], grid["x"][-1] + result = ns["sar_compute_volume"](center, None, config, timestamps, grid=grid, verbose=False) volume, grid = result["volume"], result["grid"] iy, iz, ix = ns["sar_global_max_index"](volume, grid) + from scipy.signal import peak_widths + widths = [] + for axis, peak_index, profile in (("x", ix, volume[iy, iz, :]), ("y", iy, volume[:, iz, ix]), ("z", iz, volume[iy, :, ix])): + profile = np.asarray(profile, dtype=np.float64) + width, *_ = peak_widths(profile, [peak_index], prominence_data=(np.array([profile[peak_index]]), + np.array([0], dtype=np.intp), np.array([len(profile)-1], dtype=np.intp))) + widths.append(max(float(width[0]), 1.) * (grid[axis][1]-grid[axis][0])) expected = {axis: grid[axis] for axis in "xyz"} + expected["target_width"] = np.asarray(widths) expected.update(front=volume.max(axis=1), side=volume.max(axis=2), top=volume.max(axis=0), peak=np.array([result["level"], grid["x"][ix], grid["y"][iy], grid["z"][iz]])) return expected, hashes +def previous_mean(request: dict) -> dict: + samples = {} + if request["subtract_mean"] and request["mean_count"]: + for member in request["previous"][-request["mean_count"]:]: + for trace in member["traces"]: + values = np.asarray(trace["s21"]) + samples.setdefault((trace["output"], trace["input"]), []).append(values[:, 0] + 1j * values[:, 1]) + return {key: np.mean(values, axis=0) for key, values in samples.items()} + + +def gpr_reference(notebook: Path, script: Path, request: dict) -> tuple[np.ndarray, dict]: + """Execute Horns coherent BP unchanged, adding the notebook's equalizer before its scalar norm.""" + parameters = dict(F_START=request["f_start"], F_STOP=request["f_stop"], MIN_DEPTH=request["min_depth"], + MAX_DEPTH=request["max_depth"], LOOK_ANGLE_DEG=request["pitch"], + BP_GROUND_DEPTH_M=request["ground"], COMP_RANGE_POWER=request["power"], v=3e8) + ns, hashes = reference_functions(script, parameters, ()) + ns["N_tx"], ns["N_rx"] = len(request["tx"]), len(request["rx"]) + ns["OUTPUT_IDX"], ns["INPUT_IDX"] = [int(row[0]) for row in request["tx"]], [int(row[0]) for row in request["rx"]] + for kind in ("tx", "rx"): + for axis, values in zip("xyz", np.asarray(request[kind])[:, 1:].T): + ns[f"{axis}_{kind}"] = values + xs = np.concatenate((ns["x_tx"], ns["x_rx"])) + ns["x_ant"] = xs + ns["x_min"], ns["x_max"] = xs.min()-2, xs.max()+2 + ns["x_grid_bp"] = np.linspace(ns["x_min"], ns["x_max"], ns["NX_BP"]) + ns["z_grid_bp"] = np.linspace(.1, request["max_depth"], ns["NZ_BP"]) + ns["XX_bp"], ns["ZZ_bp"] = np.meshgrid(ns["x_grid_bp"], ns["z_grid_bp"]) + ns["depth_gate"] = (ns["ZZ_bp"] >= request["min_depth"]) & (ns["ZZ_bp"] <= request["max_depth"]) + member = request["latest"] + duration = (max(t["end"] for t in member["traces"]) - min(t["start"] for t in member["traces"])) * 1e-9 + config = ns["MotionConfig"](speed_m_s=request["speed"], look_angle_deg=request["pitch"], + tx_per_segment=len(request["tx"]), tx_order_full=tuple(ns["OUTPUT_IDX"]), + frame_acq_time_s=duration, n_freq_hint=len(member["traces"][0]["frequency"]), + intra_phase_focus_only=request["motion"] == "int_focus", motion_phase_sign=-1., + pair_order_phys=[(t["output"], t["input"]) for t in member["traces"]]) + clocks, _ = ns["compute_pair_timestamps"](config, ns["INPUT_IDX"], ns["OUTPUT_IDX"]) + mean = previous_mean(request) + H, T = {}, {} + for trace in member["traces"]: + key = (trace["output"], trace["input"]) + values = np.asarray(trace["s21"], dtype=np.float32) + samples = values[:, 0].astype(float) + 1j * values[:, 1].astype(float) + samples -= mean.get(key, 0.) + frequency = np.asarray(trace["frequency"], dtype=np.float32) + samples, _ = ns["apply_intra_sweep_phase_correction"](samples, frequency, clocks[key], config, + ns["v"], ns["F_START"], ns["F_STOP"]) + t, _, h, _, _ = ns["compute_ascan_bp"](samples, frequency, ns["F_START"], ns["F_STOP"]) + H[key], T[key] = h, t + if request["equalize"]: + eq, notebook_hashes = reference_functions(notebook, {**parameters, "PAIR_DEPTH_EQUALIZE": True}) + hashes.update({f"notebook_cell_{cell}": value for cell, value in notebook_hashes.items()}) + eq.update(x_ant=xs, x_min=ns["x_min"], x_max=ns["x_max"]) + H = eq["equalize_pairs_by_depth"](H, T) + for key, h in H.items(): + depth = T[key] * ns["v"] / 2 + mask = (depth >= request["min_depth"]) & (depth <= request["max_depth"]) + scale = np.percentile(np.abs(h)[mask] if mask.any() else np.abs(h), 50.) + H[key] = h / ((scale if np.isfinite(scale) and scale > 1e-15 else 1.) + 1e-15) + raw, _ = ns["backproject_coherent"](H, T) + return ns["normalize_bp_map"](raw).astype(np.float32), hashes + + +def generate_gpr(notebook: Path, script: Path, output: Path) -> None: + cases = [dict(equalize=False), dict(equalize=True), + dict(equalize=True, subtract_mean=True, mean_count=1, speed=-.2, pitch=-7., motion="int_focus", power=2.)] + arrays = {} + for i, settings in enumerate(cases): + request = scene(1.01e9, 4.01e9, 151) if i == 1 else scene() + request.update(settings, mode="gpr") + arrays[f"case_{i}_image"], hashes = gpr_reference(notebook, script, request) + arrays[f"case_{i}_request"] = np.array(json.dumps(request)) + arrays["provenance"] = np.array(json.dumps(dict(script=script.name, source_sha256=hashes, + notebook=notebook.name, numpy=np.__version__, scipy=scipy.__version__))) + np.savez_compressed(output, **arrays) + print(f"Wrote {output}: {output.stat().st_size} bytes, {len(cases)} GPR references") + + def generate(notebook: Path, output: Path) -> None: # Read up front so a missing notebook cannot produce an incomplete artifact. notebook.read_text(encoding="utf-8") cases = [ dict(compensation=True, normalize=True), - dict(compensation=True, normalize=False, speed=-.2, pitch=-7., ground=.5), + dict(compensation=True, normalize=False, equalize=False, speed=-.2, pitch=-7., ground=.5), dict(compensation=False, normalize=True, f_start=1.3e9, f_stop=3.7e9), dict(compensation=False, normalize=False, speed=0., pitch=0., radius=0), + dict(compensation=True, normalize=True, subtract_mean=True, mean_count=1), + dict(compensation=True, normalize=True, x_size=3., z_size=2., y_min=-.2, y_max=1.4), ] arrays, configurations = {}, [] for index, settings in enumerate(cases): @@ -187,9 +289,9 @@ def generate(notebook: Path, output: Path) -> None: arrays.update({f"case_{index}_{key}": value for key, value in expected.items()}) arrays["cases"] = np.array(json.dumps(configurations, separators=(",", ":"))) arrays["provenance"] = np.array(json.dumps(dict( - notebook=notebook.name, cell_source_sha256=hashes, cells=[0, 4], cases=len(cases), + notebook=notebook.name, cell_source_sha256=hashes, cells=[1, 14], cases=len(cases), numpy=np.__version__, scipy=scipy.__version__, - method="Unchanged notebook AST functions; max projections of sar_compute_volume output", + method="Unchanged notebook AST functions; per-frame NumPy previous-mean subtraction; max projections of sar_compute_volume output", ), sort_keys=True)) output.parent.mkdir(parents=True, exist_ok=True) np.savez_compressed(output, **arrays) @@ -198,7 +300,10 @@ def generate(notebook: Path, output: Path) -> None: if __name__ == "__main__": parser = argparse.ArgumentParser(description=__doc__) - parser.add_argument("--notebook", type=Path, default=ROOT / "Multi_antenna_3D_clean_v1.ipynb") + parser.add_argument("--notebook", type=Path, default=ROOT / "Multi_antenna_demonstrate_depthnorm.ipynb") parser.add_argument("--output", type=Path, default=REFERENCE) + parser.add_argument("--gpr-script", type=Path, default=ROOT / "Horns_motion_4x4_ground.py") + parser.add_argument("--gpr-output", type=Path, default=REFERENCE.with_name("gpr_reference.npz")) args = parser.parse_args() generate(args.notebook, args.output) + generate_gpr(args.notebook, args.gpr_script, args.gpr_output) diff --git a/data_acq_and_processing/tests/gpr_reference.npz b/data_acq_and_processing/tests/gpr_reference.npz new file mode 100644 index 0000000..c137e87 Binary files /dev/null and b/data_acq_and_processing/tests/gpr_reference.npz differ diff --git a/data_acq_and_processing/tests/native_regressions.cpp b/data_acq_and_processing/tests/native_regressions.cpp index 09ce4e0..9a0dfad 100644 --- a/data_acq_and_processing/tests/native_regressions.cpp +++ b/data_acq_and_processing/tests/native_regressions.cpp @@ -209,12 +209,14 @@ void test_configuration(const std::string& seed_path, const Workspace& workspace expect_error([&] { (void)config::load_run_config(config_path.string()); }); } const auto live_path = workspace.path / "live.json"; - std::ofstream(live_path) << Json{{"bscan_gain", 2}, {"settings_revision", 11}}; + std::ofstream(live_path) << Json{{"bscan_gain", 2}, {"settings_revision", 11}, + {"legacy_gpr_mode", false}, {"relative_permittivity", "obsolete"}}; processing::ProcessingLiveConfigLoader loader(live_path.string(), config); require(loader.refresh_if_needed().bscan_gain == 2, "valid live config did not load"); const auto revision = loader.revision(); for (const auto& invalid : {Json{{"bscan_gain", 1e100}}, - Json{{"gpr_input_positions", {0.0000001}}}, Json{{"bscan_gain", true}}}) { + Json{{"gpr_input_positions", {0.0000001}}}, Json{{"bscan_gain", true}}, + Json{{"processor_mode", "legacy_gpr"}}}) { std::ofstream(live_path) << invalid; std::this_thread::sleep_for(std::chrono::milliseconds(110)); const auto current = loader.refresh_if_needed(); diff --git a/data_acq_and_processing/tests/sar3d_fixture.cpp b/data_acq_and_processing/tests/sar3d_fixture.cpp deleted file mode 100644 index 36a3be4..0000000 --- a/data_acq_and_processing/tests/sar3d_fixture.cpp +++ /dev/null @@ -1,89 +0,0 @@ -#include "sar3d_processor.hpp" - -#include -#include -#include -#include -#include - -using Json = nlohmann::json; -using namespace radar; - -namespace { -auto complex_samples(const Json& samples) -> std::vector> { - std::vector> values; - for (const auto& sample : samples) { - values.emplace_back(sample.at(0).get(), sample.at(1).get()); - } - return values; -} - -auto collection(const Json& source) -> ipc::PreprocessedCollection { - ipc::PreprocessedCollection result; - result.collection_id = source.at("id"); - result.monotonic_ns = source.at("timestamp"); - for (const auto& entry : source.at("traces")) { - ipc::SweepTraceBlock trace; - trace.combo = {entry.at("input"), entry.at("output")}; - trace.capture_start_ns = entry.at("start"); - trace.capture_end_ns = entry.at("end"); - trace.frequency_hz = entry.at("frequency").get>(); - for (const auto& sample : complex_samples(entry.at("s21"))) { - trace.s21.push_back({sample.real(), sample.imag()}); - } - if (result.traces.empty() || trace.capture_start_ns < result.capture_start_ns) { - result.capture_start_ns = trace.capture_start_ns; - } - result.capture_end_ns = std::max(result.capture_end_ns, trace.capture_end_ns); - result.traces.push_back(std::move(trace)); - } - return result; -} -} - -// Runs the sar3d processor once: the JSON request on stdin holds the antenna -// geometry, the live settings and the latest and previous collections; the -// published payloads go to stdout as JSON. -int main() { - try { - const Json input = Json::parse(std::cin); - config::RunConfig config; - config.rings.results.slot_size_bytes = input.value("slot", config.rings.results.slot_size_bytes); - for (const auto& row : input.at("tx")) { - config.gpr.tx_geometry.push_back({row.at(0), row.at(1), row.at(2), row.at(3)}); - } - for (const auto& row : input.at("rx")) { - config.gpr.rx_geometry.push_back({row.at(0), row.at(1), row.at(2), row.at(3)}); - } - processing::ProcessingLiveConfig live; - live.gpr_start_freq_mhz = input.at("f_start").get() / 1e6F; - live.gpr_stop_freq_mhz = input.at("f_stop").get() / 1e6F; - live.gpr_min_depth_m = input.at("min_depth"); - live.gpr_max_depth_m = input.at("max_depth"); - live.gpr_speed_m_s = input.at("speed"); - live.gpr_look_angle_deg = input.at("pitch"); - live.sar3d_voxel_m = input.at("voxel"); - live.sar3d_bp_ground_depth_m = input.at("ground"); - live.sar3d_frame_radius = input.at("radius"); - live.sar3d_comp_range_power = input.at("power"); - live.sar3d_use_compensation = input.at("compensation"); - live.sar3d_normalize_each_frame = input.at("normalize"); - std::vector previous; - for (const auto& member : input.at("previous")) { - previous.push_back(collection(member)); - } - processing::Sar3dProcessor processor; - const auto result = processor.process_collection(config, collection(input.at("latest")), previous, live); - Json payloads = Json::array(); - for (const auto& payload : result.collection_payloads) { - payloads.push_back({{"name", payload.processing_name}, - {"x", payload.image_x_axis}, {"y", payload.image_y_axis}, - {"image", payload.image_values}, {"table", payload.table_values}, - {"columns", payload.table_columns}}); - } - std::cout << Json{{"payloads", payloads}} << '\n'; - } catch (const std::exception& error) { - std::cout << Json{{"error", error.what()}} << '\n'; - return 1; - } -} diff --git a/data_acq_and_processing/tests/sar3d_reference.npz b/data_acq_and_processing/tests/sar3d_reference.npz index d6af539..841b989 100644 Binary files a/data_acq_and_processing/tests/sar3d_reference.npz and b/data_acq_and_processing/tests/sar3d_reference.npz differ diff --git a/data_acq_and_processing/tests/test_sar3d.py b/data_acq_and_processing/tests/test_backprojection.py similarity index 67% rename from data_acq_and_processing/tests/test_sar3d.py rename to data_acq_and_processing/tests/test_backprojection.py index 2ffc645..47de655 100644 --- a/data_acq_and_processing/tests/test_sar3d.py +++ b/data_acq_and_processing/tests/test_backprojection.py @@ -1,9 +1,9 @@ """The sar3d processor on the real CUDA backend: notebook reference and contracts. -Every case runs build/native_tests/sar3d_fixture (SAR3D_FIXTURE overrides it). +Every case runs build/native_tests/backprojection_fixture (BACKPROJECTION_FIXTURE overrides it). Without a usable GPU the processor publishes sar3d_status 3: then only that warning is checked and the remaining tests are skipped. The notebook is needed -only to regenerate sar3d_reference.npz with generate_sar3d_reference.py. +only to regenerate sar3d_reference.npz with generate_backprojection_reference.py. """ from __future__ import annotations @@ -17,15 +17,18 @@ import unittest import numpy as np ROOT = Path(__file__).resolve().parents[2] -FIXTURE = Path(os.environ.get("SAR3D_FIXTURE", ROOT / "build/native_tests/sar3d_fixture")) +FIXTURE = Path(os.environ.get("BACKPROJECTION_FIXTURE", ROOT / "build/native_tests/backprojection_fixture")) REFERENCE = Path(__file__).with_name("sar3d_reference.npz") VIEWS = ("sar3d_front", "sar3d_side", "sar3d_top", "sar3d_oblique") -def run(request: dict) -> dict[str, dict]: +def run(request: dict | list[dict]) -> dict[str, dict]: output = subprocess.run([str(FIXTURE)], input=json.dumps(request), text=True, capture_output=True, timeout=30, check=True) - return {item["name"]: item for item in json.loads(output.stdout)["payloads"]} + result = json.loads(output.stdout) + if isinstance(request, list): + result = result[-1] + return {item["name"]: item for item in result["payloads"]} def status(payloads: dict[str, dict]) -> tuple[int, str]: @@ -87,6 +90,8 @@ class Sar3dTest(unittest.TestCase): rtol=2e-6, atol=1e-6) np.testing.assert_allclose(np.asarray(payload["image"]).reshape(expected.shape), expected, rtol=2e-3, atol=3e-4) + np.testing.assert_allclose(payloads["sar3d_target_width"]["table"], + self.reference[f"case_{index}_target_width"], rtol=2e-3, atol=2e-4) peak = payloads["sar3d_peak"]["table"] expected_peak = self.reference[f"case_{index}_peak"] np.testing.assert_allclose(peak[1:], expected_peak[1:], rtol=0, atol=2e-6) @@ -112,6 +117,15 @@ class Sar3dTest(unittest.TestCase): self.assertEqual(code, 2, message) self.assertTrue(message) + def test_cached_aperture_updates_when_previous_mean_context_changes(self) -> None: + request = self.scene(subtract_mean=True, mean_count=1) + shorter = copy.deepcopy(request) + shorter["previous"] = shorter["previous"][-1:] + cached, fresh = run([request, shorter]), run(shorter) + self.assertEqual(status(cached)[0], 0, status(cached)) + for name in VIEWS: + np.testing.assert_allclose(cached[name]["image"], fresh[name]["image"], rtol=2e-3, atol=3e-4) + def test_empty_signal_is_reported_without_a_peak(self) -> None: request = self.scene(previous=[], radius=0) for trace in request["latest"]["traces"]: @@ -135,5 +149,32 @@ class Sar3dTest(unittest.TestCase): self.assertLess(values.max(), 1e-10) +class GprTest(unittest.TestCase): + def test_changing_tx_rx_counts_rebuilds_the_cached_grid(self) -> None: + with np.load(REFERENCE.with_name("gpr_reference.npz"), allow_pickle=False) as reference: + base = json.loads(reference["case_0_request"].item()) + requests = [] + for tx, rx in ((2, 1), (1, 2)): + request = copy.deepcopy(base) + request.update(tx=[[i, 0, 0, 0] for i in range(tx)], rx=[[i, 0, 0, 0] for i in range(rx)]) + request["latest"]["traces"] = [trace for trace in request["latest"]["traces"] + if trace["output"] < tx and trace["input"] < rx] + requests.append(request) + np.testing.assert_allclose(run(requests)["gpr_accumulator"]["image"], + run(requests[-1])["gpr_accumulator"]["image"]) + + def test_coherent_ground_bp_matches_horns_and_notebook_equalizer(self) -> None: + with np.load(REFERENCE.with_name("gpr_reference.npz"), allow_pickle=False) as reference: + for i in range(3): + with self.subTest(case=i): + request = json.loads(reference[f"case_{i}_request"].item()) + payload = run(request)["gpr_accumulator"] + expected = reference[f"case_{i}_image"] + np.testing.assert_allclose(np.asarray(payload["image"]).reshape(expected.shape), expected, + rtol=2e-4, atol=3e-5) + np.testing.assert_allclose(payload["y"], np.linspace(.1, request["max_depth"], expected.shape[0]), + rtol=1e-6, atol=1e-7) + + if __name__ == "__main__": unittest.main() diff --git a/device_firmware/librevna_protocol.md b/device_firmware/librevna_protocol.md index 78019b8..dceac3d 100644 --- a/device_firmware/librevna_protocol.md +++ b/device_firmware/librevna_protocol.md @@ -273,6 +273,26 @@ runtime-бандлы пересобираются из NPZ при Start. Ста не меняет. RSL2 отвергается; Python и C++ обновляются одновременно после остановки процессов. +### Живые настройки backprojection + +GUI публикует настройки через `processing_live.json`; формат кадров колец от +изменения настроек не меняется. Для coherent GPR и SAR `pair_depth_equalize` +(boolean, по умолчанию true) включает `equalize_pairs_by_depth` перед медианной +нормировкой пары. `gpr_background_subtract_enabled` и `gpr_background_mean_count` +действуют в обоих режимах: среднее предыдущих коллекций вычитается до IFFT. +Плоскость coherent GPR — земля с высотой `gpr_bp_ground_depth_m` и наклоном +`gpr_look_angle_deg`. SAR использует `sar3d_frame_radius` предыдущих кадров и +текущий. Размеры сцены: `sar3d_x_size_m`, `sar3d_z_size_m` (0 — автоматически), +`sar3d_y_min_m`, `sar3d_y_max_m` — наземные координаты, независимо от depth gate. + +SAR публикует `sar3d_target_width` как TableF32 с тремя колонками [X,Y,Z]: +ширины непрерывного отклика на половине максимума по осевым сечениям через +совместный максимум сглаженного объёма (метры, не меньше одного вокселя). +`sar3d_peak` по-прежнему содержит [уровень,X,Y,Z]. GUI использует ширины для +квадратного контура цели. Зум и `SAR_CLOUD_THRESHOLD` — настройки отображения +desktop/веба; в live-файл не входят, ревизию и расчёт не меняют. RSL3 не меняется: +новая таблица использует существующий тип payload. + ## 7. Почему поле мощности читается, но не используется В парсере стоит: diff --git a/python_app/gui/controllers/config_controller.py b/python_app/gui/controllers/config_controller.py index 59f300c..c65b459 100644 --- a/python_app/gui/controllers/config_controller.py +++ b/python_app/gui/controllers/config_controller.py @@ -10,7 +10,7 @@ import math from pathlib import Path from PyQt6.QtCore import QSignalBlocker -from PyQt6.QtWidgets import QCheckBox, QComboBox, QDoubleSpinBox, QFormLayout, QLabel, QSpinBox +from PyQt6.QtWidgets import QCheckBox, QComboBox, QDoubleSpinBox, QFormLayout, QLabel, QSlider, QSpinBox from python_app.storage.npz.paths import radar_config_filename_prefix from python_app.gui.antenna_geometry_editor import AntennaGeometryEditor @@ -61,7 +61,6 @@ _ADC_FIELDS = ( ("radar.laser_control.variation.delay_time", "_optical_delay_time_input"), ) _GEOMETRY_FIELDS = { - "relative_permittivity": "_gpr_relative_permittivity", "tx_geometry": "_gpr_tx_geometry_input", "rx_geometry": "_gpr_rx_geometry_input", } @@ -97,15 +96,15 @@ def _set_web_widget(widget: SettingWidget, value: object) -> None: if isinstance(widget, QCheckBox): if not isinstance(value, bool): raise ValueError("Checkbox value must be true or false") - elif isinstance(widget, (QSpinBox, QDoubleSpinBox)): + elif isinstance(widget, (QSpinBox, QDoubleSpinBox, QSlider)): if isinstance(value, bool) or not isinstance(value, (int, float, str)): raise ValueError("Numeric setting requires a number") numeric = float(value) if not math.isfinite(numeric) or not widget.minimum() <= numeric <= widget.maximum(): raise ValueError("Numeric setting is outside its allowed range") - if isinstance(widget, QSpinBox) and numeric != int(numeric): + if isinstance(widget, (QSpinBox, QSlider)) and numeric != int(numeric): raise ValueError("Setting requires an integer") - value = int(numeric) if isinstance(widget, QSpinBox) else numeric + value = int(numeric) if isinstance(widget, (QSpinBox, QSlider)) else numeric elif isinstance(widget, QComboBox): if not isinstance(value, str) or widget.findText(value) < 0: raise ValueError(f"Unknown setting option: {value}") @@ -192,7 +191,7 @@ class ConfigController: config.input_switch.positions <= 1 and config.output_switch.positions <= 1 or config.input_switch.driver_mode == config.output_switch.driver_mode == "mock" ) - if self._scalar("_processing_mode") not in {"gpr", "legacy_gpr", "sar3d"} and static: + if self._scalar("_processing_mode") not in {"gpr", "sar3d"} and static: config.combos = [ComboModel(0, 0)] for key in PREPROCESS_ASSET_NAMES: asset = preprocess_asset_model(config, key) @@ -224,7 +223,6 @@ class ConfigController: def geometry(self) -> GprModel: return GprModel( - relative_permittivity=float(self._scalar("_gpr_relative_permittivity")), tx_geometry=[GprTxGeometryModel(**entry) for entry in self._geometry("tx_geometry").records()], rx_geometry=[GprRxGeometryModel(**entry) for entry in self._geometry("rx_geometry").records()], ) @@ -242,9 +240,10 @@ class ConfigController: if binding.spec.live in {"gpr_input_positions", "gpr_output_positions"}: value = parse_position_csv(str(value)) setattr(config, binding.spec.live, value) - config.relative_permittivity = float(self._scalar("_gpr_relative_permittivity")) config.tx_geometry = self._geometry("tx_geometry").records() config.rx_geometry = self._geometry("rx_geometry").records() + if config.sar3d_y_max_m <= config.sar3d_y_min_m: + raise ValueError("SAR_Y_RANGE_M max must be above min") for axis in ("x", "z"): lower_name, upper_name = f"gpr_visible_{axis}_min_m", f"gpr_visible_{axis}_max_m" lower, upper = sorted((getattr(config, lower_name), getattr(config, upper_name))) @@ -279,7 +278,6 @@ class ConfigController: if binding.spec.widget: write_setting(_scalar_widget(self.fields[binding.spec.widget]), binding.read(gui)) write_setting(_scalar_widget(self.fields["_log_level_combo"]), config.logging.level.capitalize()) - write_setting(_scalar_widget(self.fields["_gpr_relative_permittivity"]), config.gpr.relative_permittivity) for name, count, records in ( ("tx_geometry", config.output_switch.positions, config.gpr.tx_geometry), ("rx_geometry", config.input_switch.positions, config.gpr.rx_geometry), @@ -292,8 +290,8 @@ class ConfigController: widgets = {binding.spec.web: (binding.spec.group, self.fields[binding.spec.widget]) for binding in active_web_fields(mode)} for name, widget_id in _GEOMETRY_FIELDS.items(): - if mode in ({"legacy_gpr"} if name == "relative_permittivity" else {"gpr", "legacy_gpr", "sar3d"}): - widgets[name] = ("Geometry & medium", self.fields[widget_id]) + if mode in {"gpr", "sar3d"}: + widgets[name] = ("Geometry", self.fields[widget_id]) return widgets def web_schema(self) -> list[dict[str, object]]: @@ -315,9 +313,10 @@ class ConfigController: entry.update(kind="select", options=[widget.itemText(i) for i in range(widget.count())]) elif isinstance(widget, QCheckBox): entry.update(kind="bool") - elif isinstance(widget, (QSpinBox, QDoubleSpinBox)): - entry.update(kind="int" if isinstance(widget, QSpinBox) else "float", - min=widget.minimum(), max=widget.maximum(), step=widget.singleStep()) + elif isinstance(widget, (QSpinBox, QDoubleSpinBox, QSlider)): + entry.update(kind="int" if isinstance(widget, (QSpinBox, QSlider)) else "float", + min=widget.minimum(), max=widget.maximum(), step=widget.singleStep(), + slider=isinstance(widget, QSlider)) if isinstance(widget, QDoubleSpinBox): entry["decimals"] = widget.decimals() elif isinstance(widget, AntennaGeometryEditor): @@ -357,7 +356,7 @@ class ConfigController: geometry_snapshots.append((widget, widget.records())) else: scalar_snapshots.append((widget, read_setting(widget))) - if isinstance(widget, (QSpinBox, QDoubleSpinBox)): + if isinstance(widget, (QSpinBox, QDoubleSpinBox, QSlider)): ranges.append((widget, widget.minimum(), widget.maximum())) try: if "processor_mode" in fields: diff --git a/python_app/gui/controllers/live_settings_controller.py b/python_app/gui/controllers/live_settings_controller.py index ff37d49..64bd391 100644 --- a/python_app/gui/controllers/live_settings_controller.py +++ b/python_app/gui/controllers/live_settings_controller.py @@ -17,8 +17,7 @@ from python_app.models.run_config_schema import RunConfigModel from python_app.models.setting_fields import active_live_fields from python_app.orchestration.live_processing_config import ProcessingLiveConfig, ProcessingLiveConfigWriter -_GPR_MODES = {"gpr", "legacy_gpr"} -_GEOMETRY_MODES = _GPR_MODES | {"sar3d"} +_GEOMETRY_MODES = {"gpr", "sar3d"} class LiveSettingsController: @@ -84,7 +83,7 @@ class LiveSettingsController: def commit() -> None: self.validate_mode(self._mode.currentText()) live_fields = {binding.spec.web for binding in active_live_fields(self._mode.currentText())} - if settings.keys() & (live_fields | {"tx_geometry", "rx_geometry", "relative_permittivity"}): + if settings.keys() & (live_fields | {"tx_geometry", "rx_geometry"}): self.publish() self.config.apply_web_fields(settings, commit) @@ -171,12 +170,12 @@ class LiveSettingsController: self._log.show_exception("Failed to update plot display settings", exc) def locator_changed(self, *_args: object) -> None: - if self._mode.currentText() in _GPR_MODES and self._publish_changes(): + if self._mode.currentText() == "gpr" and self._publish_changes(): self.visual_changed() def set_mode_page(self, mode: str) -> None: pages = self._processing.pages - pages.setCurrentIndex({"pass_through": 0, "bscan": 1, "gpr": 2, "legacy_gpr": 3, "sar3d": 4}.get(mode, 0)) + pages.setCurrentIndex({"pass_through": 0, "bscan": 1, "gpr": 2, "sar3d": 3}.get(mode, 0)) self._processing.geometry_page.setVisible(mode in _GEOMETRY_MODES) page = pages.currentWidget() if page is not None: diff --git a/python_app/gui/controllers/plot_controller.py b/python_app/gui/controllers/plot_controller.py index d08ed04..22ba0ed 100644 --- a/python_app/gui/controllers/plot_controller.py +++ b/python_app/gui/controllers/plot_controller.py @@ -55,7 +55,7 @@ class PlotController: return self.sar3d.grid if current is self.sar3d.widget else current def set_mode(self, mode: str) -> None: - renderer = {"bscan": self.bscan, "gpr": self.gpr, "legacy_gpr": self.gpr, "sar3d": self.sar3d}.get(mode, self.trace) + renderer = {"bscan": self.bscan, "gpr": self.gpr, "sar3d": self.sar3d}.get(mode, self.trace) self.stack.setCurrentWidget(renderer.widget) self.sync_display_combos() @@ -114,15 +114,15 @@ class PlotController: def draw(self, collection: ResultCollection) -> bool: settings = self._settings() mode = settings.selected_mode - if mode in {"gpr", "legacy_gpr", "sar3d"} and self.required_revision is not None: + if mode in {"gpr", "sar3d"} and self.required_revision is not None: if collection.settings_revision != self.required_revision: return False if mode == "bscan": return self.draw_bscan() - if mode in {"gpr", "legacy_gpr"}: - return self.gpr.draw(collection, settings.legacy_gpr if mode == "legacy_gpr" else settings.gpr, self._geometry()) + if mode == "gpr": + return self.gpr.draw(collection, settings.gpr, self._geometry()) if mode == "sar3d": - return self.sar3d.draw(collection) + return self.sar3d.draw(collection, settings.sar3d, self._geometry()) return self.trace.draw(collection, settings.pass_through) def redraw(self, *_args: object) -> None: @@ -144,14 +144,13 @@ class PlotController: self.sync_display_combos() def clear_gpr(self) -> None: - settings = self._settings() - self.gpr.clear(settings.legacy_gpr if settings.selected_mode == "legacy_gpr" else settings.gpr) + self.gpr.clear() def clear_current(self) -> None: mode = self._settings().selected_mode if mode == "bscan": self.clear_bscan() - elif mode in {"gpr", "legacy_gpr"}: + elif mode == "gpr": self.clear_gpr() elif mode == "sar3d": self.sar3d.clear() diff --git a/python_app/gui/controllers/renderers/common.py b/python_app/gui/controllers/renderers/common.py index 893fe41..bab5a46 100644 --- a/python_app/gui/controllers/renderers/common.py +++ b/python_app/gui/controllers/renderers/common.py @@ -3,6 +3,8 @@ import numpy as np import pyqtgraph as pg +from python_app.models.run_config_schema import GprModel +from python_app.models.run_config_validation import parse_position_csv from python_app.gui.theme import PLOT_BACKGROUND @@ -22,3 +24,16 @@ def build_lut(stops: list[str], *, size: int = 256) -> np.ndarray: for channel in range(3): lut[:, channel] = np.interp(sample_positions, stop_positions, stop_colors[:, channel]).astype(np.uint8) return lut + + +# Same Tx/Rx markers on ground BP and the three planar SAR views. +ANTENNA_STYLES = (("t", "#ff595e", "#ffca3a"), ("t1", "#4cc9f0", "#e0fbfc")) + + +def selected_antenna_coordinates(geometry: GprModel, inputs: str, outputs: str) -> tuple[np.ndarray, np.ndarray]: + selected_inputs, selected_outputs = parse_position_csv(inputs), parse_position_csv(outputs) + tx = [[row.x_m, row.y_m, row.z_m] for row in geometry.tx_geometry + if not selected_outputs or row.output_pos in selected_outputs] + rx = [[row.x_m, row.y_m, row.z_m] for row in geometry.rx_geometry + if not selected_inputs or row.input_pos in selected_inputs] + return np.asarray(tx, dtype=float).reshape(-1, 3), np.asarray(rx, dtype=float).reshape(-1, 3) diff --git a/python_app/gui/controllers/renderers/gpr.py b/python_app/gui/controllers/renderers/gpr.py index b03fccf..8e94f8e 100644 --- a/python_app/gui/controllers/renderers/gpr.py +++ b/python_app/gui/controllers/renderers/gpr.py @@ -6,18 +6,14 @@ from PyQt6.QtCore import QRectF import numpy as np import pyqtgraph as pg -from python_app.gui.controllers.renderers.common import build_lut, create_plot +from python_app.gui.controllers.renderers.common import ANTENNA_STYLES, build_lut, create_plot, selected_antenna_coordinates from python_app.gui.result_payloads import ( - collection_payload_by_name, collection_payloads_by_prefix, - filter_object_rows, gpr_object_rows, + collection_payload_by_name, gpr_object_rows, ) -from python_app.models.run_config_validation import parse_position_csv from python_app.models.dataset_model import ResultCollection, ResultPayload -from python_app.models.gui_profile_schema import GuiGprStateModel, GuiLegacyGprStateModel +from python_app.models.gui_profile_schema import GuiGprStateModel from python_app.models.run_config_schema import GprModel -GprSettings = GuiGprStateModel | GuiLegacyGprStateModel - def normalized_display_range(start: float, stop: float, *, minimum_span: float = 0.1) -> tuple[float, float]: """Order display bounds and expand an empty interval around its center.""" @@ -54,8 +50,8 @@ def object_label_candidates( ] -def _title(settings: GprSettings, *, objects_only: bool = False) -> str: - title = f"Legacy GPR {settings.mode}" if isinstance(settings, GuiLegacyGprStateModel) else "GPR coherent BP" +def _title(*, objects_only: bool = False) -> str: + title = "GPR coherent BP — ground plane" return f"{title} Objects Only" if objects_only else title @@ -63,7 +59,7 @@ class GprRenderer: """Draw one collection from explicit display settings and antenna geometry.""" def __init__(self) -> None: - self.widget = create_plot("GPR coherent BP") + self.widget = create_plot("GPR coherent BP — ground plane") plot = self.widget.getPlotItem() assert plot is not None self._plot = plot @@ -75,17 +71,14 @@ class GprRenderer: self._tx = pg.ScatterPlotItem() self._rx = pg.ScatterPlotItem() self._points = pg.ScatterPlotItem() - self._region_centers = pg.ScatterPlotItem() - self._markers = (self._tx, self._rx, self._points, self._region_centers) - for item, z_value in zip((self._image, *self._markers), (0, 20, 20, 30, 30)): + self._markers = (self._tx, self._rx, self._points) + for item, z_value in zip((self._image, *self._markers), (0, 20, 20, 30)): item.setZValue(z_value) item.hide() self._plot.addItem(item) self._labels: list[pg.TextItem] = [] - self._masks: list[pg.ImageItem] = [] - self._contours: list[pg.IsocurveItem] = [] self._geometry_signature: tuple[object, ...] | None = None - self._selected_geometry = (np.zeros(0, dtype=np.float32), np.zeros(0, dtype=np.float32)) + self._selected_geometry = (np.empty((0, 3)), np.empty((0, 3))) def _configure_axes(self) -> None: plot = self._plot @@ -103,13 +96,11 @@ class GprRenderer: view.setMenuEnabled(False) def _clear_annotations(self) -> None: - for item in (*self._labels, *self._masks): + for item in self._labels: self._plot.removeItem(item) self._labels.clear() - self._masks.clear() - self._contours.clear() - def clear(self, settings: GprSettings) -> None: + def clear(self) -> None: """Clear displayed data while retaining reusable graphics items.""" self._configure_axes() self._clear_annotations() @@ -117,41 +108,25 @@ class GprRenderer: for marker in self._markers: marker.setData(x=[], y=[]) marker.hide() - self._plot.setTitle(_title(settings)) + self._plot.setTitle(_title()) - def _geometry(self, settings: GprSettings, geometry: GprModel) -> tuple[np.ndarray, np.ndarray]: - inputs = tuple(parse_position_csv(settings.input_positions)) - outputs = tuple(parse_position_csv(settings.output_positions)) + def _geometry(self, settings: GuiGprStateModel, geometry: GprModel) -> tuple[np.ndarray, np.ndarray]: signature = ( tuple((entry.output_pos, entry.x_m, entry.y_m, entry.z_m) for entry in geometry.tx_geometry), tuple((entry.input_pos, entry.x_m, entry.y_m, entry.z_m) for entry in geometry.rx_geometry), - inputs, outputs, + settings.input_positions, settings.output_positions, ) if signature != self._geometry_signature: - self._selected_geometry = ( - np.asarray([ - entry.x_m for entry in sorted(geometry.tx_geometry, key=lambda entry: entry.output_pos) - if not outputs or entry.output_pos in outputs - ], dtype=np.float32), - np.asarray([ - entry.x_m for entry in sorted(geometry.rx_geometry, key=lambda entry: entry.input_pos) - if not inputs or entry.input_pos in inputs - ], dtype=np.float32), - ) + self._selected_geometry = selected_antenna_coordinates(geometry, settings.input_positions, settings.output_positions) self._geometry_signature = signature return self._selected_geometry - def _draw_geometry(self, settings: GprSettings, geometry: GprModel) -> None: + def _draw_geometry(self, settings: GuiGprStateModel, geometry: GprModel) -> None: tx, rx = self._geometry(settings, geometry) - for marker, xs, symbol, brush, pen in ( - (self._tx, tx, "t", "#ff595e", "#ffca3a"), - (self._rx, rx, "t1", "#4cc9f0", "#e0fbfc"), - ): - marker.setData( - x=xs, y=np.zeros_like(xs), symbol=symbol, size=13, - brush=pg.mkBrush(brush), pen=pg.mkPen(pen, width=1.0), - ) - marker.setVisible(xs.size > 0) + for marker, coordinates, (symbol, brush, pen) in zip((self._tx, self._rx), (tx, rx), ANTENNA_STYLES): + marker.setData(x=coordinates[:, 0], y=coordinates[:, 2], symbol=symbol, size=13, + brush=pg.mkBrush(brush), pen=pg.mkPen(pen, width=1.0)) + marker.setVisible(len(coordinates) > 0) def _label(self, text: str, x: float, y: float, color: str) -> pg.TextItem: label = pg.TextItem(text=text, color=color, anchor=(0.0, 1.0)) @@ -161,34 +136,6 @@ class GprRenderer: self._labels.append(label) return label - def _draw_regions(self, collection: ResultCollection, rect: QRectF) -> None: - centers_payload = collection_payload_by_name(collection, "gpr_region_centers", kind=4) - if centers_payload is not None and np.asarray(centers_payload.table).size > 0: - centers = np.asarray(centers_payload.table, dtype=np.float32) - self._region_centers.setData( - x=centers[:, 0], y=centers[:, 1], symbol="o", size=10, - brush=pg.mkBrush("#80ed99"), pen=pg.mkPen("#081c15", width=1.1), - ) - self._region_centers.show() - for row in centers: - self._label(f"{float(row[2]):.0f}", float(row[0]), float(row[1]), "#d8f3dc") - for payload in collection_payloads_by_prefix(collection, "gpr_region_mask_", kind=3): - mask = np.asarray(payload.image, dtype=np.float32) - if mask.ndim != 2 or mask.size == 0: - continue - mask_image = pg.ImageItem(axisOrder="row-major") - mask_image.setZValue(5) - mask_image.setImage(mask, autoLevels=False) - mask_image.setRect(rect) - mask_image.setOpacity(0.0) - self._plot.addItem(mask_image) - contour = pg.IsocurveItem(data=mask) - contour.setLevel(0.5) - contour.setPen("#4cc9f0", width=1.3) - contour.setParentItem(mask_image) - self._masks.append(mask_image) - self._contours.append(contour) - @staticmethod def _place_object_label( label: pg.TextItem, x: float, z: float, x_span: float, z_span: float, occupied: list[QRectF], @@ -203,7 +150,7 @@ class GprRenderer: occupied.append(label.sceneBoundingRect().adjusted(-4.0, -4.0, 4.0, 4.0)) def _draw_points( - self, rows: np.ndarray, settings: GprSettings, objects_only: bool, x_span: float, z_span: float, + self, rows: np.ndarray, objects_only: bool, x_span: float, z_span: float, ) -> None: if rows.size == 0: return @@ -217,39 +164,33 @@ class GprRenderer: occupied: list[QRectF] = [] for x_value, z_value, score in rows: x, z = float(x_value), float(z_value) - text = f"{float(score):.0f}" if isinstance(settings, GuiLegacyGprStateModel) else f"{float(score):.2f}" + text = f"{float(score):.2f}" if objects_only: text += f" | x={x:.1f} | z={z:.1f}" label = self._label(text, x, z, "#ffd6d9" if objects_only else "#ffffff") if objects_only: self._place_object_label(label, x, z, x_span, z_span, occupied) - def draw(self, collection: ResultCollection, settings: GprSettings, geometry: GprModel) -> bool: + def draw(self, collection: ResultCollection, settings: GuiGprStateModel, geometry: GprModel) -> bool: """Render a heatmap or object-only view, returning whether data is present.""" objects_only = settings.render_mode == "objects_only" accumulator = collection_payload_by_name(collection, "gpr_accumulator", kind=3) x_min, x_max = normalized_display_range(settings.visible_x_min_m, settings.visible_x_max_m) z_min, z_max = normalized_display_range(settings.visible_z_min_m, settings.visible_z_max_m) rows = gpr_object_rows(collection) - # Coherent BP rows are finalized in C++; only legacy GPR uses GUI filtering. - if isinstance(settings, GuiLegacyGprStateModel): - rows = filter_object_rows( - rows, min_score=float(settings.min_visible_pair_count), - x_bounds=(x_min, x_max), z_bounds=(z_min, z_max), - ) if objects_only: if accumulator is None and rows.size == 0: - self.clear(settings) + self.clear() return False elif accumulator is None or not self._valid_accumulator(accumulator): - self.clear(settings) + self.clear() return False self.widget.setUpdatesEnabled(False) try: self._clear_annotations() self._image.hide() - for marker in (self._points, self._region_centers): + for marker in (self._points,): marker.setData(x=[], y=[]) marker.hide() self._view.invertY(False) @@ -274,13 +215,11 @@ class GprRenderer: points = collection_payload_by_name(collection, "gpr_points", kind=4) if points is None or np.asarray(points.table).size == 0: rows = np.zeros((0, 3), dtype=np.float32) - self._draw_points(rows, settings, objects_only, x_max - x_min, z_max - z_min) - if rect is not None: - self._draw_regions(collection, rect) + self._draw_points(rows, objects_only, x_max - x_min, z_max - z_min) if objects_only: self._view.setXRange(x_min, x_max, padding=0.0) self._view.setYRange(display_y_min(z_min, z_max), z_max, padding=0.0) - self._plot.setTitle(_title(settings, objects_only=objects_only)) + self._plot.setTitle(_title(objects_only=objects_only)) finally: self.widget.setUpdatesEnabled(True) return True diff --git a/python_app/gui/controllers/renderers/sar3d.py b/python_app/gui/controllers/renderers/sar3d.py index ebb4606..a4425c8 100644 --- a/python_app/gui/controllers/renderers/sar3d.py +++ b/python_app/gui/controllers/renderers/sar3d.py @@ -9,6 +9,7 @@ away from the Y boundaries: - ``sar3d_oblique``: orthographic projection along the camera of the notebook's 3D scene; ``sar3d_oblique_basis`` (2x3) gives its image axes u and v as scene directions, which places the scene box, the ground and the maximum on it; +- ``sar3d_target_width``: axial FWHM [X, Y, Z] in metres through the joint peak. - ``sar3d_peak``: [level, X, Y, Z] of that maximum before normalization. - ``sar3d_status``: [code, UTF-8 bytes...] with ready CUDA (0), no signal (1), processing error (2), or CUDA unavailable (3), including the explanation. @@ -23,13 +24,11 @@ import pyqtgraph as pg from PyQt6.QtCore import QRectF, Qt from PyQt6.QtWidgets import QGridLayout, QLabel, QVBoxLayout, QWidget -from python_app.gui.controllers.renderers.common import build_lut, create_plot +from python_app.gui.controllers.renderers.common import ANTENNA_STYLES, build_lut, create_plot, selected_antenna_coordinates from python_app.gui.result_payloads import collection_payload_by_name from python_app.models.dataset_model import ResultCollection, ResultPayload - -# Colours start at this fraction of the normalization reference; weaker voxels stay transparent, -# like those under the notebook's point-cloud threshold (SAR_CLOUD_THRESHOLD). -DISPLAY_THRESHOLD = 0.70 +from python_app.models.gui_profile_schema import GuiSar3dStateModel +from python_app.models.run_config_schema import GprModel # Status codes published by the sar3d processor in its "sar3d_status" table. SAR3D_STATES = {0: "ready", 1: "no signal", 2: "error", 3: "warning"} @@ -93,11 +92,19 @@ class Sar3dRenderer: self.plots[name] = plot lut = build_lut(["#081c15", "#1b4332", "#ffd166", "#f94144"]) # the GPR palette lut = np.column_stack([lut, np.full(len(lut), 255, dtype=np.uint8)]) - lut[0, 3] = 0 # below DISPLAY_THRESHOLD + lut[0, 3] = 0 # below SAR_CLOUD_THRESHOLD + self.antennas: dict[str, tuple[pg.ScatterPlotItem, pg.ScatterPlotItem]] = {} self.images: dict[str, pg.ImageItem] = {} - self.peaks: dict[str, pg.ScatterPlotItem] = {} + self.targets: dict[str, pg.PlotCurveItem] = {} for name, (_title, x_label, y_label, _axes) in SAR3D_VIEWS.items(): plot = self.plots[name] + plot.setAspectLocked(True) + if _axes is not None: + markers = (pg.ScatterPlotItem(), pg.ScatterPlotItem()) + for marker in markers: + marker.setZValue(4) + plot.addItem(marker) + self.antennas[name] = markers item = plot.getPlotItem() assert item is not None view = item.getViewBox() @@ -108,11 +115,11 @@ class Sar3dRenderer: image = pg.ImageItem(axisOrder="row-major", lut=lut) plot.addItem(image) # Hollow, so the response at the maximum stays visible. - peak = pg.ScatterPlotItem(symbol="star", size=18, pen=pg.mkPen("w", width=1.5), brush=None) + peak = pg.PlotCurveItem(pen=pg.mkPen("w", width=1.5)) peak.setZValue(3) plot.addItem(peak) self.images[name] = image - self.peaks[name] = peak + self.targets[name] = peak if x_label: plot.showGrid(x=True, y=True, alpha=0.2) plot.setLabel("bottom", x_label, units="m") @@ -146,15 +153,18 @@ class Sar3dRenderer: """Clear all four 3D SAR views.""" for name, (title, *_rest) in SAR3D_VIEWS.items(): self.images[name].hide() - self.peaks[name].setData([], []) + self.targets[name].setData([], []) self.plots[name].setTitle(title) + for markers in self.antennas.values(): + for marker in markers: + marker.setData([], []) self.box.setData([], []) self.ground.setData([], []) for label in self.axis_labels.values(): label.hide() self.status_label.setText(message) - def draw(self, collection: ResultCollection) -> bool: + def draw(self, collection: ResultCollection, settings: GuiSar3dStateModel, geometry: GprModel) -> bool: """Draw four maximum projections of the volume; return success.""" status = sar3d_status(collection) if status is not None and status[0] != 0: @@ -174,20 +184,44 @@ class Sar3dRenderer: peak = collection_payload_by_name(collection, "sar3d_peak", kind=4) level, *peak_xyz = peak.table[0] if peak is not None else (None,) + antennas = selected_antenna_coordinates(geometry, settings.input_positions, settings.output_positions) + pitch = np.radians(settings.look_angle_deg) + for coordinates in antennas: + y, z = coordinates[:, 1].copy(), coordinates[:, 2].copy() + coordinates[:, 1] = settings.bp_ground_depth_m + y*np.cos(pitch) - z*np.sin(pitch) + coordinates[:, 2] = y*np.sin(pitch) + z*np.cos(pitch) + width_payload = collection_payload_by_name(collection, "sar3d_target_width", kind=4) + widths = width_payload.table[0] if width_payload is not None else np.full(3, settings.voxel_m) + zoom = settings.zoom_percent / 100.0 for name, payload in payloads.items(): title, _x_label, _y_label, axes = SAR3D_VIEWS[name] rect = _image_rect(payload.image_x_axis, payload.image_y_axis) image = self.images[name] - image.setImage(payload.image, autoLevels=False, levels=(DISPLAY_THRESHOLD, 1.0)) + image.setImage(np.where(payload.image >= settings.cloud_threshold, payload.image, 0), + autoLevels=False, levels=(0, 1)) image.setRect(rect) image.show() plot = self.plots[name] - plot.setRange(rect, padding=0.0) - if peak_xyz: - position = basis.table @ np.asarray(peak_xyz) if axes is None else [peak_xyz[axis] for axis in axes] - self.peaks[name].setData([position[0]], [position[1]]) + position = (basis.table @ np.asarray(peak_xyz) if axes is None else np.asarray(peak_xyz)[list(axes)]) if peak_xyz else None + if axes is not None: + for marker, coordinates, (symbol, brush, pen) in zip(self.antennas[name], antennas, ANTENNA_STYLES): + xx, yy = coordinates[:, axes[0]], coordinates[:, axes[1]] + marker.setData(x=xx, y=yy, symbol=symbol, size=13, brush=pg.mkBrush(brush), pen=pg.mkPen(pen)) + if len(xx): + pad = settings.voxel_m + rect = rect.united(QRectF(xx.min()-pad, yy.min()-pad, np.ptp(xx)+2*pad, np.ptp(yy)+2*pad)) + if position is not None: + extent = np.abs(basis.table) @ widths if axes is None else widths[list(axes)] + half = .5 * max(extent) + self.targets[name].setData(position[0]+np.array([-half, half, half, -half, -half]), + position[1]+np.array([-half, -half, half, half, -half])) + if zoom > 1.0: + center = position + (np.array([rect.center().x(), rect.center().y()])-position) / zoom**2 + rect = QRectF(center[0]-rect.width()/(2*zoom), center[1]-rect.height()/(2*zoom), + rect.width()/zoom, rect.height()/zoom) else: - self.peaks[name].setData([], []) + self.targets[name].setData([], []) + plot.setRange(rect, padding=0.04) plot.setTitle(title) # Scene box and ground on the oblique view, from the extents of the other views. @@ -213,5 +247,5 @@ class Sar3dRenderer: if level is not None: x, y, z = peak_xyz summary += f" | max X {x:+.2f}, Y {y:+.2f}, Z {z:.2f} m | level {level:.3g}" - self.status_label.setText(summary) + self.status_label.setText(f"{summary} | zoom {settings.zoom_percent}%") return True diff --git a/python_app/gui/controllers/sections/processing_section.py b/python_app/gui/controllers/sections/processing_section.py index 081c28e..8c4342a 100644 --- a/python_app/gui/controllers/sections/processing_section.py +++ b/python_app/gui/controllers/sections/processing_section.py @@ -2,7 +2,7 @@ from __future__ import annotations -from PyQt6.QtCore import QRegularExpression +from PyQt6.QtCore import QRegularExpression, Qt from PyQt6.QtGui import QRegularExpressionValidator from PyQt6.QtWidgets import ( QCheckBox, @@ -10,6 +10,7 @@ from PyQt6.QtWidgets import ( QFormLayout, QGroupBox, QLineEdit, + QSlider, QSizePolicy, QStackedWidget, QVBoxLayout, @@ -51,7 +52,7 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces group = view.group form = QFormLayout(group) form.setFieldGrowthPolicy(QFormLayout.FieldGrowthPolicy.AllNonFixedFieldsGrow) - processing_mode = view.add("_processing_mode", choice(['pass_through', 'bscan', 'gpr', 'legacy_gpr', 'sar3d'])) + processing_mode = view.add("_processing_mode", choice(['pass_through', 'bscan', 'gpr', 'sar3d'])) processing_mode_pages = QStackedWidget(group) processing_mode_pages.setSizePolicy(QSizePolicy.Policy.Preferred, QSizePolicy.Policy.Fixed) show_magnitude_checkbox = view.add("_show_magnitude_checkbox", QCheckBox("Show magnitude")) @@ -112,10 +113,6 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces gpr_defaults = config.gpr - # Medium permittivity is a legacy-GPR-only knob (shown on the legacy page below). - # Coherent BP fixes the medium to eps_r = 1, so it is not offered there. - gpr_relative_permittivity = view.add("_gpr_relative_permittivity", double_spin(0.0001, 1000.0, 4, 0.05)) - gpr_relative_permittivity.setValue(float(gpr_defaults.relative_permittivity)) gpr_tx_geometry_input = view.add("_gpr_tx_geometry_input", AntennaGeometryEditor( "output_pos", gpr_defaults.tx_geometry, position_count=config.output_switch.positions, )) @@ -178,10 +175,9 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces gpr_visible_x_max_m = view.add("_gpr_visible_x_max_m", double_spin(-100.0, 100.0, 2, 0.1)) gpr_visible_z_min_m = view.add("_gpr_visible_z_min_m", double_spin(0.0, 100.0, 2, 0.1)) gpr_visible_z_max_m = view.add("_gpr_visible_z_max_m", double_spin(0.0, 100.0, 2, 0.1)) - gpr_imaging_plane_y_m = view.add("_gpr_imaging_plane_y_m", double_spin(-50.0, 50.0, 3, 0.05)) - gpr_imaging_plane_y_m.setToolTip( - "Y coordinate of the BP imaging slice (m). Use 0 for legacy 1D antenna layouts." - ) + gpr_bp_ground_depth_m = view.add("_gpr_bp_ground_depth_m", double_spin(0.0, 50.0, 3, 0.05)) + gpr_bp_ground_depth_m.setToolTip("Height of the radar origin above the ground; the BP slice follows the ground and radar tilt.") + gpr_pair_depth_equalize = view.add("_gpr_pair_depth_equalize", QCheckBox("PAIR_DEPTH_EQUALIZE")) gpr_page = _build_processing_mode_page( processing_mode_pages, @@ -204,89 +200,22 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces ("Approach min frames", gpr_object_approach_min_frames), ("Start (MHz)", gpr_start_freq_mhz), ("Stop (MHz)", gpr_stop_freq_mhz), - ("Imaging plane Y (m)", gpr_imaging_plane_y_m), + ("Radar height (m)", gpr_bp_ground_depth_m), ("Visible X min (m)", gpr_visible_x_min_m), ("Visible X max (m)", gpr_visible_x_max_m), ("Visible Z min (m)", gpr_visible_z_min_m), ("Visible Z max (m)", gpr_visible_z_max_m), gpr_background_subtract_enabled, ("Mean count", gpr_background_mean_count), + gpr_pair_depth_equalize, gpr_remove_sidelobe_objects_enabled, ], split_index=8, ) processing_mode_pages.addWidget(gpr_page) - legacy_gpr_config_mode = view.add("_legacy_gpr_config_mode", choice(['point', 'extended'])) - legacy_gpr_input_positions_input = view.add("_legacy_gpr_input_positions_input", QLineEdit()) - legacy_gpr_input_positions_input.setPlaceholderText("0,1,2") - legacy_gpr_output_positions_input = view.add("_legacy_gpr_output_positions_input", QLineEdit()) - legacy_gpr_output_positions_input.setPlaceholderText("0,1") - legacy_gpr_min_depth_m = view.add("_legacy_gpr_min_depth_m", double_spin(0.0, 50.0, 2, 0.1)) - legacy_gpr_max_depth_m = view.add("_legacy_gpr_max_depth_m", double_spin(0.1, 50.0, 2, 0.1)) - legacy_gpr_comp_power = view.add("_legacy_gpr_comp_power", double_spin(0.0, 5.0, 3, 0.05)) - legacy_gpr_snr_thresh = view.add("_legacy_gpr_snr_thresh", double_spin(0.0, 1000.0, 2, 0.1)) - legacy_gpr_snr_comp_max = view.add("_legacy_gpr_snr_comp_max", double_spin(0.0, 1000.0, 2, 0.5)) - legacy_gpr_start_freq_mhz = view.add("_legacy_gpr_start_freq_mhz", double_spin(100.0, 8800.0, 1, 10.0)) - legacy_gpr_stop_freq_mhz = view.add("_legacy_gpr_stop_freq_mhz", double_spin(100.0, 8800.0, 1, 10.0)) - legacy_gpr_speed_m_s = view.add("_legacy_gpr_speed_m_s", double_spin(-100.0, 100.0, 3, 0.01)) - legacy_gpr_ignore_socket_speed_enabled = view.add("_legacy_gpr_ignore_socket_speed_enabled", QCheckBox("Ignore socket speed")) - legacy_gpr_ignore_socket_speed_enabled.setToolTip(_IGNORE_SOCKET_SPEED_TOOLTIP) - legacy_gpr_look_angle_deg = view.add("_legacy_gpr_look_angle_deg", double_spin(-90.0, 90.0, 2, 0.1)) - - # Intra-sweep phase compensation toggle. When checked, S21 phase is - # corrected per frequency point before the IFFT so the radar's motion - # during one sweep is removed at the spectrum stage. Only matters when - # speed != 0; otherwise has no effect by construction. - legacy_gpr_apply_freq_phase_correction = view.add("_legacy_gpr_apply_freq_phase_correction", QCheckBox("Apply intra-sweep phase correction")) - - # Anchor point for the motion model's per-event `dt_ref`. `frame_center` - # spreads motion symmetrically around the frame midpoint; `first_tx_event` - # accumulates it forward from the first event center. - legacy_gpr_reference_mode = view.add("_legacy_gpr_reference_mode", choice(['frame_center', 'first_tx_event'])) - legacy_gpr_background_subtract_enabled = view.add("_legacy_gpr_background_subtract_enabled", QCheckBox("Subtract mean of previous collections")) - legacy_gpr_background_mean_count = view.add("_legacy_gpr_background_mean_count", integer_spin(0, 10000)) - legacy_gpr_render_mode = view.add("_legacy_gpr_render_mode", choice(['heatmap', 'objects_only'])) - legacy_gpr_min_visible_pair_count = view.add("_legacy_gpr_min_visible_pair_count", integer_spin(1, 10000)) - legacy_gpr_visible_x_min_m = view.add("_legacy_gpr_visible_x_min_m", double_spin(-100.0, 100.0, 2, 0.1)) - legacy_gpr_visible_x_max_m = view.add("_legacy_gpr_visible_x_max_m", double_spin(-100.0, 100.0, 2, 0.1)) - legacy_gpr_visible_z_min_m = view.add("_legacy_gpr_visible_z_min_m", double_spin(0.0, 100.0, 2, 0.1)) - legacy_gpr_visible_z_max_m = view.add("_legacy_gpr_visible_z_max_m", double_spin(0.0, 100.0, 2, 0.1)) - - legacy_gpr_page = _build_processing_mode_page( - processing_mode_pages, - [ - ("Config mode", legacy_gpr_config_mode), - ("Relative permittivity", gpr_relative_permittivity), - ("Input positions", legacy_gpr_input_positions_input), - ("Output positions", legacy_gpr_output_positions_input), - ("Min depth (m)", legacy_gpr_min_depth_m), - ("Max depth (m)", legacy_gpr_max_depth_m), - ("Compensation power", legacy_gpr_comp_power), - ("SNR threshold", legacy_gpr_snr_thresh), - ("Max SNR compensation", legacy_gpr_snr_comp_max), - ("Speed (m/s)", legacy_gpr_speed_m_s), - legacy_gpr_ignore_socket_speed_enabled, - ("Render mode", legacy_gpr_render_mode), - ("Min visible pairs", legacy_gpr_min_visible_pair_count), - ("Start (MHz)", legacy_gpr_start_freq_mhz), - ("Stop (MHz)", legacy_gpr_stop_freq_mhz), - ("Look angle (°)", legacy_gpr_look_angle_deg), - ("Reference mode", legacy_gpr_reference_mode), - legacy_gpr_apply_freq_phase_correction, - ("Visible X min (m)", legacy_gpr_visible_x_min_m), - ("Visible X max (m)", legacy_gpr_visible_x_max_m), - ("Visible Z min (m)", legacy_gpr_visible_z_min_m), - ("Visible Z max (m)", legacy_gpr_visible_z_max_m), - legacy_gpr_background_subtract_enabled, - ("Mean count", legacy_gpr_background_mean_count), - ], - split_index=10, - ) - processing_mode_pages.addWidget(legacy_gpr_page) - # 3D SAR page. The labels are the names of the reference notebook - # Multi_antenna_3D_clean_v1.ipynb; the shared gpr_* fields get their own - # widgets here, like on the legacy page. + # Multi_antenna_demonstrate_depthnorm.ipynb; the shared gpr_* fields get their own + # widgets here, like on the GPR page. sar3d_input_positions_input = view.add("_sar3d_input_positions_input", QLineEdit()) sar3d_input_positions_input.setPlaceholderText("0,1,2") sar3d_output_positions_input = view.add("_sar3d_output_positions_input", QLineEdit()) @@ -311,8 +240,8 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces sar3d_speed_m_s = view.add("_sar3d_speed_m_s", double_spin(-100.0, 100.0, 3, 0.05)) sar3d_frame_radius = view.add("_sar3d_frame_radius", integer_spin(0, 24)) sar3d_frame_radius.setToolTip( - "Aperture of 2*r+1 consecutive frames. 0: the current frame only (needs a vertical antenna " - "baseline); >0 adds aperture from motion and needs an accurate speed. The image lags by r frames." + "Causal aperture: current and r previous frames. 0: the current frame only (needs a vertical antenna " + "baseline); >0 adds aperture from motion and needs an accurate speed. The scene follows the current frame." ) sar3d_bp_ground_depth_m = view.add("_sar3d_bp_ground_depth_m", double_spin(0.0, 10.0, 2, 0.05)) sar3d_bp_ground_depth_m.setToolTip("Height of the lower antenna row above the ground (m); scene Y = 0 is the ground.") @@ -321,6 +250,17 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces sar3d_use_compensation = view.add("_sar3d_use_compensation", QCheckBox("SAR_USE_COMPENSATION")) sar3d_normalize_each_frame = view.add("_sar3d_normalize_each_frame", QCheckBox("SAR_NORMALIZE_EACH_FRAME")) sar3d_normalize_each_frame.setToolTip("Scale every Tx/Rx A-scan by its median inside the depth gate.") + sar3d_pair_depth_equalize = view.add("_sar3d_pair_depth_equalize", QCheckBox("PAIR_DEPTH_EQUALIZE")) + sar3d_background_subtract_enabled = view.add("_sar3d_background_subtract_enabled", QCheckBox("Subtract mean of previous collections")) + sar3d_background_mean_count = view.add("_sar3d_background_mean_count", integer_spin(0, 10000)) + sar3d_x_size_m = view.add("_sar3d_x_size_m", double_spin(0.0, 100.0, 2, 0.1)) + sar3d_y_min_m = view.add("_sar3d_y_min_m", double_spin(-50.0, 50.0, 2, 0.1)) + sar3d_y_max_m = view.add("_sar3d_y_max_m", double_spin(-50.0, 50.0, 2, 0.1)) + sar3d_z_size_m = view.add("_sar3d_z_size_m", double_spin(0.0, 100.0, 2, 0.1)) + sar3d_cloud_threshold = view.add("_sar3d_cloud_threshold", double_spin(0.0, 1.0, 2, 0.05)) + sar3d_zoom_percent = view.add("_sar3d_zoom_percent", QSlider(Qt.Orientation.Horizontal)) + sar3d_zoom_percent.setRange(100, 1000) + sar3d_zoom_percent.setToolTip("100% shows the full scene; higher values zoom towards its brightest object.") sar3d_page = _build_processing_mode_page( processing_mode_pages, @@ -340,6 +280,15 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces ("COMP_RANGE_POWER", sar3d_comp_range_power), sar3d_use_compensation, sar3d_normalize_each_frame, + sar3d_pair_depth_equalize, + sar3d_background_subtract_enabled, + ("Mean count", sar3d_background_mean_count), + ("SAR_X_SIZE_M (0 = auto)", sar3d_x_size_m), + ("SAR_Y_RANGE_M min", sar3d_y_min_m), + ("SAR_Y_RANGE_M max", sar3d_y_max_m), + ("SAR_Z_SIZE_M (0 = auto)", sar3d_z_size_m), + ("SAR_CLOUD_THRESHOLD", sar3d_cloud_threshold), + ("Zoom (%)", sar3d_zoom_percent), ], split_index=8, ) @@ -367,17 +316,14 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces geometry_page=gpr_common_page, geometries=(gpr_tx_geometry_input, gpr_rx_geometry_input), fixed_y=pass_through_fixed_y_enabled, y_bounds=(pass_through_y_min_db, pass_through_y_max_db), speed_controls=((gpr_ignore_socket_speed_enabled, gpr_speed_m_s), - (legacy_gpr_ignore_socket_speed_enabled, legacy_gpr_speed_m_s), (sar3d_ignore_socket_speed_enabled, sar3d_speed_m_s)), trace_controls=(show_magnitude_checkbox, show_phase_checkbox, unwrap_phase_checkbox, pass_through_combo_filter_input, pass_through_fixed_y_enabled, pass_through_y_min_db, pass_through_y_max_db), - visual_controls=(gpr_render_mode, legacy_gpr_render_mode), + visual_controls=(gpr_render_mode, sar3d_cloud_threshold, sar3d_zoom_percent), locator_controls=(gpr_max_detected_objects_to_draw, gpr_draw_top_m_objects, - legacy_gpr_min_visible_pair_count, gpr_visible_x_min_m, gpr_visible_x_max_m, - gpr_visible_z_min_m, gpr_visible_z_max_m, legacy_gpr_visible_x_min_m, - legacy_gpr_visible_x_max_m, legacy_gpr_visible_z_min_m, legacy_gpr_visible_z_max_m), + gpr_visible_x_min_m, gpr_visible_x_max_m, + gpr_visible_z_min_m, gpr_visible_z_max_m), frequency_controls=(bscan_start_freq_mhz, bscan_stop_freq_mhz, gpr_start_freq_mhz, gpr_stop_freq_mhz, - legacy_gpr_start_freq_mhz, legacy_gpr_stop_freq_mhz, sar3d_start_freq_mhz, sar3d_stop_freq_mhz), ) diff --git a/python_app/gui/controllers/sections/views.py b/python_app/gui/controllers/sections/views.py index 1e8f5ef..ee5ebff 100644 --- a/python_app/gui/controllers/sections/views.py +++ b/python_app/gui/controllers/sections/views.py @@ -6,13 +6,13 @@ from typing import TypeAlias, TypeVar from PyQt6.QtWidgets import ( QCheckBox, QComboBox, QDoubleSpinBox, QGroupBox, QLineEdit, QPlainTextEdit, - QSpinBox, QWidget, QLabel, QPushButton, QStackedWidget, + QSpinBox, QSlider, QWidget, QLabel, QPushButton, QStackedWidget, ) from python_app.gui.antenna_geometry_editor import AntennaGeometryEditor from python_app.models.setting_fields import Scalar -ScalarWidget: TypeAlias = QLineEdit | QPlainTextEdit | QComboBox | QCheckBox | QSpinBox | QDoubleSpinBox +ScalarWidget: TypeAlias = QLineEdit | QPlainTextEdit | QComboBox | QCheckBox | QSpinBox | QDoubleSpinBox | QSlider SettingWidget: TypeAlias = ScalarWidget | AntennaGeometryEditor _Widget = TypeVar("_Widget", bound=SettingWidget) @@ -102,7 +102,7 @@ class RadarView(SectionView): def read_setting(widget: ScalarWidget) -> Scalar: if isinstance(widget, QCheckBox): return widget.isChecked() - if isinstance(widget, (QSpinBox, QDoubleSpinBox)): + if isinstance(widget, (QSpinBox, QDoubleSpinBox, QSlider)): return widget.value() if isinstance(widget, QComboBox): return widget.currentText() @@ -114,7 +114,7 @@ def read_setting(widget: ScalarWidget) -> Scalar: def write_setting(widget: ScalarWidget, value: Scalar) -> None: if isinstance(widget, QCheckBox): widget.setChecked(bool(value)) - elif isinstance(widget, QSpinBox): + elif isinstance(widget, (QSpinBox, QSlider)): widget.setValue(int(value)) elif isinstance(widget, QDoubleSpinBox): widget.setValue(float(value)) @@ -134,7 +134,7 @@ def write_setting(widget: ScalarWidget, value: Scalar) -> None: def connect_setting(widget: ScalarWidget, callback: Callable[..., object]) -> None: if isinstance(widget, QCheckBox): signal = widget.toggled - elif isinstance(widget, (QSpinBox, QDoubleSpinBox)): + elif isinstance(widget, (QSpinBox, QDoubleSpinBox, QSlider)): signal = widget.valueChanged elif isinstance(widget, QComboBox): signal = widget.currentTextChanged diff --git a/python_app/gui/result_payloads.py b/python_app/gui/result_payloads.py index b9c8677..acfd295 100644 --- a/python_app/gui/result_payloads.py +++ b/python_app/gui/result_payloads.py @@ -1,4 +1,4 @@ -"""Result payload lookup and GPR display filtering.""" +"""Result payload lookup and coherent-GPR objects.""" from __future__ import annotations @@ -22,22 +22,6 @@ def collection_payload_by_name( return None -def collection_payloads_by_prefix( - collection: ResultCollection, - prefix: str, - kind: int | None = None, -) -> list[ResultPayload]: - """Return collection payloads matching a processing-name prefix.""" - payloads: list[ResultPayload] = [] - for payload in collection.collection_payloads: - if not str(payload.processing_name).startswith(prefix): - continue - if kind is not None and int(payload.kind) != int(kind): - continue - payloads.append(payload) - return payloads - - def gpr_object_rows(collection: ResultCollection) -> np.ndarray: """Return object rows as `[x_m, z_m, score]` from a GPR collection.""" points_payload = collection_payload_by_name(collection, "gpr_points", kind=4) @@ -46,33 +30,4 @@ def gpr_object_rows(collection: ResultCollection) -> np.ndarray: if points.ndim == 2 and points.shape[1] >= 3: return points[:, :3] - centers_payload = collection_payload_by_name(collection, "gpr_region_centers", kind=4) - if centers_payload is not None: - centers = np.asarray(centers_payload.table, dtype=np.float32) - if centers.ndim == 2 and centers.shape[1] >= 3: - return centers[:, :3] - return np.zeros((0, 3), dtype=np.float32) - - -def filter_object_rows( - rows: np.ndarray, - *, - min_score: float, - x_bounds: tuple[float, float], - z_bounds: tuple[float, float], -) -> np.ndarray: - """Keep finite objects inside the visible window and above the score threshold.""" - if rows.size == 0: - return rows - x_min, x_max = x_bounds - z_min, z_max = z_bounds - visible_mask = ( - np.all(np.isfinite(rows[:, :3]), axis=1) - & (rows[:, 2] >= min_score) - & (rows[:, 0] >= x_min) - & (rows[:, 0] <= x_max) - & (rows[:, 1] >= z_min) - & (rows[:, 1] <= z_max) - ) - return rows[visible_mask] diff --git a/python_app/gui/runtime/constraints.py b/python_app/gui/runtime/constraints.py index 41c28b4..019e0b0 100644 --- a/python_app/gui/runtime/constraints.py +++ b/python_app/gui/runtime/constraints.py @@ -26,7 +26,7 @@ def validate_processing_mode_constraints( # sar3d shares the GPR combo/geometry requirements: it images from the same # multi-Tx/Rx pair set selected through the gpr_* live-config fields. - if processing_mode not in {"gpr", "legacy_gpr", "sar3d"}: + if processing_mode not in {"gpr", "sar3d"}: return configured_combos = {(int(combo.input), int(combo.output)) for combo in config.combos} diff --git a/python_app/gui/settings_defaults.py b/python_app/gui/settings_defaults.py index a01b3bd..6549c2e 100644 --- a/python_app/gui/settings_defaults.py +++ b/python_app/gui/settings_defaults.py @@ -5,7 +5,7 @@ from pathlib import Path from python_app.models.run_config_schema import RunConfigModel, ComboModel from python_app.models.gui_profile_schema import ( GuiStateModel, GuiSwitchStateModel, GuiProcessingStateModel, GuiGprStateModel, - GuiLegacyGprStateModel, GuiSar3dStateModel, GuiDataActionsStateModel, + GuiSar3dStateModel, GuiDataActionsStateModel, ) GUI_SAVE_HISTORY_LIMIT = 1000 @@ -91,12 +91,6 @@ def default_gui_state_for_config(config: RunConfigModel, project_root: Path) -> visible_x_min_m=default_gpr_x_min_m, visible_x_max_m=default_gpr_x_max_m, ), - legacy_gpr=GuiLegacyGprStateModel( - input_positions=inputs, - output_positions=outputs, - visible_x_min_m=default_gpr_x_min_m, - visible_x_max_m=default_gpr_x_max_m, - ), sar3d=GuiSar3dStateModel( input_positions=inputs, output_positions=outputs, diff --git a/python_app/models/gui_profile_codec.py b/python_app/models/gui_profile_codec.py index 8d1288e..fe0bdcb 100644 --- a/python_app/models/gui_profile_codec.py +++ b/python_app/models/gui_profile_codec.py @@ -10,7 +10,7 @@ from python_app.models.json_fields import finite_number, json_object from python_app.models.run_config_schema import RunConfigModel -def _decode_state(model: object, payload: dict[str, Any], context: str, *, check_constraints: bool = True) -> None: +def _decode_state(model: object, payload: dict[str, Any], context: str) -> None: """Read known fields recursively; unknown fields do not survive re-encoding.""" if not is_dataclass(model) or isinstance(model, type): raise TypeError("GUI state must be a dataclass instance") @@ -19,14 +19,14 @@ def _decode_state(model: object, payload: dict[str, Any], context: str, *, check value = payload.get(member.name, default) name = f"{context}.{member.name}" if is_dataclass(default): - _decode_state(default, json_object(payload.get(member.name), name), name, check_constraints=check_constraints) + _decode_state(default, json_object(payload.get(member.name), name), name) continue if isinstance(default, float): value = finite_number(value, name) elif type(value) is not type(default): raise ValueError(f"{name} must be a JSON {type(default).__name__}") spec = member.metadata.get("setting") - if check_constraints and isinstance(spec, Setting): + if isinstance(spec, Setting): if spec.choices and value not in spec.choices: raise ValueError(f"{name} must be one of: {', '.join(spec.choices)}") if spec.minimum is not None: @@ -39,33 +39,16 @@ def _decode_state(model: object, payload: dict[str, Any], context: str, *, check def gui_profile_from_dict(payload: dict[str, Any]) -> GuiProfileModel: - """Read the pipeline config and optional GUI state, preserving old GPR profiles.""" + """Read the pipeline config and optional GUI state; ignore unknown fields.""" profile = GuiProfileModel(run_config=RunConfigModel.from_dict(payload)) if payload.get("gui") is None: return profile gui_object = json_object(payload["gui"], "gui") - processing = json_object(gui_object.get("processing"), "gui.processing") - old_gpr = json_object(processing.get("gpr"), "gui.processing.gpr") gui = GuiStateModel() - # Older profiles kept both algorithms under gpr. Root mode and the short-lived - # algorithm selector choose the legacy page; explicit legacy values win. - root_gpr = payload.get("gpr") - root_mode = root_gpr.get("mode") if isinstance(root_gpr, dict) else None - algorithm_mode = {"legacy_point": "point", "legacy_extended": "extended"}.get( - str(old_gpr.get("algorithm", ""))) - if root_mode in ("point", "extended"): - gui.processing.legacy_gpr.mode = root_mode - if algorithm_mode is not None: - gui.processing.legacy_gpr.mode = algorithm_mode - _decode_state(gui.processing.legacy_gpr, old_gpr, "gui.processing.gpr", check_constraints=False) _decode_state(gui, gui_object, "gui") if gui.version != 1: raise ValueError(f"Unsupported gui.version: {gui.version}") - if gui.processing.selected_mode == "gpr" and ( - algorithm_mode is not None or root_mode in ("point", "extended") - ): - gui.processing.selected_mode = "legacy_gpr" gui.preprocess_dialog.median_sweep_count = max(1, gui.preprocess_dialog.median_sweep_count) profile.gui = gui return profile diff --git a/python_app/models/gui_profile_schema.py b/python_app/models/gui_profile_schema.py index 888430f..3c4ebb2 100644 --- a/python_app/models/gui_profile_schema.py +++ b/python_app/models/gui_profile_schema.py @@ -121,7 +121,8 @@ class GuiGprStateModel: remove_sidelobe_objects_enabled: bool = setting( True, '_gpr_remove_sidelobe_objects_enabled', live='gpr_remove_sidelobe_objects_enabled', group='Imaging', ) - imaging_plane_y_m: float = setting(0.0, '_gpr_imaging_plane_y_m', live='gpr_imaging_plane_y_m', group='Geometry & depth') + bp_ground_depth_m: float = setting(1.25, '_gpr_bp_ground_depth_m', live='gpr_bp_ground_depth_m', group='Geometry & depth', minimum=0) + pair_depth_equalize: bool = setting(True, '_gpr_pair_depth_equalize', live='pair_depth_equalize', group='Imaging') render_mode: str = setting("heatmap", '_gpr_render_mode', web='render_mode', group='Display', choices=('heatmap', 'objects_only')) visible_x_min_m: float = setting(-2.0, '_gpr_visible_x_min_m', live='gpr_visible_x_min_m', group='Visible window') visible_x_max_m: float = setting(2.0, '_gpr_visible_x_max_m', live='gpr_visible_x_max_m', group='Visible window') @@ -129,48 +130,11 @@ class GuiGprStateModel: visible_z_max_m: float = setting(14.0, '_gpr_visible_z_max_m', live='gpr_visible_z_max_m', group='Visible window') -@dataclass(slots=True) -class GuiLegacyGprStateModel: - """UI-only defaults for legacy point/extended GPR live settings.""" - - mode: str = setting("point", '_legacy_gpr_config_mode', live='legacy_gpr_mode', group='Mode', choices=('point', 'extended')) - input_positions: str = setting("", '_legacy_gpr_input_positions_input', live='gpr_input_positions', group='Geometry & depth') - output_positions: str = setting("", '_legacy_gpr_output_positions_input', live='gpr_output_positions', group='Geometry & depth') - min_depth_m: float = setting(2.0, '_legacy_gpr_min_depth_m', live='gpr_min_depth_m', group='Geometry & depth') - max_depth_m: float = setting(14.0, '_legacy_gpr_max_depth_m', live='gpr_max_depth_m', group='Geometry & depth') - comp_power: float = setting(0.2, '_legacy_gpr_comp_power', live='gpr_comp_power', group='Detection', minimum=0) - start_freq_mhz: float = setting(3000.0, '_legacy_gpr_start_freq_mhz', live='gpr_start_freq_mhz', group='Geometry & depth') - stop_freq_mhz: float = setting(6000.0, '_legacy_gpr_stop_freq_mhz', live='gpr_stop_freq_mhz', group='Geometry & depth') - speed_m_s: float = setting(0.0, '_legacy_gpr_speed_m_s', live='gpr_speed_m_s', group='Motion') - ignore_socket_speed_enabled: bool = setting(False, '_legacy_gpr_ignore_socket_speed_enabled', live='ignore_socket_speed', group='Motion') - look_angle_deg: float = setting(0.0, '_legacy_gpr_look_angle_deg', live='gpr_look_angle_deg', group='Motion') - apply_freq_phase_correction: bool = setting( - True, '_legacy_gpr_apply_freq_phase_correction', live='gpr_apply_freq_phase_correction', group='Motion', - ) - reference_mode: str = setting( - "frame_center", '_legacy_gpr_reference_mode', live='gpr_reference_mode', group='Motion', choices=('frame_center', 'first_tx_event'), - ) - snr_thresh: float = setting(4.5, '_legacy_gpr_snr_thresh', live='gpr_snr_thresh', group='Detection', minimum=0) - snr_comp_max: float = setting(25.0, '_legacy_gpr_snr_comp_max', live='gpr_snr_comp_max', group='Detection', minimum=0) - background_subtract_enabled: bool = setting( - True, '_legacy_gpr_background_subtract_enabled', live='gpr_background_subtract_enabled', group='Imaging', - ) - background_mean_count: int = setting(10, '_legacy_gpr_background_mean_count', live='gpr_background_mean_count', group='Imaging') - render_mode: str = setting("heatmap", '_legacy_gpr_render_mode', web='render_mode', group='Display', choices=('heatmap', 'objects_only')) - min_visible_pair_count: int = setting( - 1, '_legacy_gpr_min_visible_pair_count', live='legacy_gpr_min_visible_pair_count', group='Detection', minimum=1, - ) - visible_x_min_m: float = setting(-2.0, '_legacy_gpr_visible_x_min_m', live='gpr_visible_x_min_m', group='Visible window') - visible_x_max_m: float = setting(2.0, '_legacy_gpr_visible_x_max_m', live='gpr_visible_x_max_m', group='Visible window') - visible_z_min_m: float = setting(0.0, '_legacy_gpr_visible_z_min_m', live='gpr_visible_z_min_m', group='Visible window') - visible_z_max_m: float = setting(14.0, '_legacy_gpr_visible_z_max_m', live='gpr_visible_z_max_m', group='Visible window') - - @dataclass(slots=True) class GuiSar3dStateModel: """UI-only defaults for the 3D SAR (sar3d) live settings. - Values are the reference notebook's (Multi_antenna_3D_clean_v1.ipynb), whose + Values are the reference notebook's (Multi_antenna_demonstrate_depthnorm.ipynb), whose names the GUI shows: INPUT_IDX, F_START, MIN_DEPTH, LOOK_ANGLE_DEG, ... """ @@ -189,6 +153,15 @@ class GuiSar3dStateModel: comp_range_power: float = setting(0.3, '_sar3d_comp_range_power', live='sar3d_comp_range_power', group='3D SAR') use_compensation: bool = setting(True, '_sar3d_use_compensation', live='sar3d_use_compensation', group='3D SAR') normalize_each_frame: bool = setting(True, '_sar3d_normalize_each_frame', live='sar3d_normalize_each_frame', group='3D SAR') + pair_depth_equalize: bool = setting(True, '_sar3d_pair_depth_equalize', live='pair_depth_equalize', group='3D SAR') + background_subtract_enabled: bool = setting(True, '_sar3d_background_subtract_enabled', live='gpr_background_subtract_enabled', group='3D SAR') + background_mean_count: int = setting(10, '_sar3d_background_mean_count', live='gpr_background_mean_count', group='3D SAR', minimum=0) + x_size_m: float = setting(0.0, '_sar3d_x_size_m', live='sar3d_x_size_m', group='Scene', minimum=0) + y_min_m: float = setting(-0.5, '_sar3d_y_min_m', live='sar3d_y_min_m', group='Scene') + y_max_m: float = setting(2.0, '_sar3d_y_max_m', live='sar3d_y_max_m', group='Scene') + z_size_m: float = setting(0.0, '_sar3d_z_size_m', live='sar3d_z_size_m', group='Scene', minimum=0) + cloud_threshold: float = setting(0.65, '_sar3d_cloud_threshold', web='sar3d_cloud_threshold', group='Display', minimum=0, maximum=1) + zoom_percent: int = setting(100, '_sar3d_zoom_percent', web='sar3d_zoom_percent', group='Display', minimum=100, maximum=1000) @dataclass(slots=True) @@ -196,12 +169,11 @@ class GuiProcessingStateModel: """UI-only processing-section defaults.""" selected_mode: str = setting( - "pass_through", '_processing_mode', live='processor_mode', group='Processor', choices=('pass_through', 'bscan', 'gpr', 'legacy_gpr', 'sar3d'), + "pass_through", '_processing_mode', live='processor_mode', group='Processor', choices=('pass_through', 'bscan', 'gpr', 'sar3d'), ) pass_through: GuiPassThroughStateModel = field(default_factory=GuiPassThroughStateModel) bscan: GuiBscanStateModel = field(default_factory=GuiBscanStateModel) gpr: GuiGprStateModel = field(default_factory=GuiGprStateModel) - legacy_gpr: GuiLegacyGprStateModel = field(default_factory=GuiLegacyGprStateModel) sar3d: GuiSar3dStateModel = field(default_factory=GuiSar3dStateModel) diff --git a/python_app/models/run_config_codec.py b/python_app/models/run_config_codec.py index cf9496e..8b67b4e 100644 --- a/python_app/models/run_config_codec.py +++ b/python_app/models/run_config_codec.py @@ -343,9 +343,6 @@ def run_config_from_dict(payload: object) -> RunConfigModel: if model.preprocess.notch.taper_type not in ("cosine", "hard"): raise ValueError("preprocess.notch.taper_type must be cosine or hard") - model.gpr.relative_permittivity = _read_float( - gpr_payload, "relative_permittivity", model.gpr.relative_permittivity - ) model.gpr.tx_geometry = [] for entry in _as_list(gpr_payload.get("tx_geometry"), "gpr.tx_geometry"): entry_payload = _as_dict(entry, "gpr.tx_geometry[]") diff --git a/python_app/models/run_config_schema.py b/python_app/models/run_config_schema.py index cb86731..9a3745d 100644 --- a/python_app/models/run_config_schema.py +++ b/python_app/models/run_config_schema.py @@ -304,7 +304,6 @@ class GprRxGeometryModel: class GprModel: """Stable GPR configuration saved in run_config.json.""" - relative_permittivity: float = 1.0 tx_geometry: list[GprTxGeometryModel] = field(default_factory=list) rx_geometry: list[GprRxGeometryModel] = field(default_factory=list) diff --git a/python_app/models/run_config_validation.py b/python_app/models/run_config_validation.py index 7d67ee1..044c50c 100644 --- a/python_app/models/run_config_validation.py +++ b/python_app/models/run_config_validation.py @@ -288,8 +288,7 @@ def validate_gpr_model( ) -> None: """Validate stable GPR config against current switch dimensions. - Checks that the relative permittivity is positive and that every tx/rx - geometry entry indexes a valid, non-duplicate switch position. When ``sweep`` + Checks that every tx/rx geometry entry indexes a valid, non-duplicate switch position. When ``sweep`` is supplied (the load and GUI/save paths share this chokepoint), its bounds are validated here too so sweep failures surface alongside the GPR checks instead of as a C++ boot crash. @@ -297,9 +296,6 @@ def validate_gpr_model( if sweep is not None: validate_sweep_model(sweep) - if not math.isfinite(gpr.relative_permittivity) or float(gpr.relative_permittivity) <= 0.0: - raise ValueError("gpr.relative_permittivity must be > 0") - normalize_antenna_geometry(gpr.tx_geometry, "output_pos", output_switch_positions) normalize_antenna_geometry(gpr.rx_geometry, "input_pos", input_switch_positions) diff --git a/python_app/orchestration/live_processing_config.py b/python_app/orchestration/live_processing_config.py index 59972f3..f7c8ef0 100644 --- a/python_app/orchestration/live_processing_config.py +++ b/python_app/orchestration/live_processing_config.py @@ -21,8 +21,6 @@ class ProcessingLiveConfig: bscan_gain: float = 1.0 bscan_start_freq_mhz: float = 100.0 bscan_stop_freq_mhz: float = 8800.0 - legacy_gpr_mode: str = "point" - relative_permittivity: float | None = None tx_geometry: list[dict] | None = None rx_geometry: list[dict] | None = None gpr_input_positions: list[int] | None = None @@ -30,7 +28,6 @@ class ProcessingLiveConfig: gpr_min_depth_m: float = 2.0 gpr_max_depth_m: float = 14.0 gpr_range_comp_power: float = 0.1 - gpr_comp_power: float = 0.2 gpr_object_min_frac: float = 0.7 gpr_score_mode: str = "combined" # Backprojection intra-sweep speed-correction mode: "int_minus" (full @@ -41,28 +38,14 @@ class ProcessingLiveConfig: gpr_object_approach_min_frames: int = 3 gpr_speed_m_s: float = 0.0 gpr_look_angle_deg: float = 0.0 - # Motion-model knobs for the legacy GPR pipeline. `direction_sign` flips - # which way later events appear deeper (+1) vs shallower (-1) along Z. - # `apply_freq_phase_correction` enables intra-sweep frequency-domain phase - # compensation that fires *before* the IFFT — needed when the radar moves - # appreciably during one sweep. gpr_direction_sign: float = 1.0 - gpr_apply_freq_phase_correction: bool = True - # Anchor for the motion model's per-event `dt_ref`: 'frame_center' (default) - # or 'first_tx_event'. Mirrors Python `MOTION_CONFIG.reference_mode`. - gpr_reference_mode: str = "frame_center" - gpr_snr_thresh: float = 4.5 - gpr_snr_comp_max: float = 25.0 gpr_start_freq_mhz: float = 3000.0 gpr_stop_freq_mhz: float = 6000.0 gpr_background_subtract_enabled: bool = True gpr_background_mean_count: int = 10 gpr_remove_sidelobe_objects_enabled: bool = True - gpr_imaging_plane_y_m: float = 0.0 - # Locator filter parameter consumed by the C++ TCP locator server. Coherent BP - # objects are already finalized in the processor (no score threshold); only legacy - # GPR still thresholds, on a pair count. - legacy_gpr_min_visible_pair_count: float = 0.0 + gpr_bp_ground_depth_m: float = 1.25 + pair_depth_equalize: bool = True # Visible X/Z window (metres). The locator and the desktop plot both clip # detected objects to this window, so the socket broadcasts only what is shown. gpr_visible_x_min_m: float = -2.0 @@ -70,7 +53,7 @@ class ProcessingLiveConfig: gpr_visible_z_min_m: float = 0.0 gpr_visible_z_max_m: float = 14.0 # 3D SAR (sar3d processor) parameters, named after the reference notebook - # Multi_antenna_3D_clean_v1.ipynb (SAR_FRAME_RADIUS, BP_GROUND_DEPTH_M, ...). + # Multi_antenna_demonstrate_depthnorm.ipynb (SAR_FRAME_RADIUS, BP_GROUND_DEPTH_M, ...). # The band, depth gate, antenna selection, speed and look angle are the # shared gpr_* fields. sar3d_frame_radius: int = 0 @@ -79,6 +62,10 @@ class ProcessingLiveConfig: sar3d_comp_range_power: float = 0.3 sar3d_use_compensation: bool = True sar3d_normalize_each_frame: bool = True + sar3d_x_size_m: float = 0.0 + sar3d_y_min_m: float = -0.5 + sar3d_y_max_m: float = 2.0 + sar3d_z_size_m: float = 0.0 # When true, the C++ data_processor ignores socket-supplied vlc updates # and keeps using `gpr_speed_m_s` from this file. ignore_socket_speed: bool = False @@ -100,7 +87,7 @@ class ProcessingLiveConfig: def to_dict(self) -> dict[str, object]: """Convert live config to JSON-serializable dictionary.""" - optional_geometry = {"relative_permittivity", "tx_geometry", "rx_geometry"} + optional_geometry = {"tx_geometry", "rx_geometry"} return { name: value for name, value in asdict(self).items() if name not in optional_geometry or value is not None diff --git a/python_app/tests/test_gui_action_lifecycle.py b/python_app/tests/test_gui_action_lifecycle.py index 83c0a6e..58e8fce 100644 --- a/python_app/tests/test_gui_action_lifecycle.py +++ b/python_app/tests/test_gui_action_lifecycle.py @@ -15,7 +15,7 @@ import unittest from unittest import mock os.environ.setdefault("QT_QPA_PLATFORM", "offscreen") -from PyQt6.QtWidgets import QApplication, QCheckBox, QComboBox, QDoubleSpinBox, QLineEdit, QSpinBox # noqa: E402 +from PyQt6.QtWidgets import QApplication, QCheckBox, QComboBox, QDoubleSpinBox, QLineEdit, QSlider, QSpinBox # noqa: E402 from PyQt6.QtCore import QSignalBlocker, QTimer # noqa: E402 from python_app.gui.app_window import AppWindow # noqa: E402 @@ -532,11 +532,11 @@ class GuiProfileWriteTest(unittest.TestCase): value = not widget.isChecked() elif isinstance(widget, QComboBox): value = widget.itemText((widget.currentIndex() + 1) % widget.count()) - elif isinstance(widget, (QSpinBox, QDoubleSpinBox)): + elif isinstance(widget, (QSpinBox, QDoubleSpinBox, QSlider)): value = widget.value() + widget.singleStep() if value > widget.maximum(): value = widget.value() - widget.singleStep() - value = int(value) if isinstance(widget, QSpinBox) else round(value, widget.decimals()) + value = int(value) if isinstance(widget, (QSpinBox, QSlider)) else round(value, widget.decimals()) elif isinstance(widget, QLineEdit): value = "0:0" if binding.spec.web == "combo_filter" else "0" if "positions" in binding.path[-1] else "1" else: diff --git a/python_app/tests/test_gui_profile_codec.py b/python_app/tests/test_gui_profile_codec.py index 6df1678..bfa7eb9 100644 --- a/python_app/tests/test_gui_profile_codec.py +++ b/python_app/tests/test_gui_profile_codec.py @@ -110,14 +110,20 @@ class GuiProfileCodecTest(unittest.TestCase): with self.assertRaises(ValueError): # combos_text must be a JSON string, not a number GuiProfileModel.from_dict({"gui": {"switches": {"combos_text": 123}}}) - def test_legacy_gpr_payload_migrates_selected_mode(self) -> None: - # An old payload that selects 'gpr' but carries a legacy root gpr.mode is migrated. + def test_obsolete_inactive_settings_are_ignored(self) -> None: decoded = GuiProfileModel.from_dict({ - "gui": {"processing": {"selected_mode": "gpr"}}, - "gpr": {"relative_permittivity": 4.0, "mode": "point"}, + "gui": {"processing": {"selected_mode": "gpr", "legacy_gpr": "obsolete", + "gpr": {"algorithm": "legacy_point"}}}, + "gpr": {"relative_permittivity": "obsolete", "mode": "point"}, }) assert decoded.gui is not None - self.assertEqual(decoded.gui.processing.selected_mode, "legacy_gpr") + self.assertEqual(decoded.gui.processing.selected_mode, "gpr") + self.assertNotIn("legacy_gpr", decoded.to_dict()["gui"]["processing"]) + self.assertNotIn("relative_permittivity", decoded.to_dict()["gpr"]) + + def test_removed_processing_mode_is_rejected(self) -> None: + with self.assertRaisesRegex(ValueError, "selected_mode must be one of"): + GuiProfileModel.from_dict({"gui": {"processing": {"selected_mode": "legacy_gpr"}}}) def test_modern_gpr_selection_is_not_migrated(self) -> None: decoded = GuiProfileModel.from_dict({"gui": {"processing": {"selected_mode": "gpr"}}}) diff --git a/python_app/tests/test_live_settings_regressions.py b/python_app/tests/test_live_settings_regressions.py index 6947d78..c64dc51 100644 --- a/python_app/tests/test_live_settings_regressions.py +++ b/python_app/tests/test_live_settings_regressions.py @@ -25,11 +25,11 @@ class LiveSettingsRegressionTest(unittest.TestCase): config, gui = RunConfigModel(), GuiStateModel() config.radar.sweep.start_hz, config.radar.sweep.stop_hz = 1e6, 9.5e9 for settings in (gui.processing.bscan, gui.processing.gpr, - gui.processing.legacy_gpr, gui.processing.sar3d): + gui.processing.sar3d): settings.start_freq_mhz, settings.stop_freq_mhz = 20, 9200 section = build_processing_group(config, gui) self.addCleanup(section.group.deleteLater) - for mode in ("bscan", "gpr", "legacy_gpr", "sar3d"): + for mode in ("bscan", "gpr", "sar3d"): for bound, expected in (("start", 20), ("stop", 9200)): with self.subTest(mode=mode, bound=bound): editor = widget(section, f"_{mode}_{bound}_freq_mhz", QDoubleSpinBox) diff --git a/python_app/tests/test_live_settings_schema.py b/python_app/tests/test_live_settings_schema.py index c5beb31..cae06bc 100644 --- a/python_app/tests/test_live_settings_schema.py +++ b/python_app/tests/test_live_settings_schema.py @@ -11,7 +11,7 @@ class LiveSettingsSchemaCoverageTest(unittest.TestCase): def test_bindings_cover_live_config_and_select_each_page_once(self) -> None: config_names = {field.name for field in fields(ProcessingLiveConfig)} runtime = {"reprocess_current_result", "history_clear_before_ns", "history_removed", "settings_revision"} - geometry = {"relative_permittivity", "tx_geometry", "rx_geometry"} + geometry = {"tx_geometry", "rx_geometry"} modes = next(binding.spec.choices for binding in GUI_FIELDS if binding.spec.live == "processor_mode") widget_names = [binding.spec.widget for binding in GUI_FIELDS if binding.spec.widget] self.assertEqual(len(widget_names), len(set(widget_names))) diff --git a/python_app/tests/test_processing.py b/python_app/tests/test_processing.py index 7035f48..491adfb 100644 --- a/python_app/tests/test_processing.py +++ b/python_app/tests/test_processing.py @@ -1,4 +1,4 @@ -"""GPR object filtering, B-scan transforms and live-setting normalization.""" +"""GPR payload lookup, B-scan transforms and live-setting normalization.""" from __future__ import annotations @@ -26,8 +26,6 @@ from python_app.models.dataset_model import ( ) from python_app.gui.result_payloads import ( collection_payload_by_name, - collection_payloads_by_prefix, - filter_object_rows, gpr_object_rows, ) from python_app.orchestration.live_processing_config import ProcessingLiveConfig @@ -43,32 +41,16 @@ def _collection(payloads: Iterable[ResultPayload]) -> ResultCollection: class CollectionLookupTest(unittest.TestCase): def test_payload_by_name_matches_name_and_optional_kind(self) -> None: - col = _collection([_table("gpr_points", [[0, 0, 1]]), _table("gpr_region_centers", [[1, 1, 2, 9]])]) + col = _collection([_table("gpr_points", [[0, 0, 1]])]) self.assertIs(collection_payload_by_name(col, "gpr_points"), col.collection_payloads[0]) self.assertIsNone(collection_payload_by_name(col, "missing")) self.assertIsNone(collection_payload_by_name(col, "gpr_points", kind=3)) # wrong kind - def test_payloads_by_prefix_returns_all_in_order(self) -> None: - col = _collection([ - ResultPayload(processing_name="gpr_region_mask_0", kind=3, image=np.zeros((1, 1), dtype=np.float32)), - ResultPayload(processing_name="gpr_region_mask_1", kind=3, image=np.zeros((1, 1), dtype=np.float32)), - _table("gpr_points", [[0, 0, 1]]), - ]) - masks = collection_payloads_by_prefix(col, "gpr_region_mask_", kind=3) - self.assertEqual([p.processing_name for p in masks], ["gpr_region_mask_0", "gpr_region_mask_1"]) - self.assertEqual(collection_payloads_by_prefix(col, "nope_"), []) - - class GprObjectRowsTest(unittest.TestCase): def test_extracts_first_three_columns_of_points(self) -> None: rows = gpr_object_rows(_collection([_table("gpr_points", [[1, 2, 3], [4, 5, 6]])])) self.assertTrue(np.array_equal(rows, np.array([[1, 2, 3], [4, 5, 6]], dtype=np.float32))) - def test_falls_back_to_region_centers_and_drops_extra_columns(self) -> None: - # region_centers is [x, z, score, pixel_count]; only the first 3 cols are used. - rows = gpr_object_rows(_collection([_table("gpr_region_centers", [[1, 2, 3, 99]])])) - self.assertTrue(np.array_equal(rows, np.array([[1, 2, 3]], dtype=np.float32))) - def test_empty_when_no_object_payloads(self) -> None: self.assertEqual(gpr_object_rows(_collection([_table("gpr_accumulator", [[1, 2, 3]])])).shape, (0, 3)) @@ -76,39 +58,6 @@ class GprObjectRowsTest(unittest.TestCase): self.assertEqual(gpr_object_rows(_collection([_table("gpr_points", [[1, 2]])])).shape, (0, 3)) -class FilterObjectRowsTest(unittest.TestCase): - BOUNDS = {"x_bounds": (-2.0, 2.0), "z_bounds": (0.0, 10.0)} - - def _filter(self, rows: Sequence[Sequence[float]], *, min_score: float) -> np.ndarray: - return filter_object_rows(np.asarray(rows, dtype=np.float32), min_score=min_score, - **self.BOUNDS) - - def test_keeps_in_window_and_at_or_above_threshold(self) -> None: - out = self._filter([[0.0, 5.0, 0.5], [1.0, 1.0, 0.9]], min_score=0.5) # score==threshold kept (inclusive) - self.assertEqual(out.shape[0], 2) - - def test_drops_below_threshold(self) -> None: - out = self._filter([[0.0, 5.0, 0.4]], min_score=0.5) - self.assertEqual(out.shape[0], 0) - - def test_window_bounds_are_inclusive(self) -> None: - out = self._filter([[2.0, 10.0, 1.0], [-2.0, 0.0, 1.0]], min_score=0.0) # exactly on each edge - self.assertEqual(out.shape[0], 2) - - def test_drops_outside_window(self) -> None: - out = self._filter([[2.001, 5.0, 1.0], [0.0, 10.001, 1.0]], min_score=0.0) - self.assertEqual(out.shape[0], 0) - - def test_drops_non_finite_rows(self) -> None: - out = self._filter([[np.nan, 5.0, 1.0], [0.0, np.inf, 1.0], [0.0, 5.0, 1.0]], min_score=0.0) - self.assertEqual(out.shape[0], 1) - - def test_legacy_pair_count_threshold_uses_same_compare(self) -> None: - # legacy GPR passes an integer pair-count threshold; the >= compare is identical. - out = self._filter([[0.0, 5.0, 3.0], [0.0, 5.0, 2.0]], min_score=3) - self.assertTrue(np.array_equal(out, np.array([[0.0, 5.0, 3.0]], dtype=np.float32))) - - class BscanTransformTest(unittest.TestCase): def test_mean_ascan_subtraction(self) -> None: history = deque([np.array([1.0, 1.0], dtype=np.float32), np.array([3.0, 3.0], dtype=np.float32)]) diff --git a/python_app/tests/test_processor_history.py b/python_app/tests/test_processor_history.py index db6f79d..e2a4070 100644 --- a/python_app/tests/test_processor_history.py +++ b/python_app/tests/test_processor_history.py @@ -130,6 +130,27 @@ class _ProcessorTestCase(unittest.TestCase): class FedBackSweepsTest(_ProcessorTestCase): + def test_mean_history_is_not_limited_by_transport_ring_capacity(self) -> None: + self._live.processor_mode = "gpr" + self._live.gpr_input_positions = self._live.gpr_output_positions = [0] + self._live.gpr_background_subtract_enabled = True + self._live.gpr_background_mean_count = 60 + self._live.gpr_min_depth_m, self._live.gpr_max_depth_m = .1, 1.5 + self._live.gpr_bp_ground_depth_m = 0. + self._live_writer.write(self._live) + start = time.monotonic_ns() + for cid in range(1, 61): + sweep = _sweep(cid, start + cid*1_000_000) + sweep.traces[0].s21 *= (-1.)**np.arange(len(sweep.traces[0].s21)) + self._push([sweep]) + self._result_for_revision(0, cid) + current = _sweep(61, start + 61_000_000) + current.traces[0].s21[:] = 30.5 * (-1.)**np.arange(len(current.traces[0].s21)) + self._push([current]) + result = self._result_for_revision(0, 61) + image = next(payload.image for payload in result.collection_payloads if payload.processing_name == "gpr_accumulator") + self.assertFalse(np.any(image)) # Exact cancellation of the 60 preceding spectra (1…60). + def test_remove_targets_exact_sweep_with_reused_ids_and_preserves_missing_targets(self) -> None: start_ns = time.monotonic_ns() sweeps = [_sweep(cid, start_ns + index * 1_000_000) diff --git a/python_app/tests/test_review_gui_regressions.py b/python_app/tests/test_review_gui_regressions.py index a62f609..0d05f16 100644 --- a/python_app/tests/test_review_gui_regressions.py +++ b/python_app/tests/test_review_gui_regressions.py @@ -150,7 +150,7 @@ class ReviewGuiRegressions(unittest.TestCase): deadline.assert_called_once() window = isolated_window(self.addCleanup) window._pipeline.timer.stop() - for mode in ("gpr", "legacy_gpr", "sar3d"): + for mode in ("gpr", "sar3d"): field = widget(window._view.processing, f"_{mode}_input_positions_input", QLineEdit) validator = field.validator() self.assertIsNotNone(validator) @@ -161,7 +161,8 @@ class ReviewGuiRegressions(unittest.TestCase): # A trailing comma during editing is harmless to the shared CSV parser. for table in (np.empty((0, 0)), np.array([[99.]])): result = ResultCollection(1, 2, collection_payloads=[ResultPayload(processing_name="sar3d_status", kind=4, table=table)]) - self.assertFalse(window._plots.sar3d.draw(result)) + self.assertFalse(window._plots.sar3d.draw(result, window._config.processing_settings().sar3d, + window._config.geometry())) self.assertIn("error", window._plots.sar3d.status_label.text().lower()) def test_stop_or_close_cancels_paused_actions_without_error_dialog(self) -> None: diff --git a/python_app/tests/test_run_config_codec.py b/python_app/tests/test_run_config_codec.py index f4e23a2..7e9cd0b 100644 --- a/python_app/tests/test_run_config_codec.py +++ b/python_app/tests/test_run_config_codec.py @@ -363,7 +363,7 @@ class RingValidationTest(unittest.TestCase): class GprValidationTest(unittest.TestCase): @staticmethod def _gpr(**overrides: object) -> GprModel: - model = GprModel(relative_permittivity=4.0) + model = GprModel() for key, value in overrides.items(): setattr(model, key, value) return model @@ -371,10 +371,6 @@ class GprValidationTest(unittest.TestCase): def _validate(self, gpr: GprModel) -> None: validate_gpr_model(gpr, input_switch_positions=4, output_switch_positions=2) - def test_permittivity_must_be_positive(self) -> None: - with self.assertRaisesRegex(ValueError, "relative_permittivity"): - self._validate(self._gpr(relative_permittivity=0.0)) - def test_tx_output_pos_out_of_range(self) -> None: with self.assertRaisesRegex(ValueError, "output_pos is out of range"): self._validate(self._gpr(tx_geometry=[GprTxGeometryModel(output_pos=5, x_m=0.0, y_m=0.0, z_m=0.0)])) diff --git a/python_app/tests/test_sar3d_plot.py b/python_app/tests/test_sar3d_plot.py index dce36d4..42cf572 100644 --- a/python_app/tests/test_sar3d_plot.py +++ b/python_app/tests/test_sar3d_plot.py @@ -9,11 +9,13 @@ from unittest import mock import numpy as np os.environ.setdefault("QT_QPA_PLATFORM", "offscreen") -from PyQt6.QtCore import Qt, QSignalBlocker # noqa: E402 -from PyQt6.QtWidgets import QApplication, QCheckBox, QComboBox, QDoubleSpinBox # noqa: E402 +from PyQt6.QtCore import QPointF, QRectF, Qt, QSignalBlocker # noqa: E402 +from PyQt6.QtWidgets import QApplication, QCheckBox, QComboBox, QDoubleSpinBox, QSlider # noqa: E402 from python_app.gui.controllers.renderers.sar3d import SAR3D_STATES, SAR3D_VIEWS, Sar3dRenderer # noqa: E402 from python_app.models.dataset_model import ComboKey, ResultCollection, ResultPayload, TraceData # noqa: E402 +from python_app.models.gui_profile_schema import GuiSar3dStateModel # noqa: E402 +from python_app.models.run_config_schema import GprModel, GprTxGeometryModel, GprRxGeometryModel # noqa: E402 from python_app.tests.gui_test_support import widget @@ -47,12 +49,13 @@ class Sar3dPlotTest(unittest.TestCase): def setUp(self) -> None: self.renderer = Sar3dRenderer() + self.settings, self.geometry = GuiSar3dStateModel(), GprModel() self.addCleanup(self.renderer.widget.close) def test_all_four_projections_have_visible_images_and_correct_peak_axes(self) -> None: - self.assertTrue(self.renderer.draw(_projections())) + self.assertTrue(self.renderer.draw(_projections(), self.settings, self.geometry)) self.assertEqual(self.renderer.status_label.text(), - "SAR3D: ready — CUDA | max X +0.50, Y -0.20, Z 2.00 m | level 2") + "SAR3D: ready — CUDA | max X +0.50, Y -0.20, Z 2.00 m | level 2 | zoom 100%") positions = {"sar3d_front": (0.5, -0.2), "sar3d_side": (2, -0.2), "sar3d_top": (0.5, 2), "sar3d_oblique": (0.5, -0.2)} for name, position in positions.items(): @@ -60,27 +63,83 @@ class Sar3dPlotTest(unittest.TestCase): self.assertTrue(image.isVisible(), name) self.assertGreater(image.boundingRect().width(), 0) self.assertGreater(image.boundingRect().height(), 0) - xx, yy = self.renderer.peaks[name].getData() - np.testing.assert_allclose([xx[0], yy[0]], position) + xx, yy = self.renderer.targets[name].getData() + assert xx is not None and yy is not None + np.testing.assert_allclose([.5*(xx.min()+xx.max()), .5*(yy.min()+yy.max())], position) + self.assertAlmostEqual(float(np.ptp(xx)), float(np.ptp(yy))) self.assertTrue(all(label.isVisible() for label in self.renderer.axis_labels.values())) self.renderer.clear() self.assertEqual(self.renderer.status_label.text(), "SAR3D: waiting for measurements") def test_unavailable_warning_clears_previous_images_and_labels_and_shows_reason(self) -> None: - self.renderer.draw(_projections()) + self.renderer.draw(_projections(), self.settings, self.geometry) message = "CUDA unavailable: драйвер " - self.assertFalse(self.renderer.draw(_collection(_status(3, message)))) + self.assertFalse(self.renderer.draw(_collection(_status(3, message)), self.settings, self.geometry)) self.assertEqual(self.renderer.status_label.text(), f"SAR3D: warning — {message}") self.assertEqual(self.renderer.status_label.textFormat(), Qt.TextFormat.PlainText) self.assertTrue(all(not image.isVisible() for image in self.renderer.images.values())) self.assertTrue(all(not label.isVisible() for label in self.renderer.axis_labels.values())) def test_no_signal_and_missing_projections_explain_empty_panels(self) -> None: - self.assertFalse(self.renderer.draw(_collection(_status(1, "No signal in selected range")))) + self.assertFalse(self.renderer.draw(_collection(_status(1, "No signal in selected range")), self.settings, self.geometry)) self.assertIn("no signal", self.renderer.status_label.text()) - self.assertFalse(self.renderer.draw(_collection())) + self.assertFalse(self.renderer.draw(_collection(), self.settings, self.geometry)) self.assertEqual(self.renderer.status_label.text(), "SAR3D: waiting for projection data") + def test_response_size_antenna_positions_and_zoom_use_scene_coordinates(self) -> None: + result = _projections() + result.collection_payloads.append(ResultPayload("sar3d_target_width", 4, table=np.array([[.4, .6, .8]]))) + self.geometry = GprModel(tx_geometry=[GprTxGeometryModel(0, .7, .3, .4)], + rx_geometry=[GprRxGeometryModel(0, -.7, -.3, .2)]) + self.settings.look_angle_deg = 30 + self.settings.bp_ground_depth_m = 1.25 + self.renderer.widget.resize(800, 600) + self.renderer.widget.show() + self.renderer.draw(result, self.settings, self.geometry) + self.app.processEvents() + baseline = {} + expected_antenna = [.7, 1.25 + .3*np.cos(np.pi/6) - .4*.5, .3*.5 + .4*np.cos(np.pi/6)] + for name, (_, _, _, axes) in SAR3D_VIEWS.items(): + xx, yy = self.renderer.targets[name].getData() + assert xx is not None and yy is not None + side = .6 if name in {"sar3d_front", "sar3d_oblique"} else .8 + self.assertAlmostEqual(float(np.ptp(xx)), side) + if axes is not None: + tx, _ = self.renderer.antennas[name] + ax, ay = tx.getData() + np.testing.assert_allclose([ax[0], ay[0]], np.asarray(expected_antenna)[list(axes)]) + view = self.renderer.plots[name].getViewBox() + for marker in self.renderer.antennas[name]: + ax, ay = marker.getData() + for x, y in zip(ax, ay): + point = view.mapViewToScene(QPointF(float(x), float(y))) + self.assertTrue(view.sceneBoundingRect().contains(QRectF(point.x()-7, point.y()-7, 14, 14)), name) + rect = self.renderer.plots[name].getViewBox().targetRect() + baseline[name] = np.array([[rect.left(), rect.right()], [rect.top(), rect.bottom()]]) + self.assertNotIn("sar3d_oblique", self.renderer.antennas) + self.settings.zoom_percent = 250 + self.renderer.draw(result, self.settings, self.geometry) + self.app.processEvents() + for name, ranges in baseline.items(): + rect = self.renderer.plots[name].getViewBox().targetRect() + zoomed = np.array([[rect.left(), rect.right()], [rect.top(), rect.bottom()]]) + np.testing.assert_allclose(np.ptp(zoomed, axis=1), np.ptp(ranges, axis=1)/2.5) + xx, yy = self.renderer.targets[name].getData() + assert xx is not None and yy is not None + peak = np.array([.5*(xx.min()+xx.max()), .5*(yy.min()+yy.max())]) + self.assertLess(np.linalg.norm(zoomed.mean(axis=1)-peak), np.linalg.norm(ranges.mean(axis=1)-peak)) + + def test_cloud_threshold_masks_weak_response_without_changing_levels(self) -> None: + result = _projections() + for payload in result.collection_payloads: + if payload.kind == 3: + payload.image.flat[:2] = [.2, .8] + for threshold, expected in ((.65, [0, .8, 1]), (1., [0, 0, 1])): + self.settings.cloud_threshold = threshold + self.renderer.draw(result, self.settings, self.geometry) + for image in self.renderer.images.values(): + np.testing.assert_allclose(np.asarray(image.image).flat[:3], expected) + np.testing.assert_allclose(np.asarray(image.levels), [0, 1]) class Sar3dIntegrationTest(unittest.TestCase): @@ -100,6 +159,30 @@ class Sar3dIntegrationTest(unittest.TestCase): self.view._plots.required_revision = None self.renderer = self.view._plots.sar3d + def test_display_controls_redraw_without_processing_revision_or_live_fields(self) -> None: + result = _projections() + self.view._history.record_result(result) + revision = self.view._live.revision + widget(self.view._view.processing, "_sar3d_zoom_percent", QSlider).setValue(247) + widget(self.view._view.processing, "_sar3d_cloud_threshold", QDoubleSpinBox).setValue(1.) + self.assertEqual(self.view._live.revision, revision) + self.assertIn("zoom 247%", self.renderer.status_label.text()) + live = self.view._live.build().to_dict() + self.assertNotIn("sar3d_zoom_percent", live) + self.assertNotIn("sar3d_cloud_threshold", live) + schema = {entry["name"]: entry for entry in self.view._config.web_schema()} + self.assertEqual(schema["sar3d_zoom_percent"]["value"], 247) + self.assertEqual(schema["sar3d_zoom_percent"]["step"], 1) + self.assertEqual(schema["sar3d_cloud_threshold"]["value"], 1.) + + def test_inactive_sar_scene_is_validated_before_publishing_gpr_settings(self) -> None: + with QSignalBlocker(self.mode): + self.mode.setCurrentText("gpr") + with QSignalBlocker(widget(self.view._view.processing, "_sar3d_y_min_m", QDoubleSpinBox)): + widget(self.view._view.processing, "_sar3d_y_min_m", QDoubleSpinBox).setValue(3.) + with self.assertRaisesRegex(ValueError, "SAR_Y_RANGE_M"): + self.view._config.build_live_config() + def test_capture_preview_keeps_sar_page_and_uses_saved_trace_settings(self) -> None: self.view._plots.set_mode("sar3d") fields = self.view._view.processing diff --git a/python_app/tests/test_socket_speed_controls.py b/python_app/tests/test_socket_speed_controls.py index b387850..fac0678 100644 --- a/python_app/tests/test_socket_speed_controls.py +++ b/python_app/tests/test_socket_speed_controls.py @@ -2,7 +2,7 @@ The data_processor uses the speed a socket client streams unless the live config sets `ignore_socket_speed`; then `gpr_speed_m_s` from the GUI field applies. Each -GPR-like page (gpr, legacy_gpr, sar3d) has its own checkbox and field, and the +GPR-like page (gpr, sar3d) has its own checkbox and field, and the field is editable only while it is the one in force. """ @@ -31,7 +31,7 @@ class SocketSpeedControlsTest(unittest.TestCase): def test_socket_authority_controls_editability_and_published_speed(self) -> None: assert _window is not None section = _window._view.processing - for mode, (checkbox, speed) in zip(("gpr", "legacy_gpr", "sar3d"), section.speed_controls, strict=True): + for mode, (checkbox, speed) in zip(("gpr", "sar3d"), section.speed_controls, strict=True): section.mode.setCurrentText(mode) speed.setValue(1.25) for ignore in (True, False): diff --git a/python_app/webui/static/app.js b/python_app/webui/static/app.js index 5b74add..097fe9a 100644 --- a/python_app/webui/static/app.js +++ b/python_app/webui/static/app.js @@ -218,7 +218,7 @@ function makeFieldRow(entry) { el.value = String(entry.value); } else if (entry.kind === "int" || entry.kind === "float") { el = document.createElement("input"); - el.type = "number"; + el.type = entry.slider ? "range" : "number"; if (entry.min != null) el.min = entry.min; if (entry.max != null) el.max = entry.max; el.step = entry.kind === "float" ? (entry.step || "any") : (entry.step || 1); diff --git a/run_config.json b/run_config.json index 38a241b..64adb95 100644 --- a/run_config.json +++ b/run_config.json @@ -184,7 +184,6 @@ } }, "gpr": { - "relative_permittivity": 1.0, "tx_geometry": [ { "output_pos": 0, diff --git a/run_config_examples/run_config_librevna.example.json b/run_config_examples/run_config_librevna.example.json index 9fdaf42..14e8995 100644 --- a/run_config_examples/run_config_librevna.example.json +++ b/run_config_examples/run_config_librevna.example.json @@ -142,7 +142,6 @@ } }, "gpr": { - "relative_permittivity": 1.0, "tx_geometry": [ { "output_pos": 0, diff --git a/run_config_examples/run_config_librevna_multi.example.json b/run_config_examples/run_config_librevna_multi.example.json index 703f024..2d8b8bc 100644 --- a/run_config_examples/run_config_librevna_multi.example.json +++ b/run_config_examples/run_config_librevna_multi.example.json @@ -145,7 +145,6 @@ } }, "gpr": { - "relative_permittivity": 1.0, "tx_geometry": [ { "output_pos": 0, diff --git a/run_config_examples/run_config_librevna_multi.example2.json b/run_config_examples/run_config_librevna_multi.example2.json index 9d48a62..15d0e94 100644 --- a/run_config_examples/run_config_librevna_multi.example2.json +++ b/run_config_examples/run_config_librevna_multi.example2.json @@ -145,7 +145,6 @@ } }, "gpr": { - "relative_permittivity": 1.0, "tx_geometry": [ { "output_pos": 0, diff --git a/run_config_examples/run_config_librevna_multi_per_point.example.json b/run_config_examples/run_config_librevna_multi_per_point.example.json index 4411b64..37e6c45 100644 --- a/run_config_examples/run_config_librevna_multi_per_point.example.json +++ b/run_config_examples/run_config_librevna_multi_per_point.example.json @@ -188,7 +188,6 @@ } }, "gpr": { - "relative_permittivity": 1.0, "tx_geometry": [ { "output_pos": 0, @@ -326,31 +325,6 @@ "visible_z_min_m": 0.0, "visible_z_max_m": 8.0 }, - "legacy_gpr": { - "mode": "point", - "input_positions": "0,1,2,3", - "output_positions": "0,1,2,3", - "min_depth_m": 2.0, - "max_depth_m": 14.0, - "comp_power": 0.2, - "start_freq_mhz": 3000.0, - "stop_freq_mhz": 6000.0, - "speed_m_s": 0.0, - "ignore_socket_speed_enabled": false, - "look_angle_deg": 0.0, - "apply_freq_phase_correction": true, - "reference_mode": "frame_center", - "snr_thresh": 4.5, - "snr_comp_max": 25.0, - "background_subtract_enabled": true, - "background_mean_count": 10, - "render_mode": "heatmap", - "min_visible_pair_count": 1, - "visible_x_min_m": -2.91, - "visible_x_max_m": 2.91, - "visible_z_min_m": 0.0, - "visible_z_max_m": 14.0 - }, "sar3d": { "input_positions": "", "output_positions": "", diff --git a/run_config_examples/run_config_librevna_multi_per_point_lowband.example.json b/run_config_examples/run_config_librevna_multi_per_point_lowband.example.json index e671c66..4d160cc 100644 --- a/run_config_examples/run_config_librevna_multi_per_point_lowband.example.json +++ b/run_config_examples/run_config_librevna_multi_per_point_lowband.example.json @@ -184,7 +184,6 @@ } }, "gpr": { - "relative_permittivity": 1.0, "tx_geometry": [ { "output_pos": 0, diff --git a/run_config_examples/run_config_rph_radar.example.json b/run_config_examples/run_config_rph_radar.example.json index 1891de0..aadbe72 100644 --- a/run_config_examples/run_config_rph_radar.example.json +++ b/run_config_examples/run_config_rph_radar.example.json @@ -175,7 +175,6 @@ } }, "gpr": { - "relative_permittivity": 1.0, "tx_geometry": [ { "output_pos": 0, diff --git a/run_config_examples/run_config_rph_radar_rfg.example.json b/run_config_examples/run_config_rph_radar_rfg.example.json index 20ecda8..2894727 100644 --- a/run_config_examples/run_config_rph_radar_rfg.example.json +++ b/run_config_examples/run_config_rph_radar_rfg.example.json @@ -102,7 +102,6 @@ } }, "gpr": { - "relative_permittivity": 1.0, "tx_geometry": [ { "output_pos": 0, diff --git a/run_config_examples/run_config_simulator.example.json b/run_config_examples/run_config_simulator.example.json index b7aef83..97ba7d9 100644 --- a/run_config_examples/run_config_simulator.example.json +++ b/run_config_examples/run_config_simulator.example.json @@ -142,7 +142,6 @@ } }, "gpr": { - "relative_permittivity": 1.0, "tx_geometry": [ { "output_pos": 0, @@ -265,31 +264,6 @@ "visible_x_max_m": 2.0, "visible_z_min_m": 0.0, "visible_z_max_m": 14.0 - }, - "legacy_gpr": { - "mode": "point", - "input_positions": "0,1,2,3", - "output_positions": "0,1", - "min_depth_m": 2.0, - "max_depth_m": 14.0, - "comp_power": 0.2, - "start_freq_mhz": 3000.0, - "stop_freq_mhz": 6000.0, - "speed_m_s": 0.0, - "ignore_socket_speed_enabled": false, - "look_angle_deg": 0.0, - "apply_freq_phase_correction": true, - "reference_mode": "frame_center", - "snr_thresh": 4.5, - "snr_comp_max": 25.0, - "background_subtract_enabled": true, - "background_mean_count": 10, - "render_mode": "heatmap", - "min_visible_pair_count": 1, - "visible_x_min_m": -2.0, - "visible_x_max_m": 2.0, - "visible_z_min_m": 0.0, - "visible_z_max_m": 14.0 } }, "data_actions": { diff --git a/run_configs/run_config.json b/run_configs/run_config.json index 6924e91..bc24489 100644 --- a/run_configs/run_config.json +++ b/run_configs/run_config.json @@ -142,7 +142,6 @@ } }, "gpr": { - "relative_permittivity": 1.0, "tx_geometry": [ { "output_pos": 0, @@ -265,31 +264,6 @@ "visible_x_max_m": 2.0, "visible_z_min_m": 0.0, "visible_z_max_m": 14.0 - }, - "legacy_gpr": { - "mode": "point", - "input_positions": "0,1,2,3", - "output_positions": "0,1", - "min_depth_m": 2.0, - "max_depth_m": 14.0, - "comp_power": 0.2, - "start_freq_mhz": 3000.0, - "stop_freq_mhz": 6000.0, - "speed_m_s": 0.0, - "ignore_socket_speed_enabled": false, - "look_angle_deg": 0.0, - "apply_freq_phase_correction": true, - "reference_mode": "frame_center", - "snr_thresh": 4.5, - "snr_comp_max": 25.0, - "background_subtract_enabled": true, - "background_mean_count": 10, - "render_mode": "heatmap", - "min_visible_pair_count": 1, - "visible_x_min_m": -2.0, - "visible_x_max_m": 2.0, - "visible_z_min_m": 0.0, - "visible_z_max_m": 14.0 } }, "data_actions": {