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

#include "ProfileIntegrate2D.h"

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <limits>
#include <vector>

#include "../../common/JFJochException.h"

namespace {

constexpr int N_SHELL = 6;
constexpr double STRONG_I_OVER_SIGMA = 5.0;
constexpr int MIN_STRONG_PER_SHELL = 30;

// Rough box-sum of one reflection: total intensity over the inner disk (r1) minus the mean of the
// background ring (r2..r3), like BraggIntegrate2D. Used to seed the fit, pick strong spots for
// profile learning, and provide the per-reflection background. ok=false when the inner disk is
// clipped by a gap/edge/bad pixel (same rejection as the box-sum integrator).
struct Rough {
    double I = 0.0, sigma = NAN, bkg = 0.0;
    int cx = 0, cy = 0, shell = -1;
    bool ok = false, strong = false;
};

template<class T>
Rough BoxSum(const T *image, size_t xpixel, size_t ypixel, int64_t special, int64_t saturation,
             const Reflection &r, float r1_sq, float r2_sq, float r3_sq, float r3, float min_sigma_ratio) {
    Rough out;
    const int64_t x0 = std::max<int64_t>(0, std::floor(r.predicted_x - r3 - 1.0));
    const int64_t x1 = std::min<int64_t>(xpixel - 1, std::ceil(r.predicted_x + r3 + 1.0));
    const int64_t y0 = std::max<int64_t>(0, std::floor(r.predicted_y - r3 - 1.0));
    const int64_t y1 = std::min<int64_t>(ypixel - 1, std::ceil(r.predicted_y + r3 + 1.0));

    int64_t I_sum = 0, n_inner = 0, n_inner_valid = 0;
    double bkg_sum = 0.0;
    int n_bkg = 0;
    for (int64_t y = y0; y <= y1; ++y)
        for (int64_t x = x0; x <= x1; ++x) {
            const double d2 = (x - r.predicted_x) * (x - r.predicted_x) + (y - r.predicted_y) * (y - r.predicted_y);
            const auto px = image[y * xpixel + x];
            if (d2 < r1_sq) {
                ++n_inner;
                if (px == special || px == special + 1 || px == saturation || px == saturation - 1) continue;
                I_sum += px;
                ++n_inner_valid;
            } else if (d2 >= r2_sq && d2 < r3_sq) {
                if (px == special || px == special + 1 || px == saturation || px == saturation - 1) continue;
                bkg_sum += static_cast<double>(px);
                ++n_bkg;
            }
        }
    if (n_inner_valid != n_inner || n_bkg <= 5)
        return out;
    out.bkg = bkg_sum / n_bkg;
    out.I = static_cast<double>(I_sum) - static_cast<double>(n_inner) * out.bkg;
    out.sigma = std::max(1.0, out.I * min_sigma_ratio);
    if (I_sum > 0)
        out.sigma = std::max(out.sigma, std::sqrt(static_cast<double>(I_sum)));
    out.cx = static_cast<int>(std::lround(r.predicted_x));
    out.cy = static_cast<int>(std::lround(r.predicted_y));
    out.ok = true;
    out.strong = out.sigma > 0.0 && out.I / out.sigma >= STRONG_I_OVER_SIGMA;
    return out;
}

template<class T>
std::vector<Reflection> ProfileIntegrateInternal(const DiffractionExperiment &experiment,
                                                 const CompressedImage &image,
                                                 const std::vector<Reflection> &predicted, size_t npredicted,
                                                 int64_t special, int64_t saturation, int64_t image_number) {
    const auto settings = experiment.GetBraggIntegrationSettings();
    const auto geom = experiment.GetDiffractionGeometry();
    const bool empirical = settings.GetIntegrator() == IntegratorMode::ProfileEmpirical;

    const float r1_sq = settings.GetR1() * settings.GetR1();
    const float r2 = settings.GetR2(), r2_sq = r2 * r2;
    const float r3 = settings.GetR3(), r3_sq = r3 * r3;
    const float min_sigma_ratio = settings.GetMinimumSigmaInRegardsToI();
    const int R = static_cast<int>(std::ceil(r2));   // profile grid half-size
    const int G = 2 * R + 1, GG = G * G;
    auto idx = [G, R](int dx, int dy) { return (dy + R) * G + (dx + R); };

    std::vector<uint8_t> buffer;
    const auto *ptr = reinterpret_cast<const T *>(image.GetUncompressedPtr(buffer));
    const size_t xpixel = image.GetWidth(), ypixel = image.GetHeight();

    // --- Pass A: box-sum every reflection (rough I, background, centre, strong flag). ---
    std::vector<Rough> rough(npredicted);
    double inv_d2_min = std::numeric_limits<double>::max(), inv_d2_max = 0.0;
    for (size_t i = 0; i < npredicted; ++i) {
        rough[i] = BoxSum<T>(ptr, xpixel, ypixel, special, saturation, predicted[i],
                             r1_sq, r2_sq, r3_sq, r3, min_sigma_ratio);
        if (rough[i].ok && predicted[i].d > 0.0f) {
            const double inv_d2 = 1.0 / (predicted[i].d * predicted[i].d);
            inv_d2_min = std::min(inv_d2_min, inv_d2);
            inv_d2_max = std::max(inv_d2_max, inv_d2);
        }
    }
    auto shell_of = [&](float d) {
        if (!(d > 0.0f) || inv_d2_max <= inv_d2_min) return 0;
        const double t = (1.0 / (d * d) - inv_d2_min) / (inv_d2_max - inv_d2_min);
        return std::clamp(static_cast<int>(t * N_SHELL), 0, N_SHELL - 1);
    };
    for (size_t i = 0; i < npredicted; ++i)
        if (rough[i].ok)
            rough[i].shell = shell_of(predicted[i].d);

    // --- Learn the profile per shell from the strong spots (+ a global fallback). Each strong
    //     spot contributes its background-subtracted, intensity-normalised, centroid-aligned grid. ---
    std::vector<std::vector<double>> shell_grid(N_SHELL, std::vector<double>(GG, 0.0));
    std::vector<int> shell_n(N_SHELL, 0);
    std::vector<double> global_grid(GG, 0.0);
    int global_n = 0;
    for (size_t i = 0; i < npredicted; ++i) {
        const auto &rh = rough[i];
        if (!rh.ok || !rh.strong || rh.I <= 0.0) continue;
        for (int dy = -R; dy <= R; ++dy)
            for (int dx = -R; dx <= R; ++dx) {
                const int64_t x = rh.cx + dx, y = rh.cy + dy;
                if (x < 0 || y < 0 || x >= static_cast<int64_t>(xpixel) || y >= static_cast<int64_t>(ypixel)) continue;
                const auto px = ptr[y * xpixel + x];
                if (px == special || px == special + 1 || px == saturation || px == saturation - 1) continue;
                const double v = (static_cast<double>(px) - rh.bkg) / rh.I;
                shell_grid[rh.shell][idx(dx, dy)] += v;
                global_grid[idx(dx, dy)] += v;
            }
        ++shell_n[rh.shell];
        ++global_n;
    }

    // Build a normalised profile (sum = 1): the empirical average, or an isotropic Gaussian of the
    // same (measured) second moment - the spots are ~round in the detector plane, so radial/tangential
    // anisotropy was measured and found to add nothing (the elongation is in the discarded rocking
    // direction), and a single width is kept for simplicity.
    auto build_profile = [&](const std::vector<double> &grid, int n) {
        std::vector<double> P(GG, 0.0);
        if (n <= 0) return P;
        double sum = 0.0, m2 = 0.0, m2w = 0.0;
        for (int dy = -R; dy <= R; ++dy)
            for (int dx = -R; dx <= R; ++dx) {
                const double g = std::max(0.0, grid[idx(dx, dy)]);
                m2 += g * (dx * dx + dy * dy);
                m2w += g;
                sum += g;
                if (empirical) P[idx(dx, dy)] = g;
            }
        if (sum <= 0.0) return P;
        if (empirical) {
            for (double &p : P) p /= sum;
        } else {
            const double sigma2 = std::max(0.25, (m2w > 0.0 ? m2 / m2w : 1.0) / 2.0);   // <r^2> = 2 sigma^2 (2D)
            double gsum = 0.0;
            for (int dy = -R; dy <= R; ++dy)
                for (int dx = -R; dx <= R; ++dx) {
                    const double g = std::exp(-(dx * dx + dy * dy) / (2.0 * sigma2));
                    P[idx(dx, dy)] = g;
                    gsum += g;
                }
            for (double &p : P) p /= gsum;
        }
        return P;
    };
    const std::vector<double> global_P = build_profile(global_grid, global_n);
    std::vector<std::vector<double>> shell_P(N_SHELL);
    for (int s = 0; s < N_SHELL; ++s)
        shell_P[s] = shell_n[s] >= MIN_STRONG_PER_SHELL ? build_profile(shell_grid[s], shell_n[s]) : global_P;

    // --- Pass B: profile-fit each reflection (Kabsch, de-biased variance v = B + I*P; iterate). ---
    std::vector<Reflection> out;
    out.reserve(npredicted);
    const auto pol = experiment.GetPolarizationFactor();
    for (size_t i = 0; i < npredicted; ++i) {
        const auto &rh = rough[i];
        if (!rh.ok) continue;
        const auto &P = shell_P[rh.shell < 0 ? 0 : rh.shell];
        const double B = std::max(rh.bkg, 1.0);

        double I = rh.I, den = 0.0;
        for (int iter = 0; iter < 4; ++iter) {
            double num = 0.0;
            den = 0.0;
            for (int dy = -R; dy <= R; ++dy)
                for (int dx = -R; dx <= R; ++dx) {
                    const double Pp = P[idx(dx, dy)];
                    if (Pp <= 0.0) continue;
                    const int64_t x = rh.cx + dx, y = rh.cy + dy;
                    if (x < 0 || y < 0 || x >= static_cast<int64_t>(xpixel) || y >= static_cast<int64_t>(ypixel)) continue;
                    const auto px = ptr[y * xpixel + x];
                    if (px == special || px == special + 1 || px == saturation || px == saturation - 1) continue;
                    const double v = B + std::max(0.0, I) * Pp;
                    num += Pp * (static_cast<double>(px) - rh.bkg) / v;
                    den += Pp * Pp / v;
                }
            if (den > 0.0) I = num / den; else break;
        }
        if (!(den > 0.0)) continue;

        Reflection refl = predicted[i];
        refl.I = static_cast<float>(I);
        refl.sigma = static_cast<float>(std::sqrt(1.0 / den));
        refl.bkg = static_cast<float>(rh.bkg);
        refl.observed = true;
        if (pol)
            refl.rlp /= geom.CalcAzIntPolarizationCorr(refl.predicted_x, refl.predicted_y, pol.value());
        refl.image_scale_corr = refl.rlp / refl.partiality;
        refl.image_number = static_cast<float>(image_number);
        out.push_back(refl);
    }
    return out;
}

}   // namespace

