constexpr double kPi = 3.14159265358979323846;
constexpr double kSpeedOfLightMetersPerSec = 299'792'458.0;
constexpr double kXMarginM = 2.0;
constexpr double kGridZMinM = 0.1;
constexpr std::size_t kGridWidth = 300U;
constexpr std::size_t kGridHeight = 300U;
constexpr std::size_t kAscanOversample = 8U;

constexpr double kRangeWeightMax = 5.0;
constexpr double kAngleWeightMax = 2.0;
constexpr double kTotalWeightMax = 8.0;
constexpr double kCompensationReferenceDepthM = 5.0;

constexpr bool kPairNormalize = true;
constexpr double kPairNormPercentile = 50.0;
constexpr double kPairNormEps = 1e-15;

constexpr double kSmoothSigma = 1.5;
constexpr std::size_t kMaxObjects = 10U;
constexpr double kObjectMinFrac = 0.35;
constexpr double kRegionThresholdFrac = 0.75;
constexpr double kSuppressThresholdFrac = 0.20;
constexpr double kSuppressRadiusXM = 0.80;
constexpr double kSuppressRadiusZM = 0.40;
constexpr double kMinRegionAreaCm2 = 10.0;

constexpr double kCenterRadiusXM = 0.60;
constexpr double kCenterRadiusZM = 0.25;
constexpr double kCenterThresholdFrac = 0.88;
constexpr double kCenterWeightPower = 2.0;

constexpr double kSidelobeRangeRmsToleranceM = 0.20;
constexpr double kSidelobeMinDxM = 0.35;
constexpr double kSidelobeMaxDzM = 0.70;
constexpr double kSidelobeMaxRelativePeak = 0.85;

constexpr double kLocalBgRadiusXM = 1.20;
constexpr double kLocalBgRadiusZM = 0.80;
constexpr double kLocalBgPercentile = 50.0;
constexpr double kLocalContrastEps = 1e-12;

constexpr double kScoreCohPeakWeight = 0.45;
constexpr double kScoreCoherenceFactorWeight = 0.25;
constexpr double kScoreProminenceWeight = 0.20;
constexpr double kScoreContrastWeight = 0.10;
constexpr double kScoreContrastCap = 6.0;
constexpr double kScoreCfEps = 1e-12;

using PairKey = std::uint64_t;

struct GeometrySelection {
    // Per-local-index Tx/Rx antenna coordinates in metres. y_/z_ default to 0
    // for legacy configs so the imaging plane coincides with the antennas.
    std::vector<std::uint32_t> input_positions{};
    std::vector<std::uint32_t> output_positions{};
    std::vector<double> x_tx{};
    std::vector<double> y_tx{};
    std::vector<double> z_tx{};
    std::vector<double> x_rx{};
    std::vector<double> y_rx{};
    std::vector<double> z_rx{};
    std::unordered_map<std::uint32_t, std::uint32_t> input_local_by_pos{};
    std::unordered_map<std::uint32_t, std::uint32_t> output_local_by_pos{};
};

struct BackgroundAccumulator {
    std::vector<std::complex<double>> sum{};
    std::size_t count = 0U;
};

struct SelectedTrace {
    ipc::ComboKey combo{};
    std::uint32_t tx_local_index = 0U;
    std::uint32_t rx_local_index = 0U;
    std::size_t run_order = 0U;
    std::vector<double> frequency_hz{};
    std::vector<std::complex<double>> s21{};
};

struct AscanResult {
    std::vector<double> time_s{};
    std::vector<std::complex<double>> samples{};
    double dt_s = 0.0;
    double bandwidth_hz = 0.0;
};

struct GridDefinition {
    std::vector<double> x_grid{};
    std::vector<double> z_grid{};
    std::vector<std::vector<double>> tx_distance_grids{};
    std::vector<std::vector<double>> rx_distance_grids{};
};

struct BpMap {
    std::vector<double> image{};
    std::vector<std::complex<double>> coherent{};
    std::vector<double> incoherent{};
    std::vector<double> coherence_factor{};
};

struct ObjectRecord {
    std::size_t index = 0U;
    std::size_t peak_index = 0U;
    double x_peak_m = 0.0;
    double z_peak_m = 0.0;
    double x_m = 0.0;
    double z_m = 0.0;
    double x_region_m = 0.0;
    double z_region_m = 0.0;
    double peak = 0.0;
    double area_cm2 = 0.0;
    double center_area_cm2 = 0.0;
    double mean_value = 0.0;
    double sum_value = 0.0;
    double local_bg = 0.0;
    double local_bg_p75 = 0.0;
    double prominence = 0.0;
    double contrast = 0.0;
    double incoh_peak = 0.0;
    double incoh_mean = 0.0;
    double incoh_center_mean = 0.0;
    double coherence_factor_peak = 0.0;
    double coherence_factor_center = 0.0;
    double score_old = 0.0;
    double score_new = 0.0;
    double selected_score = 0.0;
    std::vector<std::uint8_t> region_mask{};
    std::vector<std::uint8_t> center_mask{};
    bool sidelobe_candidate = false;
    std::size_t sidelobe_parent = 0U;
    double sidelobe_range_rms_m = std::numeric_limits<double>::quiet_NaN();
};

[[nodiscard]] auto make_pair_key(std::uint32_t tx_local_index, std::uint32_t rx_local_index) -> PairKey {
    return (static_cast<PairKey>(tx_local_index) << 32U) | static_cast<PairKey>(rx_local_index);
}

[[nodiscard]] auto combo_to_string(const ipc::ComboKey& combo) -> std::string {
    return "input=" + std::to_string(combo.input_pos) + " output=" + std::to_string(combo.output_pos);
}

[[nodiscard]] auto next_power_of_two(std::size_t value) -> std::size_t {
    if (value <= 1U) {
        return 1U;
    }

    std::size_t power = 1U;
    while (power < value) {
        if (power > (std::numeric_limits<std::size_t>::max() >> 1U)) {
            return value;
        }
        power <<= 1U;
    }
    return power;
}

void fft_inplace(std::vector<std::complex<double>>& values, bool inverse) {
    const std::size_t size = values.size();
    if (size <= 1U) {
        return;
    }

    for (std::size_t index = 1U, bit_reversed = 0U; index < size; ++index) {
        std::size_t bit = size >> 1U;
        while (bit_reversed & bit) {
            bit_reversed ^= bit;
            bit >>= 1U;
        }
        bit_reversed ^= bit;
        if (index < bit_reversed) {
            std::swap(values[index], values[bit_reversed]);
        }
    }

    for (std::size_t len = 2U; len <= size; len <<= 1U) {
        const double angle = 2.0 * kPi * (inverse ? 1.0 : -1.0) / static_cast<double>(len);
        const std::complex<double> twiddle_step(std::cos(angle), std::sin(angle));
        const std::size_t half_len = len >> 1U;

        for (std::size_t offset = 0U; offset < size; offset += len) {
            std::complex<double> twiddle(1.0, 0.0);
            for (std::size_t inner = 0U; inner < half_len; ++inner) {
                const auto even = values[offset + inner];
                const auto odd = values[offset + inner + half_len] * twiddle;
                values[offset + inner] = even + odd;
                values[offset + inner + half_len] = even - odd;
                twiddle *= twiddle_step;
            }
        }
    }

    if (!inverse) {
        return;
    }

    const double scale = 1.0 / static_cast<double>(size);
    for (auto& value : values) {
        value *= scale;
    }
}

[[nodiscard]] auto median_copy(std::vector<double> values) -> double {
    if (values.empty()) {
        return 0.0;
    }

    const auto middle = values.begin() + static_cast<std::ptrdiff_t>(values.size() / 2U);
    std::nth_element(values.begin(), middle, values.end());
    double median = *middle;
    if ((values.size() % 2U) == 0U) {
        const auto lower_middle = values.begin() + static_cast<std::ptrdiff_t>((values.size() / 2U) - 1U);
        std::nth_element(values.begin(), lower_middle, values.end());
        median = 0.5 * (median + *lower_middle);
    }
    return median;
}

