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

#include <cmath>
#include <limits>

#include "AzIntEngine.h"

AzIntEngine::AzIntEngine(const AzimuthalIntegrationMapping &integration)
: integration(integration),
azint_bins(integration.GetBinNumber()),
npixel(integration.GetPixelToBin().size()),
azint_sum(integration.GetBinNumber(), 0.0f),
azint_sum2(integration.GetBinNumber(), 0.0f),
azint_count(integration.GetBinNumber(), 0u),
clip_nsigma(integration.Settings().GetSigmaClip()),
clip_lo(clip_nsigma > 0.0f ? integration.GetBinNumber() : 0),
clip_hi(clip_nsigma > 0.0f ? integration.GetBinNumber() : 0) {}

void AzIntEngine::UpdateClipLimits() {
    constexpr float inf = std::numeric_limits<float>::infinity();
    for (int i = 0; i < azint_bins; ++i) {
        // A handful of pixels have no spread worth the name, and clipping on it would reject most of
        // them. Leave such a bin alone - it is the detector edge and the beam stop, not signal.
        if (azint_count[i] < 8) {
            clip_lo[i] = -inf;
            clip_hi[i] = inf;
            continue;
        }
        const double n = azint_count[i];
        const double mean = azint_sum[i] / n;
        // sum2/n - mean^2 is the exact variance of what was accumulated; it can go slightly negative
        // through float cancellation on a bin whose pixels are all but identical, hence the floor.
        const double sd = std::sqrt(std::max(0.0, azint_sum2[i] / n - mean * mean));
        clip_lo[i] = static_cast<float>(mean - clip_nsigma * sd);
        clip_hi[i] = static_cast<float>(mean + clip_nsigma * sd);
    }
}