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

Port the SAR3D background gate, CLEAN and object boxes from Multi_depthnorm_clean

Browse Source 
- BG_GATE_ENABLE: frame peak over its median background (sar3d_background); the
  GUI scales the frame with bg_gate_gain and its cloud threshold.
- SAR_CLEAN_ENABLE: CLEAN of a one-frame aperture on a 4x10x4 cm search grid over
  the whole scene; objects in their own scale (sar3d_objects).
- SAR_MARK_ENABLE: boxes by class around object groups (sar3d_marks), with the
  legend under the views.
- SAR_VOXEL_M as three steps X/Y/Z with the notebook defaults 0.025/0.02/0.02.
- GPU backend: image/download/load/project API, per-plane block reduction for the
  maximum, depth gate without FMA; sar3d host code built with -ffp-contract=off.
- References regenerated from notebook cells 1 and 5, with new CLEAN scenes and a
  synthetic box-rule case; docs updated.
This commit is contained in:
Ayzen
2026-10-16 12:08:57 +03:00
parent 9a9b8ce5c7
commit b88d704126
26 changed files with 1895 additions and 407 deletions
Download Patch File Download Diff File Expand all files Collapse all files
data_acq_and_processing/processing/processors/include/sar3d_processor.hpp
+7 -5
View File
@@ -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();
data_acq_and_processing/processing/processors/src/backprojection.hpp
+22 -5
View File
@@ -3,8 +3,10 @@
#include "shared_types.hpp"
#include <algorithm>
#include <array>
#include <cmath>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <span>
#include <map>
@@ -12,6 +14,7 @@
#include <numbers>
#include <stdexcept>
#include <unordered_map>
#include <utility>
#include <vector>
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::ComboKey, std::vector<std::complex<double>>, ipc::ComboKeyHash>;
// numpy.median: the middle value, or the mean of the middle two for an even count.
template <class Value>
[[nodiscard]] auto median(std::vector<Value> values) -> double {
if (values.empty()) { throw std::logic_error("median of no values"); }
const auto middle = values.begin() + static_cast<std::ptrdiff_t>(values.size() / 2U);
std::nth_element(values.begin(), middle, values.end());
if (values.size() % 2U == 1U) { return static_cast<double>(*middle); }
return 0.5 * (static_cast<double>(*middle) + static_cast<double>(*std::max_element(values.begin(), middle)));
}
// Mean of up to count previous collections, separately for each physical pair.
inline auto background_mean(std::span<const ipc::PreprocessedCollection> 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<double, 2> {
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 <class Real>
@@ -220,10 +240,7 @@ void normalize_pair(std::span<std::complex<Real>> samples, double dz, double min
for (const auto sample : samples) { magnitudes.push_back(std::abs(std::complex<double>(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<Real>(scale + 1e-15); }
}
data_acq_and_processing/processing/processors/src/sar3d_clean.cpp
+558
View File
@@ -0,0 +1,558 @@
#include "sar3d_clean.hpp"
#include "backprojection.hpp"
#include <algorithm>
#include <cmath>
#include <functional>
#include <iterator>
#include <limits>
#include <numeric>
#include <stdexcept>
#include <utility>
// 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<double, 3> kCleanMerge{0.35, 0.50, 0.40}; // _SAR_CLEAN_MERGE_M (X, Y, Z)
constexpr std::array<double, 3> 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<double, 2> 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<std::complex<float>>;
// numpy.interp: linear between the points of increasing `xp`, constant beyond them.
[[nodiscard]] auto interp(double x, std::span<const double> xp, std::span<const double> fp) -> double {
if (x <= xp.front()) {
return fp.front();
}
if (x >= xp.back()) {
return fp.back();
}
const auto j = static_cast<std::size_t>(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<double, 2> {
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<const double> 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<std::complex<float>> scan, std::span<const std::complex<float>> pulse, std::size_t shift,
std::complex<double> amplitude) {
const std::size_t length = scan.size();
for (std::size_t sample = 0U; sample < length; ++sample) {
const std::complex<double> response = pulse[(sample + length - shift) % length];
scan[sample] = std::complex<float>(std::complex<double>(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<double> range;
std::vector<double> centers;
std::vector<std::vector<std::size_t>> members;
};
[[nodiscard]] auto range_rings(const VolumeRequest& frame, const Axes& grid) -> Rings {
Rings rings;
std::vector<std::size_t> 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::size_t>(
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<double>(ring) + 0.5) * kCleanRing);
}
for (const auto voxel : work) {
const auto ring = std::min(static_cast<std::size_t>(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<const float> volume, const Rings& rings) -> std::vector<double> {
std::vector<double> known_rings;
std::vector<double> 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<float> values;
values.reserve(members.size());
for (const auto voxel : members) {
values.push_back(volume[voxel]);
}
known_rings.push_back(static_cast<double>(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<double> medians(count);
for (std::size_t ring = 0U; ring < count; ++ring) {
medians[ring] = interp(static_cast<double>(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<std::ptrdiff_t>(kCleanRingSmooth / 2U);
const auto last = static_cast<std::ptrdiff_t>(count) - 1;
const double weight = 1.0 / static_cast<double>(kCleanRingSmooth);
std::vector<double> 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<std::ptrdiff_t>(ring) + tap, std::ptrdiff_t{0}, last);
smoothed[ring] += medians[static_cast<std::size_t>(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<const float> volume, const Rings& rings) -> std::pair<std::size_t, double> {
std::vector<float> 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<std::size_t, double> 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<double> amplitude{};
};
struct Component {
std::array<double, 3> point{}; // search voxel, scene X, Y, Z
std::vector<Part> 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 <class Image>
[[nodiscard]] auto subtract_components(const Image& image, const CleanRequest& request, Scans& residual)
-> std::vector<Component> {
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<std::ptrdiff_t>(
1, static_cast<std::ptrdiff_t>(std::nearbyint(2.0 * kCleanPick * frame.samples_per_metre)));
std::vector<Component> 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<std::ptrdiff_t>(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<double, 3> target{search.x[ix], y_radar, z_radar}; // in radar axes
const auto range = [&target](const std::array<float, 3>& 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::ptrdiff_t>(
std::nearbyint((range(term.tx) + range(term.rx)) * frame.samples_per_metre));
const auto low = std::max<std::ptrdiff_t>(delay - pick, 0);
const auto high = std::min(delay + pick + 1, static_cast<std::ptrdiff_t>(length));
if (high <= low) {
continue;
}
const auto scan = std::span(residual).subspan(pair * length, length);
auto peak = static_cast<std::size_t>(low);
for (auto sample = peak + 1U; sample < static_cast<std::size_t>(high); ++sample) {
if (std::norm(std::complex<double>(scan[sample])) > std::norm(std::complex<double>(scan[peak]))) {
peak = sample;
}
}
const auto amplitude = std::complex<double>(scan[peak]) / std::complex<double>(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<double, 3> seed{};
std::vector<const Component*> members;
};
[[nodiscard]] auto group(const std::vector<Component>& components) -> std::vector<Group> {
std::vector<Group> 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<const float> values, std::span<const float> 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<float> 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<double>(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<float>(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<float>(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<const CleanObject> objects, const VolumeRequest& frame) -> std::vector<Mark> {
// 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<std::size_t> 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<std::size_t> roots;
std::vector<std::vector<std::size_t>> 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<std::size_t>(found - roots.begin())].push_back(index);
}
}
const auto& x = frame.grid.x;
const auto& y = frame.grid.y;
std::vector<Mark> 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<double>::infinity();
Mark mark;
mark.box = {std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::infinity(),
std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::infinity(),
std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::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
data_acq_and_processing/processing/processors/src/sar3d_clean.hpp
+91
View File
@@ -0,0 +1,91 @@
#pragma once
#include <array>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <span>
#include <vector>
#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<double, 3> 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<const float> values, std::span<const float> 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<const float> 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<const std::complex<float>> ascans; // full A-scan of each frame pair, back to back
std::span<const std::complex<float>> pulse; // zero-delay point response, one A-scan long
bool bg_gate = true; // BG_GATE_ENABLE
};
struct CleanObject {
std::array<std::size_t, 3> index{}; // peak voxel on the frame grid, [y, z, x]
std::array<double, 3> 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<CleanObject> 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<float> 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<double, 6> 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<const CleanObject> objects, const VolumeRequest& frame) -> std::vector<Mark>;
} // namespace radar::processing::sar3d
data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp
+278 -183
View File
@@ -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<float> {
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<int>(4.0 * sigma + 0.5);
std::vector<double> 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<float> taps;
for (const double weight : weights) {
taps.push_back(static_cast<float>(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<const float> levels, std::span<const unsigned int> 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<float>::max()) {
throw std::runtime_error("sar3d grid contains an unsupported coordinate");
}
}
}
}
#endif
} // namespace
auto reference_plane(std::span<const double> 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<float> {
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<int>(4.0 * sigma + 0.5);
std::vector<double> 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<float> taps;
for (const double weight : weights) {
taps.push_back(static_cast<float>(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 <class T>
[[nodiscard]] auto download(CUdeviceptr source, std::size_t count) -> std::vector<T> {
[[nodiscard]] auto copy_to_host(CUdeviceptr source, std::size_t count) -> std::vector<T> {
std::vector<T> 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 <class... Args>
void launch_blocks(CUfunction kernel, unsigned blocks, unsigned shared_bytes, Args... args) {
void* parameters[] = {static_cast<void*>(&args)...};
check(cuLaunchKernel(kernel, blocks, 1U, 1U, kBlockSize, 1U, 1U, shared_bytes, nullptr, parameters, nullptr),
"launch");
}
// One thread per item.
template <class... Args>
void launch(CUfunction kernel, int items, Args... args) {
void* parameters[] = {static_cast<void*>(&args)...};
const auto blocks = (static_cast<unsigned>(items) + kBlockSize - 1U) / kBlockSize;
check(cuLaunchKernel(kernel, blocks, 1U, 1U, kBlockSize, 1U, 1U, 0U, nullptr, parameters, nullptr), "launch");
launch_blocks(kernel, (static_cast<unsigned>(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<const float> weights) {
smooth_radius = static_cast<int>(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<std::size_t>(nx); }
[[nodiscard]] auto gz() const -> CUdeviceptr { return gy() + sizeof(double) * static_cast<std::size_t>(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<const std::complex<float>>, std::size_t, std::size_t, std::size_t) -> std::vector<std::complex<float>>;
auto project(const VolumeRequest&) -> Projections;
};
auto Gpu::Impl::inverse_dft(
std::span<const std::complex<float>> bands,
std::size_t points,
std::size_t first_bin,
std::size_t size
) -> std::vector<std::complex<float>> {
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<float>) * 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<std::complex<float>>(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<double> 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<const double>(values));
const CUdeviceptr output = volume.reserve(sizeof(float) * grid.x.size() * grid.y.size() * grid.z.size());
nx = static_cast<int>(grid.x.size());
ny = static_cast<int>(grid.y.size());
nz = static_cast<int>(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<const double> y, double margin) -> Peak {
const auto planes = static_cast<std::size_t>(ny);
const CUdeviceptr levels = plane_stats.reserve(2U * sizeof(float) * planes);
const CUdeviceptr indices = levels + sizeof(float) * planes;
launch_blocks(peak_planes, static_cast<unsigned>(ny), static_cast<unsigned>(2U * kBlockSize * sizeof(float)),
volume.pointer(), nz * nx, levels, indices);
const auto plane_levels = copy_to_host<float>(levels, planes);
const auto plane_indices = copy_to_host<unsigned int>(indices, planes);
const auto row = static_cast<unsigned int>(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<double> 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<const double>(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<std::size_t>(voxels));
const CUdeviceptr scratch = gpu.scratch.reserve(sizeof(float) * static_cast<std::size_t>(voxels));
auto inverse_dft(
std::span<const std::complex<float>> spectra,
std::size_t points,
std::size_t first_bin,
std::size_t size
) -> std::vector<std::complex<float>> {
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<float>) * 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<std::complex<float>>(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<float>(request.samples_per_metre),
static_cast<int>(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<std::size_t>(ny));
launch(gpu.peak_planes, ny, volume, plane, ny, stats, stats + sizeof(float) * static_cast<std::size_t>(ny));
Projections raw_peak;
const auto raw_levels = download<float>(stats, static_cast<std::size_t>(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<std::size_t>(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<float>(request.samples_per_metre),
static_cast<int>(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<const float>(weights));
const int radius = static_cast<int>(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<std::size_t>(view.nv) * static_cast<std::size_t>(view.nu), "view");
const auto image_floats = static_cast<std::size_t>(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<std::size_t>(front_size);
const CUdeviceptr top = side + sizeof(float) * static_cast<std::size_t>(side_size);
const CUdeviceptr oblique = top + sizeof(float) * static_cast<std::size_t>(top_size);
const CUdeviceptr levels = oblique + sizeof(float) * static_cast<std::size_t>(oblique_size);
const CUdeviceptr indices = levels + sizeof(float) * static_cast<std::size_t>(ny);
launch(gpu.peak_planes, ny, volume, plane, ny, levels, indices);
Projections result;
result.raw_level = raw_peak.level;
const auto plane_levels = download<float>(levels, static_cast<std::size_t>(ny));
const auto plane_indices = download<unsigned int>(indices, static_cast<std::size_t>(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<int>(result.peak[2]), static_cast<int>(result.peak[0]), static_cast<int>(result.peak[1]), widths);
const auto measured = download<float>(widths, 3U);
auto download() -> std::vector<float> {
use_context();
require_volume();
return copy_to_host<float>(volume.pointer(), static_cast<std::size_t>(voxels()));
}
void load(const Axes& grid, std::span<const float> values) {
use_context();
try {
check(cuMemcpyHtoD(use_grid(grid), values.data(), values.size_bytes()), "upload");
} catch (...) {
forget_volume();
throw;
}
}
auto project(const std::array<std::size_t, 3>& peak, const ObliqueView& view) -> Projections {
use_context();
require_volume();
const auto [py, pz, px] = peak;
if (py >= static_cast<std::size_t>(ny) || pz >= static_cast<std::size_t>(nz) ||
px >= static_cast<std::size_t>(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<std::size_t>(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<std::size_t>(front_size);
const CUdeviceptr top = side + sizeof(float) * static_cast<std::size_t>(side_size);
const CUdeviceptr oblique = top + sizeof(float) * static_cast<std::size_t>(top_size);
const CUdeviceptr widths = oblique + sizeof(float) * static_cast<std::size_t>(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<std::size_t>(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<int>(px), static_cast<int>(py), static_cast<int>(pz), widths);
Projections result;
result.front = copy_to_host<float>(front, static_cast<std::size_t>(front_size));
result.side = copy_to_host<float>(side, static_cast<std::size_t>(side_size));
result.top = copy_to_host<float>(top, static_cast<std::size_t>(top_size));
result.oblique = copy_to_host<float>(oblique, static_cast<std::size_t>(oblique_size));
const auto measured = copy_to_host<float>(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<std::size_t>(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<float>(front, static_cast<std::size_t>(front_size));
result.side = download<float>(side, static_cast<std::size_t>(side_size));
result.top = download<float>(top, static_cast<std::size_t>(top_size));
result.oblique = download<float>(oblique, static_cast<std::size_t>(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<const std::complex<float>>, std::size_t, std::size_t, std::size_t)
-> std::vector<std::complex<float>> { missing(); }
auto image(const VolumeRequest&) -> Peak { missing(); }
auto download() -> std::vector<float> { missing(); }
void load(const Axes&, std::span<const float>) { missing(); }
auto project(const std::array<std::size_t, 3>&, const ObliqueView&) -> Projections { missing(); }
};
#endif
Gpu::Gpu(double smooth_sigma_voxels) {
Gpu::Gpu() {
#ifdef SAR3D_WITH_CUDA
impl_ = std::make_unique<Impl>(gaussian_taps(smooth_sigma_voxels));
impl_ = std::make_unique<Impl>();
#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<const std::complex<float>> bands, std::size_t points,
std::size_t first_bin, std::size_t size) -> std::vector<std::complex<float>> {
if (points == 0U || bands.empty() || bands.size() % points != 0U || !std::has_single_bit(size) ||
@@ -585,36 +677,20 @@ auto Gpu::inverse_dft(std::span<const std::complex<float>> 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<float>::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<float>::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<float>::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<float> {
return impl_->download();
}
void Gpu::load(const Axes& grid, std::span<const float> 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<std::size_t, 3>& 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<std::size_t>(view.nu), static_cast<std::size_t>(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
data_acq_and_processing/processing/processors/src/sar3d_gpu.hpp
+39 -16
View File
@@ -9,8 +9,8 @@
#include <vector>
// 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<const double> x;
std::span<const double> y;
std::span<const double> 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<const double> 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<const Contribution> contributions;
std::span<const Contribution> contributions{};
// One A-scan per contribution, back to back; sample k is the Tx+Rx path k / samples_per_metre.
std::span<const std::complex<float>> ascans;
std::span<const std::complex<float>> 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<std::size_t, 3> index{}; // first voxel at `level`, [y, z, x]
};
// Maximum projections of the current volume (row-major), unscaled.
struct Projections {
std::vector<float> front; // [y][x], maximum along Z
std::vector<float> side; // [y][z], maximum along X
std::vector<float> top; // [z][x], maximum along Y
std::vector<float> oblique; // [v][u]
float level = 0.0F;
float raw_level = 0.0F;
std::array<std::size_t, 3> peak{}; // joint voxel maximum, [y, z, x]
std::array<float, 3> target_width{}; // contiguous FWHM through the peak, [X, Y, Z], metres
std::array<float, 3> 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<std::complex<float>>;
// 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<float>;
// Makes non-negative host values on `grid` the current volume.
void load(const Axes& grid, std::span<const float> values);
// Max projections of the current volume and its half-maximum widths through `peak` ([y, z, x]).
[[nodiscard]] auto project(const std::array<std::size_t, 3>& peak, const ObliqueView& view) -> Projections;
private:
struct Impl;
data_acq_and_processing/processing/processors/src/sar3d_processor.cpp
+231 -81
View File
@@ -9,23 +9,28 @@
#include <cstddef>
#include <cstdint>
#include <exception>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <memory>
#include <numbers>
#include <optional>
#include <set>
#include <span>
#include <sstream>
#include <stdexcept>
#include <string>
#include <string_view>
#include <tuple>
#include <utility>
#include <vector>
#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<double, 3> 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<double>(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<double> values, double q) -> double {
const double position = q / 100.0 * static_cast<double>(values.size() - 1U);
const auto lower = static_cast<std::size_t>(position);
const auto lower_it = values.begin() + static_cast<std::ptrdiff_t>(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<double> {
std::vector<double> window(count);
const auto last = static_cast<double>(count - 1U);
for (std::size_t index = 0U; index < count; ++index) {
window[index] = 0.5 - 0.5 * std::cos(2.0 * pi * static_cast<double>(index) / last);
}
const double upper = *std::min_element(lower_it + 1, values.end());
return *lower_it + (position - static_cast<double>(lower)) * (upper - *lower_it);
return window;
}
struct Band {
@@ -197,7 +203,7 @@ struct Band {
steps.push_back(static_cast<double>(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<double>(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<double>(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<std::complex<float>> ascans{}; // [pair][sample]
// A-scan of a point reflector at zero delay, the same band and window (sar_clean_unit_pulse).
std::vector<std::complex<float>> 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<float>(weight));
}
Frame frame{};
frame.center_ns = timing.center_ns;
frame.samples = first.size;
frame.dt_s = 1.0 / (static_cast<double>(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<std::ptrdiff_t>(frame.samples), frame.ascans.end());
frame.ascans.resize(frame.ascans.size() - frame.samples);
std::vector<std::span<std::complex<float>>> 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<int>::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<std::size_t>(count);
}
@@ -327,12 +340,19 @@ struct Grid {
std::vector<double> x{};
std::vector<double> y{};
std::vector<double> 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<std::array<double, 2>, 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<double, 2>& 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::size_t>(std::numeric_limits<int>::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<double, 3>& step, std::string_view remedy)
-> Grid {
std::array<std::size_t, 3> 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::size_t>(std::numeric_limits<int>::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<int>(checked_axis_size(std::round((u_max - view.u0) / view.pixel_m) + 1.0));
view.nv = static_cast<int>(checked_axis_size(std::round((v_max - view.v0) / view.pixel_m) + 1.0));
view.nu = static_cast<int>(checked_axis_size(std::round((u_max - view.u0) / view.pixel_m) + 1.0, kVoxelRemedy));
view.nv = static_cast<int>(checked_axis_size(std::round((v_max - view.v0) / view.pixel_m) + 1.0, kVoxelRemedy));
if (view.nu > std::numeric_limits<int>::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<float> 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<std::size_t, 3> peak{}; // marked voxel [y, z, x]: the maximum or the strongest object
std::optional<sar3d::FrameBackground> background; // when BG_GATE_ENABLE or CLEAN needs it
std::optional<sar3d::Clean> 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<float>(static_cast<double>(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<sar3d::Gpu>(kSmoothSigma); }
try { device = std::make_unique<sar3d::Gpu>(); }
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<double>(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<std::array<double, 3>> 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<double>(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<float>(static_cast<double>(value) / normalization);
if (!std::isfinite(value)) { throw std::runtime_error("SAR normalization exceeds the float32 range"); }
}
}
const double normalization = static_cast<double>(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<ipc::ResultPayload> 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<float> 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<float>(normalization), static_cast<float>(grid.x[peak_x]),
static_cast<float>(grid.y[peak_y]),
static_cast<float>(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<float>(normalization), static_cast<float>(grid.x[peak_x]),
static_cast<float>(grid.y[peak_y]), static_cast<float>(grid.z[peak_z]),
}));
if (settings.bg_gate_enabled) {
payloads.push_back(table_payload("sar3d_background", 3U, {
static_cast<float>(display.background.value().ratio()),
static_cast<float>(sar3d::kBgGateShow), static_cast<float>(sar3d::kBgGateFull),
}));
}
if (display.has_objects()) {
// The objects replace the peak marker, so its width is not published.
std::vector<float> object_rows{};
for (const auto& object : display.clean->objects) {
for (const double coordinate : object.position) { object_rows.push_back(static_cast<float>(coordinate)); }
object_rows.push_back(static_cast<float>(object.level));
}
payloads.push_back(table_payload("sar3d_objects", 4U, std::move(object_rows)));
std::vector<float> mark_rows{};
for (const auto& mark : sar3d::marks(display.clean->objects, request)) {
mark_rows.push_back(static_cast<float>(mark.kind));
for (const double edge : mark.box) { mark_rows.push_back(static_cast<float>(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;
}