Jungfraujoch/image_analysis/scale_merge/IceRingMask.cpp

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

#include <algorithm>
#include <cmath>

#include "IceRingMask.h"
#include "../../common/Definitions.h"        // ICE_RING_RES_A
#include "../../common/CorrelationCoefficient.h"

std::vector<char> FindDecorrelatedIceRings(const std::vector<MergedReflection> &merged,
                                           float half_width_q_recipA, Logger &logger) {
    if (merged.empty())
        return {};

    constexpr float two_pi = 6.283185307f;
    const float w = half_width_q_recipA;
    std::vector<char> mask(ICE_RING_RES_A.size(), 0);

    for (size_t i = 0; i < ICE_RING_RES_A.size(); ++i) {
        const float q_ring = two_pi / ICE_RING_RES_A[i];
        CorrelationCoefficient ring, shoulder;
        size_t n_ring = 0, n_shoulder = 0;
        for (const auto &m : merged) {
            if (!(m.d > 0.0f) || !std::isfinite(m.I_half[0]) || !std::isfinite(m.I_half[1]))
                continue;
            const float dq = std::fabs(two_pi / m.d - q_ring);
            if (dq < w) { ring.Add(m.I_half[0], m.I_half[1]); ++n_ring; }
            else if (dq < 3.0f * w) { shoulder.Add(m.I_half[0], m.I_half[1]); ++n_shoulder; }
        }
        if (n_ring >= 20 && n_shoulder >= 20 && shoulder.GetCC() > 0.5
            && ring.GetCC() < shoulder.GetCC() - 0.05) {
            mask[i] = 1;
            logger.Info("Ice-ring mask: {:.2f} A ring CC1/2 {:.3f} << shoulders {:.3f}; masked from merge",
                        ICE_RING_RES_A[i], ring.GetCC(), shoulder.GetCC());
        }
    }

    if (std::none_of(mask.begin(), mask.end(), [](char c) { return c != 0; }))
        return {};
    return mask;
}