[[nodiscard]] auto percentile_copy(std::vector<double> values, double percentile) -> double {
    if (values.empty()) {
        return 0.0;
    }

    values.erase(
        std::remove_if(values.begin(), values.end(), [](double value) { return !std::isfinite(value); }),
        values.end()
    );
    if (values.empty()) {
        return 0.0;
    }

    std::sort(values.begin(), values.end());
    const double clamped_percentile = std::clamp(percentile, 0.0, 100.0);
    const double position = (clamped_percentile / 100.0) * static_cast<double>(values.size() - 1U);
    const auto lower_index = static_cast<std::size_t>(std::floor(position));
    const auto upper_index = std::min<std::size_t>(lower_index + 1U, values.size() - 1U);
    const double fraction = position - static_cast<double>(lower_index);
    return (values[lower_index] * (1.0 - fraction)) + (values[upper_index] * fraction);
}

[[nodiscard]] auto build_axis(double min_value, double max_value, std::size_t count) -> std::vector<double> {
    std::vector<double> axis{};
    if (count == 0U) {
        return axis;
    }

    axis.resize(count, min_value);
    if (count == 1U) {
        return axis;
    }

    const double step = (max_value - min_value) / static_cast<double>(count - 1U);
    for (std::size_t index = 0U; index < count; ++index) {
        axis[index] = min_value + (step * static_cast<double>(index));
    }
    return axis;
}

[[nodiscard]] auto clamp_index(std::ptrdiff_t value, std::size_t limit) -> std::size_t {
    if (limit == 0U) {
        return 0U;
    }
    if (value < 0) {
        return 0U;
    }
    const auto max_index = static_cast<std::ptrdiff_t>(limit - 1U);
    if (value > max_index) {
        return limit - 1U;
    }
    return static_cast<std::size_t>(value);
}

[[nodiscard]] auto max_value(const std::vector<double>& values) -> double {
    if (values.empty()) {
        return 0.0;
    }
    return *std::max_element(values.begin(), values.end());
}

void normalize_in_place(std::vector<double>& values) {
    const double maximum = max_value(values);
    if (!(maximum > 0.0)) {
        return;
    }
    for (auto& value : values) {
        value /= maximum;
    }
}

[[nodiscard]] auto build_gaussian_kernel(double sigma) -> std::vector<double> {
    if (!(sigma > 0.0)) {
        return {1.0};
    }

    const auto radius = static_cast<std::ptrdiff_t>(std::ceil(sigma * 3.0));
    std::vector<double> kernel(static_cast<std::size_t>((radius * 2) + 1), 0.0);
    double sum = 0.0;
    for (std::ptrdiff_t offset = -radius; offset <= radius; ++offset) {
        const double value = std::exp(-0.5 * std::pow(static_cast<double>(offset) / sigma, 2.0));
        kernel[static_cast<std::size_t>(offset + radius)] = value;
        sum += value;
    }
    if (sum > 0.0) {
        for (auto& value : kernel) {
            value /= sum;
        }
    }
    return kernel;
}

[[nodiscard]] auto gaussian_filter_2d(
    const std::vector<double>& values,
    std::size_t width,
    std::size_t height,
    double sigma
) -> std::vector<double> {
    if (values.empty() || width == 0U || height == 0U) {
        return {};
    }

    const auto kernel = build_gaussian_kernel(sigma);
    const auto radius = static_cast<std::ptrdiff_t>((kernel.size() - 1U) / 2U);
    std::vector<double> temp(values.size(), 0.0);
    std::vector<double> output(values.size(), 0.0);

    for (std::size_t row = 0U; row < height; ++row) {
        for (std::size_t col = 0U; col < width; ++col) {
            double sum = 0.0;
            for (std::ptrdiff_t offset = -radius; offset <= radius; ++offset) {
                const auto sample_col = clamp_index(static_cast<std::ptrdiff_t>(col) + offset, width);
                sum += values[(row * width) + sample_col] * kernel[static_cast<std::size_t>(offset + radius)];
            }
            temp[(row * width) + col] = sum;
        }
    }

    for (std::size_t row = 0U; row < height; ++row) {
        for (std::size_t col = 0U; col < width; ++col) {
            double sum = 0.0;
            for (std::ptrdiff_t offset = -radius; offset <= radius; ++offset) {
                const auto sample_row = clamp_index(static_cast<std::ptrdiff_t>(row) + offset, height);
                sum += temp[(sample_row * width) + col] * kernel[static_cast<std::size_t>(offset + radius)];
            }
            output[(row * width) + col] = sum;
        }
    }

    return output;
}

[[nodiscard]] auto contains_u32(const std::vector<std::uint32_t>& values, std::uint32_t value) -> bool {
    return std::find(values.begin(), values.end(), value) != values.end();
}

[[nodiscard]] auto selected_positions(
    const std::vector<std::uint32_t>& requested,
    std::vector<std::uint32_t> available
) -> std::vector<std::uint32_t> {
    std::sort(available.begin(), available.end());
    available.erase(std::unique(available.begin(), available.end()), available.end());
    if (requested.empty()) {
        return available;
    }

    std::vector<std::uint32_t> result{};
    for (const auto value : requested) {
        if (contains_u32(available, value) && !contains_u32(result, value)) {
            result.push_back(value);
        }
    }
    std::sort(result.begin(), result.end());
    return result;
}

struct AntennaXYZ {
    double x_m = 0.0;
    double y_m = 0.0;
    double z_m = 0.0;
};

[[nodiscard]] auto build_geometry_selection(
    const config::RunConfig& run_config,
    const ProcessingLiveConfig& live_config
) -> GeometrySelection {
    std::unordered_map<std::uint32_t, AntennaXYZ> tx_by_pos{};
    for (const auto& entry : run_config.gpr.tx_geometry) {
        tx_by_pos[entry.output_pos] = AntennaXYZ{
            static_cast<double>(entry.x_m),
            static_cast<double>(entry.y_m),
            static_cast<double>(entry.z_m),
        };
    }

    std::unordered_map<std::uint32_t, AntennaXYZ> rx_by_pos{};
    for (const auto& entry : run_config.gpr.rx_geometry) {
        rx_by_pos[entry.input_pos] = AntennaXYZ{
            static_cast<double>(entry.x_m),
            static_cast<double>(entry.y_m),
            static_cast<double>(entry.z_m),
        };
    }

    std::vector<std::uint32_t> available_outputs{};
    available_outputs.reserve(tx_by_pos.size());
    for (const auto& [position, _] : tx_by_pos) {
        available_outputs.push_back(position);
    }

    std::vector<std::uint32_t> available_inputs{};
    available_inputs.reserve(rx_by_pos.size());
    for (const auto& [position, _] : rx_by_pos) {
        available_inputs.push_back(position);
    }

    GeometrySelection selection{};
    selection.output_positions = selected_positions(live_config.gpr_output_positions, std::move(available_outputs));
    selection.input_positions = selected_positions(live_config.gpr_input_positions, std::move(available_inputs));

    const auto reserve_axes = [](GeometrySelection& target, std::size_t tx_count, std::size_t rx_count) {
        target.x_tx.reserve(tx_count);
        target.y_tx.reserve(tx_count);
        target.z_tx.reserve(tx_count);
        target.x_rx.reserve(rx_count);
        target.y_rx.reserve(rx_count);
        target.z_rx.reserve(rx_count);
    };
    reserve_axes(selection, selection.output_positions.size(), selection.input_positions.size());

    for (std::size_t index = 0U; index < selection.output_positions.size(); ++index) {
        const auto position = selection.output_positions[index];
        selection.output_local_by_pos[position] = static_cast<std::uint32_t>(index);
        const auto& xyz = tx_by_pos[position];
        selection.x_tx.push_back(xyz.x_m);
        selection.y_tx.push_back(xyz.y_m);
        selection.z_tx.push_back(xyz.z_m);
    }

    for (std::size_t index = 0U; index < selection.input_positions.size(); ++index) {
        const auto position = selection.input_positions[index];
        selection.input_local_by_pos[position] = static_cast<std::uint32_t>(index);
        const auto& xyz = rx_by_pos[position];
        selection.x_rx.push_back(xyz.x_m);
        selection.y_rx.push_back(xyz.y_m);
        selection.z_rx.push_back(xyz.z_m);
    }

    return selection;
}

