sar3d: GPU port of the notebook with four views

The sar3d processor is rewritten as a checked port of Multi_antenna_3D_clean_v1.ipynb:
NVRTC kernels, timing from the trace stamps, notebook parameter names and one cubic
SAR_VOXEL_M. It adds front, side, top and 3D views, drops gpr_cpp and the parity tools,
and the GUI speed now applies only while the socket speed is ignored.

data_acq_and_processing/processing/processors/include/sar3d_processor.hpp

@@ -8,9 +8,11 @@
 
 namespace radar::processing {
 
-// 3D SAR over a synthesized aperture, computed by the external gpr_cpp
-// CUDA/NVRTC library (VNA_AlGo). Emits 2D max-projections of the volume
-// (side XZ and top XY views) as ImageF32 collection payloads.
+// 3D SAR of the reference notebook Multi_antenna_3D_clean_v1.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.
 class Sar3dProcessor final : public ProcessorInterface {
   public:
     Sar3dProcessor();
@@ -25,8 +27,6 @@ class Sar3dProcessor final : public ProcessorInterface {
     ) -> ipc::ResultCollection override;
 
   private:
-    // Pimpl keeps gpr/CUDA headers (cuda.h, nvrtc.h) out of the shared
-    // processor include path; only sar3d_processor.cpp needs them.
     struct State;
     std::unique_ptr<State> state_;
 };

data_acq_and_processing/processing/processors/src/sar3d_gpu.cpp

@@ -0,0 +1,471 @@
+#include "sar3d_gpu.hpp"
+
+#include <stdexcept>
+
+#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 {
+
+#ifdef SAR3D_WITH_CUDA
+namespace {
+
+// Compiled at run time by NVRTC for the installed GPU, so the build needs no
+// nvcc. Arithmetic is float32 apart from a few per-voxel terms: FP64 runs at
+// 1/32-1/64 of the FP32 rate on consumer GPUs, and float32 already keeps the
+// A-scan sample position to ~1e-3 of a sample of the 8x oversampled A-scan.
+// Volumes are [y][z][x] with x fastest; their sizes stay below 2^31 (checked on
+// the host), so the hot index arithmetic is 32-bit.
+constexpr const char* kKernelSource = R"CUDA(
+struct Contribution {
+    float tx[3];
+    float rx[3];
+    float weight;
+    float inv_reference_ranges;
+};
+
+// numpy.fft.ifft of a band spectrum placed at bins first_bin.. of a size-N
+// spectrum (N a power of two), summed directly: the band is a few hundred
+// bins of thousands. One thread per (pair, sample); the phase index is exact.
+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 k = (unsigned)(item % size);
+    const unsigned mask = (unsigned)size - 1u;
+    const float2* band = bands + (item / size) * points;
+    unsigned phase = (unsigned)(((unsigned long long)first_bin * k) & mask);
+    const float turns_per_index = 2.0f / (float)size;
+    float re = 0.0f, im = 0.0f;
+    for (int n = 0; n < points; ++n) {
+        float s, c;
+        sincospif((float)phase * turns_per_index, &s, &c);
+        const float2 b = band[n];
+        re += b.x * c - b.y * s;
+        im += b.x * s + b.y * c;
+        phase = (phase + k) & mask;
+    }
+    ascans[item] = make_float2(re / (float)size, im / (float)size);
+}
+
+// Coherent backprojection, one thread per voxel: |sum w*comp*h| / sum w over the
+// valid (frame, pair) terms (notebook sar_frame_contribution and the first line
+// of sar_finish_volume).
+extern "C" __global__ void backproject(
+    const float2* __restrict__ ascans, const int samples,
+    const Contribution* __restrict__ contributions, const int count,
+    const double* __restrict__ gx, const double* __restrict__ gy, const double* __restrict__ gz,
+    const int nx, const int nz, const int voxels,
+    const double cos_p, const double sin_p, const double ground_depth,
+    const double min_depth, const double max_depth,
+    const float samples_per_metre, const int use_compensation,
+    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 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.
+    const double y0 = gy[row / nz] - ground_depth;
+    const double zs = gz[row % nz];
+    const double depth = -y0 * sin_p + zs * cos_p;
+    if (depth < min_depth || depth > max_depth) {
+        volume[voxel] = 0.0f;
+        return;
+    }
+    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);
+    float re = 0.0f, im = 0.0f, weights = 0.0f;
+    for (int c = 0; c < count; ++c) {
+        const Contribution k = contributions[c];
+        const float dxt = px - k.tx[0], dyt = py - k.tx[1], dzt = pz - k.tx[2];
+        const float dxr = px - k.rx[0], dyr = py - k.rx[1], dzr = pz - k.rx[2];
+        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;
+        const int i = min((int)position, samples - 2);
+        const float frac = position - (float)i;
+        const float2* row_ascan = ascans + (long long)c * samples;
+        const float2 a = row_ascan[i];
+        const float2 b = row_ascan[i + 1];
+        float w = k.weight;
+        if (use_compensation) {
+            // sar_compensation_weight: 1 / (geo / geo_ref)^power, clipped.
+            w *= fminf(__powf(rt * rr * k.inv_reference_ranges, comp_power), comp_max);
+        }
+        re += w * (a.x + frac * (b.x - a.x));
+        im += w * (a.y + frac * (b.y - a.y));
+        weights += k.weight;
+    }
+    volume[voxel] = sqrtf(re * re + im * im) / (weights + 1e-12f);
+}
+
+// One pass of scipy.ndimage.gaussian_filter1d (mode 'reflect') along Z or X.
+extern "C" __global__ void smooth(
+    const float* __restrict__ src, float* __restrict__ dst,
+    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 int length = along_z ? nz : nx;
+    const int index = along_z ? (voxel / nx) % nz : voxel % nx;
+    const int stride = along_z ? nx : 1;
+    float sum = 0.0f;
+    for (int k = -radius; k <= radius; ++k) {
+        int j = index + k;
+        while (j < 0 || j >= length) j = j < 0 ? -j - 1 : 2 * length - j - 1;   // d c b a | a b c d | d c b a
+        sum += taps[k + radius] * src[voxel + (j - index) * stride];
+    }
+    dst[voxel] = sum;
+}
+
+// dst[o] = max over r of src[(o / inner) * outer_stride + (o % inner) * inner_stride + r * reduce_stride].
+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;
+}
+
+// Orthographic max projection. The volume is non-negative, so float bit
+// patterns order like the values and an integer atomicMax does the job. The
+// pixel is found in double, so it does not hinge on float rounding.
+extern "C" __global__ void project_oblique(
+    const float* __restrict__ volume,
+    const double* __restrict__ gx, const double* __restrict__ gy, const double* __restrict__ gz,
+    const int nx, const int nz, const int voxels,
+    const double ux, const double uy, const double uz, const double vx, const double vy, const double vz,
+    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 float value = volume[voxel];
+    if (!(value > 0.0f)) return;
+    const int row = voxel / nx;
+    const double x = gx[voxel - row * nx];
+    const double y = gy[row / nz];
+    const double z = gz[row % nz];
+    const int pu = __double2int_rn((ux * x + uy * y + uz * z - u0) * inv_pixel);
+    const int pv = __double2int_rn((vx * x + vy * y + vz * z - v0) * inv_pixel);
+    if (pu < 0 || pu >= nu || pv < 0 || pv >= nv) return;
+    atomicMax(image + pv * nu + pu, __float_as_uint(value));
+}
+)CUDA";
+
+constexpr int kBlockSize = 256;
+
+void check(CUresult result, const char* what) {
+    if (result != CUDA_SUCCESS) {
+        const char* text = nullptr;
+        cuGetErrorString(result, &text);
+        throw std::runtime_error(std::string{"CUDA "} + what + ": " + (text != nullptr ? text : "unknown error"));
+    }
+}
+
+void check(nvrtcResult result, const char* what) {
+    if (result != NVRTC_SUCCESS) {
+        throw std::runtime_error(std::string{"NVRTC "} + what + ": " + nvrtcGetErrorString(result));
+    }
+}
+
+[[nodiscard]] auto checked_int(std::size_t value, const char* what) -> int {
+    if (value > static_cast<std::size_t>(INT_MAX)) {
+        throw std::runtime_error(std::string{"sar3d "} + what + " exceeds the 32-bit kernel indexing");
+    }
+    return static_cast<int>(value);
+}
+
+// Native code when this NVRTC knows the GPU's architecture, otherwise PTX for
+// the newest architecture it knows, which the driver compiles on load.
+[[nodiscard]] auto compile_kernels(CUdevice device) -> std::string {
+    int major = 0;
+    int minor = 0;
+    check(cuDeviceGetAttribute(&major, CU_DEVICE_ATTRIBUTE_COMPUTE_CAPABILITY_MAJOR, device), "compute capability");
+    check(cuDeviceGetAttribute(&minor, CU_DEVICE_ATTRIBUTE_COMPUTE_CAPABILITY_MINOR, device), "compute capability");
+    const int device_arch = 10 * major + minor;
+    int arch_count = 0;
+    check(nvrtcGetNumSupportedArchs(&arch_count), "supported architectures");
+    std::vector<int> archs(static_cast<std::size_t>(arch_count));
+    check(nvrtcGetSupportedArchs(archs.data()), "supported architectures");
+    int target = 0;
+    for (const int arch : archs) {
+        if (arch <= device_arch) {
+            target = std::max(target, arch);
+        }
+    }
+    if (target == 0) {
+        throw std::runtime_error(
+            "GPU sm_" + std::to_string(device_arch) + " is older than every architecture this NVRTC supports"
+        );
+    }
+    const bool native = target == device_arch;
+    const std::string arch_option =
+        std::string{native ? "--gpu-architecture=sm_" : "--gpu-architecture=compute_"} + std::to_string(target);
+
+    nvrtcProgram raw_program = nullptr;
+    check(nvrtcCreateProgram(&raw_program, kKernelSource, "sar3d.cu", 0, nullptr, nullptr), "create program");
+    const auto destroy = [](nvrtcProgram program) { nvrtcDestroyProgram(&program); };
+    const std::unique_ptr<std::remove_pointer_t<nvrtcProgram>, decltype(destroy)> program(raw_program, destroy);
+
+    const char* options[] = {arch_option.c_str()};
+    if (nvrtcCompileProgram(program.get(), 1, options) != NVRTC_SUCCESS) {
+        std::size_t log_size = 0U;
+        nvrtcGetProgramLogSize(program.get(), &log_size);
+        std::string log(log_size, '\0');
+        nvrtcGetProgramLog(program.get(), log.data());
+        throw std::runtime_error("NVRTC failed to compile the sar3d kernels:\n" + log);
+    }
+    const auto image_size = native ? nvrtcGetCUBINSize : nvrtcGetPTXSize;
+    const auto image_data = native ? nvrtcGetCUBIN : nvrtcGetPTX;
+    std::size_t size = 0U;
+    check(image_size(program.get(), &size), "kernel image size");
+    std::string image(size, '\0');
+    check(image_data(program.get(), image.data()), "kernel image");
+    return image;
+}
+
+// Device allocation that only grows, so steady-state frames never reallocate.
+class DeviceBuffer {
+  public:
+    DeviceBuffer() = default;
+    DeviceBuffer(const DeviceBuffer&) = delete;
+    auto operator=(const DeviceBuffer&) -> DeviceBuffer& = delete;
+    ~DeviceBuffer() {
+        if (pointer_ != 0U) {
+            cuMemFree(pointer_);
+        }
+    }
+
+    auto reserve(std::size_t bytes) -> CUdeviceptr {
+        if (bytes > capacity_) {
+            if (pointer_ != 0U) {
+                check(cuMemFree(pointer_), "free");
+                pointer_ = 0U;
+                capacity_ = 0U;
+            }
+            check(cuMemAlloc(&pointer_, bytes), "allocate");
+            capacity_ = bytes;
+        }
+        return pointer_;
+    }
+
+    template <class T>
+    auto upload(std::span<const T> values) -> CUdeviceptr {
+        reserve(values.size_bytes());
+        if (!values.empty()) {
+            check(cuMemcpyHtoD(pointer_, values.data(), values.size_bytes()), "upload");
+        }
+        return pointer_;
+    }
+
+    [[nodiscard]] auto pointer() const -> CUdeviceptr {
+        return pointer_;
+    }
+
+  private:
+    CUdeviceptr pointer_ = 0U;
+    std::size_t capacity_ = 0U;
+};
+
+template <class T>
+[[nodiscard]] auto download(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.
+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);
+    check(cuLaunchKernel(kernel, blocks, 1U, 1U, kBlockSize, 1U, 1U, 0U, nullptr, parameters, nullptr), "launch");
+}
+
+}  // namespace
+
+struct Gpu::Impl {
+    // 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(cuDevicePrimaryCtxRetain(&handle, device), "context");
+            if (const CUresult result = cuCtxSetCurrent(handle); result != CUDA_SUCCESS) {
+                cuDevicePrimaryCtxRelease(device);
+                check(result, "context");
+            }
+        }
+        ~Context() { cuDevicePrimaryCtxRelease(device); }
+    } context;
+    struct Module {
+        CUmodule handle = nullptr;
+        explicit Module(const std::string& image) { check(cuModuleLoadData(&handle, image.data()), "load kernels"); }
+        ~Module() { cuModuleUnload(handle); }
+        [[nodiscard]] auto function(const char* name) const -> CUfunction {
+            CUfunction kernel = nullptr;
+            check(cuModuleGetFunction(&kernel, handle, name), name);
+            return kernel;
+        }
+    } module{compile_kernels(context.device)};
+    CUfunction inverse_dft = module.function("inverse_dft");
+    CUfunction backproject = module.function("backproject");
+    CUfunction smooth = module.function("smooth");
+    CUfunction max_project = module.function("max_project");
+    CUfunction project_oblique = module.function("project_oblique");
+    int smooth_radius = 0;
+    DeviceBuffer taps;
+    DeviceBuffer bands;
+    DeviceBuffer axes;
+    DeviceBuffer contributions;
+    DeviceBuffer ascans;
+    DeviceBuffer volume;
+    DeviceBuffer scratch;
+    DeviceBuffer images;
+
+    ~Impl() { cuCtxSetCurrent(context.handle); }
+};
+
+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(
+    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_;
+    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"),
+           checked_int(size, "A-scan"), items, output);
+    return download<std::complex<float>>(output, samples);
+}
+
+auto Gpu::project(const VolumeRequest& request) -> Projections {
+    Impl& gpu = *impl_;
+    check(cuCtxSetCurrent(gpu.context.handle), "context");
+
+    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");
+
+    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));
+
+    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);
+    // 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);
+
+    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);
+    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;
+    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 {};
+}
+
+#endif
+
+}  // namespace radar::processing::sar3d

