diff --git a/Makefile b/Makefile index 0578cf9..80eba5c 100644 --- a/Makefile +++ b/Makefile @@ -25,6 +25,7 @@ endif # Rewritten only when the CUDA setup changes, so installing (or removing) CUDA # rebuilds the sar3d GPU code instead of keeping a stale build. SAR3D_CUDA_STAMP := $(BUILD_DIR)/sar3d_cuda.flags +SAR3D_HOST_FLAGS := -ffp-contract=off $(shell mkdir -p $(BUILD_DIR) && echo '$(SAR3D_CXXFLAGS) $(SAR3D_LDFLAGS)' | cmp -s - $(SAR3D_CUDA_STAMP) || echo '$(SAR3D_CXXFLAGS) $(SAR3D_LDFLAGS)' > $(SAR3D_CUDA_STAMP)) INCLUDES := \ @@ -57,6 +58,7 @@ PROCESSOR_SOURCES := \ data_acq_and_processing/processing/processors/src/gpr_processor.cpp \ data_acq_and_processing/processing/processors/src/passthrough_processor.cpp \ data_acq_and_processing/processing/processors/src/sar3d_processor.cpp \ + data_acq_and_processing/processing/processors/src/sar3d_clean.cpp \ data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp \ data_acq_and_processing/processing/data_processor/src/processing_live_config.cpp \ data_acq_and_processing/processing/data_processor/src/data_processor.cpp \ @@ -122,13 +124,14 @@ check-native: $(NATIVE_BUILD_DIR)/native_regressions $(BIN_DIR)/collector_regres 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_clean.cpp \ 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. +# The sar3d processor on the real CUDA backend against the notebook reference; without a +# usable GPU the test checks the warning and the box rules and skips the rest. check-sar3d: $(NATIVE_BUILD_DIR)/backprojection_fixture BACKPROJECTION_FIXTURE="$(abspath $<)" $(PYTHON) data_acq_and_processing/tests/test_backprojection.py Sar3dTest @@ -139,6 +142,7 @@ $(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) +$(NATIVE_BUILD_DIR)/data_acq_and_processing/processing/processors/src/sar3d_%.o: NATIVE_CXXFLAGS += $(SAR3D_HOST_FLAGS) $(NATIVE_BUILD_DIR)/data_acq_and_processing/processing/processors/src/sar3d_gpu.o: $(SAR3D_CUDA_STAMP) $(NATIVE_BUILD_DIR)/native_regressions: $(NATIVE_TEST_OBJS) @@ -159,6 +163,9 @@ $(BIN_DIR)/data_processor: $(DATA_PROCESSOR_OBJS) $(SAR3D_CUDA_STAMP) $(BUILD_DIR)/data_acq_and_processing/processing/processors/src/sar3d_gpu.o: CXXFLAGS += $(SAR3D_CXXFLAGS) $(BUILD_DIR)/data_acq_and_processing/processing/processors/src/sar3d_gpu.o: $(SAR3D_CUDA_STAMP) +# The sar3d host code breaks the same ties as NumPy (grid corners on the depth gate, range +# rings), so it is not contracted into fused multiply-adds where the target has them. +$(BUILD_DIR)/data_acq_and_processing/processing/processors/src/sar3d_%.o: CXXFLAGS += $(SAR3D_HOST_FLAGS) # Keep source-relative paths so moving a source also replaces its dependencies. $(RPH_BUILD_DIR)/%.o: %.cpp diff --git a/README.md b/README.md index 8b5cf2a..d7c47d5 100644 --- a/README.md +++ b/README.md @@ -71,13 +71,13 @@ controls, processing settings and a copy of the live plot through an embedded we - **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` +- **sar3d** — 3D SAR of the reference notebook `Multi_depthnorm_clean.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 set by the selected antennas with a 2 m X margin and height −0.5…2 m above the ground; its range extent follows - `MIN_DEPTH`/`MAX_DEPTH`, and the voxels are cubes with edge `SAR_VOXEL_M` - (default 2 cm). The radar motion is corrected from each trace's capture + `MIN_DEPTH`/`MAX_DEPTH`, and the voxel steps `SAR_VOXEL_M` X, Y and Z default to the + notebook's 2.5, 2 and 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 @@ -100,7 +100,27 @@ 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. +reflector's physical dimensions. The square is left out when CLEAN objects are shown +or when the background gate empties the frame. + +**BG_GATE_ENABLE** shows a frame relative to its background, the median of the volume +over its working voxels: below 6 backgrounds the frame is empty, at 6 its maximum +appears at the cloud threshold, and from 8.5 the scale is full. **SAR_CLEAN_ENABLE** +runs the notebook's CLEAN on a one-frame aperture: always with `SAR_FRAME_RADIUS` 0, +otherwise while no previous frame is available, as after Start. It subtracts the strongest point response +from the A-scans of every pair, step by step on a 4 × 10 × 4 cm grid over the whole scene +region, until the residual is background. The search ignores `SAR_X_SIZE_M` and +`SAR_Z_SIZE_M`, so reflectors just outside a narrowed scene are subtracted too. Each object is then imaged apart and drawn in its own 0–1 +scale over the residual, so a weak object is not hidden by the sidelobes of a strong +one. **SAR_MARK_ENABLE** draws boxes around groups of CLEAN objects: red for very low +objects (like a mine), orange on the ground, blue above it, dashed grey at the scene +boundary, where the response may lie partly outside. The legend is in the title of the +3D view. Without boxes the objects are numbered with their level in backgrounds. All +three switches are on by default, like in the notebook, and the status line reports the +background ratio and the CLEAN result. The gate and CLEAN are processing settings; boxes +only redraw. Their thresholds and rules are the notebook's fixed values. Both cost time: +on a GTX 1060 with the default scene of 10 million voxels, the gate adds about 0.1 s per +frame and CLEAN about 0.4 s more. The desktop uses a light theme with white plot canvases. **Mode** stays below the **Actions** buttons while the settings sections scroll; the settings panel and @@ -222,7 +242,7 @@ The preprocessor still requires S21 calibration/reference sets. Start the GUI: Open **Preprocessing…**, use **Create Neutral S21 Sets**, then **Start**. This makes identity calibration + zero reference, preserving the simulated S21 samples. -Finer voxels increase SAR3D work as 1/`SAR_VOXEL_M`³. Increase the voxel size when +SAR3D work grows inversely with the product of the three `SAR_VOXEL_M` steps. Increase them when a faster preview matters more than spatial resolution. Processing failures appear above the plots and a failed single capture finishes with its explanation while retaining the acquired data. @@ -469,19 +489,24 @@ physical timing. Real devices, GPIO and the CUDA backend need their own hardware `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, 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 +maximum projections of the notebook's displayed volume, its peak/level, response +widths and background ratio for ten cases. Six cover compensation/normalization, +equalization, background, scene and timing. Four are frames of point reflectors in noise +for CLEAN, which also compare the objects and their boxes: with and without the background +gate, with a scene where the stop at each range's background adds steps, and with a +narrowed scene. The box rules are also checked +one by one on synthetic objects. The expected values come from unchanged functions in cells 1 and 5 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): +status for invalid input, a volume beyond float32, a result too large for the ring slot, +an empty and a very weak signal, and that the views are normalized at their peak on +finer grids. Without a usable GPU it checks only the CUDA warning and the box rules and +skips the rest; the notebook itself is not needed to run it. To regenerate the references, +pass the original notebook and the Horns script (both are kept outside the repository): ```bash .venv/bin/python data_acq_and_processing/tests/generate_backprojection_reference.py \ - --notebook /path/to/Multi_antenna_demonstrate_depthnorm.ipynb \ + --notebook /path/to/Multi_depthnorm_clean.ipynb \ --gpr-script /path/to/Horns_motion_4x4_ground.py ``` 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 e3d2def..ab18c85 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 @@ -60,12 +60,14 @@ 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_demonstrate_depthnorm.ipynb (SAR_FRAME_RADIUS, BP_GROUND_DEPTH_M, ...). + // Multi_depthnorm_clean.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; - // Edge of the cubic voxel: the notebook's SAR_VOXEL_M with one value for X, Y, Z. - float sar3d_voxel_m = 0.02F; + // SAR_VOXEL_M: the voxel steps along X, Y and Z. + float sar3d_voxel_x_m = 0.025F; + float sar3d_voxel_y_m = 0.02F; + float sar3d_voxel_z_m = 0.02F; float sar3d_bp_ground_depth_m = 1.25F; float sar3d_comp_range_power = 0.3F; bool sar3d_use_compensation = true; @@ -74,6 +76,8 @@ struct ProcessingLiveConfig { float sar3d_y_min_m = -0.5F; float sar3d_y_max_m = 2.0F; float sar3d_z_size_m = 0.0F; + bool sar3d_bg_gate_enabled = true; // BG_GATE_ENABLE: peak-to-background ratio for the display + bool sar3d_clean_enabled = true; // SAR_CLEAN_ENABLE: CLEAN of a one-frame aperture into objects // 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/processing_live_config.cpp b/data_acq_and_processing/processing/data_processor/src/processing_live_config.cpp index c783972..398d82e 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 @@ -231,7 +231,9 @@ auto parse_geometry(const Json& value, const char* position_field, std::uint32_t config.sar3d_frame_radius = parse_u32_number(*found, "processing.sar3d_frame_radius"); } for (const auto& [key, target] : { - std::pair{"sar3d_voxel_m", &config.sar3d_voxel_m}, + std::pair{"sar3d_voxel_x_m", &config.sar3d_voxel_x_m}, + std::pair{"sar3d_voxel_y_m", &config.sar3d_voxel_y_m}, + std::pair{"sar3d_voxel_z_m", &config.sar3d_voxel_z_m}, std::pair{"sar3d_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}, @@ -244,6 +246,8 @@ auto parse_geometry(const Json& value, const char* position_field, std::uint32_t for (const auto& [key, target] : { std::pair{"sar3d_use_compensation", &config.sar3d_use_compensation}, std::pair{"sar3d_normalize_each_frame", &config.sar3d_normalize_each_frame}, + std::pair{"sar3d_bg_gate_enabled", &config.sar3d_bg_gate_enabled}, + std::pair{"sar3d_clean_enabled", &config.sar3d_clean_enabled}, }) { assign_optional(root, key, *target); } 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 48c3b6d..4b2f198 100644 --- a/data_acq_and_processing/processing/processors/include/sar3d_processor.hpp +++ b/data_acq_and_processing/processing/processors/include/sar3d_processor.hpp @@ -8,11 +8,13 @@ namespace radar::processing { -// 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 -// sar3d_peak = [level, X, Y, Z] of the volume maximum. +// 3D SAR of the reference notebook Multi_depthnorm_clean.ipynb on an NVIDIA GPU. +// Publishes max-projections of the normalized volume, or of the CLEAN display, 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); sar3d_peak = [level, X, Y, Z] of the shown peak; the frame's +// sar3d_background ratio and the CLEAN sar3d_objects and sar3d_marks tables (see +// python_app/gui/controllers/renderers/sar3d.py). class Sar3dProcessor final : public ProcessorInterface { public: Sar3dProcessor(); diff --git a/data_acq_and_processing/processing/processors/src/backprojection.hpp b/data_acq_and_processing/processing/processors/src/backprojection.hpp index 77882fe..7cd6ded 100644 --- a/data_acq_and_processing/processing/processors/src/backprojection.hpp +++ b/data_acq_and_processing/processing/processors/src/backprojection.hpp @@ -3,8 +3,10 @@ #include "shared_types.hpp" #include +#include #include #include +#include #include #include #include @@ -12,6 +14,7 @@ #include #include #include +#include #include namespace radar::processing::bp { @@ -19,6 +22,16 @@ namespace radar::processing::bp { constexpr double kWaveSpeed = 3e8; // v in both reference algorithms using Background = std::unordered_map>, ipc::ComboKeyHash>; +// numpy.median: the middle value, or the mean of the middle two for an even count. +template +[[nodiscard]] auto median(std::vector values) -> double { + if (values.empty()) { throw std::logic_error("median of no values"); } + const auto middle = values.begin() + static_cast(values.size() / 2U); + std::nth_element(values.begin(), middle, values.end()); + if (values.size() % 2U == 1U) { return static_cast(*middle); } + return 0.5 * (static_cast(*middle) + static_cast(*std::max_element(values.begin(), middle))); +} + // 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 { @@ -159,13 +172,20 @@ inline auto equalize_depth_limit(double max_depth, double x_min, double x_max, d return std::sqrt(max_depth * max_depth + x_far * x_far + 2.5 * 2.5) + 1.5 + 0.75; } +// Radar axes (Y, Z) of a scene point at height y above the ground and z along it, for a radar +// pitched by `pitch` with its origin `height` above the ground (sar_scene_to_radar_yz). +inline auto radar_yz(double y, double z, double height, double pitch) -> std::array { + const double y0 = y - height; + return {y0 * std::cos(pitch) + z * std::sin(pitch), -y0 * std::sin(pitch) + z * std::cos(pitch)}; +} + // 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. +// equalize_pairs_by_depth, Multi_depthnorm_clean.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 @@ -220,10 +240,7 @@ void normalize_pair(std::span> samples, double dz, double min 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)); } + double scale = median(std::move(magnitudes)); if (!std::isfinite(scale) || scale <= 1e-15) { scale = 1.0; } for (auto& sample : samples) { sample /= static_cast(scale + 1e-15); } } diff --git a/data_acq_and_processing/processing/processors/src/sar3d_clean.cpp b/data_acq_and_processing/processing/processors/src/sar3d_clean.cpp new file mode 100644 index 0000000..e689417 --- /dev/null +++ b/data_acq_and_processing/processing/processors/src/sar3d_clean.cpp @@ -0,0 +1,558 @@ +#include "sar3d_clean.hpp" +#include "backprojection.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +// Port of the CLEAN and frame parts of notebook cell 5. Comments name the notebook's +// functions and constants; the constants are values the notebook checked on scenes with +// several objects, not settings. +namespace radar::processing::sar3d { +namespace { + +// CLEAN (sar_clean_components, sar_clean_frame). +constexpr double kCleanStop = 5.5; // SAR_CLEAN_STOP_BG, backgrounds of the whole volume +constexpr double kCleanRangeStop = 4.5; // SAR_CLEAN_STOP_RANGE_BG, backgrounds at the voxel range +constexpr std::size_t kCleanMaxSteps = 15U; // SAR_CLEAN_MAX_STEPS +constexpr std::array kCleanMerge{0.35, 0.50, 0.40}; // _SAR_CLEAN_MERGE_M (X, Y, Z) +constexpr std::array kCleanBox{0.40, 0.50, 0.45}; // _SAR_CLEAN_BOX_M (X, Y, Z), half-size +constexpr double kCleanPick = 0.05; // _SAR_CLEAN_PICK_M, pair peak window around the delay +constexpr double kCleanSearchSigma = 1.0; // _SAR_CLEAN_SEARCH_SIGMA, voxels +constexpr double kCleanRing = 0.25; // _SAR_CLEAN_RING_M, range ring width +constexpr std::size_t kCleanRingSmooth = 5U; // _SAR_CLEAN_RING_SMOOTH, rings +constexpr std::size_t kCleanRingMinVoxels = 300U; // _SAR_CLEAN_RING_MIN_VOXELS + +// Boxes (sar_mark_groups). +constexpr double kMarkLevel = 0.6; // _SAR_MARK_LEVEL, response height level of the peak +constexpr double kMarkMemberRel = 0.4; // _SAR_MARK_MEMBER_REL +constexpr double kMarkPad = 0.25; // _SAR_MARK_PAD_M, along X and Z +constexpr double kMarkEdgeFar = 0.25; // _SAR_MARK_EDGE_FAR_M +constexpr double kMarkEdgeNear = 0.10; // _SAR_MARK_EDGE_NEAR_M +constexpr double kMarkEdgeX = 0.20; // _SAR_MARK_EDGE_X_M +constexpr double kMarkEdgeTop = 0.40; // _SAR_MARK_EDGE_TOP_M +constexpr std::array kMarkGroup{0.8, 0.5}; // SAR_MARK_GROUP_M (X, Z) +constexpr double kMarkGround = 0.35; // SAR_MARK_GROUND_M +constexpr double kMarkLow = 0.15; // SAR_MARK_LOW_M +constexpr double kMarkLowTop = 0.45; // SAR_MARK_LOW_TOP_M + +using Scans = std::vector>; + +// numpy.interp: linear between the points of increasing `xp`, constant beyond them. +[[nodiscard]] auto interp(double x, std::span xp, std::span fp) -> double { + if (x <= xp.front()) { + return fp.front(); + } + if (x >= xp.back()) { + return fp.back(); + } + const auto j = static_cast(std::ranges::upper_bound(xp, x) - xp.begin()) - 1U; + if (xp[j] == x) { + return fp[j]; + } + return (fp[j + 1U] - fp[j]) / (xp[j + 1U] - xp[j]) * (x - xp[j]) + fp[j]; +} + +// Radar axes (Y, Z) of a scene point of the frame. +[[nodiscard]] auto radar_yz(const VolumeRequest& frame, double y, double z) -> std::array { + return bp::radar_yz(y, z, frame.ground_depth_m, frame.pitch_rad); +} + +// Axis indices [first, first + count). +struct Range { + std::size_t first = 0U; + std::size_t count = 0U; +}; + +// The axis points within `half_width` of `center`; they are contiguous on a monotonic axis. +[[nodiscard]] auto within(std::span axis, double center, double half_width) -> Range { + Range range; + for (std::size_t index = 0U; index < axis.size(); ++index) { + if (std::abs(axis[index] - center) <= half_width) { + if (range.count == 0U) { + range.first = index; + } + range.count = index - range.first + 1U; + } + } + return range; +} + +// scan += amplitude * numpy.roll(pulse, shift) +void add_pulse(std::span> scan, std::span> pulse, std::size_t shift, + std::complex amplitude) { + const std::size_t length = scan.size(); + for (std::size_t sample = 0U; sample < length; ++sample) { + const std::complex response = pulse[(sample + length - shift) % length]; + scan[sample] = std::complex(std::complex(scan[sample]) + amplitude * response); + } +} + +// Range rings of the search grid (sar_clean_range_rings): every voxel's range from the +// radar, the ring centres and the working voxels (depth gate, reference planes) of each ring. +struct Rings { + std::vector range; + std::vector centers; + std::vector> members; +}; + +[[nodiscard]] auto range_rings(const VolumeRequest& frame, const Axes& grid) -> Rings { + Rings rings; + std::vector work; + for (std::size_t iy = 0U; iy < grid.y.size(); ++iy) { + const bool reference = reference_plane(grid.y, frame.norm_edge_margin_m, iy); + for (const double z : grid.z) { + const auto [y_radar, z_radar] = radar_yz(frame, grid.y[iy], z); + const bool working = reference && z_radar >= frame.min_depth_m && z_radar <= frame.max_depth_m; + for (const double x : grid.x) { + if (working) { + work.push_back(rings.range.size()); + } + rings.range.push_back(std::sqrt(x * x + y_radar * y_radar + z_radar * z_radar)); + } + } + } + // Ring edges numpy.arange(0, R.max() + 2 ring, ring); the voxel ring is numpy.digitize - 1, clipped. + const auto edges = static_cast( + std::ceil((*std::ranges::max_element(rings.range) + 2.0 * kCleanRing) / kCleanRing)); + rings.members.resize(edges - 1U); + for (std::size_t ring = 0U; ring + 1U < edges; ++ring) { + rings.centers.push_back((static_cast(ring) + 0.5) * kCleanRing); + } + for (const auto voxel : work) { + const auto ring = std::min(static_cast(rings.range[voxel] / kCleanRing), edges - 2U); + rings.members[ring].push_back(voxel); + } + return rings; +} + +// Background at the centre of every ring (sar_clean_range_background): the median of the +// ring, filled in over sparse rings and averaged over neighbouring ones; empty when no ring +// has enough voxels. +[[nodiscard]] auto range_background(std::span volume, const Rings& rings) -> std::vector { + std::vector known_rings; + std::vector known_medians; + for (std::size_t ring = 0U; ring < rings.members.size(); ++ring) { + const auto& members = rings.members[ring]; + if (members.size() >= kCleanRingMinVoxels) { + std::vector values; + values.reserve(members.size()); + for (const auto voxel : members) { + values.push_back(volume[voxel]); + } + known_rings.push_back(static_cast(ring)); + known_medians.push_back(bp::median(std::move(values))); + } + } + if (known_rings.empty()) { + return {}; + } + const std::size_t count = rings.members.size(); + std::vector medians(count); + for (std::size_t ring = 0U; ring < count; ++ring) { + medians[ring] = interp(static_cast(ring), known_rings, known_medians); + } + // Moving average over kCleanRingSmooth rings with edge padding (numpy.convolve 'valid' + // with numpy.ones(k) / k). + const auto half = static_cast(kCleanRingSmooth / 2U); + const auto last = static_cast(count) - 1; + const double weight = 1.0 / static_cast(kCleanRingSmooth); + std::vector smoothed(count); + for (std::size_t ring = 0U; ring < count; ++ring) { + for (std::ptrdiff_t tap = -half; tap <= half; ++tap) { + const auto index = std::clamp(static_cast(ring) + tap, std::ptrdiff_t{0}, last); + smoothed[ring] += medians[static_cast(index)] * weight; + } + } + return smoothed; +} + +// The search voxel highest above its stop threshold, the lower of kCleanStop backgrounds +// of the whole volume and kCleanRangeStop backgrounds at its range; its value in units of +// that threshold. Empty voxels score 0. +[[nodiscard]] auto strongest(std::span volume, const Rings& rings) -> std::pair { + std::vector positive; + std::ranges::copy_if(volume, std::back_inserter(positive), [](float value) { return value > 0.0F; }); + const double overall = positive.empty() ? 0.0 : kCleanStop * bp::median(std::move(positive)); + const auto local = range_background(volume, rings); + std::pair best{0U, 0.0}; + for (std::size_t voxel = 0U; voxel < volume.size(); ++voxel) { + if (!(volume[voxel] > 0.0F)) { + continue; + } + double threshold = overall; + if (!local.empty()) { + threshold = std::min(overall, kCleanRangeStop * interp(rings.range[voxel], rings.centers, local)); + } + const double score = volume[voxel] / (threshold + 1e-30); + if (score > best.second) { + best = {voxel, score}; + } + } + return best; +} + +// One subtracted point response: the pair, its peak sample and the pulse amplitude there. +struct Part { + std::size_t pair = 0U; + std::size_t sample = 0U; + std::complex amplitude{}; +}; + +struct Component { + std::array point{}; // search voxel, scene X, Y, Z + std::vector parts; +}; + +// CLEAN steps on the search grid (sar_clean_components). Each step takes the search voxel +// furthest above its stop threshold and subtracts a point response from the A-scan of every +// pair at the pair's own peak near the voxel delay; `residual` keeps what remains. +template +[[nodiscard]] auto subtract_components(const Image& image, const CleanRequest& request, Scans& residual) + -> std::vector { + const VolumeRequest& frame = request.frame; + const Axes& search = request.search; + const std::size_t length = request.pulse.size(); + const std::size_t plane = search.x.size() * search.z.size(); + const auto rings = range_rings(frame, search); + // Half-width of the pair peak window in samples: _SAR_CLEAN_PICK_M of range, which is + // twice that of Tx+Rx path. + const auto pick = std::max( + 1, static_cast(std::nearbyint(2.0 * kCleanPick * frame.samples_per_metre))); + std::vector components; + while (true) { + auto volume = image(residual, search, kCleanSearchSigma); + for (std::size_t iy = 0U; iy < search.y.size(); ++iy) { + if (!reference_plane(search.y, frame.norm_edge_margin_m, iy)) { + std::fill_n(volume.begin() + static_cast(iy * plane), plane, 0.0F); + } + } + const auto [voxel, level] = strongest(volume, rings); + if (level < 1.0 || components.size() >= kCleanMaxSteps) { + return components; + } + const std::size_t ix = voxel % search.x.size(); + const std::size_t iz = voxel / search.x.size() % search.z.size(); + const std::size_t iy = voxel / plane; + Component component{.point = {search.x[ix], search.y[iy], search.z[iz]}, .parts = {}}; + const auto [y_radar, z_radar] = radar_yz(frame, search.y[iy], search.z[iz]); + const std::array target{search.x[ix], y_radar, z_radar}; // in radar axes + const auto range = [&target](const std::array& antenna) { + const double dx = target[0] - antenna[0]; + const double dy = target[1] - antenna[1]; + const double dz = target[2] - antenna[2]; + return std::sqrt(dx * dx + dy * dy + dz * dz); + }; + for (std::size_t pair = 0U; pair < frame.contributions.size(); ++pair) { + const auto& term = frame.contributions[pair]; + // Sample of the voxel's Tx+Rx path (sar_clean_pair_delays). + const auto delay = static_cast( + std::nearbyint((range(term.tx) + range(term.rx)) * frame.samples_per_metre)); + const auto low = std::max(delay - pick, 0); + const auto high = std::min(delay + pick + 1, static_cast(length)); + if (high <= low) { + continue; + } + const auto scan = std::span(residual).subspan(pair * length, length); + auto peak = static_cast(low); + for (auto sample = peak + 1U; sample < static_cast(high); ++sample) { + if (std::norm(std::complex(scan[sample])) > std::norm(std::complex(scan[peak]))) { + peak = sample; + } + } + const auto amplitude = std::complex(scan[peak]) / std::complex(request.pulse[0]); + add_pulse(scan, request.pulse, peak, -amplitude); + component.parts.push_back({.pair = pair, .sample = peak, .amplitude = amplitude}); + } + components.push_back(std::move(component)); + } +} + +// Subtractions within _SAR_CLEAN_MERGE_M of an object's first subtraction belong to that +// object (sar_clean_group). +struct Group { + std::array seed{}; + std::vector members; +}; + +[[nodiscard]] auto group(const std::vector& components) -> std::vector { + std::vector groups; + for (const auto& component : components) { + const auto near_seed = [&component](const Group& candidate) { + for (std::size_t axis = 0U; axis < 3U; ++axis) { + if (std::abs(component.point[axis] - candidate.seed[axis]) > kCleanMerge[axis]) { + return false; + } + } + return true; + }; + if (const auto found = std::ranges::find_if(groups, near_seed); found != groups.end()) { + found->members.push_back(&component); + } else { + groups.push_back({.seed = component.point, .members = {&component}}); + } + } + return groups; +} + +} // namespace + +auto frame_background(std::span values, std::span volume, const VolumeRequest& frame) + -> FrameBackground { + const Axes& grid = frame.grid; + const std::size_t plane = grid.x.size() * grid.z.size(); + if (values.size() != volume.size() || volume.size() != plane * grid.y.size()) { + throw std::runtime_error("sar3d background volumes do not match the frame grid"); + } + std::vector work; + work.reserve(volume.size()); + FrameBackground result; + for (std::size_t iy = 0U; iy < grid.y.size(); ++iy) { + if (!reference_plane(grid.y, frame.norm_edge_margin_m, iy)) { + continue; + } + for (std::size_t voxel = iy * plane; voxel < (iy + 1U) * plane; ++voxel) { + if (volume[voxel] > 0.0F) { + work.push_back(values[voxel]); + result.peak = std::max(result.peak, static_cast(values[voxel])); + } + } + } + if (!work.empty()) { + result.median = bp::median(std::move(work)); + } + return result; +} + +auto clean(Gpu& gpu, const CleanRequest& request) -> Clean { + const VolumeRequest& frame = request.frame; + const Axes& grid = frame.grid; + const std::size_t pairs = frame.contributions.size(); + const std::size_t length = request.pulse.size(); + if (pairs == 0U || request.ascans.size() != pairs * length || + request.volume.size() != grid.x.size() * grid.y.size() * grid.z.size() || !(request.level > 0.0)) { + throw std::runtime_error("sar3d CLEAN input does not match its frame"); + } + // Smoothed |BP| of the frame pairs with other A-scans on another grid (sar_abs_volume). + // The A-scans are full-length: the search grid may reach past the frame grid that the + // frame's shortened A-scans cover. + const auto image = [&gpu, &frame, length](const Scans& scans, const Axes& axes, double sigma) { + VolumeRequest imaged = frame; + imaged.grid = axes; + imaged.ascans = scans; + imaged.samples_per_ascan = length; + imaged.smooth_sigma_voxels = sigma; + gpu.image(imaged); + return gpu.download(); + }; + + Scans residual(request.ascans.begin(), request.ascans.end()); + const auto components = subtract_components(image, request, residual); + Clean result; + result.steps = components.size(); + if (components.empty()) { + result.residual_ratio = request.frame_background.ratio(); + return result; + } + // The residual in units of the frame maximum, hidden when it is background. + const auto residual_image = image(residual, grid, frame.smooth_sigma_voxels); + result.residual_ratio = frame_background(residual_image, request.volume, frame).ratio(); + result.display.assign(residual_image.size(), 0.0F); + if (!request.bg_gate || result.residual_ratio >= kBgGateShow) { + std::ranges::transform(residual_image, result.display.begin(), + [&request](float value) { return static_cast(value / request.level); }); + } + + const std::size_t grid_nx = grid.x.size(); + const std::size_t grid_nz = grid.z.size(); + for (const auto& object : group(components)) { + const auto& seed = object.seed; + // The object apart from the others: the residual plus its own subtractions, imaged on + // the frame grid around it only (sar_clean_subgrid). + const Range window_x = within(grid.x, seed[0], 2.0 * kCleanBox[0]); + const Range window_y = within(grid.y, seed[1], 2.0 * kCleanBox[1]); + const Range window_z = within(grid.z, seed[2], 2.0 * kCleanBox[2]); + if (window_x.count < 2U || window_y.count < 2U || window_z.count < 2U) { + continue; + } + const Axes local{ + .x = grid.x.subspan(window_x.first, window_x.count), + .y = grid.y.subspan(window_y.first, window_y.count), + .z = grid.z.subspan(window_z.first, window_z.count), + }; + Scans scans = residual; + for (const Component* member : object.members) { + for (const Part& part : member->parts) { + add_pulse(std::span(scans).subspan(part.pair * length, length), request.pulse, part.sample, + part.amplitude); + } + } + const auto values = image(scans, local, frame.smooth_sigma_voxels); + const std::size_t nx = local.x.size(); + const std::size_t nz = local.z.size(); + const auto at = [&values, nx, nz](std::size_t iy, std::size_t iz, std::size_t ix) { + return values[(iy * nz + iz) * nx + ix]; + }; + + // Peak within _SAR_CLEAN_BOX_M of the first subtraction. + const Range near_x = within(local.x, seed[0], kCleanBox[0]); + const Range near_y = within(local.y, seed[1], kCleanBox[1]); + const Range near_z = within(local.z, seed[2], kCleanBox[2]); + if (near_x.count == 0U || near_y.count == 0U || near_z.count == 0U) { + continue; + } + std::size_t py = near_y.first; + std::size_t pz = near_z.first; + std::size_t px = near_x.first; + for (std::size_t iy = near_y.first; iy < near_y.first + near_y.count; ++iy) { + for (std::size_t iz = near_z.first; iz < near_z.first + near_z.count; ++iz) { + for (std::size_t ix = near_x.first; ix < near_x.first + near_x.count; ++ix) { + if (at(iy, iz, ix) > at(py, pz, px)) { + py = iy; + pz = iz; + px = ix; + } + } + } + } + const double peak = at(py, pz, px); + if (!(peak > 0.0)) { + continue; + } + + // Height of the response: where the maximum over 5 x 5 voxels in Z and X around the + // peak stays at _SAR_MARK_LEVEL of the peak. + const auto column = [&](std::size_t iy) { + float maximum = 0.0F; + for (std::size_t iz = pz < 2U ? 0U : pz - 2U; iz < std::min(pz + 3U, nz); ++iz) { + for (std::size_t ix = px < 2U ? 0U : px - 2U; ix < std::min(px + 3U, nx); ++ix) { + maximum = std::max(maximum, at(iy, iz, ix)); + } + } + return maximum / peak; + }; + std::size_t y_low = py; + std::size_t y_high = py; + while (y_low > 0U && column(y_low - 1U) >= kMarkLevel) { + --y_low; + } + while (y_high + 1U < local.y.size() && column(y_high + 1U) >= kMarkLevel) { + ++y_high; + } + + // The object in its own 0-1 scale over the box around its peak. + const Range box_x = within(local.x, local.x[px], kCleanBox[0]); + const Range box_y = within(local.y, local.y[py], kCleanBox[1]); + const Range box_z = within(local.z, local.z[pz], kCleanBox[2]); + for (std::size_t iy = box_y.first; iy < box_y.first + box_y.count; ++iy) { + for (std::size_t iz = box_z.first; iz < box_z.first + box_z.count; ++iz) { + for (std::size_t ix = box_x.first; ix < box_x.first + box_x.count; ++ix) { + float& shown = result.display[((window_y.first + iy) * grid_nz + window_z.first + iz) * grid_nx + + window_x.first + ix]; + shown = std::max(shown, static_cast(at(iy, iz, ix) / peak)); + } + } + } + result.objects.push_back({ + .index = {window_y.first + py, window_z.first + pz, window_x.first + px}, + .position = {local.x[px], local.y[py], local.z[pz]}, + .y_low = local.y[y_low], + .y_high = local.y[y_high], + .level = peak / (request.frame_background.median + 1e-30), + }); + } + std::ranges::stable_sort(result.objects, std::ranges::greater{}, &CleanObject::level); + if (result.objects.empty()) { + result.display.clear(); + } + return result; +} + +auto marks(std::span objects, const VolumeRequest& frame) -> std::vector { + // Objects within SAR_MARK_GROUP_M of each other along X and Z, at any height and through + // chains of them, are one item. + std::vector parent(objects.size()); + std::iota(parent.begin(), parent.end(), 0U); + const auto root = [&parent](std::size_t index) { + while (parent[index] != index) { + index = parent[index]; + } + return index; + }; + for (std::size_t i = 0U; i < objects.size(); ++i) { + for (std::size_t j = i + 1U; j < objects.size(); ++j) { + if (std::abs(objects[i].position[0] - objects[j].position[0]) <= kMarkGroup[0] && + std::abs(objects[i].position[2] - objects[j].position[2]) <= kMarkGroup[1]) { + parent[root(i)] = root(j); + } + } + } + // Members in object order, so the strongest leads each group and the groups come in the + // order of their strongest members. + std::vector roots; + std::vector> groups; + for (std::size_t index = 0U; index < objects.size(); ++index) { + const auto found = std::ranges::find(roots, root(index)); + if (found == roots.end()) { + roots.push_back(root(index)); + groups.push_back({index}); + } else { + groups[static_cast(found - roots.begin())].push_back(index); + } + } + + const auto& x = frame.grid.x; + const auto& y = frame.grid.y; + std::vector result; + for (const auto& members : groups) { + const CleanObject& top = objects[members.front()]; + // Height of the centre: the mean height of the strong members weighted by level. + double weight = 0.0; + double weighted_height = 0.0; + double top_y = -std::numeric_limits::infinity(); + Mark mark; + mark.box = {std::numeric_limits::infinity(), -std::numeric_limits::infinity(), + std::numeric_limits::infinity(), -std::numeric_limits::infinity(), + std::numeric_limits::infinity(), -std::numeric_limits::infinity()}; + for (const auto index : members) { + const CleanObject& member = objects[index]; + if (member.level >= kMarkMemberRel * top.level) { + weight += member.level; + weighted_height += member.level * member.position[1]; + top_y = std::max(top_y, member.y_high); + } + auto& box = mark.box; + box[0] = std::min(box[0], member.position[0] - kMarkPad); + box[1] = std::max(box[1], member.position[0] + kMarkPad); + box[2] = std::min(box[2], member.y_low); + box[3] = std::max(box[3], member.y_high); + box[4] = std::min(box[4], member.position[2] - kMarkPad); + box[5] = std::max(box[5], member.position[2] + kMarkPad); + } + const double centre_y = weighted_height / weight; + // The strongest point near the scene boundary may belong to a reflector partly + // outside it, so its position and height are not yet reliable. + const double depth = radar_yz(frame, top.position[1], top.position[2])[1]; + const bool at_edge = depth >= frame.max_depth_m - kMarkEdgeFar || depth <= frame.min_depth_m + kMarkEdgeNear || + top.position[0] <= x.front() + kMarkEdgeX || top.position[0] >= x.back() - kMarkEdgeX || + top.position[1] >= y.back() - kMarkEdgeTop; + if (at_edge) { + mark.kind = MarkKind::edge; + } else if (centre_y <= kMarkLow && top_y <= kMarkLowTop) { + mark.kind = MarkKind::low; + } else if (centre_y <= kMarkGround) { + mark.kind = MarkKind::ground; + } else { + mark.kind = MarkKind::up; + } + result.push_back(mark); + } + return result; +} + +} // namespace radar::processing::sar3d diff --git a/data_acq_and_processing/processing/processors/src/sar3d_clean.hpp b/data_acq_and_processing/processing/processors/src/sar3d_clean.hpp new file mode 100644 index 0000000..5ff9ce9 --- /dev/null +++ b/data_acq_and_processing/processing/processors/src/sar3d_clean.hpp @@ -0,0 +1,91 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "sar3d_gpu.hpp" + +// Display stage of one 3D SAR frame (notebook Multi_depthnorm_clean.ipynb, cell 5): the +// frame's peak over its background, CLEAN of a one-frame image into separate objects, and +// the boxes drawn around them. +namespace radar::processing::sar3d { + +// BG_GATE_3D_SHOW and BG_GATE_3D_FULL: peak-to-background ratios at which a frame appears +// and from which it is shown at full scale. +inline constexpr double kBgGateShow = 6.0; +inline constexpr double kBgGateFull = 8.5; + +// SAR_CLEAN_GRID_M along X and Z, SAR_CLEAN_GRID_Y_M along Y: steps of the CLEAN search grid [m]. +inline constexpr std::array kCleanSearchStep{0.04, 0.10, 0.04}; + +// Background of a frame: the median over its working voxels and the maximum there +// (frame_peak_over_background). +struct FrameBackground { + double median = 0.0; + double peak = 0.0; + + [[nodiscard]] auto ratio() const -> double { return median > 0.0 ? peak / median : 0.0; } +}; + +// Background of `values` over the working voxels of `volume`, an image of `frame`: those it +// reaches from inside the depth gate (non-zero) on the reference planes (sar_work_mask). +[[nodiscard]] auto frame_background(std::span values, std::span volume, + const VolumeRequest& frame) -> FrameBackground; + +struct CleanRequest { + // Request of the frame image, a one-frame aperture; CLEAN images other A-scans with it. + VolumeRequest frame; + std::span volume; // the frame image, unscaled + double level = 0.0; // its normalization (sar_finish_volume) + FrameBackground frame_background; // of that image + // Search grid over the full scene region around the frame grid (sar_clean_search_grid). + Axes search; + std::span> ascans; // full A-scan of each frame pair, back to back + std::span> pulse; // zero-delay point response, one A-scan long + bool bg_gate = true; // BG_GATE_ENABLE +}; + +struct CleanObject { + std::array index{}; // peak voxel on the frame grid, [y, z, x] + std::array position{}; // peak, scene X, Y, Z [m] + double y_low = 0.0; // heights where the response stays at + double y_high = 0.0; // _SAR_MARK_LEVEL of the peak [m] + double level = 0.0; // peak over the frame background (level_bg) +}; + +struct Clean { + std::vector objects; // strongest first + // The residual in units of the frame maximum with every object over it in its own + // 0-1 scale, on the frame grid; empty without objects. + std::vector display; + std::size_t steps = 0U; // subtractions + double residual_ratio = 0.0; // residual peak over its background +}; + +// CLEAN of a one-frame image (sar_clean_frame): subtracts the strongest point responses +// from the A-scans until the residual is background, then images every object apart. +// Imaging leaves one of these images, not the frame, as the current volume of `gpu`. +[[nodiscard]] auto clean(Gpu& gpu, const CleanRequest& request) -> Clean; + +// Classes of the boxes; the values are their codes in the results. +enum class MarkKind : std::uint8_t { + low = 0U, // centre and top of the response near the ground (like a mine) + ground = 1U, // on the ground + up = 2U, // above the ground + edge = 3U, // strongest point near the scene boundary: not yet decidable +}; + +struct Mark { + MarkKind kind = MarkKind::edge; + std::array box{}; // X0, X1, Y0, Y1, Z0, Z1 [m] +}; + +// Boxes around groups of the CLEAN objects of an image of `frame` (sar_mark_groups); +// `objects` come strongest first, and so do the boxes. +[[nodiscard]] auto marks(std::span objects, const VolumeRequest& frame) -> std::vector; + +} // namespace radar::processing::sar3d 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 dec6acb..7468557 100644 --- a/data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp +++ b/data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp @@ -27,51 +27,31 @@ auto checked_product(std::size_t a, std::size_t b, std::size_t limit, const char return a * b; } -#ifdef SAR3D_WITH_CUDA -auto gaussian_taps(double sigma) -> std::vector { - if (!(sigma > 0.0) || !std::isfinite(sigma) || sigma > 100.0) { - throw std::runtime_error("sar3d smoothing sigma must be between 0 and 100"); +// Two points per axis, 32-bit voxel indices and coordinates representable in float. +void check_grid(const Axes& grid) { + const auto plane = checked_product(grid.x.size(), grid.z.size(), kMaxIndexedItems, "grid plane"); + (void)checked_product(plane, grid.y.size(), kMaxIndexedItems, "grid (32-bit indexing; increase SAR_VOXEL_M)"); + if (grid.x.size() < 2U || grid.y.size() < 2U || grid.z.size() < 2U) { + throw std::runtime_error("sar3d grid needs two points on every axis"); } - const int radius = static_cast(4.0 * sigma + 0.5); - std::vector weights; - double total = 0.0; - for (int k = -radius; k <= radius; ++k) { - weights.push_back(std::exp(-0.5 * k * k / (sigma * sigma))); - total += weights.back(); - } - std::vector taps; - for (const double weight : weights) { - taps.push_back(static_cast(weight / total)); - } - return taps; -} - -auto reference_plane(const VolumeRequest& request, std::size_t iy) -> bool { - const double margin = request.norm_edge_margin_m; - if (request.y.size() == 2U || request.y.back() - request.y.front() <= 2.5 * margin) { - return true; - } - return request.y[iy] >= request.y.front() + margin && request.y[iy] <= request.y.back() - margin; -} - -void choose_peak(Projections& result, const VolumeRequest& request, - std::span levels, std::span indices = {}) { - for (std::size_t iy = 0U; iy < request.y.size(); ++iy) { - if (!std::isfinite(levels[iy])) { - throw std::runtime_error("sar3d volume contains a non-finite value"); - } - if (reference_plane(request, iy) && levels[iy] > result.level) { - result.level = levels[iy]; - if (!indices.empty()) { - result.peak = {iy, indices[iy] / request.x.size(), indices[iy] % request.x.size()}; + for (const auto axis : {grid.x, grid.y, grid.z}) { + for (const double coordinate : axis) { + if (!std::isfinite(coordinate) || std::abs(coordinate) > std::numeric_limits::max()) { + throw std::runtime_error("sar3d grid contains an unsupported coordinate"); } } } } -#endif } // namespace +auto reference_plane(std::span y, double margin, std::size_t iy) -> bool { + if (y.size() == 2U || y.back() - y.front() <= 2.5 * margin) { + return true; + } + return y[iy] >= y.front() + margin && y[iy] <= y.back() - margin; +} + #ifdef SAR3D_WITH_CUDA namespace { @@ -134,11 +114,12 @@ extern "C" __global__ void backproject( if (index >= (unsigned)voxels) return; const int voxel = (int)index; const int row = voxel / nx; - // Scene -> radar axes (sar_scene_to_radar_yz). The depth gate is decided in - // double so voxels on its edge are classified like in the reference. + // Scene -> radar axes (sar_scene_to_radar_yz). The depth gate is decided in double + // and rounded like NumPy, without a fused multiply-add: the grid corners lie on its + // edge by construction and must be classified like in the reference. const double y0 = gy[row / nz] - ground_depth; const double zs = gz[row % nz]; - const double depth = -y0 * sin_p + zs * cos_p; + const double depth = __dadd_rn(__dmul_rn(-y0, sin_p), __dmul_rn(zs, cos_p)); if (depth < min_depth || depth > max_depth) { volume[voxel] = 0.0f; return; @@ -208,25 +189,48 @@ extern "C" __global__ void max_project( dst[o] = maximum; } -// Keep the first joint (y,z,x) maximum. Separate projections cannot locate it. +// Maximum of each Y plane and its first in-plane index, one block per plane, so the +// host finds the first joint (y,z,x) maximum; separate projections cannot locate it. +// A non-finite voxel makes its plane level NaN. Shared memory holds blockDim.x +// (a power of two) levels and as many indices. extern "C" __global__ void peak_planes( - const float* volume, const int plane, const int ny, float* levels, unsigned int* indices) + const float* __restrict__ volume, const int plane, float* __restrict__ levels, + unsigned int* __restrict__ indices) { - const unsigned thread = blockIdx.x * blockDim.x + threadIdx.x; - if (thread >= (unsigned)ny) return; - const int y = (int)thread; - const float* row = volume + y * plane; + extern __shared__ float best[]; + unsigned int* where = (unsigned int*)(best + blockDim.x); + __shared__ int non_finite; + if (threadIdx.x == 0) non_finite = 0; + __syncthreads(); + const float* row = volume + blockIdx.x * plane; float level = 0.0f; unsigned int index = 0u; - for (int i = 0; i < plane; ++i) { - if (!isfinite(row[i])) { levels[y] = row[i]; indices[y] = (unsigned int)i; return; } - if (row[i] > level) { level = row[i]; index = (unsigned int)i; } + for (int i = threadIdx.x; i < plane; i += blockDim.x) { + const float value = row[i]; + if (!isfinite(value)) non_finite = 1; + if (value > level) { level = value; index = (unsigned int)i; } + } + best[threadIdx.x] = level; + where[threadIdx.x] = index; + __syncthreads(); + for (unsigned int stride = blockDim.x / 2; stride > 0; stride /= 2) { + if (threadIdx.x < stride) { + const float other = best[threadIdx.x + stride]; + const unsigned int other_index = where[threadIdx.x + stride]; + if (other > best[threadIdx.x] || (other == best[threadIdx.x] && other_index < where[threadIdx.x])) { + best[threadIdx.x] = other; + where[threadIdx.x] = other_index; + } + } + __syncthreads(); + } + if (threadIdx.x == 0) { + levels[blockIdx.x] = non_finite ? __int_as_float(0x7fc00000) : best[0]; + indices[blockIdx.x] = where[0]; } - levels[y] = level; - indices[y] = index; } -// Contiguous half-maximum interval through the joint peak on each scene axis. +// Contiguous half-maximum interval through the given peak voxel 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, @@ -281,7 +285,26 @@ extern "C" __global__ void project_oblique( } )CUDA"; -constexpr int kBlockSize = 256; +constexpr unsigned kBlockSize = 256U; + +// scipy.ndimage.gaussian_filter1d weights (truncate = 4). +auto gaussian_taps(double sigma) -> std::vector { + if (!(sigma > 0.0) || !std::isfinite(sigma) || sigma > 100.0) { + throw std::runtime_error("sar3d smoothing sigma must be between 0 and 100"); + } + const int radius = static_cast(4.0 * sigma + 0.5); + std::vector weights; + double total = 0.0; + for (int k = -radius; k <= radius; ++k) { + weights.push_back(std::exp(-0.5 * k * k / (sigma * sigma))); + total += weights.back(); + } + std::vector taps; + for (const double weight : weights) { + taps.push_back(static_cast(weight / total)); + } + return taps; +} void check(CUresult result, const char* what, bool device_required = false) { if (result != CUDA_SUCCESS) { @@ -398,18 +421,24 @@ class DeviceBuffer { }; template -[[nodiscard]] auto download(CUdeviceptr source, std::size_t count) -> std::vector { +[[nodiscard]] auto copy_to_host(CUdeviceptr source, std::size_t count) -> std::vector { std::vector values(count); check(cuMemcpyDtoH(values.data(), source, sizeof(T) * count), "download"); return values; } -// One thread per item; argument types must match the kernel signature exactly. +// `blocks` blocks of kBlockSize threads; argument types must match the kernel signature exactly. +template +void launch_blocks(CUfunction kernel, unsigned blocks, unsigned shared_bytes, Args... args) { + void* parameters[] = {static_cast(&args)...}; + check(cuLaunchKernel(kernel, blocks, 1U, 1U, kBlockSize, 1U, 1U, shared_bytes, nullptr, parameters, nullptr), + "launch"); +} + +// One thread per item. template void launch(CUfunction kernel, int items, Args... args) { - void* parameters[] = {static_cast(&args)...}; - const auto blocks = (static_cast(items) + kBlockSize - 1U) / kBlockSize; - check(cuLaunchKernel(kernel, blocks, 1U, 1U, kBlockSize, 1U, 1U, 0U, nullptr, parameters, nullptr), "launch"); + launch_blocks(kernel, (static_cast(items) + kBlockSize - 1U) / kBlockSize, 0U, args...); } } // namespace @@ -447,129 +476,193 @@ struct Gpu::Impl { 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; DeviceBuffer bands; - DeviceBuffer axes; + DeviceBuffer axes; // of the current volume: x, then y, then z DeviceBuffer contributions; DeviceBuffer ascans; - DeviceBuffer volume; + DeviceBuffer volume; // the current volume DeviceBuffer scratch; DeviceBuffer images; - DeviceBuffer peak_stats; + DeviceBuffer plane_stats; + int nx = 0; // current volume size; 0 before the first + int ny = 0; + int nz = 0; ~Impl() { cuCtxSetCurrent(context.handle); } - explicit Impl(std::span weights) { - smooth_radius = static_cast(weights.size() / 2U); - taps.upload(weights); + + [[nodiscard]] auto gx() const -> CUdeviceptr { return axes.pointer(); } + [[nodiscard]] auto gy() const -> CUdeviceptr { return gx() + sizeof(double) * static_cast(nx); } + [[nodiscard]] auto gz() const -> CUdeviceptr { return gy() + sizeof(double) * static_cast(ny); } + [[nodiscard]] auto voxels() const -> int { return nx * ny * nz; } + + void use_context() const { check(cuCtxSetCurrent(context.handle), "context"); } + + void require_volume() const { + if (voxels() == 0) { + throw std::runtime_error("sar3d has no current volume"); + } } - auto inverse_dft(std::span>, std::size_t, std::size_t, std::size_t) -> std::vector>; - auto project(const VolumeRequest&) -> Projections; -}; -auto Gpu::Impl::inverse_dft( - std::span> bands, - std::size_t points, - std::size_t first_bin, - std::size_t size -) -> std::vector> { - Impl& gpu = *this; - check(cuCtxSetCurrent(gpu.context.handle), "context"); - const std::size_t samples = bands.size() / points * size; - const int items = checked_int(samples, "A-scan stack"); - const CUdeviceptr input = gpu.bands.upload(bands); - const CUdeviceptr output = gpu.ascans.reserve(sizeof(std::complex) * samples); - launch(gpu.inverse_dft_kernel, items, input, checked_int(points, "band"), checked_int(first_bin, "first bin"), - checked_int(size, "A-scan"), items, output); - return download>(output, samples); -} + // There is no current volume until an update of it has succeeded. + void forget_volume() { nx = ny = nz = 0; } -auto Gpu::Impl::project(const VolumeRequest& request) -> Projections { - Impl& gpu = *this; - check(cuCtxSetCurrent(gpu.context.handle), "context"); + // Makes `grid` the grid of the current volume; returns the volume, sized for it. + auto use_grid(const Axes& grid) -> CUdeviceptr { + forget_volume(); + std::vector values(grid.x.begin(), grid.x.end()); + values.insert(values.end(), grid.y.begin(), grid.y.end()); + values.insert(values.end(), grid.z.begin(), grid.z.end()); + axes.upload(std::span(values)); + const CUdeviceptr output = volume.reserve(sizeof(float) * grid.x.size() * grid.y.size() * grid.z.size()); + nx = static_cast(grid.x.size()); + ny = static_cast(grid.y.size()); + nz = static_cast(grid.z.size()); + return output; + } - const int nx = checked_int(request.x.size(), "grid"); - const int ny = checked_int(request.y.size(), "grid"); - const int nz = checked_int(request.z.size(), "grid"); - const int voxels = checked_int(request.x.size() * request.y.size() * request.z.size(), "grid"); - const int samples = checked_int(request.samples_per_ascan, "A-scan"); + // Maximum of the current volume over the reference planes of `y`, and its first voxel. + auto plane_peak(std::span y, double margin) -> Peak { + const auto planes = static_cast(ny); + const CUdeviceptr levels = plane_stats.reserve(2U * sizeof(float) * planes); + const CUdeviceptr indices = levels + sizeof(float) * planes; + launch_blocks(peak_planes, static_cast(ny), static_cast(2U * kBlockSize * sizeof(float)), + volume.pointer(), nz * nx, levels, indices); + const auto plane_levels = copy_to_host(levels, planes); + const auto plane_indices = copy_to_host(indices, planes); + const auto row = static_cast(nx); + Peak peak; + for (std::size_t iy = 0U; iy < planes; ++iy) { + if (!std::isfinite(plane_levels[iy])) { + throw std::runtime_error("sar3d volume contains a non-finite value"); + } + if (reference_plane(y, margin, iy) && plane_levels[iy] > peak.level) { + peak.level = plane_levels[iy]; + peak.index = {iy, plane_indices[iy] / row, plane_indices[iy] % row}; + } + } + return peak; + } - std::vector axes(request.x.begin(), request.x.end()); - axes.insert(axes.end(), request.y.begin(), request.y.end()); - axes.insert(axes.end(), request.z.begin(), request.z.end()); - const CUdeviceptr gx = gpu.axes.upload(std::span(axes)); - const CUdeviceptr gy = gx + sizeof(double) * request.x.size(); - const CUdeviceptr gz = gy + sizeof(double) * request.y.size(); - const CUdeviceptr contributions = gpu.contributions.upload(request.contributions); - const CUdeviceptr ascans = gpu.ascans.upload(request.ascans); - const CUdeviceptr volume = gpu.volume.reserve(sizeof(float) * static_cast(voxels)); - const CUdeviceptr scratch = gpu.scratch.reserve(sizeof(float) * static_cast(voxels)); + auto inverse_dft( + std::span> spectra, + std::size_t points, + std::size_t first_bin, + std::size_t size + ) -> std::vector> { + use_context(); + const std::size_t samples = spectra.size() / points * size; + const int items = checked_int(samples, "A-scan stack"); + const CUdeviceptr input = bands.upload(spectra); + const CUdeviceptr output = ascans.reserve(sizeof(std::complex) * samples); + launch(inverse_dft_kernel, items, input, checked_int(points, "band"), checked_int(first_bin, "first bin"), + checked_int(size, "A-scan"), items, output); + return copy_to_host>(output, samples); + } - launch(gpu.backproject, voxels, ascans, samples, contributions, - checked_int(request.contributions.size(), "aperture"), gx, gy, gz, nx, nz, voxels, - std::cos(request.pitch_rad), std::sin(request.pitch_rad), request.ground_depth_m, - request.min_depth_m, request.max_depth_m, static_cast(request.samples_per_metre), - static_cast(request.use_compensation), request.comp_range_power, request.comp_weight_max, volume); - const int plane = nz * nx; - const auto stats = gpu.peak_stats.reserve(2U * sizeof(float) * static_cast(ny)); - launch(gpu.peak_planes, ny, volume, plane, ny, stats, stats + sizeof(float) * static_cast(ny)); - Projections raw_peak; - const auto raw_levels = download(stats, static_cast(ny)); - choose_peak(raw_peak, request, raw_levels); - // Separable Gaussian with sigma (0, s, s) over [y][z][x]: Z first, then X, like scipy. - const CUdeviceptr taps = gpu.taps.pointer(); - launch(gpu.smooth, voxels, volume, scratch, taps, gpu.smooth_radius, nx, nz, voxels, 1); - launch(gpu.smooth, voxels, scratch, volume, taps, gpu.smooth_radius, nx, nz, voxels, 0); + auto image(const VolumeRequest& request) -> Peak { + const auto weights = gaussian_taps(request.smooth_sigma_voxels); + use_context(); + try { + const CUdeviceptr output = use_grid(request.grid); + const CUdeviceptr work = scratch.reserve(sizeof(float) * static_cast(voxels())); + const CUdeviceptr terms = contributions.upload(request.contributions); + const CUdeviceptr scans = ascans.upload(request.ascans); + launch(backproject, voxels(), scans, checked_int(request.samples_per_ascan, "A-scan"), terms, + checked_int(request.contributions.size(), "aperture"), gx(), gy(), gz(), nx, nz, voxels(), + std::cos(request.pitch_rad), std::sin(request.pitch_rad), request.ground_depth_m, + request.min_depth_m, request.max_depth_m, static_cast(request.samples_per_metre), + static_cast(request.use_compensation), request.comp_range_power, request.comp_weight_max, + output); + const float raw_level = plane_peak(request.grid.y, request.norm_edge_margin_m).level; + // Separable Gaussian with sigma (0, s, s) over [y][z][x]: Z first, then X, like scipy. + const CUdeviceptr kernel = taps.upload(std::span(weights)); + const int radius = static_cast(weights.size() / 2U); + launch(smooth, voxels(), output, work, kernel, radius, nx, nz, voxels(), 1); + launch(smooth, voxels(), work, output, kernel, radius, nx, nz, voxels(), 0); + Peak peak = plane_peak(request.grid.y, request.norm_edge_margin_m); + peak.raw_level = raw_level; + return peak; + } catch (...) { + forget_volume(); + throw; + } + } - const ObliqueView& view = request.oblique; - const int front_size = ny * nx; - const int side_size = ny * nz; - const int top_size = nz * nx; - const int oblique_size = checked_int(static_cast(view.nv) * static_cast(view.nu), "view"); - const auto image_floats = static_cast(front_size) + side_size + top_size + oblique_size + 2U * ny; - const CUdeviceptr front = gpu.images.reserve(sizeof(float) * image_floats); - const CUdeviceptr side = front + sizeof(float) * static_cast(front_size); - const CUdeviceptr top = side + sizeof(float) * static_cast(side_size); - const CUdeviceptr oblique = top + sizeof(float) * static_cast(top_size); - const CUdeviceptr levels = oblique + sizeof(float) * static_cast(oblique_size); - const CUdeviceptr indices = levels + sizeof(float) * static_cast(ny); - launch(gpu.peak_planes, ny, volume, plane, ny, levels, indices); - Projections result; - result.raw_level = raw_peak.level; - 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); + auto download() -> std::vector { + use_context(); + require_volume(); + return copy_to_host(volume.pointer(), static_cast(voxels())); + } + + void load(const Axes& grid, std::span values) { + use_context(); + try { + check(cuMemcpyHtoD(use_grid(grid), values.data(), values.size_bytes()), "upload"); + } catch (...) { + forget_volume(); + throw; + } + } + + auto project(const std::array& peak, const ObliqueView& view) -> Projections { + use_context(); + require_volume(); + const auto [py, pz, px] = peak; + if (py >= static_cast(ny) || pz >= static_cast(nz) || + px >= static_cast(nx)) { + throw std::runtime_error("sar3d peak lies outside the volume"); + } + const int plane = nz * nx; + const int front_size = ny * nx; + const int side_size = ny * nz; + const int top_size = plane; + const int oblique_size = view.nv * view.nu; + const auto floats = static_cast(front_size) + side_size + top_size + oblique_size + 3U; + const CUdeviceptr front = images.reserve(sizeof(float) * floats); + const CUdeviceptr side = front + sizeof(float) * static_cast(front_size); + const CUdeviceptr top = side + sizeof(float) * static_cast(side_size); + const CUdeviceptr oblique = top + sizeof(float) * static_cast(top_size); + const CUdeviceptr widths = oblique + sizeof(float) * static_cast(oblique_size); + const CUdeviceptr values = volume.pointer(); + launch(max_project, front_size, values, front, front_size, nx, plane, 1, nz, nx); + launch(max_project, side_size, values, side, side_size, nz, plane, nx, nx, 1); + launch(max_project, top_size, values, top, top_size, nx, nx, 1, ny, plane); + check(cuMemsetD32(oblique, 0U, static_cast(oblique_size)), "clear"); + launch(project_oblique, voxels(), values, gx(), gy(), gz(), nx, nz, voxels(), view.u[0], view.u[1], view.u[2], + view.v[0], view.v[1], view.v[2], view.u0, view.v0, 1.0 / view.pixel_m, view.nu, view.nv, oblique); + launch(peak_widths, 3, values, gx(), gy(), gz(), nx, ny, nz, + static_cast(px), static_cast(py), static_cast(pz), widths); + Projections result; + result.front = copy_to_host(front, static_cast(front_size)); + result.side = copy_to_host(side, static_cast(side_size)); + result.top = copy_to_host(top, static_cast(top_size)); + result.oblique = copy_to_host(oblique, static_cast(oblique_size)); + const auto measured = copy_to_host(widths, 3U); std::copy(measured.begin(), measured.end(), result.target_width.begin()); + return result; } - 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); - check(cuMemsetD32(oblique, 0U, static_cast(oblique_size)), "clear"); - launch(gpu.project_oblique, voxels, volume, gx, gy, gz, nx, nz, voxels, view.u[0], view.u[1], view.u[2], - view.v[0], view.v[1], view.v[2], view.u0, view.v0, 1.0 / view.pixel_m, view.nu, view.nv, oblique); - result.front = download(front, static_cast(front_size)); - result.side = download(side, static_cast(side_size)); - result.top = download(top, static_cast(top_size)); - result.oblique = download(oblique, static_cast(oblique_size)); - return result; -} +}; #else -struct Gpu::Impl {}; +// Built without CUDA the constructor throws, so these are never reached. +struct Gpu::Impl { + [[noreturn]] static void missing() { throw CudaUnavailable("SAR3D requires an NVIDIA GPU with CUDA"); } + auto inverse_dft(std::span>, std::size_t, std::size_t, std::size_t) + -> std::vector> { missing(); } + auto image(const VolumeRequest&) -> Peak { missing(); } + auto download() -> std::vector { missing(); } + void load(const Axes&, std::span) { missing(); } + auto project(const std::array&, const ObliqueView&) -> Projections { missing(); } +}; #endif -Gpu::Gpu(double smooth_sigma_voxels) { +Gpu::Gpu() { #ifdef SAR3D_WITH_CUDA - impl_ = std::make_unique(gaussian_taps(smooth_sigma_voxels)); + impl_ = std::make_unique(); #else - (void)smooth_sigma_voxels; throw CudaUnavailable("SAR3D requires an NVIDIA GPU with CUDA; this program was built without " "nvrtc.h, cuda.h, libnvrtc or libcuda. Install CUDA and rebuild."); #endif @@ -577,7 +670,6 @@ Gpu::Gpu(double smooth_sigma_voxels) { Gpu::~Gpu() = default; -// Without CUDA the constructor throws, so the #else branches below only satisfy the compiler. auto Gpu::inverse_dft(std::span> bands, std::size_t points, std::size_t first_bin, std::size_t size) -> std::vector> { if (points == 0U || bands.empty() || bands.size() % points != 0U || !std::has_single_bit(size) || @@ -585,36 +677,20 @@ auto Gpu::inverse_dft(std::span> bands, std::size_t po throw std::runtime_error("sar3d invalid inverse FFT dimensions"); } (void)checked_product(bands.size() / points, size, kMaxIndexedItems, "A-scan stack"); -#ifdef SAR3D_WITH_CUDA return impl_->inverse_dft(bands, points, first_bin, size); -#else - throw CudaUnavailable("SAR3D requires an NVIDIA GPU with CUDA"); -#endif } -auto Gpu::project(const VolumeRequest& request) -> Projections { - const auto nx = request.x.size(), ny = request.y.size(), nz = request.z.size(); - const auto plane = checked_product(nx, nz, kMaxIndexedItems, "grid plane"); - (void)checked_product(plane, ny, kMaxIndexedItems, "grid (32-bit indexing; increase SAR_VOXEL_M)"); - const auto& view = request.oblique; - if (nx < 2U || ny < 2U || nz < 2U || request.samples_per_ascan < 2U || - request.contributions.empty() || request.samples_per_metre <= 0.0 || - !std::isfinite(request.samples_per_metre) || request.samples_per_metre > std::numeric_limits::max() || - view.nu <= 0 || view.nv <= 0 || !(view.pixel_m > 0.0) || !std::isfinite(view.pixel_m)) { - throw std::runtime_error("sar3d invalid grid, aperture, A-scan or oblique dimensions"); +auto Gpu::image(const VolumeRequest& request) -> Peak { + check_grid(request.grid); + if (request.samples_per_ascan < 2U || request.contributions.empty() || request.samples_per_metre <= 0.0 || + !std::isfinite(request.samples_per_metre) || request.samples_per_metre > std::numeric_limits::max()) { + throw std::runtime_error("sar3d invalid aperture or A-scan dimensions"); } const auto scan_count = checked_product(request.contributions.size(), request.samples_per_ascan, kMaxIndexedItems, "A-scan stack"); if (request.ascans.size() != scan_count) { throw std::runtime_error("sar3d A-scan stack does not match the aperture"); } - for (const auto axis : {request.x, request.y, request.z}) { - for (const double coordinate : axis) { - if (!std::isfinite(coordinate) || std::abs(coordinate) > std::numeric_limits::max()) { - throw std::runtime_error("sar3d grid contains an unsupported coordinate"); - } - } - } for (const auto& term : request.contributions) { if (!std::isfinite(term.weight) || term.weight < 0.0F || !std::isfinite(term.inv_reference_ranges) || term.inv_reference_ranges <= 0.0F) { @@ -626,13 +702,32 @@ auto Gpu::project(const VolumeRequest& request) -> Projections { } } } + return impl_->image(request); +} + +auto Gpu::download() -> std::vector { + return impl_->download(); +} + +void Gpu::load(const Axes& grid, std::span values) { + check_grid(grid); + if (values.size() != grid.x.size() * grid.y.size() * grid.z.size()) { + throw std::runtime_error("sar3d volume does not match its grid"); + } + // The oblique projection orders values by their bit patterns, which needs them non-negative. + if (!std::ranges::all_of(values, [](float value) { return std::isfinite(value) && value >= 0.0F; })) { + throw std::runtime_error("sar3d volume must be finite and non-negative"); + } + impl_->load(grid, values); +} + +auto Gpu::project(const std::array& peak, const ObliqueView& view) -> Projections { + if (view.nu <= 0 || view.nv <= 0 || !(view.pixel_m > 0.0) || !std::isfinite(view.pixel_m)) { + throw std::runtime_error("sar3d invalid oblique view"); + } (void)checked_product(static_cast(view.nu), static_cast(view.nv), kMaxIndexedItems, "oblique image"); -#ifdef SAR3D_WITH_CUDA - return impl_->project(request); -#else - throw CudaUnavailable("SAR3D requires an NVIDIA GPU with CUDA"); -#endif + return impl_->project(peak, view); } } // namespace radar::processing::sar3d 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 623c15c..5783eae 100644 --- a/data_acq_and_processing/processing/processors/src/sar3d_gpu.hpp +++ b/data_acq_and_processing/processing/processors/src/sar3d_gpu.hpp @@ -9,8 +9,8 @@ #include // 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, peak and its axial half-maximum widths leave the backend. +// backprojection of an aperture into a smoothed 3D volume. The volume stays on the +// device for its projections and is downloaded only when its voxels are needed. namespace radar::processing::sar3d { class CudaUnavailable : public std::runtime_error { @@ -40,44 +40,57 @@ struct ObliqueView { int nv = 0; }; -struct VolumeRequest { - // Scene axes [m]: X across, Y height above the ground, Z along the ground. +// Scene axes of a volume [m]: X across, Y height above the ground, Z along the ground. +// Volumes are [y][z][x] with x fastest. +struct Axes { std::span x; std::span y; std::span z; +}; + +// Whether Y plane `iy` is clear of the grid's top and bottom `margin`: a constant-range +// ridge leaves the grid there and would set a false maximum (sar_norm_reference_mask). +[[nodiscard]] auto reference_plane(std::span y, double margin, std::size_t iy) -> bool; + +struct VolumeRequest { + Axes grid{}; double pitch_rad = 0.0; // radar tilt that maps scene axes onto radar axes double ground_depth_m = 0.0; // height of the radar axes origin above the ground double min_depth_m = 0.0; // radar-depth gate; voxels outside it stay empty double max_depth_m = 0.0; - std::span contributions; + std::span contributions{}; // One A-scan per contribution, back to back; sample k is the Tx+Rx path k / samples_per_metre. - std::span> ascans; + std::span> ascans{}; std::size_t samples_per_ascan = 0U; double samples_per_metre = 0.0; bool use_compensation = true; float comp_range_power = 0.0F; float comp_weight_max = 0.0F; - double norm_edge_margin_m = 0.0; - ObliqueView oblique{}; + double norm_edge_margin_m = 0.0; // SAR_NORM_EDGE_MARGIN_M of reference_plane + double smooth_sigma_voxels = 0.0; // Gaussian over Z and X }; -// Maximum projections of the smoothed volume (row-major), before normalization. +// Maximum of a volume over its reference planes. +struct Peak { + float level = 0.0F; // of the smoothed volume + float raw_level = 0.0F; // before smoothing + std::array index{}; // first voxel at `level`, [y, z, x] +}; + +// Maximum projections of the current volume (row-major), unscaled. struct Projections { std::vector front; // [y][x], maximum along Z std::vector side; // [y][z], maximum along X std::vector top; // [z][x], maximum along Y std::vector oblique; // [v][u] - 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 + std::array target_width{}; // contiguous FWHM through the given peak, [X, Y, Z], metres }; // CUDA context and runtime-compiled kernels. A missing device raises // CudaUnavailable; computation and compilation failures remain ordinary errors. class Gpu { public: - explicit Gpu(double smooth_sigma_voxels); + Gpu(); ~Gpu(); Gpu(const Gpu&) = delete; auto operator=(const Gpu&) -> Gpu& = delete; @@ -91,8 +104,18 @@ class Gpu { std::size_t size ) -> std::vector>; - // Backprojects, applies the Gaussian smoothing over (Z, X) and projects. - [[nodiscard]] auto project(const VolumeRequest& request) -> Projections; + // Backprojects and smooths over (Z, X) into the current volume; returns its maximum over + // the reference planes. + auto image(const VolumeRequest& request) -> Peak; + + // The current volume, [y][z][x]. + [[nodiscard]] auto download() -> std::vector; + + // Makes non-negative host values on `grid` the current volume. + void load(const Axes& grid, std::span values); + + // Max projections of the current volume and its half-maximum widths through `peak` ([y, z, x]). + [[nodiscard]] auto project(const std::array& peak, const ObliqueView& view) -> Projections; private: struct Impl; 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 09a7457..4f62d83 100644 --- a/data_acq_and_processing/processing/processors/src/sar3d_processor.cpp +++ b/data_acq_and_processing/processing/processors/src/sar3d_processor.cpp @@ -9,23 +9,28 @@ #include #include #include +#include #include #include #include #include #include +#include #include #include +#include #include #include +#include #include #include #include +#include "sar3d_clean.hpp" #include "sar3d_gpu.hpp" -// 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 +// Port of the 3D SAR part of the reference notebook Multi_depthnorm_clean.ipynb +// (cells 1 and 5). Comments name the notebook functions and constants each step // follows; the parameters keep the notebook names in the GUI. namespace radar::processing { namespace { @@ -62,7 +67,7 @@ struct Settings { double look_angle_rad = 0.0; // LOOK_ANGLE_DEG double ground_depth_m = 0.0; // BP_GROUND_DEPTH_M std::size_t frame_radius = 0U; // SAR_FRAME_RADIUS - double voxel_m = 0.0; // SAR_VOXEL_M, one edge for X, Y and Z + std::array voxel_m{}; // SAR_VOXEL_M, steps along X, Y and Z float comp_range_power = 0.0F; // COMP_RANGE_POWER bool use_compensation = true; // SAR_USE_COMPENSATION bool normalize_each_frame = true; // SAR_NORMALIZE_EACH_FRAME @@ -72,6 +77,8 @@ struct Settings { 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 + bool bg_gate_enabled = true; // BG_GATE_ENABLE + bool clean_enabled = true; // SAR_CLEAN_ENABLE // Everything the cached A-scans depend on. [[nodiscard]] auto ascan_key() const { @@ -126,7 +133,7 @@ void select_antennas( settings.look_angle_rad = static_cast(live.gpr_look_angle_deg) * (pi / 180.0); settings.ground_depth_m = live.sar3d_bp_ground_depth_m; settings.frame_radius = live.sar3d_frame_radius; - settings.voxel_m = live.sar3d_voxel_m; + settings.voxel_m = {live.sar3d_voxel_x_m, live.sar3d_voxel_y_m, live.sar3d_voxel_z_m}; 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; @@ -137,6 +144,8 @@ void select_antennas( 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; + settings.bg_gate_enabled = live.sar3d_bg_gate_enabled; + settings.clean_enabled = live.sar3d_clean_enabled; 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"); } @@ -146,23 +155,20 @@ void select_antennas( if (!(settings.max_depth_m > settings.min_depth_m)) { throw std::runtime_error("MAX_DEPTH must be above MIN_DEPTH"); } - if (!(settings.voxel_m > 0.0)) { - throw std::runtime_error("SAR_VOXEL_M must be positive"); + if (!std::ranges::all_of(settings.voxel_m, [](double step) { return step > 0.0; })) { + throw std::runtime_error("SAR_VOXEL_M must be positive along X, Y and Z"); } return settings; } -// numpy.percentile with linear interpolation. -[[nodiscard]] auto percentile(std::vector values, double q) -> double { - const double position = q / 100.0 * static_cast(values.size() - 1U); - const auto lower = static_cast(position); - const auto lower_it = values.begin() + static_cast(lower); - std::nth_element(values.begin(), lower_it, values.end()); - if (lower + 1U == values.size()) { - return *lower_it; +// numpy.hanning (BP_WINDOW). +[[nodiscard]] auto hann(std::size_t count) -> std::vector { + std::vector window(count); + const auto last = static_cast(count - 1U); + for (std::size_t index = 0U; index < count; ++index) { + window[index] = 0.5 - 0.5 * std::cos(2.0 * pi * static_cast(index) / last); } - const double upper = *std::min_element(lower_it + 1, values.end()); - return *lower_it + (position - static_cast(lower)) * (upper - *lower_it); + return window; } struct Band { @@ -197,7 +203,7 @@ struct Band { steps.push_back(static_cast(frequency[points[index + 1U]]) - frequency[points[index]]); } Band band{}; - band.df_hz = percentile(std::move(steps), 50.0); + band.df_hz = bp::median(std::move(steps)); if (!(band.df_hz > 0.0)) { throw std::runtime_error("trace frequencies do not increase"); } @@ -219,17 +225,13 @@ struct Band { 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); + const auto window = hann(points.size()); 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 auto sample = samples[point]; + const auto sample = samples[points[index]]; if (!std::isfinite(sample.real()) || !std::isfinite(sample.imag())) { throw std::runtime_error("S21 contains a non-finite sample"); } - band.spectrum.emplace_back( - sample * hann - ); + band.spectrum.emplace_back(sample * window[index]); } return band; } @@ -240,6 +242,8 @@ struct Frame { double dt_s = 0.0; std::size_t samples = 0U; // per A-scan std::vector> ascans{}; // [pair][sample] + // A-scan of a point reflector at zero delay, the same band and window (sar_clean_unit_pulse). + std::vector> pulse{}; }; [[nodiscard]] auto compute_frame( @@ -275,11 +279,17 @@ struct Frame { } spectra.insert(spectra.end(), band.spectrum.begin(), band.spectrum.end()); } + // The pulse is the A-scan of a flat unit spectrum, computed after the pairs. + for (const double weight : hann(first.spectrum.size())) { + spectra.emplace_back(static_cast(weight)); + } Frame frame{}; frame.center_ns = timing.center_ns; 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); + frame.pulse.assign(frame.ascans.end() - static_cast(frame.samples), frame.ascans.end()); + frame.ascans.resize(frame.ascans.size() - frame.samples); 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)); @@ -316,9 +326,12 @@ struct Frame { return axis; } -[[nodiscard]] auto checked_axis_size(double count) -> std::size_t { +// Every kernel uses signed 32-bit indices; `remedy` tells how to stay within them. +constexpr std::string_view kVoxelRemedy = "increase SAR_VOXEL_M"; + +[[nodiscard]] auto checked_axis_size(double count, std::string_view remedy) -> std::size_t { if (!std::isfinite(count) || count < 1.0 || count > std::numeric_limits::max()) { - throw std::runtime_error("SAR axis exceeds 32-bit kernel indexing; increase SAR_VOXEL_M"); + throw std::runtime_error("SAR axis exceeds 32-bit kernel indexing; " + std::string(remedy)); } return static_cast(count); } @@ -327,12 +340,19 @@ struct Grid { std::vector x{}; std::vector y{}; std::vector z{}; + + [[nodiscard]] auto axes() const -> sar3d::Axes { + return {x, y, z}; + } }; -// Scene grid (sar_build_full_region_grid, 'ground' axes): Y = 0 is the ground, +// Bounds [low, high] of the scene along X, Y and Z. +using Region = std::array, 3>; + +// Scene region (sar_build_full_region_grid, 'ground' axes): Y = 0 is the ground, // Z runs along the ground and spans every point of the window whose radar depth // is inside [MIN_DEPTH, MAX_DEPTH] at some height of the window. -[[nodiscard]] auto scene_grid(const Settings& settings) -> Grid { +[[nodiscard]] auto scene_region(const Settings& settings) -> Region { const double cos_p = std::cos(settings.look_angle_rad); const double sin_p = std::sin(settings.look_angle_rad); if (std::abs(cos_p) < 1e-6) { @@ -357,24 +377,47 @@ 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; } + return {{{x_min, x_max}, {settings.y_min_m, settings.y_max_m}, {z_min, z_max}}}; +} + +// Region of the frame grid: the scene region with X and Z resized about their centres to +// SAR_X_SIZE_M and SAR_Z_SIZE_M where those are set. +[[nodiscard]] auto frame_region(Region region, const Settings& settings) -> Region { + const auto resize = [](std::array& bounds, double size) { + if (size > 0.0) { + const double center = (bounds[0] + bounds[1]) / 2.0; + bounds = {center - size / 2.0, 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(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"); + resize(region[0], settings.x_size_m); + resize(region[2], settings.z_size_m); + return region; +} + +// Region of the CLEAN search: the scene region around the frame grid (sar_clean_search_grid), +// grown to a frame that SAR_X_SIZE_M or SAR_Z_SIZE_M widen past it. +[[nodiscard]] auto search_region(const Region& scene, const Region& frame) -> Region { + Region region{}; + for (std::size_t axis = 0U; axis < 3U; ++axis) { + region[axis] = {std::min(scene[axis][0], frame[axis][0]), std::max(scene[axis][1], frame[axis][1])}; } - return {linspace(x_min, x_max, nx), linspace(settings.y_min_m, settings.y_max_m, ny), - linspace(z_min, z_max, nz)}; + return region; +} + +// Grid over the region with the given steps along X, Y and Z (sar_build_full_region_grid). +// Oversized products are rejected before the axes are allocated. +[[nodiscard]] auto region_grid(const Region& region, const std::array& step, std::string_view remedy) + -> Grid { + std::array count{}; + for (std::size_t axis = 0U; axis < 3U; ++axis) { + const auto [low, high] = region[axis]; + count[axis] = checked_axis_size(std::max(std::nearbyint((high - low) / step[axis]) + 1.0, 2.0), remedy); + } + if (count[0] > static_cast(std::numeric_limits::max()) / count[1] / count[2]) { + throw std::runtime_error("SAR grid exceeds 32-bit kernel indexing; " + std::string(remedy)); + } + return {linspace(region[0][0], region[0][1], count[0]), linspace(region[1][0], region[1][1], count[1]), + linspace(region[2][0], region[2][1], count[2])}; } // Cosine taper over the aperture: 1 at the centre frame, SAR_EDGE_FRAME_WEIGHT at @@ -456,10 +499,10 @@ struct Oblique { } } } - view.nu = static_cast(checked_axis_size(std::round((u_max - view.u0) / view.pixel_m) + 1.0)); - view.nv = static_cast(checked_axis_size(std::round((v_max - view.v0) / view.pixel_m) + 1.0)); + view.nu = static_cast(checked_axis_size(std::round((u_max - view.u0) / view.pixel_m) + 1.0, kVoxelRemedy)); + view.nv = static_cast(checked_axis_size(std::round((v_max - view.v0) / view.pixel_m) + 1.0, kVoxelRemedy)); if (view.nu > std::numeric_limits::max() / view.nv) { - throw std::runtime_error("SAR oblique image exceeds 32-bit kernel indexing; increase SAR_VOXEL_M"); + throw std::runtime_error("SAR oblique image exceeds 32-bit kernel indexing; " + std::string(kVoxelRemedy)); } for (int index = 0; index < view.nu; ++index) { oblique.u.push_back(view.u0 + index * view.pixel_m); @@ -512,6 +555,21 @@ struct Oblique { return payload; } +// CLEAN part of the status line (sar_display_note, sar_clean_report). +[[nodiscard]] auto clean_note(const sar3d::Clean& clean) -> std::string { + std::ostringstream note; + note << std::fixed << std::setprecision(1) << "CLEAN: "; + if (clean.objects.empty()) { + note << "no objects"; + } else { + note << clean.objects.size() << (clean.objects.size() == 1U ? " object" : " objects") + << ", each in its own scale"; + } + note << " (" << clean.steps << (clean.steps == 1U ? " subtraction" : " subtractions") + << ", residual maximum " << clean.residual_ratio << " backgrounds)"; + return note.str(); +} + [[nodiscard]] auto status_payload(float code, const std::string& message) -> ipc::ResultPayload { std::vector values{code}; for (const unsigned char byte : message) { @@ -521,6 +579,71 @@ struct Oblique { return table_payload("sar3d_status", columns, std::move(values)); } +// What the views show of one frame image (sar_prepare_display). +struct Display { + sar3d::Projections views; // the image in units of its maximum, or the CLEAN display + std::array peak{}; // marked voxel [y, z, x]: the maximum or the strongest object + std::optional background; // when BG_GATE_ENABLE or CLEAN needs it + std::optional clean; // when CLEAN ran + + [[nodiscard]] auto has_objects() const -> bool { return clean.has_value() && !clean->objects.empty(); } +}; + +// Views of the frame image of `request`, the current volume of `backend` with maximum +// `peak`. Measures the frame background when BG_GATE_ENABLE or CLEAN needs it and, when +// the aperture is the one frame `single` (nullptr for a longer aperture, as sar_clean_frame), +// runs CLEAN, whose objects then replace the image on the views. +[[nodiscard]] auto display_frame( + sar3d::Gpu& backend, + const sar3d::VolumeRequest& request, + const sar3d::Peak& peak, + double normalization, + const sar3d::ObliqueView& view, + const Settings& settings, + const Region& search_area, + const Frame* single +) -> Display { + Display display{.views = backend.project(peak.index, view), .peak = peak.index, .background = {}, .clean = {}}; + const bool run_clean = settings.clean_enabled && single != nullptr; + if (peak.level > 0.0F && (settings.bg_gate_enabled || run_clean)) { + const auto volume = backend.download(); + display.background = sar3d::frame_background(volume, volume, request); + if (run_clean) { + // The search spans the whole scene region, so responses just outside a narrowed frame + // are subtracted too. + const Grid search = region_grid(search_area, sar3d::kCleanSearchStep, + "narrow the scene or turn SAR_CLEAN_ENABLE off"); + display.clean = sar3d::clean(backend, { + .frame = request, + .volume = volume, + .level = normalization, + .frame_background = *display.background, + .search = search.axes(), + .ascans = single->ascans, + .pulse = single->pulse, + .bg_gate = settings.bg_gate_enabled, + }); + if (display.has_objects()) { + // CLEAN left one of its own images current; the display replaces it. + display.peak = display.clean->objects.front().index; + backend.load(request.grid, display.clean->display); + display.views = backend.project(display.peak, view); + return display; + } + } + } + auto& views = display.views; + for (auto* image : {&views.front, &views.side, &views.top, &views.oblique}) { + for (auto& value : *image) { + value = static_cast(static_cast(value) / normalization); + if (!std::isfinite(value)) { + throw std::runtime_error("SAR normalization exceeds the float32 range"); + } + } + } + return display; +} + } // namespace struct Sar3dProcessor::State { @@ -570,7 +693,7 @@ struct Sar3dProcessor::State { auto gpu() -> sar3d::Gpu& { if (device == nullptr && device_error == nullptr) { - try { device = std::make_unique(kSmoothSigma); } + try { device = std::make_unique(); } catch (...) { device_error = std::current_exception(); } } if (device == nullptr) { std::rethrow_exception(device_error); } @@ -617,15 +740,17 @@ struct Sar3dProcessor::State { offsets_s.push_back(1e-9 * static_cast(delta_ns)); } const auto weights = frame_weights(offsets_s); - const Grid grid = scene_grid(settings); + const Region scene = scene_region(settings); + const Region frame_area = frame_region(scene, settings); + const Grid grid = region_grid(frame_area, settings.voxel_m, kVoxelRemedy); const double cos_p = std::cos(settings.look_angle_rad); const double sin_p = std::sin(settings.look_angle_rad); std::vector> corners{}; // of the scene box, in radar axes for (const double x : {grid.x.front(), grid.x.back()}) { for (const double y : {grid.y.front(), grid.y.back()}) { for (const double z : {grid.z.front(), grid.z.back()}) { - const double y0 = y - settings.ground_depth_m; - corners.push_back({x, y0 * cos_p + z * sin_p, -y0 * sin_p + z * cos_p}); + const auto [y_radar, z_radar] = bp::radar_yz(y, z, settings.ground_depth_m, settings.look_angle_rad); + corners.push_back({x, y_radar, z_radar}); } } } @@ -666,11 +791,8 @@ struct Sar3dProcessor::State { } } - const Oblique oblique = oblique_view(grid, settings.voxel_m); - auto images = backend.project({ - .x = grid.x, - .y = grid.y, - .z = grid.z, + const sar3d::VolumeRequest request{ + .grid = grid.axes(), .pitch_rad = settings.look_angle_rad, .ground_depth_m = settings.ground_depth_m, .min_depth_m = settings.min_depth_m, @@ -683,44 +805,72 @@ struct Sar3dProcessor::State { .comp_range_power = settings.comp_range_power, .comp_weight_max = kCompRangeWeightMax, .norm_edge_margin_m = kNormEdgeMargin, - .oblique = oblique.view, - }); - + .smooth_sigma_voxels = kSmoothSigma, + }; // Normalization level and maximum over the Y planes clear of the window // edges, where a constant-range ridge leaves the grid and would set a // false maximum (sar_norm_reference_mask, sar_finish_volume, // sar_global_max_index). The level is the volume maximum, so the - // Values at excluded Y edges can exceed 1, as in sar_finish_volume. - const float level = images.level; - const auto [peak_y, peak_z, peak_x] = images.peak; + // values at excluded Y edges can exceed 1, as in sar_finish_volume. + const auto peak = backend.image(request); // The notebook normalizes before and after the linear Gaussian pass. // Combining both factors in double preserves its epsilons and avoids // float32 reciprocal overflow for a weak, nonzero signal. - const double normalization = static_cast(level) + (images.raw_level + 1e-12) * 1e-12; - for (auto* image : {&images.front, &images.side, &images.top, &images.oblique}) { - for (auto& value : *image) { - value = static_cast(static_cast(value) / normalization); - if (!std::isfinite(value)) { throw std::runtime_error("SAR normalization exceeds the float32 range"); } - } - } + const double normalization = static_cast(peak.level) + (peak.raw_level + 1e-12) * 1e-12; + // Pixels of the largest voxel step leave no gaps between projected voxels. + const Oblique oblique = oblique_view(grid, std::ranges::max(settings.voxel_m)); + auto display = display_frame(backend, request, peak, normalization, oblique.view, settings, + search_region(scene, frame_area), aperture.size() == 1U ? newest.get() : nullptr); std::vector payloads{}; - payloads.push_back(image_payload("sar3d_front", grid.x, grid.y, std::move(images.front))); - payloads.push_back(image_payload("sar3d_side", grid.z, grid.y, std::move(images.side))); - payloads.push_back(image_payload("sar3d_top", grid.x, grid.z, std::move(images.top))); - payloads.push_back(image_payload("sar3d_oblique", oblique.u, oblique.v, std::move(images.oblique))); + payloads.push_back(image_payload("sar3d_front", grid.x, grid.y, std::move(display.views.front))); + payloads.push_back(image_payload("sar3d_side", grid.z, grid.y, std::move(display.views.side))); + payloads.push_back(image_payload("sar3d_top", grid.x, grid.z, std::move(display.views.top))); + payloads.push_back(image_payload("sar3d_oblique", oblique.u, oblique.v, std::move(display.views.oblique))); std::vector basis(oblique.view.u.begin(), oblique.view.u.end()); 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])})); + if (!(peak.level > 0.0F)) { + payloads.push_back(status_payload(1.0F, "No signal inside the selected depth range")); + return payloads; } - payloads.push_back(status_payload(level > 0.0F ? 0.0F : 1.0F, - level > 0.0F ? "CUDA" : "No signal inside the selected depth range")); + const auto [peak_y, peak_z, peak_x] = display.peak; + 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]), + })); + if (settings.bg_gate_enabled) { + payloads.push_back(table_payload("sar3d_background", 3U, { + static_cast(display.background.value().ratio()), + static_cast(sar3d::kBgGateShow), static_cast(sar3d::kBgGateFull), + })); + } + if (display.has_objects()) { + // The objects replace the peak marker, so its width is not published. + std::vector object_rows{}; + for (const auto& object : display.clean->objects) { + for (const double coordinate : object.position) { object_rows.push_back(static_cast(coordinate)); } + object_rows.push_back(static_cast(object.level)); + } + payloads.push_back(table_payload("sar3d_objects", 4U, std::move(object_rows))); + std::vector mark_rows{}; + for (const auto& mark : sar3d::marks(display.clean->objects, request)) { + mark_rows.push_back(static_cast(mark.kind)); + for (const double edge : mark.box) { mark_rows.push_back(static_cast(edge)); } + } + payloads.push_back(table_payload("sar3d_marks", 7U, std::move(mark_rows))); + } else { + const auto& widths = display.views.target_width; + payloads.push_back(table_payload("sar3d_target_width", 3U, {widths[0], widths[1], widths[2]})); + } + std::string status = "CUDA"; + if (display.clean.has_value()) { + status += "; " + clean_note(*display.clean); + } else if (settings.clean_enabled) { + status += "; CLEAN skipped: it needs a one-frame aperture, this one has " + + std::to_string(aperture.size()) + " frames"; + } + payloads.push_back(status_payload(0.0F, status)); return payloads; } diff --git a/data_acq_and_processing/tests/backprojection_fixture.cpp b/data_acq_and_processing/tests/backprojection_fixture.cpp index 2d5d8fa..7a1ddd7 100644 --- a/data_acq_and_processing/tests/backprojection_fixture.cpp +++ b/data_acq_and_processing/tests/backprojection_fixture.cpp @@ -1,9 +1,11 @@ #include "sar3d_processor.hpp" #include "gpr_processor.hpp" +#include "../processing/processors/src/sar3d_clean.hpp" #include #include #include +#include #include #include @@ -40,11 +42,37 @@ auto collection(const Json& source) -> ipc::PreprocessedCollection { } return result; } + +// Boxes around given CLEAN objects, rows [x, y, z, y_low, y_high, level], in a scene with the +// given X and Y bounds and radar placement (sar3d::marks). +auto marks(const Json& input) -> Json { + std::vector objects; + for (const auto& row : input.at("objects")) { + objects.push_back({.index = {}, .position = {row.at(0), row.at(1), row.at(2)}, + .y_low = row.at(3), .y_high = row.at(4), .level = row.at(5)}); + } + const auto x = input.at("x").get>(); + const auto y = input.at("y").get>(); + const processing::sar3d::VolumeRequest frame{ + .grid = {x, y, {}}, + .pitch_rad = input.at("pitch").get() * (std::numbers::pi / 180.0), + .ground_depth_m = input.at("ground"), + .min_depth_m = input.at("min_depth"), + .max_depth_m = input.at("max_depth"), + }; + Json rows = Json::array(); + for (const auto& mark : processing::sar3d::marks(objects, frame)) { + Json row = Json::array({static_cast(mark.kind)}); + for (const double edge : mark.box) { row.push_back(edge); } + rows.push_back(std::move(row)); + } + return Json{{"marks", rows}}; +} } // 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. +// published payloads go to stdout as JSON. A request with mode "marks" boxes given objects instead. int main() { try { const Json source = Json::parse(std::cin); @@ -54,6 +82,10 @@ int main() { processing::Sar3dProcessor sar; processing::GprProcessor gpr; for (const auto& input : requests) { + if (input.value("mode", "sar3d") == "marks") { + results.push_back(marks(input)); + continue; + } 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")) { @@ -62,6 +94,7 @@ int main() { 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)}); } + const bool gpr_mode = input.value("mode", "sar3d") == "gpr"; 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; @@ -69,31 +102,39 @@ int main() { 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; + if (gpr_mode) { + 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; + } else { + const auto& voxel = input.at("voxel"); + live.sar3d_voxel_x_m = voxel.at(0); + live.sar3d_voxel_y_m = voxel.at(1); + live.sar3d_voxel_z_m = voxel.at(2); + 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.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.sar3d_bg_gate_enabled = input.value("bg_gate", live.sar3d_bg_gate_enabled); + live.sar3d_clean_enabled = input.value("clean", live.sar3d_clean_enabled); + } 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" + processing::ProcessorInterface& processor = gpr_mode ? static_cast(gpr) : static_cast(sar); const auto result = processor.process_collection(config, latest, previous, live); Json payloads = Json::array(); diff --git a/data_acq_and_processing/tests/generate_backprojection_reference.py b/data_acq_and_processing/tests/generate_backprojection_reference.py index db90345..1b4d13d 100644 --- a/data_acq_and_processing/tests/generate_backprojection_reference.py +++ b/data_acq_and_processing/tests/generate_backprojection_reference.py @@ -1,6 +1,6 @@ -"""Generate SAR3D/GPR references from notebook cells 1/14 and Horns coherent BP. +"""Generate SAR3D/GPR references from notebook cells 1/5 and Horns coherent BP. -Requires the original, untracked Multi_antenna_demonstrate_depthnorm.ipynb; a missing notebook +Requires the original, untracked Multi_depthnorm_clean.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. @@ -16,6 +16,7 @@ import hashlib import json from pathlib import Path import tempfile +from types import SimpleNamespace import numpy as np import scipy @@ -24,7 +25,7 @@ ROOT = Path(__file__).resolve().parents[2] REFERENCE = Path(__file__).with_name("sar3d_reference.npz") -def reference_functions(source_path: Path, parameters: dict, cells: tuple[int, ...] = (1, 14)) -> tuple[dict, dict]: +def reference_functions(source_path: Path, parameters: dict, cells: tuple[int, ...] = (1, 5)) -> 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__} @@ -56,8 +57,9 @@ def scene(f_start=1e9, f_stop=4e9, points=201) -> dict: tx=[[0, -.4, 0, 0], [1, .45, .2, .05]], 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=2, power=1.5, compensation=True, normalize=True, + voxel=[.4, .4, .4], radius=2, power=1.5, compensation=True, normalize=True, equalize=True, subtract_mean=False, mean_count=1, motion="int_minus", + bg_gate=True, clean=False, ) # Match the float32 settings/geometry accepted by the native fixture. for key in ("tx", "rx"): @@ -79,6 +81,45 @@ def scene(f_start=1e9, f_stop=4e9, points=201) -> dict: return request +# Point reflectors (X, height Y, Z, amplitude) at heights of 0, 0.6 and 0.25 m. +REFLECTORS = [(.3, 0., 2.2, 1.), (-.8, .6, 2.6, .5), (1.3, .25, 1.9, .7)] + + +def point_scene(targets: list[tuple[float, float, float, float]], noise: float, seed: int, **settings) -> dict: + """Settings of one static frame of point reflectors (X, height Y, Z, amplitude) in noise, for CLEAN. + + Three Tx and three Rx at different heights resolve Y in a single frame. The voxel steps + keep the notebook's proportions, 2.5 : 2 : 2. + """ + request = dict( + tx=[[0, -.6, 0., 0.], [1, 0., .35, .02], [2, .6, .1, 0.]], + rx=[[0, -.3, .2, 0.], [1, .3, .45, .03], [2, .9, 0., .01]], + min_depth=1., max_depth=3.5, speed=0., pitch=13., ground=1., voxel=[.0625, .05, .05], radius=0, + clean=True, previous=[], + ) + request.update(settings) + 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) + cos_p, sin_p = np.cos(np.radians(request["pitch"])), np.sin(np.radians(request["pitch"])) + random = np.random.default_rng(seed) + traces = [] + for tx in request["tx"]: + for rx in request["rx"]: + samples = noise * (random.standard_normal(frequency.size) + 1j * random.standard_normal(frequency.size)) + for x, height, z, amplitude in targets: + # Scene to radar axes, as sar_scene_to_radar_yz. + y0 = height - request["ground"] + point = np.array([x, y0 * cos_p + z * sin_p, -y0 * sin_p + z * cos_p]) + path = np.linalg.norm(point - tx[1:]) + np.linalg.norm(point - rx[1:]) + samples += amplitude * np.exp(-2j * np.pi * frequency.astype(float) * path / 3e8) + samples = samples.astype(np.complex64) + traces.append(dict(input=int(rx[0]), output=int(tx[0]), start=10_000_000_000, end=10_165_000_000, + frequency=frequency.tolist(), s21=np.column_stack((samples.real, samples.imag)).tolist())) + request["latest"] = dict(id=1, timestamp=10_165_000_000, traces=traces) + return request + + def write_dataset(directory: Path, request: dict) -> list[Path]: """Supply the notebook's original config_profile/meta/NPY input contract.""" profile = {"gpr": {}, "run": {"combos": []}} @@ -115,6 +156,8 @@ def reference(notebook: Path, request: dict, directory: Path) -> tuple[dict, dic snapshots = write_dataset(directory, request) center = snapshots[-1] parameters = dict( + SAR_Y_RANGE_M=(request.get("y_min", -.5), request.get("y_max", 2.)), + BG_GATE_ENABLE=request["bg_gate"], SAR_CLEAN_ENABLE=request["clean"], SAR_MARK_ENABLE=True, INPUT_IDX=[int(row[0]) for row in request["rx"]], OUTPUT_IDX=[int(row[0]) for row in request["tx"]], F_START=request["f_start"], F_STOP=request["f_stop"], @@ -122,7 +165,7 @@ def reference(notebook: Path, request: dict, directory: Path) -> tuple[dict, dic PAIR_NORM_DEPTH_MIN=request["min_depth"], PAIR_NORM_DEPTH_MAX=request["max_depth"], SPEED_M_S=request["speed"], LOOK_ANGLE_DEG=request["pitch"], BP_GROUND_DEPTH_M=request["ground"], COMP_RANGE_POWER=request["power"], - SAR_VOXEL_M=(request["voxel"],) * 3, SAR_FRAME_RADIUS=request["radius"], + SAR_VOXEL_M=tuple(request["voxel"]), 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, @@ -150,31 +193,76 @@ 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) - 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.) + grid = ns["sar_build_full_region_grid"](config) + for index, axis in ((0, "x"), (2, "z")): + size, step = request.get(f"{axis}_size", 0.), request["voxel"][index] 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[axis] = np.linspace(center_axis - size/2, center_axis + size/2, max(round(size/step)+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])) + # What the processor publishes: the frame volume, or the CLEAN display when CLEAN found + # objects, and the peak it marks; the background gain stays to the GUI. + ns["sar_prepare_display"](result, None, timestamps) + objects = result["clean_objects"] or [] + volume, grid = result["display_volume"] if objects else result["volume"], result["grid"] + iy, iz, ix = result["display_max_idx"] 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]])) + if not objects: + # Width of the response through the marked maximum; objects are marked instead. + 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["target_width"] = np.asarray(widths) + if request["bg_gate"]: + expected["background"] = np.array([result["bg_ratio"], ns["BG_GATE_3D_SHOW"], ns["BG_GATE_3D_FULL"]]) + if result["clean_info"] is not None: + expected["clean"] = np.array([result["clean_info"]["steps"], result["clean_info"]["residual_ratio"]]) + if objects: + kinds = {"low": 0, "ground": 1, "up": 2, "edge": 3} + expected["objects"] = np.array([[item["x"], item["y"], item["z"], item["level_bg"]] for item in objects]) + expected["marks"] = np.array([[kinds[group["kind"]], *group["box"]] for group in result["mark_groups"]]) return expected, hashes +def marks_reference(notebook: Path) -> tuple[dict, np.ndarray]: + """sar_mark_groups on synthetic CLEAN objects, each of which tests one rule of the boxes. + + Returns the fixture's "marks" request and the expected rows [class, X0, X1, Y0, Y1, Z0, Z1]. + """ + scene = dict(x=[-2.6, 2.9], y=[-.5, 2.], pitch=13., ground=1., min_depth=1., max_depth=3.5) + objects = [ # [x, y, z, y_low, y_high, level], strongest first + [0., 0., 2., -.1, .3, 60.], # low + [1., .05, 2.2, -.1, .6, 55.], # centre low, top above SAR_MARK_LOW_TOP_M: ground + [-1., .3, 2.4, .1, .5, 50.], # ground + [-.2, .45, 3., .25, .65, 45.], # up while its weak member (last) stays out of the centre + [-1.6, 0., 1.3, -.1, .3, 40.], # a chain of three: the ends are 1.4 m apart, still one box + [-.9, .1, 1.3, -.1, .35, 35.], + [-.2, .05, 1.3, -.15, .3, 30.], + [-2.5, .2, 2., 0., .4, 28.], # near the X minimum: edge + [2.75, .2, 2., 0., .4, 26.], # near the X maximum: edge + [2.2, .1, 3.2, -.1, .3, 24.], # near MAX_DEPTH: edge + [2.2, 0., .85, -.2, .2, 22.], # near MIN_DEPTH: edge + [1.2, 1.7, 2.8, 1.5, 1.9, 20.], # near the top: edge + [.2, -.4, 3.1, -.5, -.25, 15.], # weak member of the "up" box, below 40 % of its top + ] + ns, _ = reference_functions(notebook, dict(BP_GROUND_DEPTH_M=scene["ground"], MIN_DEPTH=scene["min_depth"], + MAX_DEPTH=scene["max_depth"])) + groups = ns["sar_mark_groups"]( + [dict(zip(("x", "y", "z", "y_lo", "y_hi", "level_bg"), row)) for row in objects], + {axis: np.asarray(scene[axis]) for axis in "xy"}, SimpleNamespace(look_angle_deg=scene["pitch"])) + kinds = {"low": 0, "ground": 1, "up": 2, "edge": 3} + return dict(mode="marks", objects=objects, **scene), np.array([[kinds[g["kind"]], *g["box"]] for g in groups]) + + def previous_mean(request: dict) -> dict: samples = {} if request["subtract_mean"] and request["mean_count"]: @@ -258,17 +346,32 @@ def generate(notebook: Path, output: Path) -> None: cases = [ dict(compensation=True, normalize=True), 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), + # Three different voxel steps: each axis keeps its own. + dict(compensation=False, normalize=True, f_start=1.3e9, f_stop=3.7e9, voxel=[.45, .4, .35]), 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), + # CLEAN stops below its thresholds; the background gate hides the residual. + point_scene(REFLECTORS, noise=.02, seed=1), + # CLEAN stops at SAR_CLEAN_MAX_STEPS; without the gate the residual stays visible. + point_scene(REFLECTORS, noise=.01, seed=3, bg_gate=False), + # The stop at the background of each range (SAR_CLEAN_STOP_RANGE_BG) adds two steps to + # those of the stop at the background of the whole volume. + point_scene([(.2, 0., 2., 1.), (-.5, .1, 4.4, .08)], noise=.01, seed=4, max_depth=5.), + # A narrowed frame: CLEAN still searches the whole scene, so the reflector at X = -0.8 m + # just outside the frame is subtracted instead of its sidelobe becoming an object. + point_scene(REFLECTORS, noise=.02, seed=19, x_size=1.6, z_size=.8), ] arrays, configurations = {}, [] for index, settings in enumerate(cases): request = scene() request.update(settings) - for key in ("min_depth", "max_depth", "speed", "pitch", "ground", "voxel", "power"): - request[key] = float(np.float32(request[key])) + # Match the float32 live settings of the native fixture. + for key in ("min_depth", "max_depth", "speed", "pitch", "ground", "power", + "x_size", "z_size", "y_min", "y_max"): + if key in request: + request[key] = float(np.float32(request[key])) + request["voxel"] = [float(np.float32(step)) for step in request["voxel"]] if index == 1: # A separate, unequal capture window per output instead of interleaving. for member in [*request["previous"], request["latest"]]: @@ -288,10 +391,13 @@ 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=(",", ":"))) + marks_request, arrays["marks"] = marks_reference(notebook) + arrays["marks_request"] = np.array(json.dumps(marks_request)) arrays["provenance"] = np.array(json.dumps(dict( - notebook=notebook.name, cell_source_sha256=hashes, cells=[1, 14], cases=len(cases), + notebook=notebook.name, cell_source_sha256=hashes, cells=[1, 5], cases=len(cases), numpy=np.__version__, scipy=scipy.__version__, - method="Unchanged notebook AST functions; per-frame NumPy previous-mean subtraction; max projections of sar_compute_volume output", + method="Unchanged notebook AST functions; per-frame NumPy previous-mean subtraction; " + "max projections of the sar_prepare_display volume, its background, CLEAN objects and marks", ), sort_keys=True)) output.parent.mkdir(parents=True, exist_ok=True) np.savez_compressed(output, **arrays) @@ -300,7 +406,7 @@ 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_demonstrate_depthnorm.ipynb") + parser.add_argument("--notebook", type=Path, default=ROOT / "Multi_depthnorm_clean.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")) diff --git a/data_acq_and_processing/tests/gpr_reference.npz b/data_acq_and_processing/tests/gpr_reference.npz index c137e87..63546b3 100644 Binary files a/data_acq_and_processing/tests/gpr_reference.npz 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 9a0dfad..287bd07 100644 --- a/data_acq_and_processing/tests/native_regressions.cpp +++ b/data_acq_and_processing/tests/native_regressions.cpp @@ -210,9 +210,12 @@ void test_configuration(const std::string& seed_path, const Workspace& workspace } const auto live_path = workspace.path / "live.json"; std::ofstream(live_path) << Json{{"bscan_gain", 2}, {"settings_revision", 11}, + {"sar3d_bg_gate_enabled", false}, {"sar3d_clean_enabled", false}, {"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 loaded = loader.refresh_if_needed(); + require(loaded.bscan_gain == 2 && !loaded.sar3d_bg_gate_enabled && !loaded.sar3d_clean_enabled, + "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}}, diff --git a/data_acq_and_processing/tests/sar3d_reference.npz b/data_acq_and_processing/tests/sar3d_reference.npz index 841b989..8a5edc4 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_backprojection.py b/data_acq_and_processing/tests/test_backprojection.py index 47de655..b23e5c7 100644 --- a/data_acq_and_processing/tests/test_backprojection.py +++ b/data_acq_and_processing/tests/test_backprojection.py @@ -2,8 +2,9 @@ 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_backprojection_reference.py. +warning and the box rules, which are host code, are checked, and the remaining tests +are skipped. The notebook is needed only to regenerate the references with +generate_backprojection_reference.py. """ from __future__ import annotations @@ -11,6 +12,7 @@ import copy import json import os from pathlib import Path +import re import subprocess import unittest @@ -90,22 +92,80 @@ 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) + widths = self.reference.get(f"case_{index}_target_width") + if widths is None: + self.assertNotIn("sar3d_target_width", payloads) + else: + np.testing.assert_allclose(payloads["sar3d_target_width"]["table"], widths, 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) np.testing.assert_allclose(peak[0], expected_peak[0], rtol=2e-3, atol=0) + self.check_display(index, payloads) + + def check_display(self, index: int, payloads: dict[str, dict]) -> None: + """Background, CLEAN objects and boxes against sar_prepare_display. + + Steps at the noise level may pick a neighbouring A-scan sample where two samples are + equal to float32 precision; the reference scenes were checked to keep their decisions + with the kernels compiled without fused multiply-adds. + """ + background = self.reference.get(f"case_{index}_background") + if background is None: + self.assertNotIn("sar3d_background", payloads) + else: + np.testing.assert_allclose(payloads["sar3d_background"]["table"], background, rtol=2e-3, atol=0) + objects = self.reference.get(f"case_{index}_objects") + if objects is None: + self.assertNotIn("sar3d_objects", payloads) + self.assertNotIn("sar3d_marks", payloads) + else: + table = np.reshape(payloads["sar3d_objects"]["table"], (-1, 4)) + np.testing.assert_allclose(table[:, :3], objects[:, :3], rtol=0, atol=2e-6) + np.testing.assert_allclose(table[:, 3], objects[:, 3], rtol=2e-3, atol=0) + np.testing.assert_allclose(np.reshape(payloads["sar3d_marks"]["table"], (-1, 7)), + self.reference[f"case_{index}_marks"], rtol=0, atol=2e-6) + clean = self.reference.get(f"case_{index}_clean") + if clean is not None: + message = status(payloads)[1] + match = re.search(r"\((\d+) subtractions?, residual maximum ([\d.]+) backgrounds\)", message) + assert match is not None, message + self.assertEqual(int(match[1]), clean[0]) + self.assertAlmostEqual(float(match[2]), clean[1], delta=0.051) + + def test_clean_needs_a_one_frame_aperture(self) -> None: + # The actual aperture decides, as in sar_clean_frame: three frames skip CLEAN, while the + # first frame of a run is a one-frame aperture even with SAR_FRAME_RADIUS = 2. + skipped = run(self.scene(clean=True)) + code, message = status(skipped) + self.assertEqual(code, 0, message) + self.assertIn("CLEAN skipped: it needs a one-frame aperture, this one has 3 frames", message) + self.assertNotIn("sar3d_objects", skipped) + code, message = status(run(self.scene(clean=True, previous=[]))) + self.assertEqual(code, 0, message) + self.assertIn("; CLEAN: ", message) + + def test_boxes_follow_the_notebook_rules(self) -> None: + # Host code only, so it runs without a GPU as well. + request = json.loads(self.reference["marks_request"].item()) + output = subprocess.run([str(FIXTURE)], input=json.dumps(request), text=True, + capture_output=True, timeout=30, check=True) + np.testing.assert_allclose(json.loads(output.stdout)["marks"], self.reference["marks"], rtol=0, atol=1e-9) def test_invalid_input_is_an_error_with_a_reason(self) -> None: def no_traces(request): request["latest"]["traces"] = [] def empty_window(request): request["latest"]["traces"][0]["end"] = request["latest"]["traces"][0]["start"] def falling_frequency(request): request["latest"]["traces"][0]["frequency"].reverse() + def overflow(request): + for trace in request["latest"]["traces"]: + samples = np.asarray(trace["s21"]) + trace["s21"] = (samples / np.abs(samples).max() * 1e38).tolist() for name, settings, damage in ( ("no traces", {}, no_traces), ("empty capture window", {}, empty_window), ("falling frequency", {}, falling_frequency), - ("grid beyond 32-bit indexing", {"voxel": 1e-8}, None), + ("volume beyond float32", {"previous": [], "radius": 0, "normalize": False, "equalize": False}, overflow), + ("grid beyond 32-bit indexing", {"voxel": [1e-8] * 3}, None), ("vertical look angle", {"pitch": 90}, None), ("result larger than the ring slot", {"slot": 1024}, None), ): @@ -117,6 +177,20 @@ class Sar3dTest(unittest.TestCase): self.assertEqual(code, 2, message) self.assertTrue(message) + def test_views_are_normalized_at_the_peak_of_the_reference_planes(self) -> None: + # Fine grids spread every Y plane over all threads of the per-plane maximum. + for voxel in (.05, .06, .07, .08, .09, .11): + with self.subTest(voxel=voxel): + payloads = run(self.scene(voxel=[voxel] * 3)) + self.assertEqual(status(payloads)[0], 0, status(payloads)) + front = payloads["sar3d_front"] + x, y = np.asarray(front["x"]), np.asarray(front["y"]) + image = np.where(((y >= y[0] + .3) & (y <= y[-1] - .3))[:, None], + np.reshape(front["image"], (len(y), len(x))), -1.) + row, column = np.unravel_index(np.argmax(image), image.shape) + self.assertAlmostEqual(float(image[row, column]), 1., places=5) + np.testing.assert_allclose(payloads["sar3d_peak"]["table"][1:3], [x[column], y[row]], rtol=0, atol=1e-6) + 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) @@ -133,6 +207,7 @@ class Sar3dTest(unittest.TestCase): payloads = run(request) self.assertEqual(status(payloads)[0], 1) self.assertNotIn("sar3d_peak", payloads) + self.assertNotIn("sar3d_background", payloads) for name in VIEWS: self.assertFalse(np.any(payloads[name]["image"])) diff --git a/device_firmware/librevna_protocol.md b/device_firmware/librevna_protocol.md index 25c4da8..a3efb54 100644 --- a/device_firmware/librevna_protocol.md +++ b/device_firmware/librevna_protocol.md @@ -287,16 +287,30 @@ GUI публикует настройки через `processing_live.json`; ф действуют в обоих режимах: среднее предыдущих коллекций вычитается до 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 — автоматически), +текущий; шаги вокселя по X, Y и Z — `sar3d_voxel_x_m`, `sar3d_voxel_y_m`, `sar3d_voxel_z_m`. Размеры сцены: `sar3d_x_size_m`, `sar3d_z_size_m` (0 — автоматически), `sar3d_y_min_m`, `sar3d_y_max_m` — наземные координаты, независимо от depth gate. +`sar3d_bg_gate_enabled` (BG_GATE_ENABLE) и `sar3d_clean_enabled` (SAR_CLEAN_ENABLE), +boolean, по умолчанию true: отношение максимума кадра к его фону и CLEAN апертуры +из одного кадра. -SAR публикует `sar3d_target_width` как TableF32 с тремя колонками [X,Y,Z]: -ширины непрерывного отклика на половине максимума по осевым сечениям через -совместный максимум сглаженного объёма (метры, не меньше одного вокселя). -`sar3d_peak` по-прежнему содержит [уровень,X,Y,Z]. GUI использует ширины для -квадратного контура цели. Зум и `SAR_CLOUD_THRESHOLD` — настройки отображения -desktop/веба; в live-файл не входят, ревизию и расчёт не меняют. RSL3 не меняется: -новая таблица использует существующий тип payload. +SAR публикует `sar3d_peak` как TableF32 [уровень нормировки,X,Y,Z] отмеченного +максимума и, если объектов CLEAN нет, `sar3d_target_width` с тремя колонками [X,Y,Z]: +ширины непрерывного отклика на половине максимума по осевым сечениям через этот +максимум сглаженного объёма (метры, не меньше одного вокселя). GUI использует ширины +для квадратного контура цели. Без сигнала (статус 1) этих таблиц нет. Остальные +таблицы TableF32: + +- `sar3d_background` — [максимум/фон, BG_GATE_3D_SHOW, BG_GATE_3D_FULL], только при + `sar3d_bg_gate_enabled` и сигнале; множитель шкалы считает GUI с порогом облака. +- `sar3d_objects` — объекты CLEAN, сильнейший первым, строки [X,Y,Z,уровень в фонах]. + Если они есть, проекции — остаток в долях максимума кадра с каждым объектом в своей + шкале 0–1, а `sar3d_peak` указывает на первый объект. +- `sar3d_marks` — рамки групп объектов, строки [класс,X0,X1,Y0,Y1,Z0,Z1]; класс 0 — + очень низкий, 1 — на земле, 2 — выше земли, 3 — у границы сцены. + +Зум, `SAR_CLOUD_THRESHOLD` и `SAR_MARK_ENABLE` — настройки отображения desktop/веба; +в live-файл не входят, ревизию и расчёт не меняют. RSL3 не меняется: таблицы +используют существующий тип payload. ## 7. Почему поле мощности читается, но не используется diff --git a/python_app/gui/controllers/renderers/sar3d.py b/python_app/gui/controllers/renderers/sar3d.py index 2d93508..69a294a 100644 --- a/python_app/gui/controllers/renderers/sar3d.py +++ b/python_app/gui/controllers/renderers/sar3d.py @@ -9,23 +9,33 @@ 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_peak``: [level, X, Y, Z] of the marked maximum; level is the normalization. +- ``sar3d_target_width``: axial FWHM [X, Y, Z] in metres through that maximum, + without CLEAN objects. +- ``sar3d_background``: [peak over background, BG_GATE_3D_SHOW, BG_GATE_3D_FULL] + while BG_GATE_ENABLE is on; the frame is scaled by ``bg_gate_gain``. +- ``sar3d_objects``: CLEAN objects, strongest first, rows [X, Y, Z, level in + backgrounds]. The projections then show the residual with every object in its + own 0-1 scale, and the marked maximum is the first object. +- ``sar3d_marks``: boxes around groups of the objects, rows [class, X0, X1, Y0, + Y1, Z0, Z1], classes as in ``SAR3D_MARKS``. - ``sar3d_status``: [code, UTF-8 bytes...] with ready CUDA (0), no signal (1), processing error (2), or CUDA unavailable (3), including the explanation. """ from __future__ import annotations +from collections.abc import Callable import itertools import numpy as np +from numpy.typing import ArrayLike 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 ( - ANTENNA_STYLES, build_lut, create_plot, image_rect, selected_antenna_coordinates, + ANTENNA_STYLES, OVERLAY_FILL, build_lut, create_plot, image_rect, selected_antenna_coordinates, ) from python_app.gui.result_payloads import collection_payload_by_name from python_app.gui.theme import PLOT_FOREGROUND @@ -36,6 +46,19 @@ 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"} +# Box classes of "sar3d_marks" (notebook _SAR_MARK_STYLE): name, color, dashed. A box at the +# boundary is not yet decidable: part of its response may lie outside the scene. +SAR3D_MARKS: dict[int, tuple[str, tuple[int, int, int], bool]] = { + 0: ("very low", (229, 38, 43), False), + 1: ("on the ground", (255, 154, 31), False), + 2: ("above the ground", (47, 111, 222), False), + 3: ("at the boundary", (138, 138, 138), True), +} + +# Under the views and inside the captured grid, so the web view shows it too. +_BOX_LEGEND = "Boxes: " + "   ".join( + f"■ {name}" for name, color, _dashed in SAR3D_MARKS.values()) + _GROUND_COLOR = (115, 105, 88) _SCENE_BOX_COLOR = "#898781" @@ -47,7 +70,7 @@ SAR3D_VIEWS: dict[str, tuple[str, str, str, tuple[int, int] | None]] = { "sar3d_oblique": ("3D Projection", "", "", None), } -# Corners of the scene box indexed by bits (x, y, z); an edge joins corners one bit apart. +# Corners of a box indexed by bits (x, y, z); an edge joins corners one bit apart. _BOX_EDGES = [(a, b) for a, b in itertools.combinations(range(8), 2) if bin(a ^ b).count("1") == 1] # Axis names go on the edges meeting at the bottom corner nearest the camera (X max, Y min, Z min). _AXIS_EDGES = {"X": (0, 1), "Y": (1, 3), "Z": (1, 5)} @@ -67,6 +90,45 @@ def sar3d_status(collection: ResultCollection) -> tuple[int, str] | None: return code, bytes(row[1:].astype(np.uint8)).decode("utf-8", errors="replace") +def bg_gate_gain(ratio: float, show: float, full: float, floor: float) -> float: + """Scale of a frame by its peak over background (notebook bg_gate_gain). + + Below ``show`` the frame is empty; at ``show`` its maximum lands on ``floor``, the + cloud threshold, so an object appears as a point; from ``full`` the scale is full. + """ + if ratio < show: + return 0.0 + if full <= show: + return 1.0 + return float(np.clip(floor + (1.0 - floor) * (ratio - show) / (full - show), floor, 1.0)) + + +Projection = Callable[[ArrayLike], np.ndarray] + + +def _projection(axes: tuple[int, int] | None, basis: np.ndarray) -> Projection: + """Scene points -> plot coordinates of a view: two scene axes, or the oblique image axes.""" + if axes is None: + return lambda points: np.asarray(points, dtype=float) @ basis.T + return lambda points: np.asarray(points, dtype=float)[..., list(axes)] + + +def _box_corners(bounds: np.ndarray) -> np.ndarray: + """Eight corners of the box [X0, X1, Y0, Y1, Z0, Z1], indexed by bits (x, y, z).""" + return np.array([[bounds[2*axis + ((index >> axis) & 1)] for axis in range(3)] for index in range(8)]) + + +def _table(collection: ResultCollection, name: str, columns: int) -> np.ndarray: + """Rows of a table payload; none when it is absent.""" + payload = collection_payload_by_name(collection, name, kind=4) + if payload is None: + return np.empty((0, columns)) + table = np.asarray(payload.table, dtype=float) + if table.ndim != 2 or table.shape[1] != columns: + raise ValueError(f"{name} must have {columns} columns") + return table + + class Sar3dRenderer: """Draws the 3D SAR views.""" @@ -91,6 +153,9 @@ class Sar3dRenderer: self.antennas: dict[str, tuple[pg.ScatterPlotItem, pg.ScatterPlotItem]] = {} self.images: dict[str, pg.ImageItem] = {} self.targets: dict[str, pg.PlotCurveItem] = {} + self.marks: dict[str, dict[int, pg.PlotCurveItem]] = {} + self.objects: dict[str, pg.ScatterPlotItem] = {} + self.object_labels: dict[str, list[pg.TextItem]] = {} for index, (name, (title, x_label, y_label, axes)) in enumerate(SAR3D_VIEWS.items()): plot = create_plot(title) grid_layout.addWidget(plot, index // 2, index % 2) @@ -116,6 +181,21 @@ class Sar3dRenderer: plot.addItem(peak) self.images[name] = image self.targets[name] = peak + self.marks[name] = {} + for kind, (_name, color, dashed) in SAR3D_MARKS.items(): + style = Qt.PenStyle.DashLine if dashed else Qt.PenStyle.SolidLine + # A dark outline keeps the box visible over pixels of its own colour. + curve = pg.PlotCurveItem(pen=pg.mkPen(color, width=2, style=style), connect="pairs", + shadowPen=pg.mkPen(PLOT_FOREGROUND, width=4, style=style)) + curve.setZValue(3) + plot.addItem(curve) + self.marks[name][kind] = curve + objects = pg.ScatterPlotItem(symbol="d", size=11, brush=pg.mkBrush("#ffffff"), + pen=pg.mkPen("#111111", width=1.5)) + objects.setZValue(5) + plot.addItem(objects) + self.objects[name] = objects + self.object_labels[name] = [] if x_label: plot.setLabel("bottom", x_label, units="m") plot.setLabel("left", y_label, units="m") @@ -142,6 +222,11 @@ class Sar3dRenderer: label = pg.TextItem(axis_name, color=PLOT_FOREGROUND, anchor=(0.5, 0.5)) oblique.addItem(label) self.axis_labels[axis_name] = label + self.box_legend = QLabel(_BOX_LEGEND) + self.box_legend.setTextFormat(Qt.TextFormat.RichText) + self.box_legend.setWordWrap(True) + self.box_legend.hide() + grid_layout.addWidget(self.box_legend, 2, 0, 1, 2) self.clear() def clear(self, message: str = "SAR3D: waiting for measurements") -> None: @@ -153,12 +238,45 @@ class Sar3dRenderer: for markers in self.antennas.values(): for marker in markers: marker.setData([], []) + self._draw_marks(np.empty((0, 7)), {}) + self._draw_objects(np.empty((0, 4)), {}) self.box.setData([], []) self.ground.setData([], []) for label in self.axis_labels.values(): label.hide() self.status_label.setText(message) + def _draw_marks(self, boxes: np.ndarray, projections: dict[str, Projection]) -> None: + """Boxes around the CLEAN object groups, one curve per class and view.""" + for view, curves in self.marks.items(): + for kind, curve in curves.items(): + rows = boxes[boxes[:, 0] == kind, 1:] + if not len(rows) or view not in projections: + curve.setData([], []) + continue + points = projections[view](np.concatenate([_box_corners(row)[np.ravel(_BOX_EDGES)] for row in rows])) + curve.setData(points[:, 0], points[:, 1]) + + def _draw_objects(self, objects: np.ndarray, projections: dict[str, Projection]) -> None: + """Numbered CLEAN objects with their level in frame backgrounds.""" + for name, markers in self.objects.items(): + plot = self.plots[name] + for label in self.object_labels[name]: + plot.removeItem(label) + self.object_labels[name].clear() + if not len(objects) or name not in projections: + markers.setData([], []) + continue + points = projections[name](objects[:, :3]) + markers.setData(points[:, 0], points[:, 1]) + for number, ((x, y), level) in enumerate(zip(points, objects[:, 3]), start=1): + label = pg.TextItem(f"{number}: {level:.0f} bg", color=PLOT_FOREGROUND, + fill=pg.mkBrush(*OVERLAY_FILL), anchor=(0.5, 1.2)) + label.setZValue(6) + label.setPos(float(x), float(y)) + plot.addItem(label) + self.object_labels[name].append(label) + def draw(self, collection: ResultCollection, settings: GuiSar3dStateModel, geometry: GprModel) -> bool: """Draw four maximum projections of the volume; return success.""" status = sar3d_status(collection) @@ -178,6 +296,22 @@ class Sar3dRenderer: return False peak = collection_payload_by_name(collection, "sar3d_peak", kind=4) level, *peak_xyz = peak.table[0] if peak is not None else (None,) + objects = _table(collection, "sar3d_objects", 4) + marks = _table(collection, "sar3d_marks", 7) + # Frame relative to its background; CLEAN objects are each in their own scale already. + background = _table(collection, "sar3d_background", 3) + gain = 1.0 + background_note = "" + if len(background): + ratio, appear, full = (float(value) for value in background[0]) + background_note = f"peak {ratio:.1f} backgrounds" + if not len(objects): + gain = bg_gate_gain(ratio, appear, full, settings.cloud_threshold) + state = ("below the appearance level, frame empty" if gain <= 0.0 + else "full scale" if gain >= 1.0 else f"scale ×{gain:.2f}") + background_note += f" (appears from {appear:g}, full scale from {full:g}): {state}" + # Objects replace the peak marker, and so does a frame below the appearance level. + show_peak = not len(objects) and gain > 0.0 antennas = selected_antenna_coordinates(geometry, settings.input_positions, settings.output_positions) pitch = np.radians(settings.look_angle_deg) @@ -186,48 +320,51 @@ class Sar3dRenderer: 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) + steps = np.array([settings.voxel_x_m, settings.voxel_y_m, settings.voxel_z_m]) + widths = width_payload.table[0] if width_payload is not None else steps zoom = settings.zoom_percent / 100.0 + projections = {name: _projection(axes, basis.table) for name, (*_labels, axes) in SAR3D_VIEWS.items()} 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(np.where(payload.image >= settings.cloud_threshold, payload.image, 0), - autoLevels=False, levels=(0, 1)) + shown = payload.image * gain + image.setImage(np.where(shown >= settings.cloud_threshold, shown, 0), autoLevels=False, levels=(0, 1)) image.setRect(rect) image.show() plot = self.plots[name] - position = (basis.table @ np.asarray(peak_xyz) if axes is None else np.asarray(peak_xyz)[list(axes)]) if peak_xyz else None + position = projections[name](peak_xyz) 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 + pad = steps.max() 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: + if position is not None and show_peak: 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.targets[name].setData([], []) + if position is not None and 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) plot.setRange(rect, padding=0.04) plot.setTitle(title) + boxes = marks if settings.mark_enabled else np.empty((0, 7)) + self._draw_marks(boxes, projections) + # The legend follows the settings rather than each frame, so the views keep their size. + self.box_legend.setVisible(settings.mark_enabled and settings.clean_enabled) + self._draw_objects(np.empty((0, 4)) if settings.mark_enabled else objects, projections) # Scene box and ground on the oblique view, from the extents of the other views. front, top = payloads["sar3d_front"], payloads["sar3d_top"] - bounds = [ - (front.image_x_axis[0], front.image_x_axis[-1]), - (front.image_y_axis[0], front.image_y_axis[-1]), - (top.image_y_axis[0], top.image_y_axis[-1]), - ] - corners = np.array([[bounds[axis][(index >> axis) & 1] for axis in range(3)] for index in range(8)]) - projected = corners @ basis.table.T + bounds = np.array([front.image_x_axis[0], front.image_x_axis[-1], front.image_y_axis[0], + front.image_y_axis[-1], top.image_y_axis[0], top.image_y_axis[-1]]) + projected = _box_corners(bounds) @ basis.table.T box = projected[np.ravel(_BOX_EDGES)] self.box.setData(box[:, 0], box[:, 1]) centre = projected.mean(axis=0) @@ -235,12 +372,20 @@ class Sar3dRenderer: middle = (projected[start] + projected[end]) / 2 self.axis_labels[axis_name].setPos(*(middle + 0.08 * (middle - centre))) # just outside self.axis_labels[axis_name].show() - (x0, x1), _, (z0, z1) = bounds + x0, x1, _y0, _y1, z0, z1 = bounds ground = (np.array([[x0, 0.0, z0], [x1, 0.0, z0], [x1, 0.0, z1], [x0, 0.0, z1], [x0, 0.0, z0]]) @ basis.table.T).T self.ground.setData(ground[0], ground[1]) - summary = f"SAR3D: ready — {status[1]}" if status else "SAR3D: ready" + + notes = [f"SAR3D: ready — {status[1]}" if status else "SAR3D: ready"] + if background_note: + notes.append(background_note) + if len(boxes): + counts = {kind: int(np.sum(boxes[:, 0] == kind)) for kind in SAR3D_MARKS} + notes.append("boxes: " + ", ".join(f"{counts[kind]} {name}" for kind, (name, *_style) + in SAR3D_MARKS.items() if kind != 3 or counts[kind])) 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(f"{summary} | zoom {settings.zoom_percent}%") + notes.append(f"max X {x:+.2f}, Y {y:+.2f}, Z {z:.2f} m | level {level:.3g}") + notes.append(f"zoom {settings.zoom_percent}%") + self.status_label.setText(" | ".join(notes)) return True diff --git a/python_app/gui/controllers/sections/processing_section.py b/python_app/gui/controllers/sections/processing_section.py index 8c4342a..c76e2fb 100644 --- a/python_app/gui/controllers/sections/processing_section.py +++ b/python_app/gui/controllers/sections/processing_section.py @@ -214,7 +214,7 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces ) processing_mode_pages.addWidget(gpr_page) # 3D SAR page. The labels are the names of the reference notebook - # Multi_antenna_demonstrate_depthnorm.ipynb; the shared gpr_* fields get their own + # Multi_depthnorm_clean.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") @@ -226,11 +226,12 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces sar3d_min_depth_m.setToolTip("Radar depth gate start (m); the scene grid along Z follows the gate.") sar3d_max_depth_m = view.add("_sar3d_max_depth_m", double_spin(0.1, 50.0, 2, 0.1)) sar3d_max_depth_m.setToolTip("Radar depth gate end (m); the scene grid along Z follows the gate.") - sar3d_voxel_m = view.add("_sar3d_voxel_m", double_spin(0.01, 0.1, 3, 0.005)) - sar3d_voxel_m.setToolTip( - "Edge of the cubic voxel (m), the same along X, Y and Z. Smaller is finer and slower: " - "halving it makes 8 times more voxels." - ) + sar3d_voxel_x_m = view.add("_sar3d_voxel_x_m", double_spin(0.01, 0.1, 3, 0.005)) + sar3d_voxel_y_m = view.add("_sar3d_voxel_y_m", double_spin(0.01, 0.1, 3, 0.005)) + sar3d_voxel_z_m = view.add("_sar3d_voxel_z_m", double_spin(0.01, 0.1, 3, 0.005)) + for axis, step in zip("XYZ", (sar3d_voxel_x_m, sar3d_voxel_y_m, sar3d_voxel_z_m)): + step.setToolTip(f"Voxel step along {axis} (m). Smaller is finer and slower: halving all three " + "steps makes 8 times more voxels.") sar3d_look_angle_deg = view.add("_sar3d_look_angle_deg", double_spin(-90.0, 90.0, 2, 0.5)) sar3d_look_angle_deg.setToolTip( "Downward tilt of the radar (deg): motion correction along range and the ground scene axes." @@ -253,6 +254,17 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces 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_bg_gate_enabled = view.add("_sar3d_bg_gate_enabled", QCheckBox("BG_GATE_ENABLE")) + sar3d_bg_gate_enabled.setToolTip( + "Show the frame relative to its background, the median of the volume: below 6 backgrounds the " + "frame is empty, at 6 its maximum appears at SAR_CLOUD_THRESHOLD, from 8.5 the scale is full." + ) + sar3d_clean_enabled = view.add("_sar3d_clean_enabled", QCheckBox("SAR_CLEAN_ENABLE")) + sar3d_clean_enabled.setToolTip( + "CLEAN: subtract the strongest responses one by one, so weak objects show beside strong ones; " + "each object is drawn in its own 0-1 scale. Runs on a one-frame aperture: always with " + "SAR_FRAME_RADIUS = 0, otherwise while no previous frame is available, as after Start." + ) 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)) @@ -261,6 +273,11 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces 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_mark_enabled = view.add("_sar3d_mark_enabled", QCheckBox("SAR_MARK_ENABLE")) + sar3d_mark_enabled.setToolTip( + "Box the CLEAN objects by height: red very low (like a mine), orange on the ground, blue above it, " + "dashed grey at the scene boundary (not yet decidable). Off: numbered objects instead." + ) sar3d_page = _build_processing_mode_page( processing_mode_pages, @@ -271,7 +288,9 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces ("F_STOP MHz", sar3d_stop_freq_mhz), ("MIN_DEPTH", sar3d_min_depth_m), ("MAX_DEPTH", sar3d_max_depth_m), - ("SAR_VOXEL_M", sar3d_voxel_m), + ("SAR_VOXEL_M X", sar3d_voxel_x_m), + ("SAR_VOXEL_M Y", sar3d_voxel_y_m), + ("SAR_VOXEL_M Z", sar3d_voxel_z_m), ("LOOK_ANGLE_DEG", sar3d_look_angle_deg), sar3d_ignore_socket_speed_enabled, ("SPEED_M_S", sar3d_speed_m_s), @@ -283,14 +302,16 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces sar3d_pair_depth_equalize, sar3d_background_subtract_enabled, ("Mean count", sar3d_background_mean_count), + sar3d_bg_gate_enabled, + sar3d_clean_enabled, ("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), + sar3d_mark_enabled, ], - split_index=8, ) processing_mode_pages.addWidget(sar3d_page) @@ -320,7 +341,7 @@ def build_processing_group(config: RunConfigModel, gui: GuiStateModel) -> Proces 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, sar3d_cloud_threshold, sar3d_zoom_percent), + visual_controls=(gpr_render_mode, sar3d_cloud_threshold, sar3d_zoom_percent, sar3d_mark_enabled), locator_controls=(gpr_max_detected_objects_to_draw, gpr_draw_top_m_objects, gpr_visible_x_min_m, gpr_visible_x_max_m, gpr_visible_z_min_m, gpr_visible_z_max_m), diff --git a/python_app/models/gui_profile_schema.py b/python_app/models/gui_profile_schema.py index 3c4ebb2..70448ee 100644 --- a/python_app/models/gui_profile_schema.py +++ b/python_app/models/gui_profile_schema.py @@ -134,7 +134,7 @@ class GuiGprStateModel: class GuiSar3dStateModel: """UI-only defaults for the 3D SAR (sar3d) live settings. - Values are the reference notebook's (Multi_antenna_demonstrate_depthnorm.ipynb), whose + Values are the reference notebook's (Multi_depthnorm_clean.ipynb), whose names the GUI shows: INPUT_IDX, F_START, MIN_DEPTH, LOOK_ANGLE_DEG, ... """ @@ -148,7 +148,9 @@ class GuiSar3dStateModel: ignore_socket_speed_enabled: bool = setting(False, '_sar3d_ignore_socket_speed_enabled', live='ignore_socket_speed', group='Motion') look_angle_deg: float = setting(19.5, '_sar3d_look_angle_deg', live='gpr_look_angle_deg', group='Motion') frame_radius: int = setting(0, '_sar3d_frame_radius', live='sar3d_frame_radius', group='3D SAR', minimum=0) - voxel_m: float = setting(0.02, '_sar3d_voxel_m', live='sar3d_voxel_m', group='3D SAR', minimum=0, exclusive_minimum=True) + voxel_x_m: float = setting(0.025, '_sar3d_voxel_x_m', live='sar3d_voxel_x_m', group='3D SAR', minimum=0, exclusive_minimum=True) + voxel_y_m: float = setting(0.02, '_sar3d_voxel_y_m', live='sar3d_voxel_y_m', group='3D SAR', minimum=0, exclusive_minimum=True) + voxel_z_m: float = setting(0.02, '_sar3d_voxel_z_m', live='sar3d_voxel_z_m', group='3D SAR', minimum=0, exclusive_minimum=True) bp_ground_depth_m: float = setting(1.25, '_sar3d_bp_ground_depth_m', live='sar3d_bp_ground_depth_m', group='3D SAR') 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') @@ -156,12 +158,15 @@ class GuiSar3dStateModel: 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) + bg_gate_enabled: bool = setting(True, '_sar3d_bg_gate_enabled', live='sar3d_bg_gate_enabled', group='3D SAR') + clean_enabled: bool = setting(True, '_sar3d_clean_enabled', live='sar3d_clean_enabled', group='3D SAR') 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) + mark_enabled: bool = setting(True, '_sar3d_mark_enabled', web='sar3d_mark_enabled', group='Display') @dataclass(slots=True) diff --git a/python_app/orchestration/live_processing_config.py b/python_app/orchestration/live_processing_config.py index f7c8ef0..fb6a856 100644 --- a/python_app/orchestration/live_processing_config.py +++ b/python_app/orchestration/live_processing_config.py @@ -53,11 +53,13 @@ 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_demonstrate_depthnorm.ipynb (SAR_FRAME_RADIUS, BP_GROUND_DEPTH_M, ...). + # Multi_depthnorm_clean.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 - sar3d_voxel_m: float = 0.02 # SAR_VOXEL_M, one edge for X, Y and Z + sar3d_voxel_x_m: float = 0.025 # SAR_VOXEL_M, the voxel steps along X, Y and Z + sar3d_voxel_y_m: float = 0.02 + sar3d_voxel_z_m: float = 0.02 sar3d_bp_ground_depth_m: float = 1.25 sar3d_comp_range_power: float = 0.3 sar3d_use_compensation: bool = True @@ -66,6 +68,8 @@ class ProcessingLiveConfig: sar3d_y_min_m: float = -0.5 sar3d_y_max_m: float = 2.0 sar3d_z_size_m: float = 0.0 + sar3d_bg_gate_enabled: bool = True # BG_GATE_ENABLE: peak-to-background ratio for the display + sar3d_clean_enabled: bool = True # SAR_CLEAN_ENABLE: CLEAN of a one-frame aperture into objects # 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 diff --git a/python_app/tests/test_gui_profile_codec.py b/python_app/tests/test_gui_profile_codec.py index bfa7eb9..1ca5fa9 100644 --- a/python_app/tests/test_gui_profile_codec.py +++ b/python_app/tests/test_gui_profile_codec.py @@ -82,7 +82,9 @@ class GuiProfileCodecTest(unittest.TestCase): ignore_socket_speed_enabled=True, look_angle_deg=6.5, frame_radius=3, - voxel_m=0.035, + voxel_x_m=0.035, + voxel_y_m=0.03, + voxel_z_m=0.025, bp_ground_depth_m=1.1, comp_range_power=0.5, use_compensation=False, diff --git a/python_app/tests/test_processor_history.py b/python_app/tests/test_processor_history.py index 17fe4a0..cc480b3 100644 --- a/python_app/tests/test_processor_history.py +++ b/python_app/tests/test_processor_history.py @@ -382,7 +382,7 @@ class Sar3dRecoveryTest(_ProcessorTestCase): self._live.gpr_speed_m_s = 0.0 self._live.gpr_start_freq_mhz = 1500.0 self._live.gpr_stop_freq_mhz = 6000.0 - self._live.sar3d_voxel_m = 0.4 + self._live.sar3d_voxel_x_m = self._live.sar3d_voxel_y_m = self._live.sar3d_voxel_z_m = 0.4 self._live.sar3d_frame_radius = 0 self._live.tx_geometry = [{"output_pos": 0, "x_m": -0.3, "y_m": 0.0, "z_m": 0.0}] self._live.rx_geometry = [{"input_pos": 0, "x_m": 0.2, "y_m": 0.4, "z_m": 0.0}] diff --git a/python_app/tests/test_sar3d_plot.py b/python_app/tests/test_sar3d_plot.py index cc868eb..71e5bcd 100644 --- a/python_app/tests/test_sar3d_plot.py +++ b/python_app/tests/test_sar3d_plot.py @@ -2,17 +2,19 @@ from __future__ import annotations +import json import os import unittest from unittest import mock import numpy as np +import pyqtgraph as pg os.environ.setdefault("QT_QPA_PLATFORM", "offscreen") 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.gui.controllers.renderers.sar3d import SAR3D_MARKS, 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 @@ -29,6 +31,17 @@ def _collection(*payloads: ResultPayload) -> ResultCollection: capture_start_ns=1200, capture_end_ns=1400) +def _points(item: pg.PlotCurveItem | pg.ScatterPlotItem) -> tuple[np.ndarray, np.ndarray]: + """Drawn x and y of a curve or scatter item, empty when nothing is drawn.""" + xx, yy = item.getData() + assert xx is not None and yy is not None + return np.asarray(xx), np.asarray(yy) + + +def _table(name: str, *rows: list[float]) -> ResultPayload: + return ResultPayload(processing_name=name, kind=4, table=np.array(rows, dtype=np.float32)) + + def _projections() -> ResultCollection: x, y, z = np.array([-1, 0, 1]), np.array([-0.5, 0.5]), np.array([1, 2, 3]) axes = [(x, y), (z, y), (x, z), (x, y)] @@ -129,6 +142,68 @@ class Sar3dPlotTest(unittest.TestCase): 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_bg_gate_scales_the_frame_and_hides_it_below_appearance(self) -> None: + self.settings.cloud_threshold = .5 + for ratio, gain, state in ((9., 1., "full scale"), (7.5, .8, "scale ×0.80"), + (5., 0., "below the appearance level, frame empty")): + with self.subTest(ratio=ratio): + result = _projections() + result.collection_payloads.append(_table("sar3d_background", [ratio, 6, 8.5])) + self.assertTrue(self.renderer.draw(result, self.settings, self.geometry)) + for name, image in self.renderer.images.items(): + np.testing.assert_allclose(np.asarray(image.image), gain, err_msg=name) + self.assertEqual(len(_points(self.renderer.targets[name])[0]) > 0, gain > 0, name) + self.assertIn(f"peak {ratio:.1f} backgrounds (appears from 6, full scale from 8.5): {state}", + self.renderer.status_label.text()) + + def test_clean_objects_are_boxed_by_class_or_numbered(self) -> None: + result = _projections() + result.collection_payloads += [ + _table("sar3d_background", [7.5, 6, 8.5]), + _table("sar3d_objects", [.5, -.2, 2., 45.], [-.5, .3, 1.5, 12.]), + _table("sar3d_marks", [0, .25, .75, -.3, .1, 1.75, 2.25], [2, -.75, -.25, .1, .5, 1.25, 1.75]), + ] + self.assertTrue(self.renderer.draw(result, self.settings, self.geometry)) + text = self.renderer.status_label.text() + self.assertIn("peak 7.5 backgrounds | boxes: 1 very low, 0 on the ground, 1 above the ground |", text) + self.assertNotIn("at the boundary", text) + for name, (_title, _x, _y, axes) in SAR3D_VIEWS.items(): + # Objects are in their own scale: no background gain and no peak marker. + np.testing.assert_allclose(np.asarray(self.renderer.images[name].image), 1) + self.assertFalse(len(_points(self.renderer.targets[name])[0])) + self.assertFalse(len(_points(self.renderer.objects[name])[0])) + curves = self.renderer.marks[name] + self.assertFalse(len(_points(curves[1])[0]) or len(_points(curves[3])[0])) + plane = axes or (0, 1) # the test basis shows the oblique view along X and Y + for kind, low, high in ((0, (.25, -.3, 1.75), (.75, .1, 2.25)), (2, (-.75, .1, 1.25), (-.25, .5, 1.75))): + xx, yy = _points(curves[kind]) + self.assertEqual(len(xx), 24) + np.testing.assert_allclose([xx.min(), yy.min()], np.asarray(low)[list(plane)]) + np.testing.assert_allclose([xx.max(), yy.max()], np.asarray(high)[list(plane)]) + # The legend lists every class under the views while boxes are on. + self.assertFalse(self.renderer.box_legend.isHidden()) + self.assertTrue(all(name in self.renderer.box_legend.text() for name, *_style in SAR3D_MARKS.values())) + + self.settings.mark_enabled = False + self.renderer.draw(result, self.settings, self.geometry) + self.assertTrue(self.renderer.box_legend.isHidden()) + self.assertNotIn("boxes:", self.renderer.status_label.text()) + for name, (_title, _x, _y, axes) in SAR3D_VIEWS.items(): + self.assertTrue(all(not len(_points(curve)[0]) for curve in self.renderer.marks[name].values())) + xx, yy = _points(self.renderer.objects[name]) + plane = list(axes or (0, 1)) + np.testing.assert_allclose(np.column_stack([xx, yy]), np.array([[.5, -.2, 2.], [-.5, .3, 1.5]])[:, plane]) + self.assertEqual([label.textItem.toPlainText() for label in self.renderer.object_labels[name]], + ["1: 45 bg", "2: 12 bg"]) + self.renderer.clear() + self.assertFalse(any(self.renderer.object_labels.values())) + + def test_table_with_another_width_is_an_error(self) -> None: + result = _projections() + result.collection_payloads.append(_table("sar3d_objects", [.5, -.2, 2.])) + with self.assertRaisesRegex(ValueError, "sar3d_objects must have 4 columns"): + self.renderer.draw(result, self.settings, self.geometry) + def test_cloud_threshold_masks_weak_response_without_changing_levels(self) -> None: result = _projections() for payload in result.collection_payloads: @@ -167,13 +242,35 @@ class Sar3dIntegrationTest(unittest.TestCase): 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()) + widget(self.view._view.processing, "_sar3d_mark_enabled", QCheckBox).setChecked(False) + self.assertEqual(self.view._live.revision, revision) live = self.view._live.build().to_dict() self.assertNotIn("sar3d_zoom_percent", live) self.assertNotIn("sar3d_cloud_threshold", live) + self.assertNotIn("sar3d_mark_enabled", 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.) + self.assertIs(schema["sar3d_mark_enabled"]["value"], False) + + def test_bg_gate_and_clean_publish_a_processing_revision(self) -> None: + self.view._pipeline.drain_results = mock.Mock() + live = self.view._live.build().to_dict() + self.assertIs(live["sar3d_bg_gate_enabled"], True) + self.assertIs(live["sar3d_clean_enabled"], True) + for name, key in (("_sar3d_bg_gate_enabled", "sar3d_bg_gate_enabled"), + ("_sar3d_clean_enabled", "sar3d_clean_enabled")): + with self.subTest(key): + revision = self.view._live.revision + widget(self.view._view.processing, name, QCheckBox).setChecked(False) + self.assertNotEqual(self.view._live.revision, revision) + published = json.loads(self.view._live.writer.path.read_text(encoding="utf-8")) + self.assertIs(published[key], False) + self.assertEqual(published["settings_revision"], self.view._live.revision) + schema = {entry["name"]: entry for entry in self.view._config.web_schema()} + self.assertEqual(schema["sar3d_bg_gate_enabled"]["label"], "BG_GATE_ENABLE") + self.assertEqual(schema["sar3d_clean_enabled"]["label"], "SAR_CLEAN_ENABLE") def test_inactive_sar_scene_is_validated_before_publishing_gpr_settings(self) -> None: with QSignalBlocker(self.mode): 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 37e6c45..82d2d04 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 @@ -336,7 +336,6 @@ "ignore_socket_speed_enabled": false, "look_angle_deg": 10.0, "frame_radius": 0, - "voxel_m": 0.02, "bp_ground_depth_m": 1.25, "comp_range_power": 0.3, "use_compensation": true,