void validate_collection_trace_order(
    const config::RunConfig& run_config,
    const ipc::PreprocessedCollection& collection
) {
    if (collection.traces.size() != run_config.run_combos.size()) {
        throw std::runtime_error(
            "GPR requires collection trace order to match run.combos exactly: trace_count=" +
            std::to_string(collection.traces.size()) +
            ", run_combo_count=" + std::to_string(run_config.run_combos.size())
        );
    }

    for (std::size_t index = 0U; index < collection.traces.size(); ++index) {
        const auto& actual = collection.traces[index].combo;
        const auto& expected = run_config.run_combos[index];
        if (actual.input_pos == expected.input_pos && actual.output_pos == expected.output_pos) {
            continue;
        }

        throw std::runtime_error(
            "GPR requires preprocessed trace order to match run.combos: index=" + std::to_string(index) +
            ", expected=(" + combo_to_string(expected) + "), actual=(" + combo_to_string(actual) + ")"
        );
    }
}

[[nodiscard]] auto build_background_mean(
    std::span<const ipc::PreprocessedCollection> previous_collections,
    const GeometrySelection& selection,
    const ProcessingLiveConfig& live_config
) -> std::unordered_map<PairKey, std::vector<std::complex<double>>> {
    std::unordered_map<PairKey, std::vector<std::complex<double>>> result{};
    if (!live_config.gpr_background_subtract_enabled || live_config.gpr_background_mean_count == 0U) {
        return result;
    }

    const std::size_t mean_count = static_cast<std::size_t>(live_config.gpr_background_mean_count);
    const std::size_t start_index =
        previous_collections.size() > mean_count ? previous_collections.size() - mean_count : 0U;

    std::unordered_map<PairKey, BackgroundAccumulator> accumulators{};
    for (std::size_t collection_index = start_index; collection_index < previous_collections.size(); ++collection_index) {
        const auto& collection = previous_collections[collection_index];
        for (const auto& trace : collection.traces) {
            const auto output_it = selection.output_local_by_pos.find(trace.combo.output_pos);
            const auto input_it = selection.input_local_by_pos.find(trace.combo.input_pos);
            if (output_it == selection.output_local_by_pos.end() || input_it == selection.input_local_by_pos.end()) {
                continue;
            }

            const auto key = make_pair_key(output_it->second, input_it->second);
            auto& accumulator = accumulators[key];
            if (accumulator.sum.empty()) {
                accumulator.sum.assign(trace.s21.size(), std::complex<double>(0.0, 0.0));
            }
            if (accumulator.sum.size() != trace.s21.size()) {
                continue;
            }

            for (std::size_t sample_index = 0U; sample_index < trace.s21.size(); ++sample_index) {
                const auto& sample = trace.s21[sample_index];
                accumulator.sum[sample_index] += std::complex<double>(sample.re, sample.im);
            }
            accumulator.count += 1U;
        }
    }

    for (auto& [key, accumulator] : accumulators) {
        if (accumulator.count == 0U) {
            continue;
        }

        auto& mean_trace = result[key];
        mean_trace = std::move(accumulator.sum);
        const double inverse_count = 1.0 / static_cast<double>(accumulator.count);
        for (auto& sample : mean_trace) {
            sample *= inverse_count;
        }
    }

    return result;
}

[[nodiscard]] auto collect_selected_traces(
    const ipc::PreprocessedCollection& collection,
    const GeometrySelection& selection,
    const std::unordered_map<PairKey, std::vector<std::complex<double>>>& background_mean
) -> std::vector<SelectedTrace> {
    std::vector<SelectedTrace> traces{};
    traces.reserve(collection.traces.size());
    std::unordered_set<PairKey> seen_keys{};

    for (std::size_t trace_index = 0U; trace_index < collection.traces.size(); ++trace_index) {
        const auto& trace = collection.traces[trace_index];
        const auto output_it = selection.output_local_by_pos.find(trace.combo.output_pos);
        const auto input_it = selection.input_local_by_pos.find(trace.combo.input_pos);
        if (output_it == selection.output_local_by_pos.end() || input_it == selection.input_local_by_pos.end()) {
            continue;
        }
        if (trace.frequency_hz.size() != trace.s21.size()) {
            continue;
        }

        SelectedTrace selected{};
        selected.combo = trace.combo;
        selected.tx_local_index = output_it->second;
        selected.rx_local_index = input_it->second;
        selected.run_order = trace_index;
        selected.frequency_hz.reserve(trace.frequency_hz.size());
        for (const auto value : trace.frequency_hz) {
            selected.frequency_hz.push_back(static_cast<double>(value));
        }

        const auto key = make_pair_key(selected.tx_local_index, selected.rx_local_index);
        if (!seen_keys.insert(key).second) {
            throw std::runtime_error("GPR requires unique selected combos; duplicate combo: " + combo_to_string(trace.combo));
        }

        const auto background_it = background_mean.find(key);
        selected.s21.reserve(trace.s21.size());
        for (std::size_t sample_index = 0U; sample_index < trace.s21.size(); ++sample_index) {
            std::complex<double> sample(trace.s21[sample_index].re, trace.s21[sample_index].im);
            if (background_it != background_mean.end() && background_it->second.size() == trace.s21.size()) {
                sample -= background_it->second[sample_index];
            }
            selected.s21.push_back(sample);
        }

        traces.push_back(std::move(selected));
    }

    return traces;
}

[[nodiscard]] auto compute_ascan(
    const SelectedTrace& trace,
    double start_hz,
    double stop_hz
) -> AscanResult {
    AscanResult result{};
    if (trace.frequency_hz.size() != trace.s21.size()) {
        return result;
    }

    const double low_hz = std::min(start_hz, stop_hz);
    const double high_hz = std::max(start_hz, stop_hz);
    std::vector<double> frequency_hz{};
    std::vector<std::complex<double>> s21{};
    frequency_hz.reserve(trace.frequency_hz.size());
    s21.reserve(trace.s21.size());

    for (std::size_t index = 0U; index < trace.frequency_hz.size(); ++index) {
        const double frequency_value = trace.frequency_hz[index];
        if (frequency_value < low_hz || frequency_value > high_hz) {
            continue;
        }
        frequency_hz.push_back(frequency_value);
        s21.push_back(trace.s21[index]);
    }

    if (frequency_hz.size() < 2U) {
        return result;
    }

    std::vector<double> df_values{};
    df_values.reserve(frequency_hz.size() - 1U);
    for (std::size_t index = 1U; index < frequency_hz.size(); ++index) {
        const double df = frequency_hz[index] - frequency_hz[index - 1U];
        if (!(df > 0.0)) {
            return result;
        }
        df_values.push_back(df);
    }

    const double df_hz = median_copy(std::move(df_values));
    if (!(df_hz > 0.0)) {
        return result;
    }

    const auto start_bin = static_cast<std::int64_t>(std::llround(frequency_hz.front() / df_hz));
    if (start_bin < 0) {
        return result;
    }

    const std::size_t point_count = frequency_hz.size();
    const auto start_index = static_cast<std::size_t>(start_bin);
    const std::size_t min_fft_len = 2U * (start_index + point_count - 1U);
    const std::size_t base_fft_len = next_power_of_two(min_fft_len);
    if (base_fft_len < min_fft_len || base_fft_len > (std::numeric_limits<std::size_t>::max() / kAscanOversample)) {
        return result;
    }

    const std::size_t fft_len = base_fft_len * kAscanOversample;
    if (start_index > fft_len || point_count > (fft_len - start_index)) {
        return result;
    }

    std::vector<std::complex<double>> spectrum(fft_len, std::complex<double>(0.0, 0.0));
    for (std::size_t index = 0U; index < point_count; ++index) {
        const double window = point_count > 1U
            ? 0.5 - (0.5 * std::cos((2.0 * kPi * static_cast<double>(index)) / static_cast<double>(point_count - 1U)))
            : 1.0;
        spectrum[start_index + index] = s21[index] * window;
    }

    fft_inplace(spectrum, true);

    result.bandwidth_hz = frequency_hz.back() - frequency_hz.front();
    result.dt_s = 1.0 / (static_cast<double>(fft_len) * df_hz);
    result.samples = std::move(spectrum);
    result.time_s.resize(fft_len, 0.0);
    for (std::size_t index = 0U; index < fft_len; ++index) {
        result.time_s[index] = static_cast<double>(index) * result.dt_s;
    }
    return result;
}

