// 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"   // RingMatchWindow
#include "PowderAutoSeed.h"        // ProfileRingTrack
#include "../../common/JFJochMath.h"

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; };

    // A bin no pixel fell in is NaN, not zero (AzimuthalIntegrationProfile::GetResult), and the four
    // background bins are where a module gap or the beam stop shows up first. Say so rather than
    // relying on NaN comparisons to fail the peak test further down: a sector whose background cannot
    // be measured has no measurable peak either.
    for (int i : {lo_bin, lo_bin + 1, hi_bin - 1, hi_bin}) {
        if (!std::isfinite(value(i)))
            return NAN;
    }

    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);
}

std::vector<RingOptimizerInput> RingsFromAzimuthalProfile(const std::vector<float> &profile,
                                                          const AzimuthalIntegrationMapping &mapping,
                                                          const DiffractionGeometry &geom,
                                                          const std::vector<float> &calibrant_ring_q,
                                                          float q_window_recipA,
                                                          float min_peak_over_noise,
                                                          const DiffractionGeometry *seeded) {
    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;

    // Where to LOOK, ring by ring and sector by sector. Without a seed a ring is looked for at its own
    // q in every sector, which is the right answer only when the geometry that binned the profile was
    // already close; with one, each ring is tracked through the profile it really made.
    const size_t rings = calibrant_ring_q.size();
    std::vector<std::vector<float>> track(rings);
    for (size_t i = 0; i < rings; ++i) {
        if (seeded)
            track[i] = ProfileRingTrack(calibrant_ring_q[i], *seeded, geom, azim_bins);
        else
            track[i].assign(azim_bins, calibrant_ring_q[i]);
    }

    for (size_t i = 0; i < rings; ++i) {
        for (int phi_bin = 0; phi_bin < azim_bins; ++phi_bin) {
            const float q_ring = track[i][phi_bin];
            if (!std::isfinite(q_ring))
                continue;

            // Never let the window reach into the neighbouring ring. SectorPeakQ takes the background
            // under the peak from the two bins at each end of the window, so a window wider than half
            // the gap to the next ring measures that ring's flank as this one's background. Hexagonal
            // ice has three rings within 0.06 1/A of one another, which a fixed window merges into one
            // peak. Measured against the neighbours IN THIS SECTOR, since that is where they are here.
            float window = q_window_recipA;
            if (i > 0 && std::isfinite(track[i - 1][phi_bin]))
                window = std::min(window, 0.5f * std::abs(q_ring - track[i - 1][phi_bin]));
            if (i + 1 < rings && std::isfinite(track[i + 1][phi_bin]))
                window = std::min(window, 0.5f * std::abs(track[i + 1][phi_bin] - q_ring));

            if (!(q_ring - window > low_q) || !(q_ring + window < high_q))
                continue;
            const int window_bins = static_cast<int>(std::lround(window / q_spacing));
            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);
            // Two background bins at each end and a peak between them is the least this can work with;
            // a ring whose window is narrower than that is not resolved at this q spacing.
            if (hi_bin - lo_bin < 6)
                continue;

            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, calibrant_ring_q[i]});
        }
    }
    return out;
}