data_acq_and_processing/processing/processors/src/sar3d_gpu.hpp

@@ -0,0 +1,91 @@
+#pragma once
+
+#include <array>
+#include <complex>
+#include <cstddef>
+#include <memory>
+#include <span>
+#include <vector>
+
+// GPU half 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.
+namespace radar::processing::sar3d {
+
+// One (aperture frame, Tx/Rx pair) term of the backprojection sum. The layout
+// is shared with the CUDA kernel.
+struct Contribution {
+    std::array<float, 3> tx{};             // Tx in radar axes, shifted by the frame offset [m]
+    std::array<float, 3> rx{};             // Rx in radar axes, shifted by the frame offset [m]
+    float weight = 0.0F;                   // aperture frame weight
+    float inv_reference_ranges = 0.0F;     // 1 / (Rtx * Rrx) at the compensation reference point
+};
+static_assert(sizeof(Contribution) == 8U * sizeof(float));
+
+// Orthographic view of the volume: a voxel at scene point p lands in pixel
+// (round((u.p - u0) / pixel), round((v.p - v0) / pixel)).
+struct ObliqueView {
+    std::array<double, 3> u{};   // screen right, scene axes (X, Y, Z)
+    std::array<double, 3> v{};   // screen up
+    double u0 = 0.0;
+    double v0 = 0.0;
+    double pixel_m = 0.0;
+    int nu = 0;
+    int nv = 0;
+};
+
+struct VolumeRequest {
+    // Scene axes [m]: X across, Y height above the ground, Z along the ground.
+    std::span<const double> x;
+    std::span<const double> y;
+    std::span<const double> z;
+    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;
+    // 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::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;
+    ObliqueView oblique{};
+};
+
+// Max-projections of the smoothed, not yet normalized volume (row-major).
+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]
+};
+
+// 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 {
+  public:
+    explicit Gpu(double smooth_sigma_voxels);
+    ~Gpu();
+    Gpu(const Gpu&) = delete;
+    auto operator=(const Gpu&) -> Gpu& = 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.
+    [[nodiscard]] auto 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>>;
+
+    // Backprojects, applies the Gaussian smoothing over (Z, X) and projects.
+    [[nodiscard]] auto project(const VolumeRequest& request) -> Projections;
+
+  private:
+    struct Impl;
+    std::unique_ptr<Impl> impl_;
+};
+
+}  // namespace radar::processing::sar3d