void normalize_pair_ascans(
    std::unordered_map<PairKey, AscanResult>& ascans_by_pair,
    double velocity_mps,
    double min_depth_m,
    double max_depth_m
) {
    if (!kPairNormalize) {
        return;
    }

    for (auto& item : ascans_by_pair) {
        auto& ascan = item.second;
        if (ascan.samples.empty() || ascan.time_s.size() != ascan.samples.size()) {
            continue;
        }

        std::vector<double> amplitudes{};
        amplitudes.reserve(ascan.samples.size());
        for (std::size_t index = 0U; index < ascan.samples.size(); ++index) {
            const double depth_m = 0.5 * velocity_mps * ascan.time_s[index];
            if (depth_m < min_depth_m || depth_m > max_depth_m) {
                continue;
            }
            amplitudes.push_back(std::abs(ascan.samples[index]));
        }
        if (amplitudes.empty()) {
            amplitudes.reserve(ascan.samples.size());
            for (const auto& sample : ascan.samples) {
                amplitudes.push_back(std::abs(sample));
            }
        }

        double scale = percentile_copy(std::move(amplitudes), kPairNormPercentile);
        if (!std::isfinite(scale) || !(scale > kPairNormEps)) {
            scale = 1.0;
        }

        const double denominator = scale + kPairNormEps;
        for (auto& sample : ascan.samples) {
            sample /= denominator;
        }
    }
}

[[nodiscard]] auto distance_3d(double dx, double dy, double dz) -> double {
    return std::sqrt((dx * dx) + (dy * dy) + (dz * dz));
}

[[nodiscard]] auto build_grid(
    const GeometrySelection& selection,
    double max_depth_m,
    double min_z_m,
    double imaging_plane_y_m
) -> GridDefinition {
    GridDefinition grid{};
    if (selection.x_tx.empty() || selection.x_rx.empty() || !(max_depth_m > min_z_m)) {
        return grid;
    }

    const auto [tx_min_it, tx_max_it] = std::minmax_element(selection.x_tx.begin(), selection.x_tx.end());
    const auto [rx_min_it, rx_max_it] = std::minmax_element(selection.x_rx.begin(), selection.x_rx.end());
    const double x_min = std::min(*tx_min_it, *rx_min_it) - kXMarginM;
    const double x_max = std::max(*tx_max_it, *rx_max_it) + kXMarginM;

    grid.x_grid = build_axis(x_min, x_max, kGridWidth);
    grid.z_grid = build_axis(min_z_m, max_depth_m, kGridHeight);

    const std::size_t cell_count = grid.x_grid.size() * grid.z_grid.size();
    grid.tx_distance_grids.assign(selection.x_tx.size(), std::vector<double>(cell_count, 0.0));
    grid.rx_distance_grids.assign(selection.x_rx.size(), std::vector<double>(cell_count, 0.0));

    for (std::size_t row = 0U; row < grid.z_grid.size(); ++row) {
        const double z_value = grid.z_grid[row];
        for (std::size_t col = 0U; col < grid.x_grid.size(); ++col) {
            const double x_value = grid.x_grid[col];
            const auto cell_index = (row * grid.x_grid.size()) + col;

            for (std::size_t tx_index = 0U; tx_index < selection.x_tx.size(); ++tx_index) {
                grid.tx_distance_grids[tx_index][cell_index] = distance_3d(
                    x_value - selection.x_tx[tx_index],
                    imaging_plane_y_m - selection.y_tx[tx_index],
                    z_value - selection.z_tx[tx_index]
                );
            }
            for (std::size_t rx_index = 0U; rx_index < selection.x_rx.size(); ++rx_index) {
                grid.rx_distance_grids[rx_index][cell_index] = distance_3d(
                    x_value - selection.x_rx[rx_index],
                    imaging_plane_y_m - selection.y_rx[rx_index],
                    z_value - selection.z_rx[rx_index]
                );
            }
        }
    }

    return grid;
}

[[nodiscard]] auto interpolate_complex(const AscanResult& ascan, double tau_s) -> std::complex<double> {
    if (ascan.samples.empty() || !(ascan.dt_s > 0.0) || tau_s < 0.0) {
        return {0.0, 0.0};
    }

    const double position = tau_s / ascan.dt_s;
    const auto index = static_cast<std::size_t>(std::floor(position));
    if (index >= ascan.samples.size()) {
        return {0.0, 0.0};
    }
    if (index + 1U >= ascan.samples.size()) {
        return ascan.samples[index];
    }

    const double frac = position - static_cast<double>(index);
    return (ascan.samples[index] * (1.0 - frac)) + (ascan.samples[index + 1U] * frac);
}

[[nodiscard]] auto attenuation_components(
    double r_tx,
    double r_rx,
    double dz_tx,
    double dz_rx
) -> std::pair<double, double> {
    // Antennas boresight along +Z, so cos(theta) = (z_pixel - z_antenna) / R.
    const double geo = 1.0 / ((r_tx * r_rx) + 1e-12);
    const double angle =
        std::pow(dz_tx / (r_tx + 1e-12), 2.0) *
        std::pow(dz_rx / (r_rx + 1e-12), 2.0);
    return {geo + 1e-30, angle + 1e-30};
}

[[nodiscard]] auto attenuation_components_at_ref_depth(
    std::uint32_t tx_index,
    std::uint32_t rx_index,
    const GeometrySelection& selection,
    double imaging_plane_y_m
) -> std::pair<double, double> {
    const double x_center = 0.5 * (selection.x_tx[tx_index] + selection.x_rx[rx_index]);
    const double dz_tx = kCompensationReferenceDepthM - selection.z_tx[tx_index];
    const double dz_rx = kCompensationReferenceDepthM - selection.z_rx[rx_index];
    const double r_tx = distance_3d(
        x_center - selection.x_tx[tx_index],
        imaging_plane_y_m - selection.y_tx[tx_index],
        dz_tx
    );
    const double r_rx = distance_3d(
        x_center - selection.x_rx[rx_index],
        imaging_plane_y_m - selection.y_rx[rx_index],
        dz_rx
    );
    return attenuation_components(r_tx, r_rx, dz_tx, dz_rx);
}

[[nodiscard]] auto compensation_weight(
    double r_tx,
    double r_rx,
    double dz_tx,
    double dz_rx,
    double geo_ref,
    double angle_ref,
    double range_power,
    double angle_power
) -> double {
    const auto [geo, angle] = attenuation_components(r_tx, r_rx, dz_tx, dz_rx);
    const double geo_norm = geo / geo_ref;
    const double angle_norm = angle / angle_ref;

    const double range_weight = std::clamp(
        1.0 / (std::pow(geo_norm, range_power) + 1e-12),
        0.0,
        kRangeWeightMax
    );
    const double angle_weight = std::clamp(
        1.0 / (std::pow(angle_norm, angle_power) + 1e-12),
        0.0,
        kAngleWeightMax
    );
    return std::clamp(range_weight * angle_weight, 0.0, kTotalWeightMax);
}

