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

Refactor radar system

Browse Source
This commit is contained in:
Ayzen
2026-10-05 16:46:28 +03:00
parent 515049d2ed
commit 8f097dfa2e
179 changed files with 8212 additions and 6676 deletions
Download Patch File Download Diff File Expand all files Collapse all files
data_acq_and_processing/processing/processors/src/gpr_legacy_processor.ipp
+3 -3
View File
@@ -20,10 +20,10 @@
// 5. Optional extended-mode region detection for diffuse reflectors.
//
// Sweep-event model — radar topology matters:
// * Matrix radars (`librevna_multi`, `sn9000`) fire one Tx at a time and
// * Matrix radars (`librevna_multi`) fire one Tx at a time and
// receive on all Rx channels in parallel. A run of 8 traces is just 2
// Tx-events; all (tx_k, *) pairs share one timestamp.
// * Sequential radars (single librevna with switches, kamil_adc, k209…)
// * Sequential radars (single librevna with switches, kamil_adc)
// measure each pair separately. A run of 8 traces is 8 separate sweeps;
// every pair has its own timestamp.
// Both cases reduce to: there are N events in the frame, each event lasts
@@ -283,7 +283,7 @@ enum class LegacyPeakDomain {
}
[[nodiscard]] auto is_matrix_radar_model(const std::string& model) -> bool {
return model == "librevna_multi" || model == "sn9000";
return model == "librevna_multi";
}
// Assign an `event_index` to every selected pair. The mapping depends on the
data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp
+225 -87
View File
@@ -1,19 +1,76 @@
#include "sar3d_gpu.hpp"
#include <algorithm>
#include <bit>
#include <climits>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <string>
#ifdef SAR3D_WITH_CUDA
#include <cuda.h>
#include <nvrtc.h>
#include <algorithm>
#include <climits>
#include <cmath>
#include <string>
#include <type_traits>
#endif
namespace radar::processing::sar3d {
namespace {
constexpr std::size_t kMaxIndexedItems = static_cast<std::size_t>(INT_MAX);
auto checked_product(std::size_t a, std::size_t b, std::size_t limit, const char* what) -> std::size_t {
if (a == 0U || b == 0U || a > limit / b) {
throw std::runtime_error(std::string("sar3d ") + what + " is too large or empty");
}
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");
}
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()};
}
}
}
}
#endif
} // namespace
#ifdef SAR3D_WITH_CUDA
namespace {
@@ -39,8 +96,9 @@ extern "C" __global__ void inverse_dft(
const float2* __restrict__ bands, const int points, const int first_bin,
const int size, const int items, float2* __restrict__ ascans)
{
const int item = blockIdx.x * blockDim.x + threadIdx.x;
if (item >= items) return;
const unsigned index = blockIdx.x * blockDim.x + threadIdx.x;
if (index >= (unsigned)items) return;
const int item = (int)index;
const unsigned k = (unsigned)(item % size);
const unsigned mask = (unsigned)size - 1u;
const float2* band = bands + (item / size) * points;
@@ -72,8 +130,9 @@ extern "C" __global__ void backproject(
const float comp_power, const float comp_max,
float* __restrict__ volume)
{
const int voxel = blockIdx.x * blockDim.x + threadIdx.x;
if (voxel >= voxels) return;
const unsigned index = blockIdx.x * blockDim.x + threadIdx.x;
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.
@@ -87,7 +146,7 @@ extern "C" __global__ void backproject(
const float px = (float)gx[voxel - row * nx];
const float py = (float)(y0 * cos_p + zs * sin_p);
const float pz = (float)depth;
const float last = (float)(samples - 1);
const double last = (double)(samples - 1);
float re = 0.0f, im = 0.0f, weights = 0.0f;
for (int c = 0; c < count; ++c) {
const Contribution k = contributions[c];
@@ -96,7 +155,7 @@ extern "C" __global__ void backproject(
const float rt = sqrtf(dxt * dxt + dyt * dyt + dzt * dzt);
const float rr = sqrtf(dxr * dxr + dyr * dyr + dzr * dzr);
const float position = (rt + rr) * samples_per_metre;
if (position > last) continue;
if (!isfinite(position) || position < 0.0f || (double)position > last) continue;
const int i = min((int)position, samples - 2);
const float frac = position - (float)i;
const float2* row_ascan = ascans + (long long)c * samples;
@@ -111,7 +170,7 @@ extern "C" __global__ void backproject(
im += w * (a.y + frac * (b.y - a.y));
weights += k.weight;
}
volume[voxel] = sqrtf(re * re + im * im) / (weights + 1e-12f);
volume[voxel] = hypotf(re, im) / (weights + 1e-12f);
}
// One pass of scipy.ndimage.gaussian_filter1d (mode 'reflect') along Z or X.
@@ -120,8 +179,9 @@ extern "C" __global__ void smooth(
const float* __restrict__ taps, const int radius,
const int nx, const int nz, const int voxels, const int along_z)
{
const int voxel = blockIdx.x * blockDim.x + threadIdx.x;
if (voxel >= voxels) return;
const unsigned thread = blockIdx.x * blockDim.x + threadIdx.x;
if (thread >= (unsigned)voxels) return;
const int voxel = (int)thread;
const int length = along_z ? nz : nx;
const int index = along_z ? (voxel / nx) % nz : voxel % nx;
const int stride = along_z ? nx : 1;
@@ -134,17 +194,36 @@ extern "C" __global__ void smooth(
dst[voxel] = sum;
}
// dst[o] = max over r of src[(o / inner) * outer_stride + (o % inner) * inner_stride + r * reduce_stride].
// Strides select the displayed axes and the axis reduced to its maximum.
extern "C" __global__ void max_project(
const float* __restrict__ src, float* __restrict__ dst, const int outputs, const int inner,
const int outer_stride, const int inner_stride, const int reduce, const int reduce_stride)
{
const int o = blockIdx.x * blockDim.x + threadIdx.x;
if (o >= outputs) return;
const float* p = src + (o / inner) * outer_stride + (o % inner) * inner_stride;
float m = 0.0f;
for (int r = 0; r < reduce; ++r) m = fmaxf(m, p[r * reduce_stride]);
dst[o] = m;
const unsigned thread = blockIdx.x * blockDim.x + threadIdx.x;
if (thread >= (unsigned)outputs) return;
const int o = (int)thread;
const float* row = src + (o / inner) * outer_stride + (o % inner) * inner_stride;
float maximum = 0.0f;
for (int index = 0; index < reduce; ++index) maximum = fmaxf(maximum, row[index * reduce_stride]);
dst[o] = maximum;
}
// Keep the first joint (y,z,x) maximum. Separate projections cannot locate it.
extern "C" __global__ void peak_planes(
const float* volume, const int plane, const int ny, float* levels, unsigned int* 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;
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; }
}
levels[y] = level;
indices[y] = index;
}
// Orthographic max projection. The volume is non-negative, so float bit
@@ -158,8 +237,9 @@ extern "C" __global__ void project_oblique(
const double u0, const double v0, const double inv_pixel, const int nu, const int nv,
unsigned int* __restrict__ image)
{
const int voxel = blockIdx.x * blockDim.x + threadIdx.x;
if (voxel >= voxels) return;
const unsigned index = blockIdx.x * blockDim.x + threadIdx.x;
if (index >= (unsigned)voxels) return;
const int voxel = (int)index;
const float value = volume[voxel];
if (!(value > 0.0f)) return;
const int row = voxel / nx;
@@ -175,11 +255,13 @@ extern "C" __global__ void project_oblique(
constexpr int kBlockSize = 256;
void check(CUresult result, const char* what) {
void check(CUresult result, const char* what, bool device_required = false) {
if (result != CUDA_SUCCESS) {
const char* text = nullptr;
cuGetErrorString(result, &text);
throw std::runtime_error(std::string{"CUDA "} + what + ": " + (text != nullptr ? text : "unknown error"));
const auto message = std::string{"CUDA "} + what + ": " + (text != nullptr ? text : "unknown error");
if (device_required) { throw CudaUnavailable(message); }
throw std::runtime_error(message);
}
}
@@ -299,20 +381,20 @@ template <class T>
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 - 1) / kBlockSize);
const auto blocks = (static_cast<unsigned>(items) + kBlockSize - 1U) / kBlockSize;
check(cuLaunchKernel(kernel, blocks, 1U, 1U, kBlockSize, 1U, 1U, 0U, nullptr, parameters, nullptr), "launch");
}
} // namespace
struct Gpu::Impl {
struct CudaBackend {
// Declared first, destroyed last: buffers and the module need the context.
struct Context {
CUdevice device = 0;
CUcontext handle = nullptr;
Context() {
check(cuInit(0U), "init");
check(cuDeviceGet(&device, 0), "device");
check(cuInit(0U), "init", true);
check(cuDeviceGet(&device, 0), "device", true);
check(cuDevicePrimaryCtxRetain(&handle, device), "context");
if (const CUresult result = cuCtxSetCurrent(handle); result != CUDA_SUCCESS) {
cuDevicePrimaryCtxRelease(device);
@@ -331,10 +413,11 @@ struct Gpu::Impl {
return kernel;
}
} module{compile_kernels(context.device)};
CUfunction inverse_dft = module.function("inverse_dft");
CUfunction inverse_dft_kernel = module.function("inverse_dft");
CUfunction backproject = module.function("backproject");
CUfunction smooth = module.function("smooth");
CUfunction max_project = module.function("max_project");
CUfunction peak_planes = module.function("peak_planes");
CUfunction project_oblique = module.function("project_oblique");
int smooth_radius = 0;
DeviceBuffer taps;
@@ -345,50 +428,36 @@ struct Gpu::Impl {
DeviceBuffer volume;
DeviceBuffer scratch;
DeviceBuffer images;
DeviceBuffer peak_stats;
~Impl() { cuCtxSetCurrent(context.handle); }
~CudaBackend() { cuCtxSetCurrent(context.handle); }
explicit CudaBackend(std::span<const float> weights) {
smooth_radius = static_cast<int>(weights.size() / 2U);
taps.upload(weights);
}
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;
};
Gpu::Gpu(double smooth_sigma_voxels) : impl_(std::make_unique<Impl>()) {
// scipy.ndimage.gaussian_filter1d taps: radius int(4 sigma + 0.5), normalized.
const int radius = static_cast<int>(4.0 * smooth_sigma_voxels + 0.5);
std::vector<double> weights{};
for (int k = -radius; k <= radius; ++k) {
weights.push_back(std::exp(-0.5 * k * k / (smooth_sigma_voxels * smooth_sigma_voxels)));
}
double total = 0.0;
for (const double weight : weights) {
total += weight;
}
std::vector<float> taps{};
for (const double weight : weights) {
taps.push_back(static_cast<float>(weight / total));
}
impl_->smooth_radius = radius;
impl_->taps.upload(std::span<const float>(taps));
}
Gpu::~Gpu() = default;
auto Gpu::inverse_dft(
auto CudaBackend::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 = *impl_;
CudaBackend& 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, items, input, checked_int(points, "band"), checked_int(first_bin, "first bin"),
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);
}
auto Gpu::project(const VolumeRequest& request) -> Projections {
Impl& gpu = *impl_;
auto CudaBackend::project(const VolumeRequest& request) -> Projections {
CudaBackend& gpu = *this;
check(cuCtxSetCurrent(gpu.context.handle), "context");
const int nx = checked_int(request.x.size(), "grid");
@@ -413,6 +482,12 @@ auto Gpu::project(const VolumeRequest& request) -> Projections {
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);
@@ -423,49 +498,112 @@ auto Gpu::project(const VolumeRequest& request) -> Projections {
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);
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 int plane = nz * nx;
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);
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);
return Projections{
.front = download<float>(front, static_cast<std::size_t>(front_size)),
.side = download<float>(side, static_cast<std::size_t>(side_size)),
.top = download<float>(top, static_cast<std::size_t>(top_size)),
.oblique = download<float>(oblique, static_cast<std::size_t>(oblique_size)),
};
}
#else // !SAR3D_WITH_CUDA
struct Gpu::Impl {};
Gpu::Gpu(double /*smooth_sigma_voxels*/) {
throw std::runtime_error(
"data_processor was built without CUDA (nvrtc.h, cuda.h, libnvrtc or libcuda not found), "
"so sar3d cannot run; install the CUDA requirements and rebuild"
);
}
Gpu::~Gpu() = default;
auto Gpu::inverse_dft(std::span<const std::complex<float>>, std::size_t, std::size_t, std::size_t)
-> std::vector<std::complex<float>> {
return {};
}
auto Gpu::project(const VolumeRequest& /*request*/) -> Projections {
return {};
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;
}
#endif
struct Engine::Impl {
#ifdef SAR3D_WITH_CUDA
std::unique_ptr<CudaBackend> cuda;
#endif
};
Engine::Engine(double smooth_sigma_voxels) : impl_(std::make_unique<Impl>()) {
#ifdef SAR3D_WITH_CUDA
impl_->cuda = std::make_unique<CudaBackend>(gaussian_taps(smooth_sigma_voxels));
#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
}
Engine::~Engine() = default;
auto Engine::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) ||
size > kMaxIndexedItems || first_bin >= size || points > size - first_bin) {
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_->cuda->inverse_dft(bands, points, first_bin, size);
#else
throw CudaUnavailable("SAR3D requires an NVIDIA GPU with CUDA");
#endif
}
auto Engine::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");
}
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& sample : request.ascans) {
if (!std::isfinite(sample.real()) || !std::isfinite(sample.imag())) {
throw std::runtime_error("sar3d A-scan contains a non-finite sample");
}
}
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) {
throw std::runtime_error("sar3d aperture contains an invalid weight or reference range");
}
for (const auto point : {term.tx, term.rx}) {
if (!std::ranges::all_of(point, [](float value) { return std::isfinite(value); })) {
throw std::runtime_error("sar3d aperture contains a non-finite antenna position");
}
}
}
(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_->cuda->project(request);
#else
throw CudaUnavailable("SAR3D requires an NVIDIA GPU with CUDA");
#endif
}
} // namespace radar::processing::sar3d
data_acq_and_processing/processing/processors/src/sar3d_gpu.hpp
+20 -10
View File
@@ -5,13 +5,19 @@
#include <cstddef>
#include <memory>
#include <span>
#include <stdexcept>
#include <vector>
// GPU half of the sar3d processor: A-scans of the band spectra, and the
// Numerical 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 of the volume leave the GPU.
// Only the projections and peak leave the backend.
namespace radar::processing::sar3d {
class CudaUnavailable : public std::runtime_error {
public:
using std::runtime_error::runtime_error;
};
// One (aperture frame, Tx/Rx pair) term of the backprojection sum. The layout
// is shared with the CUDA kernel.
struct Contribution {
@@ -51,25 +57,29 @@ struct VolumeRequest {
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{};
};
// Max-projections of the smoothed, not yet normalized volume (row-major).
// Maximum projections of the smoothed volume (row-major), before normalization.
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]
};
// CUDA context with the kernels compiled for the installed GPU. Construction
// throws when no usable GPU is present or the processor was built without CUDA.
class Gpu {
// CUDA context and runtime-compiled kernels. A missing device raises
// CudaUnavailable; computation and compilation failures remain ordinary errors.
class Engine {
public:
explicit Gpu(double smooth_sigma_voxels);
~Gpu();
Gpu(const Gpu&) = delete;
auto operator=(const Gpu&) -> Gpu& = delete;
explicit Engine(double smooth_sigma_voxels);
~Engine();
Engine(const Engine&) = delete;
auto operator=(const Engine&) -> Engine& = delete;
// numpy.fft.ifft of each band spectrum placed at bins [first_bin, first_bin +
// points) of a `size`-point spectrum; `bands` holds the bands back to back.
data_acq_and_processing/processing/processors/src/sar3d_processor.cpp
+142 -85
View File
@@ -7,11 +7,13 @@
#include <complex>
#include <cstddef>
#include <cstdint>
#include <exception>
#include <iostream>
#include <limits>
#include <map>
#include <memory>
#include <numbers>
#include <set>
#include <span>
#include <stdexcept>
#include <string>
@@ -158,22 +160,26 @@ void select_antennas(
struct TraceClock {
double start_s = 0.0; // sweep start
double event_s = 0.0; // sweep duration / (stages * points)
long long stage = 0;
long long stages = 1;
std::size_t stage = 0U;
std::size_t stages = 1U;
[[nodiscard]] auto sample_time_s(std::size_t point) const -> double {
return start_s + static_cast<double>(stages * static_cast<long long>(point) + stage) * event_s;
return start_s + (static_cast<double>(stages) * point + stage) * event_s;
}
};
struct CollectionTiming {
double center_s = 0.0; // absolute acquisition centre, positions the frame in the aperture
std::uint64_t center_ns = 0U; // integer centre avoids overflow and loss of relative timing
double center_fraction_ns = 0.0;
std::vector<TraceClock> clocks; // per collection trace
};
// Timing from the per-trace capture windows (snapshot_acquisition_timing); every
// trace takes part, so the stages do not depend on the INPUT_IDX/OUTPUT_IDX choice.
[[nodiscard]] auto collection_timing(const ipc::PreprocessedCollection& collection) -> CollectionTiming {
if (collection.traces.empty()) {
throw std::runtime_error("collection has no traces");
}
std::map<std::pair<std::uint64_t, std::uint64_t>, std::vector<std::uint32_t>> sweep_outputs{};
std::uint64_t first = std::numeric_limits<std::uint64_t>::max();
std::uint64_t last = 0U;
@@ -181,6 +187,9 @@ struct CollectionTiming {
if (trace.capture_end_ns <= trace.capture_start_ns) {
throw std::runtime_error("collection has traces without a capture window");
}
if (trace.frequency_hz.empty()) {
throw std::runtime_error("collection has a trace without frequency points");
}
sweep_outputs[{trace.capture_start_ns, trace.capture_end_ns}].push_back(trace.combo.output_pos);
first = std::min(first, trace.capture_start_ns);
last = std::max(last, trace.capture_end_ns);
@@ -191,18 +200,19 @@ struct CollectionTiming {
}
CollectionTiming timing{};
const auto center_x2_ns = static_cast<std::int64_t>(first + last);
timing.center_s = 0.5e-9 * static_cast<double>(center_x2_ns);
const auto duration_ns = last - first;
timing.center_ns = first + duration_ns / 2U;
timing.center_fraction_ns = 0.5 * static_cast<double>(duration_ns % 2U);
for (const auto& trace : collection.traces) {
const auto& stages = sweep_outputs.at({trace.capture_start_ns, trace.capture_end_ns});
const auto stage = std::find(stages.begin(), stages.end(), trace.combo.output_pos) - stages.begin();
TraceClock clock{};
clock.start_s =
0.5e-9 * static_cast<double>(2 * static_cast<std::int64_t>(trace.capture_start_ns) - center_x2_ns);
clock.start_s = 1e-9 * (static_cast<double>(trace.capture_start_ns - first) -
0.5 * static_cast<double>(duration_ns));
clock.stage = stage;
clock.stages = static_cast<long long>(stages.size());
clock.stages = stages.size();
clock.event_s = 1e-9 * static_cast<double>(trace.capture_end_ns - trace.capture_start_ns) /
static_cast<double>(clock.stages * static_cast<long long>(trace.frequency_hz.size()));
static_cast<double>(clock.stages) / trace.frequency_hz.size();
timing.clocks.push_back(clock);
}
return timing;
@@ -226,6 +236,10 @@ struct Band {
}
std::vector<std::size_t> points{};
for (std::size_t point = 0U; point < frequency.size(); ++point) {
if (!std::isfinite(frequency[point]) || frequency[point] < 0.0F ||
(point > 0U && frequency[point] <= frequency[point - 1U])) {
throw std::runtime_error("trace frequencies must be finite, non-negative and increasing");
}
if (frequency[point] >= settings.f_start_hz && frequency[point] <= settings.f_stop_hz) {
points.push_back(point);
}
@@ -242,7 +256,15 @@ struct Band {
if (!(band.df_hz > 0.0)) {
throw std::runtime_error("trace frequencies do not increase");
}
band.first_bin = static_cast<std::size_t>(std::nearbyint(frequency[points.front()] / band.df_hz));
// Bound before converting to an integer or calling bit_ceil. CUDA uses
// 32-bit indexing; pathological sub-Hz steps must not wrap CUDA indices.
constexpr std::size_t max_fft = 1U << 30U; // largest power of two inside CUDA's signed index
const double first_bin = std::nearbyint(frequency[points.front()] / band.df_hz);
if (!std::isfinite(first_bin) || first_bin < 0.0 ||
first_bin + points.size() - 1U > max_fft / (2U * kOversample)) {
throw std::runtime_error("frequency band exceeds the 32-bit CUDA A-scan indexing");
}
band.first_bin = static_cast<std::size_t>(first_bin);
band.size = std::bit_ceil(2U * (band.first_bin + points.size() - 1U)) * kOversample;
// The radar moves SPEED_M_S * cos(LOOK_ANGLE_DEG) along range; its offset
@@ -255,6 +277,9 @@ struct Band {
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 std::complex<double> sample(trace.s21[point].re, trace.s21[point].im);
if (!std::isfinite(sample.real()) || !std::isfinite(sample.imag())) {
throw std::runtime_error("S21 contains a non-finite sample");
}
band.spectrum.emplace_back(
sample * std::polar(hann, phase_per_hz_s * frequency[point] * clock.sample_time_s(point))
);
@@ -289,7 +314,8 @@ void normalize_ascan(std::span<std::complex<float>> ascan, double dt_s, const Se
// A-scans of every selected Tx/Rx pair of one collection, (Tx, Rx) order.
struct Frame {
double center_s = 0.0;
std::uint64_t center_ns = 0U;
double center_fraction_ns = 0.0;
double dt_s = 0.0;
std::size_t samples = 0U; // per A-scan
std::vector<std::complex<float>> ascans{}; // [pair][sample]
@@ -298,7 +324,7 @@ struct Frame {
[[nodiscard]] auto compute_frame(
const ipc::PreprocessedCollection& collection,
const Settings& settings,
sar3d::Gpu& gpu
sar3d::Engine& backend
) -> Frame {
const auto timing = collection_timing(collection);
std::vector<Band> bands{};
@@ -327,10 +353,11 @@ struct Frame {
spectra.insert(spectra.end(), band.spectrum.begin(), band.spectrum.end());
}
Frame frame{};
frame.center_s = timing.center_s;
frame.center_ns = timing.center_ns;
frame.center_fraction_ns = timing.center_fraction_ns;
frame.samples = first.size;
frame.dt_s = 1.0 / (static_cast<double>(first.size) * first.df_hz);
frame.ascans = gpu.inverse_dft(spectra, first.spectrum.size(), first.first_bin, first.size);
frame.ascans = backend.inverse_dft(spectra, first.spectrum.size(), first.first_bin, first.size);
if (settings.normalize_each_frame) {
for (std::size_t offset = 0U; offset < frame.ascans.size(); offset += frame.samples) {
normalize_ascan(std::span(frame.ascans).subspan(offset, frame.samples), frame.dt_s, settings);
@@ -350,9 +377,11 @@ struct Frame {
return axis;
}
// Axis with about the given pitch (sar_build_full_region_grid: round(span/pitch)+1 points).
[[nodiscard]] auto grid_axis(double start, double stop, double pitch) -> std::vector<double> {
return linspace(start, stop, static_cast<std::size_t>(std::max(std::nearbyint((stop - start) / pitch) + 1.0, 2.0)));
[[nodiscard]] auto checked_axis_size(double count) -> 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");
}
return static_cast<std::size_t>(count);
}
struct Grid {
@@ -367,6 +396,9 @@ struct Grid {
[[nodiscard]] auto scene_grid(const Settings& settings) -> Grid {
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) {
throw std::runtime_error("SAR grid is undefined at a radar look angle of 90 degrees");
}
double z_min = std::numeric_limits<double>::infinity();
double z_max = -std::numeric_limits<double>::infinity();
for (const double height : {kSceneYMin, kSceneYMax}) {
@@ -377,11 +409,19 @@ struct Grid {
z_max = std::max(z_max, z);
}
}
return {
grid_axis(kSceneXMin, kSceneXMax, settings.voxel_m),
grid_axis(kSceneYMin, kSceneYMax, settings.voxel_m),
grid_axis(z_min, z_max, settings.voxel_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(kSceneXMax - kSceneXMin);
const auto ny = axis_size(kSceneYMax - kSceneYMin);
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");
}
return {linspace(kSceneXMin, kSceneXMax, nx), linspace(kSceneYMin, kSceneYMax, ny),
linspace(z_min, z_max, nz)};
}
// Cosine taper over the aperture: 1 at the centre frame, SAR_EDGE_FRAME_WEIGHT at
@@ -463,8 +503,11 @@ struct Oblique {
}
}
}
view.nu = static_cast<int>(std::lround((u_max - view.u0) / view.pixel_m)) + 1;
view.nv = static_cast<int>(std::lround((v_max - view.v0) / view.pixel_m)) + 1;
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));
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");
}
for (int index = 0; index < view.nu; ++index) {
oblique.u.push_back(view.u0 + index * view.pixel_m);
}
@@ -479,6 +522,11 @@ struct Oblique {
}
[[nodiscard]] auto to_float(const std::array<double, 3>& point) -> std::array<float, 3> {
for (const double coordinate : point) {
if (!std::isfinite(coordinate) || std::abs(coordinate) > std::numeric_limits<float>::max()) {
throw std::runtime_error("SAR antenna motion exceeds the supported coordinate range");
}
}
return {static_cast<float>(point[0]), static_cast<float>(point[1]), static_cast<float>(point[2])};
}
@@ -511,6 +559,15 @@ struct Oblique {
return payload;
}
[[nodiscard]] auto status_payload(float code, const std::string& message) -> ipc::ResultPayload {
std::vector<float> values{code};
for (const unsigned char byte : message) {
values.push_back(byte);
}
const auto columns = static_cast<std::uint32_t>(values.size());
return table_payload("sar3d_status", columns, std::move(values));
}
} // namespace
struct Sar3dProcessor::State {
@@ -521,44 +578,42 @@ struct Sar3dProcessor::State {
return {collection.collection_id, collection.monotonic_ns};
}
// A-scans by sweep for `frames_key`; null marks a sweep that could not be
// processed, so its error is reported once and not retried.
std::map<SweepKey, std::shared_ptr<const Frame>> frames{};
struct CachedFrame {
std::shared_ptr<const Frame> value;
std::string error;
};
std::map<SweepKey, CachedFrame> frames{};
AscanKey frames_key{};
std::unique_ptr<sar3d::Gpu> gpu{};
std::string gpu_error{};
std::unique_ptr<sar3d::Engine> engine{};
std::exception_ptr engine_error{};
std::string last_error{};
std::uint64_t repeated_errors = 0U;
auto frame(const ipc::PreprocessedCollection& collection, const Settings& settings, sar3d::Gpu& gpu)
auto frame(const ipc::PreprocessedCollection& collection, const Settings& settings, sar3d::Engine& backend)
-> std::shared_ptr<const Frame> {
if (const auto found = frames.find(sweep_key(collection)); found != frames.end()) {
return found->second;
if (found->second.value == nullptr) { report_error(found->second.error); }
return found->second.value;
}
std::shared_ptr<const Frame> frame{};
std::string failure;
try {
frame = std::make_shared<const Frame>(compute_frame(collection, settings, gpu));
frame = std::make_shared<const Frame>(compute_frame(collection, settings, backend));
} catch (const std::exception& error) {
report_error(error.what());
failure = error.what();
report_error(failure);
}
frames.emplace(sweep_key(collection), frame);
frames.emplace(sweep_key(collection), CachedFrame{frame, std::move(failure)});
return frame;
}
// Set up once: compiling the kernels takes seconds, so a failure is kept
// instead of being retried on every frame.
auto gpu_engine() -> sar3d::Gpu& {
if (gpu == nullptr && gpu_error.empty()) {
try {
gpu = std::make_unique<sar3d::Gpu>(kSmoothSigma);
} catch (const std::exception& error) {
gpu_error = error.what();
}
auto numerical_engine() -> sar3d::Engine& {
if (engine == nullptr && engine_error == nullptr) {
try { engine = std::make_unique<sar3d::Engine>(kSmoothSigma); }
catch (...) { engine_error = std::current_exception(); }
}
if (gpu == nullptr) {
throw std::runtime_error(gpu_error);
}
return *gpu;
if (engine == nullptr) { std::rethrow_exception(engine_error); }
return *engine;
}
auto render(
@@ -566,23 +621,29 @@ struct Sar3dProcessor::State {
std::span<const ipc::PreprocessedCollection> previous,
const Settings& settings
) -> std::vector<ipc::ResultPayload> {
auto& gpu = gpu_engine();
auto& backend = numerical_engine();
if (settings.ascan_key() != frames_key) {
frames.clear();
frames_key = settings.ascan_key();
}
const auto newest = frame(collection, settings, gpu);
const auto history_count = settings.frame_radius >= previous.size() / 2U + previous.size() % 2U
? previous.size() : 2U * settings.frame_radius;
const auto recent = previous.last(history_count);
std::set<SweepKey> retained{sweep_key(collection)};
for (const auto& member : recent) { retained.insert(sweep_key(member)); }
// Prune even when the newest frame is invalid, so a bad input stream
// cannot grow the failed-frame cache forever.
std::erase_if(frames, [&](const auto& entry) { return !retained.contains(entry.first); });
const auto newest = frame(collection, settings, backend);
if (newest == nullptr) {
return {};
return {status_payload(2.0F, last_error)};
}
// Aperture: the newest 2r+1 collections, centred r collections before
// the newest one, so the image trails the input by r frames. Collections
// that left the aperture are dropped from the cache.
std::vector<std::shared_ptr<const Frame>> aperture{};
std::map<SweepKey, std::shared_ptr<const Frame>> kept{{sweep_key(collection), newest}};
for (const auto& member : previous.last(std::min(previous.size(), 2U * settings.frame_radius))) {
auto member_frame = frame(member, settings, gpu);
kept.emplace(sweep_key(member), member_frame);
for (const auto& member : recent) {
auto member_frame = frame(member, settings, backend);
const bool same_axis = member_frame != nullptr && member_frame->samples == newest->samples &&
member_frame->dt_s == newest->dt_s;
if (same_axis) {
@@ -590,14 +651,16 @@ struct Sar3dProcessor::State {
}
}
aperture.push_back(newest);
frames = std::move(kept);
const std::size_t centre_index =
aperture.size() > settings.frame_radius ? aperture.size() - 1U - settings.frame_radius : 0U;
const Frame& centre = *aperture[centre_index];
std::vector<double> offsets_s{};
for (const auto& member : aperture) {
offsets_s.push_back(member->center_s - centre.center_s);
const double delta_ns = member->center_ns >= centre.center_ns
? static_cast<double>(member->center_ns - centre.center_ns)
: -static_cast<double>(centre.center_ns - member->center_ns);
offsets_s.push_back(1e-9 * (delta_ns + member->center_fraction_ns - centre.center_fraction_ns));
}
const auto weights = frame_weights(offsets_s);
const Grid grid = scene_grid(settings);
@@ -637,8 +700,9 @@ struct Sar3dProcessor::State {
// A Tx+Rx path is convex, so the box corners bound it: samples beyond the
// longest path are never read and are not sent to the GPU.
const double samples_per_metre = 1.0 / (kWaveSpeed * newest->dt_s);
const std::size_t samples =
std::min(newest->samples, static_cast<std::size_t>(std::ceil(longest_path * samples_per_metre)) + 2U);
const double needed = std::ceil(longest_path * samples_per_metre) + 2.0;
const std::size_t samples = needed >= newest->samples
? newest->samples : static_cast<std::size_t>(needed);
std::vector<std::complex<float>> ascans{};
for (const auto& member : aperture) {
const std::span<const std::complex<float>> rows = member->ascans;
@@ -649,7 +713,7 @@ struct Sar3dProcessor::State {
}
const Oblique oblique = oblique_view(grid, settings.voxel_m);
auto images = gpu.project({
auto images = backend.project({
.x = grid.x,
.y = grid.y,
.z = grid.z,
@@ -664,6 +728,7 @@ struct Sar3dProcessor::State {
.use_compensation = settings.use_compensation,
.comp_range_power = settings.comp_range_power,
.comp_weight_max = kCompRangeWeightMax,
.norm_edge_margin_m = kNormEdgeMargin,
.oblique = oblique.view,
});
@@ -671,30 +736,17 @@ struct Sar3dProcessor::State {
// 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
// published projections peak at 1.
const std::size_t nx = grid.x.size();
const std::size_t nz = grid.z.size();
float level = 0.0F;
std::size_t peak_x = 0U;
std::size_t peak_y = 0U;
for (std::size_t iy = 0U; iy < grid.y.size(); ++iy) {
if (grid.y[iy] < grid.y.front() + kNormEdgeMargin || grid.y[iy] > grid.y.back() - kNormEdgeMargin) {
continue;
}
for (std::size_t ix = 0U; ix < nx; ++ix) {
if (images.front[iy * nx + ix] > level) {
level = images.front[iy * nx + ix];
peak_x = ix;
peak_y = iy;
}
}
}
const auto side_row = std::span<const float>(images.side).subspan(peak_y * nz, nz);
const auto peak_z = static_cast<std::size_t>(std::ranges::max_element(side_row) - side_row.begin());
const float scale = level > 0.0F ? 1.0F / level : 0.0F;
// 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;
// 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 *= scale;
value = static_cast<float>(static_cast<double>(value) / normalization);
if (!std::isfinite(value)) { throw std::runtime_error("SAR normalization exceeds the float32 range"); }
}
}
@@ -707,10 +759,12 @@ struct Sar3dProcessor::State {
basis.insert(basis.end(), oblique.view.v.begin(), oblique.view.v.end());
payloads.push_back(table_payload("sar3d_oblique_basis", 3U, std::move(basis)));
if (level > 0.0F) {
payloads.push_back(table_payload("sar3d_peak", 4U, {level, static_cast<float>(grid.x[peak_x]),
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])}));
}
payloads.push_back(status_payload(level > 0.0F ? 0.0F : 1.0F,
level > 0.0F ? "CUDA" : "No signal inside the selected depth range"));
return payloads;
}
@@ -748,10 +802,13 @@ auto Sar3dProcessor::process_collection(
try {
results.collection_payloads =
state_->render(collection, previous_collections, read_settings(run_config, live_config));
} catch (const std::exception& error) {
// A data or CUDA error must not stop the data_processor loop: publish an
// empty result and keep going.
} catch (const sar3d::CudaUnavailable& error) {
state_->report_error(error.what());
results.collection_payloads = {status_payload(3.0F, error.what())};
} catch (const std::exception& error) {
// Keep the processing loop alive and make the reason visible in the GUI.
state_->report_error(error.what());
results.collection_payloads = {status_payload(2.0F, error.what())};
}
return results;
}