Jungfraujoch/common/AzimuthalIntegrationProfile.cpp

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

#include "AzimuthalIntegrationProfile.h"
#include "JFJochException.h"
#include "Definitions.h"

#include <algorithm>

inline float sum_to_count(float sum, uint64_t count) {
    if (count == 0)
        return NAN;
    return sum / (static_cast<float>(count));
}

inline float calc_std(float sum, float sum2, uint64_t count) {
    if (count == 0 || count == 1)
        return NAN;
    const auto fp_count = static_cast<float>(count);
    const float variance = (sum2 - sum * sum / fp_count) / (fp_count - 1);
    return std::sqrt(variance);
}

AzimuthalIntegrationProfile::AzimuthalIntegrationProfile(const AzimuthalIntegrationMapping &mapping)
    : sum(mapping.GetBinNumber(), 0),
      sum2(mapping.GetBinNumber(), 0),
      count(mapping.GetBinNumber(), 0),
      bin_to_q(mapping.GetBinToQ()),
      bin_to_d(mapping.GetBinToD()),
      bin_to_2theta(mapping.GetBinToTwoTheta()),
      bin_to_phi(mapping.GetBinToPhi()),
      q_bins(mapping.GetQBinCount()),
      azim_bins(mapping.GetAzimuthalBinCount()) {
}

void AzimuthalIntegrationProfile::Clear(const AzimuthalIntegrationMapping &mapping) {
    std::unique_lock ul(m);

    bin_to_d = mapping.GetBinToD();
    bin_to_q = mapping.GetBinToQ();
    bin_to_2theta = mapping.GetBinToTwoTheta();
    bin_to_phi = mapping.GetBinToPhi();
    q_bins = mapping.GetQBinCount();
    azim_bins = mapping.GetAzimuthalBinCount();

    sum = std::vector<float>(mapping.GetBinNumber(), 0);
    count = std::vector<uint64_t>(mapping.GetBinNumber(), 0);
}

void AzimuthalIntegrationProfile::Add(int64_t bin, int64_t value) {
    if (bin < 0 || bin >= sum.size())
        return;
    std::unique_lock ul(m);
    sum[bin] += static_cast<float>(value);
    sum2[bin] += static_cast<float>(value * value);
    count[bin]++;
}

void AzimuthalIntegrationProfile::Add(const std::vector<float> &in_sum,
                                      const std::vector<float> &in_sum2,
                                      const std::vector<uint32_t> &in_count) {
    std::unique_lock ul(m);
    if ((in_sum.size() == sum.size()) && (in_count.size() == count.size())) {
        for (int i = 0; i < sum.size(); i++) {
            sum[i] += in_sum[i];
            count[i] += in_count[i];
        }
        if (in_sum2.size() == sum2.size()) {
            for (int i = 0; i < sum.size(); i++) {
                sum2[i] += in_sum2[i];
            }
        }
    } else if (!in_sum.empty() && !in_count.empty())
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Mismatch in size of sum/count datasets");
}

std::vector<float> AzimuthalIntegrationProfile::GetResult() const {
    std::unique_lock ul(m);

    std::vector<float> rad_int_profile(sum.size(), 0);
    for (int i = 0; i < sum.size(); i++)
        rad_int_profile[i] = sum_to_count(sum[i], count[i]);
    return rad_int_profile;
}

std::vector<float> AzimuthalIntegrationProfile::GetStd() const {
    std::unique_lock ul(m);

    std::vector<float> rad_int_profile(sum.size(), 0);
    for (int i = 0; i < sum.size(); i++)
        rad_int_profile[i] = calc_std(sum[i], sum2[i], count[i]);
    return rad_int_profile;
}

std::vector<uint64_t> AzimuthalIntegrationProfile::GetPixelCount() const {
    std::unique_lock ul(m);
    return count;
}

std::vector<float> AzimuthalIntegrationProfile::GetResult1D() const {
    std::unique_lock ul(m);

    std::vector<float> sum_q(q_bins, 0.0f);
    std::vector<uint64_t> count_q(q_bins, 0);

    for (int i = 0; i < sum.size(); i++) {
        const int q_bin = i % q_bins;
        sum_q[q_bin] += sum[i];
        count_q[q_bin] += count[i];
    }

    std::vector<float> rad_int_profile(q_bins, 0.0f);
    for (int q = 0; q < q_bins; q++)
        rad_int_profile[q] = sum_to_count(sum_q[q], count_q[q]);

    return rad_int_profile;
}
void AzimuthalIntegrationProfile::SetTitle(const std::string &input) {
    title = input;
}