[[nodiscard]] auto backproject_coherent(
    const std::vector<SelectedTrace>& selected_traces,
    const std::unordered_map<PairKey, AscanResult>& ascans_by_pair,
    const GridDefinition& grid,
    const GeometrySelection& selection,
    double imaging_plane_y_m,
    double velocity_mps,
    double min_depth_m,
    double max_depth_m,
    double range_power,
    double angle_power
) -> BpMap {
    const std::size_t width = grid.x_grid.size();
    const std::size_t height = grid.z_grid.size();
    const std::size_t cell_count = width * height;

    BpMap result{};
    result.image.assign(cell_count, 0.0);
    result.coherent.assign(cell_count, std::complex<double>(0.0, 0.0));
    result.incoherent.assign(cell_count, 0.0);
    result.coherence_factor.assign(cell_count, 0.0);
    std::vector<double> contribution_count(cell_count, 0.0);

    for (const auto& trace : selected_traces) {
        const auto key = make_pair_key(trace.tx_local_index, trace.rx_local_index);
        const auto ascan_it = ascans_by_pair.find(key);
        if (ascan_it == ascans_by_pair.end()) {
            continue;
        }
        const auto& ascan = ascan_it->second;
        if (ascan.time_s.empty()) {
            continue;
        }

        const auto [geo_ref, angle_ref] = attenuation_components_at_ref_depth(
            trace.tx_local_index,
            trace.rx_local_index,
            selection,
            imaging_plane_y_m
        );
        const auto& tx_distances = grid.tx_distance_grids[trace.tx_local_index];
        const auto& rx_distances = grid.rx_distance_grids[trace.rx_local_index];
        const double z_tx_ant = selection.z_tx[trace.tx_local_index];
        const double z_rx_ant = selection.z_rx[trace.rx_local_index];

        for (std::size_t row = 0U; row < height; ++row) {
            const double z_m = grid.z_grid[row];
            const bool in_depth_gate = z_m >= min_depth_m && z_m <= max_depth_m;
            if (!in_depth_gate) {
                continue;
            }
            const double dz_tx = z_m - z_tx_ant;
            const double dz_rx = z_m - z_rx_ant;

            for (std::size_t col = 0U; col < width; ++col) {
                const auto cell_index = (row * width) + col;
                const double r_tx = tx_distances[cell_index];
                const double r_rx = rx_distances[cell_index];
                const double tau_s = (r_tx + r_rx) / velocity_mps;
                if (tau_s < ascan.time_s.front() || tau_s > ascan.time_s.back()) {
                    continue;
                }
                const auto sample = interpolate_complex(ascan, tau_s);

                const double weight = compensation_weight(
                    r_tx,
                    r_rx,
                    dz_tx,
                    dz_rx,
                    geo_ref,
                    angle_ref,
                    range_power,
                    angle_power
                );
                result.coherent[cell_index] += sample * weight;
                result.incoherent[cell_index] += std::abs(sample) * weight;
                contribution_count[cell_index] += 1.0;
            }
        }
    }

    for (std::size_t index = 0U; index < cell_count; ++index) {
        if (!(contribution_count[index] > 0.0)) {
            continue;
        }
        result.coherent[index] /= contribution_count[index];
        result.incoherent[index] /= contribution_count[index];
        result.image[index] = std::abs(result.coherent[index]);
        result.coherence_factor[index] =
            std::clamp(result.image[index] / (result.incoherent[index] + kScoreCfEps), 0.0, 1.0);
    }

    return result;
}

void apply_depth_gate(
    std::vector<double>& image,
    const std::vector<double>& z_grid,
    std::size_t width,
    double min_depth_m,
    double max_depth_m
) {
    for (std::size_t row = 0U; row < z_grid.size(); ++row) {
        const double z_m = z_grid[row];
        if (z_m >= min_depth_m && z_m <= max_depth_m) {
            continue;
        }
        const auto row_offset = row * width;
        std::fill(image.begin() + static_cast<std::ptrdiff_t>(row_offset),
                  image.begin() + static_cast<std::ptrdiff_t>(row_offset + width),
                  0.0);
    }
}

[[nodiscard]] auto normalize_bp_map(
    const std::vector<double>& raw_image,
    const GridDefinition& grid,
    double min_depth_m,
    double max_depth_m,
    double smooth_sigma
) -> std::vector<double> {
    auto normalized = raw_image;
    normalize_in_place(normalized);
    normalized = gaussian_filter_2d(normalized, grid.x_grid.size(), grid.z_grid.size(), smooth_sigma);
    apply_depth_gate(normalized, grid.z_grid, grid.x_grid.size(), min_depth_m, max_depth_m);
    normalize_in_place(normalized);
    return normalized;
}

[[nodiscard]] auto neighbor_indices(
    std::size_t index,
    std::size_t width,
    std::size_t height
) -> std::vector<std::size_t> {
    const auto row = index / width;
    const auto col = index % width;
    std::vector<std::size_t> neighbors{};
    neighbors.reserve(8U);

    for (std::ptrdiff_t row_offset = -1; row_offset <= 1; ++row_offset) {
        for (std::ptrdiff_t col_offset = -1; col_offset <= 1; ++col_offset) {
            if (row_offset == 0 && col_offset == 0) {
                continue;
            }
            const auto next_row = static_cast<std::ptrdiff_t>(row) + row_offset;
            const auto next_col = static_cast<std::ptrdiff_t>(col) + col_offset;
            if (next_row < 0 || next_col < 0) {
                continue;
            }
            if (next_row >= static_cast<std::ptrdiff_t>(height) || next_col >= static_cast<std::ptrdiff_t>(width)) {
                continue;
            }
            neighbors.push_back((static_cast<std::size_t>(next_row) * width) + static_cast<std::size_t>(next_col));
        }
    }

    return neighbors;
}

[[nodiscard]] auto component_containing_peak(
    const std::vector<double>& image,
    std::size_t width,
    std::size_t height,
    std::size_t peak_index,
    double threshold,
    const std::vector<std::uint8_t>* window_mask = nullptr
) -> std::vector<std::uint8_t> {
    std::vector<std::uint8_t> result(image.size(), 0U);
    if (image.empty() || peak_index >= image.size()) {
        return result;
    }
    if (image[peak_index] < threshold || (window_mask != nullptr && (*window_mask)[peak_index] == 0U)) {
        result[peak_index] = 1U;
        return result;
    }

    std::queue<std::size_t> queue{};
    result[peak_index] = 1U;
    queue.push(peak_index);

    while (!queue.empty()) {
        const auto current = queue.front();
        queue.pop();

        for (const auto neighbor : neighbor_indices(current, width, height)) {
            if (result[neighbor] != 0U) {
                continue;
            }
            if (window_mask != nullptr && (*window_mask)[neighbor] == 0U) {
                continue;
            }
            if (image[neighbor] < threshold) {
                continue;
            }
            result[neighbor] = 1U;
            queue.push(neighbor);
        }
    }

    return result;
}

[[nodiscard]] auto mask_count(const std::vector<std::uint8_t>& mask) -> std::size_t {
    return static_cast<std::size_t>(std::count(mask.begin(), mask.end(), static_cast<std::uint8_t>(1U)));
}

