radar_system/data_acq_and_processing/tests/sar3d_kernel_emulation.cpp

// Test-only execution of the exact CUDA kernel source, one logical thread at a time.
// Linked solely into sar3d_math_fixture; never used as a production backend.
#include "../processing/processors/src/sar3d_gpu.hpp"

#include <algorithm>
#include <bit>
#include <climits>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <string>

#define __global__
using std::isfinite;
using std::min;
struct Index { unsigned x; };
Index blockIdx{0}, blockDim{1}, threadIdx{0};
struct float2 { float x, y; };
auto make_float2(float x, float y) -> float2 { return {x, y}; }
void sincospif(float x, float* sine, float* cosine) {
    *sine = std::sin(x * 3.14159265358979323846F);
    *cosine = std::cos(x * 3.14159265358979323846F);
}
auto __powf(float x, float y) -> float { return std::pow(x, y); }
auto __float_as_uint(float x) -> unsigned int { return std::bit_cast<unsigned int>(x); }
auto __double2int_rn(double x) -> int { return static_cast<int>(std::nearbyint(x)); }
void atomicMax(unsigned int* destination, unsigned int value) { *destination = std::max(*destination, value); }
#include "sar3d_kernels.inc"

auto probe_ascan_upper_bound() -> float {
    blockIdx.x = 0;
    const double x[]{1073741824.0}, zero[]{0.0};
    const ::Contribution contribution{{}, {}, 1.0F, 1.0F};
    float volume = -1.0F;
    ::backproject(nullptr, INT_MAX, &contribution, 1, x, zero, zero, 1, 1, 1,
                  1.0, 0.0, 0.0, 0.0, 1.0, 1.0F, 0, 0.0F, 0.0F, &volume);
    return volume;
}

namespace {
template <class Operation>
void each(int count, Operation operation) {
    for (blockIdx.x = 0; blockIdx.x < static_cast<unsigned>(count); ++blockIdx.x) { operation(); }
}

auto float_pairs(std::span<const std::complex<float>> values) -> std::vector<float2> {
    std::vector<float2> pairs;
    pairs.reserve(values.size());
    for (const auto value : values) { pairs.push_back({value.real(), value.imag()}); }
    return pairs;
}
}

namespace radar::processing::sar3d {
namespace {
#include "sar3d_helpers.inc"
}

struct Engine::Impl { std::vector<float> taps; };
Engine::Engine(double sigma) : impl_(std::make_unique<Impl>()) { impl_->taps = gaussian_taps(sigma); }
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");
    }
    const auto count = checked_product(bands.size() / points, size, kMaxIndexedItems, "A-scan stack");
    const auto input = float_pairs(bands);
    std::vector<float2> output(count);
    each(static_cast<int>(count), [&] {
        ::inverse_dft(input.data(), points, first_bin, size, count, output.data());
    });
    std::vector<std::complex<float>> result;
    result.reserve(count);
    for (const auto value : output) { result.emplace_back(value.x, value.y); }
    return result;
}

auto Engine::project(const VolumeRequest& request) -> Projections {
    const int nx = request.x.size(), ny = request.y.size(), nz = request.z.size();
    const int plane = checked_product(nx, nz, kMaxIndexedItems, "grid plane");
    const int voxels = checked_product(plane, ny, kMaxIndexedItems, "grid");
    std::vector<::Contribution> contributions;
    for (const auto& term : request.contributions) {
        contributions.push_back(std::bit_cast<::Contribution>(term));
    }
    const auto ascans = float_pairs(request.ascans);
    std::vector<float> volume(voxels), scratch(voxels), levels(ny);
    std::vector<unsigned int> indices(ny);
    each(voxels, [&] {
        ::backproject(ascans.data(), request.samples_per_ascan, contributions.data(), contributions.size(),
            request.x.data(), request.y.data(), request.z.data(), 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, request.samples_per_metre,
            request.use_compensation, request.comp_range_power, request.comp_weight_max, volume.data());
    });
    each(ny, [&] { ::peak_planes(volume.data(), plane, ny, levels.data(), indices.data()); });
    Projections raw;
    choose_peak(raw, request, levels);
    const int radius = impl_->taps.size() / 2;
    each(voxels, [&] {
        ::smooth(volume.data(), scratch.data(), impl_->taps.data(), radius, nx, nz, voxels, 1);
    });
    each(voxels, [&] {
        ::smooth(scratch.data(), volume.data(), impl_->taps.data(), radius, nx, nz, voxels, 0);
    });
    each(ny, [&] { ::peak_planes(volume.data(), plane, ny, levels.data(), indices.data()); });
    Projections result;
    result.raw_level = raw.level;
    choose_peak(result, request, levels, indices);
    result.front.resize(ny * nx);
    result.side.resize(ny * nz);
    result.top.resize(plane);
    const auto& view = request.oblique;
    std::vector<unsigned int> oblique(checked_product(view.nu, view.nv, kMaxIndexedItems, "oblique image"));
    each(result.front.size(), [&] {
        ::max_project(volume.data(), result.front.data(), result.front.size(), nx, plane, 1, nz, nx);
    });
    each(result.side.size(), [&] {
        ::max_project(volume.data(), result.side.data(), result.side.size(), nz, plane, nx, nx, 1);
    });
    each(result.top.size(), [&] {
        ::max_project(volume.data(), result.top.data(), result.top.size(), nx, nx, 1, ny, plane);
    });
    each(voxels, [&] {
        ::project_oblique(volume.data(), request.x.data(), request.y.data(), request.z.data(), 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 / view.pixel_m, view.nu, view.nv, oblique.data());
    });
    result.oblique.reserve(oblique.size());
    for (const auto value : oblique) { result.oblique.push_back(std::bit_cast<float>(value)); }
    return result;
}
}