// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include <algorithm>
#include <cmath>

#include "RingsFromProfile.h"
#include "AssignSpotsToRings.h"   // CalculateXtalRings
#include "../../common/JFJochMath.h"

namespace {

// Peak position of one ring in one azimuthal sector, in q, or NaN if there is no peak worth using.
//
// The window is narrow and centred on where the ring is expected, so the background under it is close
// to a straight line: take it from the two bins at each end and interpolate. The position itself is the
// intensity-weighted centroid of everything above half the peak height, which is insensitive to the
// exact half-maximum crossing and needs no line-shape assumption - a powder ring is not Gaussian, it is
// the instrumental profile convolved with whatever strain and size broadening the standard has.
float SectorPeakQ(const std::vector<float> &profile, int32_t q_bins, int phi_bin,
                  int lo_bin, int hi_bin, float low_q, float q_spacing, float min_peak_over_noise) {
    const size_t row = static_cast<size_t>(phi_bin) * static_cast<size_t>(q_bins);
    const auto value = [&](int i) { return profile[row + static_cast<size_t>(i)]; };
    const auto q_of = [&](int i) { return low_q + (static_cast<float>(i) + 0.5f) * q_spacing; };

    const float bkg_lo = 0.5f * (value(lo_bin) + value(lo_bin + 1));
    const float bkg_hi = 0.5f * (value(hi_bin) + value(hi_bin - 1));
    const auto bkg_at = [&](int i) {
        const float t = static_cast<float>(i - lo_bin) / static_cast<float>(hi_bin - lo_bin);
        return bkg_lo + t * (bkg_hi - bkg_lo);
    };

    int peak = -1;
    float peak_height = 0.0f;
    for (int i = lo_bin + 2; i <= hi_bin - 2; ++i) {
        const float h = value(i) - bkg_at(i);
        if (h > peak_height) { peak_height = h; peak = i; }
    }
    if (peak < 0)
        return NAN;

    // Scatter of the background shoulders, as the noise this peak has to stand clear of. A sector with
    // no ring in it has a "peak" that is just the largest background fluctuation, and this is what
    // rejects it - the alternative, an absolute intensity cut, would need a value per detector and beam.
    float s = 0.0f;
    int n = 0;
    for (int i : {lo_bin, lo_bin + 1, hi_bin - 1, hi_bin}) {
        const float r = value(i) - bkg_at(i);
        s += r * r;
        ++n;
    }
    const float noise = std::sqrt(s / static_cast<float>(n));
    if (!(peak_height > min_peak_over_noise * noise))
        return NAN;

    const float half = 0.5f * peak_height;
    double sum_wq = 0.0, sum_w = 0.0;
    for (int i = peak; i >= lo_bin && value(i) - bkg_at(i) >= half; --i) {
        const double w = value(i) - bkg_at(i);
        sum_wq += w * q_of(i);
        sum_w += w;
    }
    for (int i = peak + 1; i <= hi_bin && value(i) - bkg_at(i) >= half; ++i) {
        const double w = value(i) - bkg_at(i);
        sum_wq += w * q_of(i);
        sum_w += w;
    }
    if (!(sum_w > 0.0))
        return NAN;
    return static_cast<float>(sum_wq / sum_w);
}

}   // namespace

std::vector<RingOptimizerInput> RingsFromAzimuthalProfile(const std::vector<float> &profile,
                                                          const AzimuthalIntegrationMapping &mapping,
                                                          const DiffractionGeometry &geom,
                                                          const UnitCell &calibrant,
                                                          float q_window_recipA,
                                                          float min_peak_over_noise) {
    std::vector<RingOptimizerInput> out;

    const int32_t q_bins = mapping.GetQBinCount();
    const int32_t azim_bins = mapping.GetAzimuthalBinCount();
    // One azimuthal bin is a plain radial profile: the ring is averaged over every direction at once, so
    // nothing remains to say where its centre is. This needs the run to have been integrated with
    // azimuthal bins (jfjoch_broker azim_int_settings.azimuthal_bins, rugnux --azim-phi-bins).
    if (azim_bins < 4 || q_bins < 8
        || profile.size() != static_cast<size_t>(q_bins) * static_cast<size_t>(azim_bins))
        return out;

    const auto &settings = mapping.Settings();
    const float low_q = settings.GetLowQ_recipA();
    const float q_spacing = settings.GetQSpacing_recipA();
    const float high_q = low_q + static_cast<float>(q_bins) * q_spacing;
    const int window_bins = std::max(3, static_cast<int>(std::lround(q_window_recipA / q_spacing)));

    for (const float q_ring : CalculateXtalRings(calibrant)) {
        if (!(q_ring - q_window_recipA > low_q) || !(q_ring + q_window_recipA < high_q))
            continue;
        const int centre_bin = static_cast<int>((q_ring - low_q) / q_spacing);
        const int lo_bin = std::max(0, centre_bin - window_bins);
        const int hi_bin = std::min(q_bins - 1, centre_bin + window_bins);
        if (hi_bin - lo_bin < 6)
            continue;

        for (int phi_bin = 0; phi_bin < azim_bins; ++phi_bin) {
            const float q_obs = SectorPeakQ(profile, q_bins, phi_bin, lo_bin, hi_bin,
                                            low_q, q_spacing, min_peak_over_noise);
            if (!std::isfinite(q_obs))
                continue;

            // The sector's CENTRE, not its lower edge: GetBin() floors phi into the sector, so a bin
            // stands for [j, j+1) and taking its edge would rotate every ring point by half a sector -
            // which is exactly the cos(phi) signal the beam centre is read from.
            const float phi_rad = static_cast<float>((static_cast<double>(phi_bin) + 0.5)
                                                     * 2.0 * PI / static_cast<double>(azim_bins));
            const auto [x, y] = geom.ResPhiToPxl(static_cast<float>(2.0 * PI) / q_obs, phi_rad);
            if (!std::isfinite(x) || !std::isfinite(y))
                continue;
            out.push_back({x, y, q_ring});
        }
    }
    return out;
}