[[nodiscard]] auto weighted_centroid(
    const std::vector<double>& image,
    const std::vector<std::uint8_t>& mask,
    const GridDefinition& grid,
    double threshold
) -> std::pair<double, double> {
    double weight_sum = 0.0;
    double x_weighted_sum = 0.0;
    double z_weighted_sum = 0.0;

    for (std::size_t index = 0U; index < image.size(); ++index) {
        if (mask[index] == 0U) {
            continue;
        }
        double weight = std::max(image[index] - threshold, 0.0);
        if (!(weight > 0.0)) {
            weight = image[index];
        }
        if (!(weight > 0.0)) {
            continue;
        }

        const auto row = index / grid.x_grid.size();
        const auto col = index % grid.x_grid.size();
        weight_sum += weight;
        x_weighted_sum += grid.x_grid[col] * weight;
        z_weighted_sum += grid.z_grid[row] * weight;
    }

    if (!(weight_sum > 0.0)) {
        const auto found = std::find(mask.begin(), mask.end(), static_cast<std::uint8_t>(1U));
        if (found == mask.end()) {
            return {0.0, 0.0};
        }
        const auto index = static_cast<std::size_t>(std::distance(mask.begin(), found));
        return {grid.x_grid[index % grid.x_grid.size()], grid.z_grid[index / grid.x_grid.size()]};
    }

    return {x_weighted_sum / weight_sum, z_weighted_sum / weight_sum};
}

[[nodiscard]] auto compact_peak_centroid(
    const std::vector<double>& image,
    const GridDefinition& grid,
    std::size_t peak_index,
    double peak
) -> std::tuple<double, double, std::vector<std::uint8_t>> {
    const std::size_t width = grid.x_grid.size();
    const auto peak_row = peak_index / width;
    const auto peak_col = peak_index % width;
    const double x_peak = grid.x_grid[peak_col];
    const double z_peak = grid.z_grid[peak_row];
    const double threshold = kCenterThresholdFrac * peak;

    std::vector<std::uint8_t> center_mask(image.size(), 0U);
    for (std::size_t row = 0U; row < grid.z_grid.size(); ++row) {
        const double z_m = grid.z_grid[row];
        if (std::abs(z_m - z_peak) > kCenterRadiusZM) {
            continue;
        }
        for (std::size_t col = 0U; col < grid.x_grid.size(); ++col) {
            const double x_m = grid.x_grid[col];
            const auto index = (row * width) + col;
            if (std::abs(x_m - x_peak) <= kCenterRadiusXM && image[index] >= threshold) {
                center_mask[index] = 1U;
            }
        }
    }

    if (mask_count(center_mask) == 0U) {
        center_mask[peak_index] = 1U;
    }

    double weight_sum = 0.0;
    double x_weighted_sum = 0.0;
    double z_weighted_sum = 0.0;
    for (std::size_t index = 0U; index < image.size(); ++index) {
        if (center_mask[index] == 0U) {
            continue;
        }
        double weight = std::pow(std::max(image[index] - threshold, 0.0), kCenterWeightPower);
        if (!(weight > 0.0)) {
            weight = image[index];
        }
        if (!(weight > 0.0)) {
            continue;
        }

        const auto row = index / width;
        const auto col = index % width;
        weight_sum += weight;
        x_weighted_sum += grid.x_grid[col] * weight;
        z_weighted_sum += grid.z_grid[row] * weight;
    }

    if (!(weight_sum > 0.0)) {
        return {x_peak, z_peak, std::move(center_mask)};
    }
    return {x_weighted_sum / weight_sum, z_weighted_sum / weight_sum, std::move(center_mask)};
}

[[nodiscard]] auto build_window_mask(
    const GridDefinition& grid,
    double x_center_m,
    double z_center_m,
    double radius_x_m,
    double radius_z_m
) -> std::vector<std::uint8_t> {
    const std::size_t width = grid.x_grid.size();
    std::vector<std::uint8_t> mask(width * grid.z_grid.size(), 0U);
    for (std::size_t row = 0U; row < grid.z_grid.size(); ++row) {
        if (std::abs(grid.z_grid[row] - z_center_m) > radius_z_m) {
            continue;
        }
        for (std::size_t col = 0U; col < grid.x_grid.size(); ++col) {
            if (std::abs(grid.x_grid[col] - x_center_m) <= radius_x_m) {
                mask[(row * width) + col] = 1U;
            }
        }
    }
    return mask;
}

[[nodiscard]] auto find_bp_objects(
    const std::vector<double>& bp_image,
    const GridDefinition& grid
) -> std::vector<ObjectRecord> {
    std::vector<ObjectRecord> objects{};
    if (bp_image.empty() || grid.x_grid.size() < 2U || grid.z_grid.size() < 2U) {
        return objects;
    }

    std::vector<double> work = bp_image;
    const double global_peak = max_value(work);
    const double stop_level = kObjectMinFrac * global_peak;
    if (!(global_peak > 0.0)) {
        return objects;
    }

    const std::size_t width = grid.x_grid.size();
    const std::size_t height = grid.z_grid.size();
    const double dx_cm = std::abs(grid.x_grid[1] - grid.x_grid[0]) * 100.0;
    const double dz_cm = std::abs(grid.z_grid[1] - grid.z_grid[0]) * 100.0;
    const double pixel_area_cm2 = dx_cm * dz_cm;

    for (std::size_t step = 0U; step < kMaxObjects; ++step) {
        const auto peak_it = std::max_element(work.begin(), work.end());
        const double peak = *peak_it;
        if (peak <= stop_level || !(peak > 0.0)) {
            break;
        }

        const auto peak_index = static_cast<std::size_t>(std::distance(work.begin(), peak_it));
        const auto peak_row = peak_index / width;
        const auto peak_col = peak_index % width;
        const double x_peak = grid.x_grid[peak_col];
        const double z_peak = grid.z_grid[peak_row];

        const double region_threshold = kRegionThresholdFrac * peak;
        auto region_mask = component_containing_peak(work, width, height, peak_index, region_threshold);
        const double region_area_cm2 = static_cast<double>(mask_count(region_mask)) * pixel_area_cm2;
        if (region_area_cm2 < kMinRegionAreaCm2) {
            work[peak_index] = 0.0;
            continue;
        }

        const auto [x_region, z_region] = weighted_centroid(work, region_mask, grid, region_threshold);
        auto [x_center, z_center, center_mask] = compact_peak_centroid(work, grid, peak_index, peak);
        const double center_area_cm2 = static_cast<double>(mask_count(center_mask)) * pixel_area_cm2;

        double sum_value = 0.0;
        std::size_t value_count = 0U;
        for (std::size_t index = 0U; index < work.size(); ++index) {
            if (region_mask[index] == 0U) {
                continue;
            }
            sum_value += work[index];
            value_count += 1U;
        }

        ObjectRecord object{};
        object.index = objects.size() + 1U;
        object.peak_index = peak_index;
        object.x_peak_m = x_peak;
        object.z_peak_m = z_peak;
        object.x_m = x_center;
        object.z_m = z_center;
        object.x_region_m = x_region;
        object.z_region_m = z_region;
        object.peak = peak;
        object.area_cm2 = region_area_cm2;
        object.center_area_cm2 = center_area_cm2;
        object.mean_value = value_count > 0U ? sum_value / static_cast<double>(value_count) : 0.0;
        object.sum_value = sum_value;
        object.region_mask = std::move(region_mask);
        object.center_mask = std::move(center_mask);
        objects.push_back(std::move(object));

        const auto suppress_window = build_window_mask(grid, x_peak, z_peak, kSuppressRadiusXM, kSuppressRadiusZM);
        auto suppress_mask = component_containing_peak(
            work,
            width,
            height,
            peak_index,
            kSuppressThresholdFrac * peak,
            &suppress_window
        );
        if (mask_count(suppress_mask) < mask_count(objects.back().region_mask)) {
            suppress_mask = objects.back().region_mask;
        }

        for (std::size_t index = 0U; index < work.size(); ++index) {
            if (suppress_mask[index] != 0U) {
                work[index] = 0.0;
            }
        }
    }

    return objects;
}