std::vector<Reflection> ProfileIntegrate2D(const DiffractionExperiment &experiment,
                                           const CompressedImage &image,
                                           const std::vector<Reflection> &predicted, size_t npredicted,
                                           int64_t image_number) {
    if (image.GetCompressedSize() == 0 || predicted.empty())
        return {};
    switch (image.GetMode()) {
        case CompressedImageMode::Int8:
            return ProfileIntegrateInternal<int8_t>(experiment, image, predicted, npredicted, INT8_MIN, INT8_MAX, image_number);
        case CompressedImageMode::Int16:
            return ProfileIntegrateInternal<int16_t>(experiment, image, predicted, npredicted, INT16_MIN, INT16_MAX, image_number);
        case CompressedImageMode::Int32:
            return ProfileIntegrateInternal<int32_t>(experiment, image, predicted, npredicted, INT32_MIN, INT32_MAX, image_number);
        case CompressedImageMode::Uint8:
            return ProfileIntegrateInternal<uint8_t>(experiment, image, predicted, npredicted, UINT8_MAX, UINT8_MAX, image_number);
        case CompressedImageMode::Uint16:
            return ProfileIntegrateInternal<uint16_t>(experiment, image, predicted, npredicted, UINT16_MAX, UINT16_MAX, image_number);
        case CompressedImageMode::Uint32:
            return ProfileIntegrateInternal<uint32_t>(experiment, image, predicted, npredicted, UINT32_MAX, UINT32_MAX, image_number);
        default:
            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Image mode not supported");
    }
}