const std::vector<float> &AzimuthalIntegrationProfile::GetXAxis(PlotAzintUnit unit) const {
    switch (unit) {
        case PlotAzintUnit::TwoTheta_deg:
            return bin_to_2theta;
        case PlotAzintUnit::d_A:
            return bin_to_d;
        default:
        case PlotAzintUnit::Q_recipA:
            return bin_to_q;
    }
}

MultiLinePlot AzimuthalIntegrationProfile::GetPlot(bool force_1d, PlotAzintUnit plot_unit) const {
    MultiLinePlot ret;

    const std::vector<float> &x_coord = GetXAxis(plot_unit);

    if (azim_bins == 1)
        ret.AddPlot(MultiLinePlotStruct{.title = title, .x = x_coord, .y = GetResult()});
    else {
        if (force_1d) {
            std::vector<float> x_shortened(q_bins);
            for (int i = 0; i < q_bins; i++)
                x_shortened[i] = x_coord[i];
            ret.AddPlot(MultiLinePlotStruct{.title = title, .x = x_shortened, .y = GetResult1D()});
        } else {
            ret.AddPlot(MultiLinePlotStruct{.title = title, .x = x_coord, .y= bin_to_phi, .z = GetResult()});
        }
    }
    return ret;
}

float AzimuthalIntegrationProfile::GetMeanValueOfBins(uint16_t min_bin, uint16_t max_bin) const {
    std::unique_lock ul(m);

    float ret_sum = 0;
    uint64_t ret_count = 0;

    for (int i = 0; i < sum.size(); i++) {
        uint16_t q_bin = i % q_bins;
        if (q_bin >= min_bin && q_bin <= max_bin) {
            ret_sum += sum[i];
            ret_count += count[i];
        }
    }

    return sum_to_count(ret_sum, ret_count);
}

float AzimuthalIntegrationProfile::GetBkgEstimate(const AzimuthalIntegrationSettings &settings) const {
    auto min_bin = settings.QToBin(settings.GetBkgEstimateLowQ_recipA());
    auto max_bin = settings.QToBin(settings.GetBkgEstimateHighQ_recipA());
    return GetMeanValueOfBins(min_bin, max_bin);
}

float AzimuthalIntegrationProfile::GetIceRingScore(const AzimuthalIntegrationSettings &settings,
                                                   float half_width_q) const {
    // For each hexagonal-ice powder ring, the mean profile intensity in the ring band (+/- half_width in
    // q = 2*pi/d) over a baseline from the two shoulders just outside it. Report the strongest ring's
    // ratio (1 = no ice); rings whose band+shoulders fall off the measured q-range are skipped.
    constexpr float two_pi = 6.283185307f;
    const float low_q = settings.GetLowQ_recipA();
    const float high_q = settings.GetHighQ_recipA();
    float score = 1.0f;
    for (const float ice_d : ICE_RING_RES_A) {
        const float q = two_pi / ice_d;
        if (q - 2 * half_width_q < low_q || q + 2 * half_width_q > high_q)
            continue;
        const float ring = GetMeanValueOfBins(settings.QToBin(q - half_width_q), settings.QToBin(q + half_width_q));
        const float lo = GetMeanValueOfBins(settings.QToBin(q - 2 * half_width_q), settings.QToBin(q - half_width_q));
        const float hi = GetMeanValueOfBins(settings.QToBin(q + half_width_q), settings.QToBin(q + 2 * half_width_q));
        const float baseline = 0.5f * (lo + hi);
        if (std::isfinite(baseline) && baseline > 0.0f && std::isfinite(ring))
            score = std::max(score, ring / baseline);
    }
    return score;
}

AzimuthalIntegrationProfile &AzimuthalIntegrationProfile::operator+=(const AzimuthalIntegrationProfile &other) {
    if ((other.bin_to_q != bin_to_q) || (sum.size() != other.sum.size())) {
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Error combining two radial integration profiles");
    }

    for (int i = 0; i < sum.size(); i++) {
        sum[i] += other.sum[i];
        sum2[i] += other.sum2[i];
        count[i] += other.count[i];
    }

    return *this;
}

void AzimuthalIntegrationProfile::Add(const DeviceOutput &result) {
    std::unique_lock ul(m);

    if (sum.size() > FPGA_INTEGRATION_BIN_COUNT )
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Error in getting result from FPGA");

    for (int i = 0; i < sum.size(); i++) {
        sum[i] += result.integration_result[i].sum;
        count[i] += result.integration_result[i].count;
    }
}