[[nodiscard]] auto bistatic_depth_signature(
    const ObjectRecord& object,
    const std::vector<SelectedTrace>& selected_traces,
    const GeometrySelection& selection,
    double imaging_plane_y_m
) -> std::vector<double> {
    std::vector<double> signature{};
    signature.reserve(selected_traces.size());
    for (const auto& trace : selected_traces) {
        const auto tx = trace.tx_local_index;
        const auto rx = trace.rx_local_index;
        const double r_tx = distance_3d(
            object.x_m - selection.x_tx[tx],
            imaging_plane_y_m - selection.y_tx[tx],
            object.z_m - selection.z_tx[tx]
        );
        const double r_rx = distance_3d(
            object.x_m - selection.x_rx[rx],
            imaging_plane_y_m - selection.y_rx[rx],
            object.z_m - selection.z_rx[rx]
        );
        signature.push_back(0.5 * (r_tx + r_rx));
    }
    return signature;
}

void mark_sidelobe_candidates(
    std::vector<ObjectRecord>& objects,
    const std::vector<SelectedTrace>& selected_traces,
    const GeometrySelection& selection,
    double imaging_plane_y_m
) {
    std::vector<std::vector<double>> signatures{};
    signatures.reserve(objects.size());
    for (const auto& object : objects) {
        signatures.push_back(bistatic_depth_signature(object, selected_traces, selection, imaging_plane_y_m));
    }

    for (std::size_t object_index = 0U; object_index < objects.size(); ++object_index) {
        auto& object = objects[object_index];
        std::size_t best_parent = 0U;
        double best_rms = std::numeric_limits<double>::infinity();

        for (std::size_t parent_index = 0U; parent_index < object_index; ++parent_index) {
            const auto& parent = objects[parent_index];
            const double relative_peak = object.peak / (parent.peak + 1e-12);
            const double dx_m = std::abs(object.x_m - parent.x_m);
            const double dz_m = std::abs(object.z_m - parent.z_m);

            double squared_sum = 0.0;
            const auto& object_signature = signatures[object_index];
            const auto& parent_signature = signatures[parent_index];
            for (std::size_t index = 0U; index < object_signature.size(); ++index) {
                const double delta = object_signature[index] - parent_signature[index];
                squared_sum += delta * delta;
            }
            const double rms = object_signature.empty()
                ? std::numeric_limits<double>::infinity()
                : std::sqrt(squared_sum / static_cast<double>(object_signature.size()));

            const bool candidate =
                relative_peak <= kSidelobeMaxRelativePeak &&
                dx_m >= kSidelobeMinDxM &&
                dz_m <= kSidelobeMaxDzM &&
                rms <= kSidelobeRangeRmsToleranceM;

            if (candidate && rms < best_rms) {
                best_parent = parent.index;
                best_rms = rms;
            }
        }

        if (best_parent != 0U) {
            object.sidelobe_candidate = true;
            object.sidelobe_parent = best_parent;
            object.sidelobe_range_rms_m = best_rms;
        }
    }
}

void add_local_prominence_metrics(
    std::vector<ObjectRecord>& objects,
    const std::vector<double>& bp_image,
    const GridDefinition& grid,
    double min_depth_m,
    double max_depth_m
) {
    const std::size_t width = grid.x_grid.size();
    if (bp_image.empty() || width == 0U) {
        return;
    }

    for (auto& object : objects) {
        std::vector<double> bg_values{};
        for (std::size_t index = 0U; index < bp_image.size(); ++index) {
            const auto row = index / width;
            const auto col = index % width;
            const double x_m = grid.x_grid[col];
            const double z_m = grid.z_grid[row];
            const bool in_local_window =
                std::abs(x_m - object.x_peak_m) <= kLocalBgRadiusXM &&
                std::abs(z_m - object.z_peak_m) <= kLocalBgRadiusZM;
            const bool in_depth_gate = z_m >= min_depth_m && z_m <= max_depth_m;
            if (in_local_window && in_depth_gate && object.region_mask[index] == 0U) {
                bg_values.push_back(bp_image[index]);
            }
        }

        if (bg_values.size() < 10U) {
            bg_values.clear();
            for (std::size_t index = 0U; index < bp_image.size(); ++index) {
                const auto row = index / width;
                const double z_m = grid.z_grid[row];
                if (z_m >= min_depth_m && z_m <= max_depth_m && object.region_mask[index] == 0U) {
                    bg_values.push_back(bp_image[index]);
                }
            }
        }

        object.local_bg = percentile_copy(bg_values, kLocalBgPercentile);
        object.local_bg_p75 = percentile_copy(std::move(bg_values), 75.0);
        object.prominence = object.peak - object.local_bg;
        object.contrast = object.peak / (object.local_bg + kLocalContrastEps);
    }
}

[[nodiscard]] auto mean_values_under_mask(
    const std::vector<double>& values,
    const std::vector<std::uint8_t>& mask
) -> double {
    if (values.empty() || values.size() != mask.size()) {
        return 0.0;
    }

    double sum = 0.0;
    std::size_t count = 0U;
    for (std::size_t index = 0U; index < values.size(); ++index) {
        if (mask[index] == 0U) {
            continue;
        }
        sum += values[index];
        count += 1U;
    }
    return count > 0U ? sum / static_cast<double>(count) : 0.0;
}

[[nodiscard]] auto max_values_under_mask(
    const std::vector<double>& values,
    const std::vector<std::uint8_t>& mask
) -> double {
    if (values.empty() || values.size() != mask.size()) {
        return 0.0;
    }

    double maximum = 0.0;
    bool has_value = false;
    for (std::size_t index = 0U; index < values.size(); ++index) {
        if (mask[index] == 0U) {
            continue;
        }
        maximum = has_value ? std::max(maximum, values[index]) : values[index];
        has_value = true;
    }
    return has_value ? maximum : 0.0;
}

void add_incoherent_support_metrics(
    std::vector<ObjectRecord>& objects,
    const std::vector<double>& bp_incoherent_image,
    const std::vector<double>& bp_coherence_factor
) {
    for (auto& object : objects) {
        object.incoh_peak = max_values_under_mask(bp_incoherent_image, object.region_mask);
        object.incoh_mean = mean_values_under_mask(bp_incoherent_image, object.region_mask);
        object.incoh_center_mean = mean_values_under_mask(bp_incoherent_image, object.center_mask);
        if (!(object.incoh_center_mean > 0.0)) {
            object.incoh_center_mean = object.incoh_mean;
        }

        if (!bp_coherence_factor.empty() && object.peak_index < bp_coherence_factor.size()) {
            object.coherence_factor_peak = std::clamp(bp_coherence_factor[object.peak_index], 0.0, 1.0);
        } else {
            object.coherence_factor_peak =
                std::clamp(object.peak / (object.incoh_peak + kScoreCfEps), 0.0, 1.0);
        }

        const double center_cf = mean_values_under_mask(bp_coherence_factor, object.center_mask);
        object.coherence_factor_center = std::clamp(
            center_cf > 0.0
                ? center_cf
                : object.mean_value / (object.incoh_center_mean + kScoreCfEps),
            0.0,
            1.0
        );
    }
}

[[nodiscard]] auto contrast_score_unit(double contrast) -> double {
    if (!std::isfinite(contrast) || kScoreContrastCap <= 1.0) {
        return 0.0;
    }
    return std::clamp((contrast - 1.0) / (kScoreContrastCap - 1.0), 0.0, 1.0);
}

[[nodiscard]] auto use_combined_gpr_score(const ProcessingLiveConfig& live_config) -> bool {
    return live_config.gpr_score_mode == "combined";
}

void add_bp_score_metrics(
    std::vector<ObjectRecord>& objects,
    const ProcessingLiveConfig& live_config
) {
    double total_weight =
        kScoreCohPeakWeight +
        kScoreCoherenceFactorWeight +
        kScoreProminenceWeight +
        kScoreContrastWeight;
    if (!(total_weight > 0.0)) {
        total_weight = 1.0;
    }

    const bool combined_score = use_combined_gpr_score(live_config);
    for (auto& object : objects) {
        const double coh_peak_score = std::clamp(object.peak, 0.0, 1.0);
        const double coherence_factor_score = std::clamp(object.coherence_factor_peak, 0.0, 1.0);
        const double prominence_score = std::clamp(object.prominence, 0.0, 1.0);
        const double contrast_score = contrast_score_unit(object.contrast);

        object.score_old = coh_peak_score;
        object.score_new = (
            (kScoreCohPeakWeight * coh_peak_score) +
            (kScoreCoherenceFactorWeight * coherence_factor_score) +
            (kScoreProminenceWeight * prominence_score) +
            (kScoreContrastWeight * contrast_score)
        ) / total_weight;
        object.selected_score = combined_score ? object.score_new : object.score_old;
    }
}

[[nodiscard]] auto output_objects_sorted(
    const std::vector<ObjectRecord>& objects,
    const ProcessingLiveConfig& live_config
) -> std::vector<const ObjectRecord*> {
    std::vector<const ObjectRecord*> visible{};
    visible.reserve(objects.size());
    for (const auto& object : objects) {
        if (live_config.gpr_remove_sidelobe_objects_enabled && object.sidelobe_candidate) {
            continue;
        }
        visible.push_back(&object);
    }

    std::sort(visible.begin(), visible.end(), [](const ObjectRecord* left, const ObjectRecord* right) {
        if (left->selected_score == right->selected_score) {
            return left->peak > right->peak;
        }
        return left->selected_score > right->selected_score;
    });
    return visible;
}

[[nodiscard]] auto flatten_table(
    const std::vector<std::vector<float>>& rows,
    std::uint32_t column_count
) -> std::vector<float> {
    std::vector<float> values{};
    values.reserve(rows.size() * column_count);
    for (const auto& row : rows) {
        values.insert(values.end(), row.begin(), row.end());
    }
    return values;
}

[[nodiscard]] auto build_table_payload(
    const std::string& processing_name,
    const std::vector<std::vector<float>>& rows,
    std::uint32_t column_count
) -> ipc::ResultPayload {
    ipc::ResultPayload payload{};
    payload.processing_name = processing_name;
    payload.kind = ipc::ResultKind::TableF32;
    payload.table_columns = column_count;
    payload.table_values = flatten_table(rows, column_count);
    return payload;
}

[[nodiscard]] auto build_image_payload(
    const std::string& processing_name,
    const std::vector<double>& x_axis,
    const std::vector<double>& y_axis,
    const std::vector<double>& values
) -> ipc::ResultPayload {
    ipc::ResultPayload payload{};
    payload.processing_name = processing_name;
    payload.kind = ipc::ResultKind::ImageF32;
    payload.image_x_axis.reserve(x_axis.size());
    for (const auto value : x_axis) {
        payload.image_x_axis.push_back(static_cast<float>(value));
    }
    payload.image_y_axis.reserve(y_axis.size());
    for (const auto value : y_axis) {
        payload.image_y_axis.push_back(static_cast<float>(value));
    }
    payload.image_values.reserve(values.size());
    for (const auto value : values) {
        payload.image_values.push_back(static_cast<float>(value));
    }
    return payload;
}

[[nodiscard]] auto process_backprojection_gpr(
    const config::RunConfig& run_config,
    const ipc::PreprocessedCollection& collection,
    std::span<const ipc::PreprocessedCollection> previous_collections,
    const ProcessingLiveConfig& live_config
) -> ipc::ResultCollection {
    ipc::ResultCollection results{};
    results.collection_id = collection.collection_id;
    results.monotonic_ns = collection.monotonic_ns;

    const auto selection = build_geometry_selection(run_config, live_config);
    if (selection.input_positions.empty() || selection.output_positions.empty()) {
        return results;
    }

    validate_collection_trace_order(run_config, collection);
    const auto background_mean = build_background_mean(previous_collections, selection, live_config);
    const auto selected_traces = collect_selected_traces(collection, selection, background_mean);
    if (selected_traces.empty()) {
        return results;
    }

    const double velocity_mps =
        kSpeedOfLightMetersPerSec / std::sqrt(std::max(1e-6, static_cast<double>(run_config.gpr.relative_permittivity)));
    const double start_hz = static_cast<double>(live_config.gpr_start_freq_mhz) * 1'000'000.0;
    const double stop_hz = static_cast<double>(live_config.gpr_stop_freq_mhz) * 1'000'000.0;
    const double min_depth_m = static_cast<double>(live_config.gpr_min_depth_m);
    const double max_depth_m = static_cast<double>(live_config.gpr_max_depth_m);
    if (!(max_depth_m > min_depth_m) || !(max_depth_m > kGridZMinM)) {
        return results;
    }

    std::unordered_map<PairKey, AscanResult> ascans_by_pair{};
    for (const auto& trace : selected_traces) {
        auto ascan = compute_ascan(trace, start_hz, stop_hz);
        if (ascan.samples.empty() || !(ascan.bandwidth_hz > 0.0) || !(ascan.dt_s > 0.0)) {
            continue;
        }
        ascans_by_pair.emplace(
            make_pair_key(trace.tx_local_index, trace.rx_local_index),
            std::move(ascan)
        );
    }
    if (ascans_by_pair.empty()) {
        return results;
    }
    normalize_pair_ascans(ascans_by_pair, velocity_mps, min_depth_m, max_depth_m);

    const double imaging_plane_y_m = static_cast<double>(live_config.gpr_imaging_plane_y_m);
    const auto grid = build_grid(selection, max_depth_m, kGridZMinM, imaging_plane_y_m);
    if (grid.x_grid.empty() || grid.z_grid.empty()) {
        return results;
    }

    auto bp = backproject_coherent(
        selected_traces,
        ascans_by_pair,
        grid,
        selection,
        imaging_plane_y_m,
        velocity_mps,
        min_depth_m,
        max_depth_m,
        std::max(0.0, static_cast<double>(live_config.gpr_range_comp_power)),
        std::max(0.0, static_cast<double>(live_config.gpr_angle_comp_power))
    );
    if (bp.image.empty()) {
        return results;
    }

    auto display_map = normalize_bp_map(bp.image, grid, min_depth_m, max_depth_m, kSmoothSigma);
    const auto incoherent_display_map =
        normalize_bp_map(bp.incoherent, grid, min_depth_m, max_depth_m, kSmoothSigma);

    auto objects = find_bp_objects(display_map, grid);
    add_local_prominence_metrics(objects, display_map, grid, min_depth_m, max_depth_m);
    add_incoherent_support_metrics(objects, incoherent_display_map, bp.coherence_factor);
    mark_sidelobe_candidates(objects, selected_traces, selection, imaging_plane_y_m);
    add_bp_score_metrics(objects, live_config);

    if (live_config.gpr_remove_sidelobe_objects_enabled) {
        for (const auto& object : objects) {
            if (!object.sidelobe_candidate) {
                continue;
            }
            for (std::size_t index = 0U; index < display_map.size(); ++index) {
                if (object.region_mask[index] != 0U) {
                    display_map[index] = 0.0;
                }
            }
        }
    }

    results.collection_payloads.push_back(build_image_payload("gpr_accumulator", grid.x_grid, grid.z_grid, display_map));

    std::vector<std::vector<float>> point_rows{};
    point_rows.reserve(objects.size());
    for (const auto* object : output_objects_sorted(objects, live_config)) {
        point_rows.push_back(
            {
                static_cast<float>(object->x_m),
                static_cast<float>(object->z_m),
                static_cast<float>(object->selected_score),
            }
        );
    }
    results.collection_payloads.push_back(build_table_payload("gpr_points", point_rows, 3U));

    return results;
}
