Jungfraujoch/tests/AdaptiveSpotFinderCPUTest.cpp

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

#include <cmath>
#include <algorithm>

#include <catch2/catch_all.hpp>

#include "../common/AzimuthalIntegrationMapping.h"
#include "../image_analysis/spot_finding/AdaptiveSpotFinderCPU.h"

namespace {

SpotFindingSettings AdaptiveSettings() {
    SpotFindingSettings s{};
    s.adaptive_threshold = true;
    s.false_pixels_per_frame = 100.0f;
    s.min_pix_per_spot = 1;
    s.max_pix_per_spot = 50;
    s.high_resolution_limit = 0.0f;
    s.low_resolution_limit = 1.0e6f;
    s.high_res_gap_Q_recipA = std::nullopt;
    return s;
}

} // namespace

// Raw (untransformed) geometry: the mapping is built over the raw module layout, which is smaller
// than the converted image. The finder has to walk the raw image, so a spot planted at a raw pixel
// comes back at that pixel.
TEST_CASE("AdaptiveSpotFinderCPU_RawGeometry", "[AdaptiveSpotFinder][portable]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080);
    x.QSpacingForAzimInt_recipA(0.05).QRangeForAzimInt_recipA(0.05, 5.0);
    x.GeometryTransformation(false);

    PixelMask pixel_mask(x);
    AzimuthalIntegrationMapping mapping(x, pixel_mask);

    const size_t w = x.GetXPixelsNum();
    const size_t h = x.GetYPixelsNum();
    REQUIRE(w * h == mapping.GetPixelToBin().size());
    REQUIRE(w * h < static_cast<size_t>(x.GetPixelsNumConv()));

    ImagePreprocessorBuffer buffer(x.GetPixelsNum());
    for (size_t i = 0; i < w * h; i++)
        buffer[i] = 8 + static_cast<int32_t>(i % 5);   // background 8..12

    // A 3x3 blob on a pixel that has a ring - with all of its neighbours on one too.
    const auto &pixel_to_bin = mapping.GetPixelToBin();
    size_t spot_row = 0, spot_col = 0;
    for (size_t row = 100; row < h - 100 && spot_row == 0; row++) {
        for (size_t col = 100; col < w - 100; col++) {
            bool all_binned = true;
            for (int dr = -1; dr <= 1; dr++)
                for (int dc = -1; dc <= 1; dc++)
                    all_binned &= pixel_to_bin[(row + dr) * w + col + dc] != UINT16_MAX;
            if (all_binned) {
                spot_row = row;
                spot_col = col;
                break;
            }
        }
    }
    REQUIRE(spot_row > 0);

    for (int dr = -1; dr <= 1; dr++)
        for (int dc = -1; dc <= 1; dc++)
            buffer[(spot_row + dr) * w + spot_col + dc] = 200;

    std::vector<bool> res_mask(x.GetPixelsNum(), false);
    AdaptiveSpotFinderCPU finder(mapping);
    finder.SetResolutionMask(res_mask);
    const auto spots = finder.Run(buffer, AdaptiveSettings());

    REQUIRE(spots.size() == 1);
    CHECK(std::lround(spots[0].RawCoord().x) == static_cast<long>(spot_col));
    CHECK(std::lround(spots[0].RawCoord().y) == static_cast<long>(spot_row));
}

// The property the whole engine exists for: the threshold comes from the image's OWN noise, so the
// same settings behave the same way on a frame whose background is ten times higher. A frame is built
// with background spread S around a mean, one pixel planted a few S above it (must stay unfound) and
// one planted far above (must be found); then the identical frame scaled by ten must give the identical
// answer. Any threshold that does not track the background - a constant, or one that drops the sigma
// term - finds the weak pixel in the scaled frame, or loses the strong one.
TEST_CASE("AdaptiveSpotFinderCPU_ThresholdTracksBackground", "[AdaptiveSpotFinder]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080);
    x.QSpacingForAzimInt_recipA(0.05).QRangeForAzimInt_recipA(0.05, 5.0);
    x.GeometryTransformation(false);

    PixelMask pixel_mask(x);
    AzimuthalIntegrationMapping mapping(x, pixel_mask);
    const auto &pixel_to_bin = mapping.GetPixelToBin();

    const size_t w = x.GetXPixelsNum();
    const size_t h = x.GetYPixelsNum();

    // Two well-separated pixels that carry a ring, so both are seen by the finder.
    std::vector<size_t> planted;
    for (size_t row = 300; row < h - 300 && planted.size() < 2; row += 137)
        for (size_t col = 300; col < w - 300; col += 149)
            if (pixel_to_bin[row * w + col] != UINT16_MAX) {
                planted.push_back(row * w + col);
                break;
            }
    REQUIRE(planted.size() == 2);

    // Background takes 5 evenly spaced levels one S apart, i.e. mean + 2S and sigma = sqrt(2) S. With
    // ~100 expected noise pixels per frame the cut lands near mean + 4.1 sigma = mean + 5.8 S.
    const auto run_at_scale = [&](int32_t scale) {
        ImagePreprocessorBuffer buffer(x.GetPixelsNum());
        for (size_t i = 0; i < w * h; i++)
            buffer[i] = scale * (10 + static_cast<int32_t>(i % 5));
        buffer[planted[0]] = scale * (10 + 5);    // mean + 3 S: below the cut
        buffer[planted[1]] = scale * (10 + 30);   // mean + 28 S: well above it
        std::vector<bool> res_mask(x.GetPixelsNum(), false);
        AdaptiveSpotFinderCPU finder(mapping);
        finder.SetResolutionMask(res_mask);
        return finder.Run(buffer, AdaptiveSettings());
    };

    const auto plain = run_at_scale(1);
    const auto scaled = run_at_scale(10);

    REQUIRE(plain.size() == 1);
    CHECK(std::lround(plain[0].RawCoord().x) == static_cast<long>(planted[1] % w));
    CHECK(std::lround(plain[0].RawCoord().y) == static_cast<long>(planted[1] / w));

    // Ten times the background, ten times the noise, ten times the signal - same answer.
    REQUIRE(scaled.size() == plain.size());
    CHECK(std::lround(scaled[0].RawCoord().x) == std::lround(plain[0].RawCoord().x));
    CHECK(std::lround(scaled[0].RawCoord().y) == std::lround(plain[0].RawCoord().y));
}

// The per-ring background the finder hands out is the peak-EXCLUDED one: a handful of very bright
// pixels planted on a ring must not move it, which is the property that lets the ice score be read
// off it instead of off the plain azimuthal profile. Rings with too few pixels to be their own
// background come back as NaN rather than as a stale value from the previous frame.
TEST_CASE("AdaptiveSpotFinderCPU_RingBackgroundExcludesPeaks", "[AdaptiveSpotFinder]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080);
    x.QSpacingForAzimInt_recipA(0.05).QRangeForAzimInt_recipA(0.05, 5.0);
    x.GeometryTransformation(false);

    PixelMask pixel_mask(x);
    AzimuthalIntegrationMapping mapping(x, pixel_mask);
    const auto &pixel_to_bin = mapping.GetPixelToBin();

    const size_t w = x.GetXPixelsNum();
    const size_t h = x.GetYPixelsNum();

    ImagePreprocessorBuffer buffer(x.GetPixelsNum());
    for (size_t i = 0; i < w * h; i++)
        buffer[i] = 8 + static_cast<int32_t>(i % 5);   // background 8..12

    std::vector<bool> res_mask(x.GetPixelsNum(), false);
    AdaptiveSpotFinderCPU finder(mapping);
    finder.SetResolutionMask(res_mask);

    finder.Detect(buffer, AdaptiveSettings());
    const std::vector<float> clean = finder.GetRingBackground();
    REQUIRE(clean.size() == mapping.GetBinNumber());

    std::vector<int64_t> pixels_per_bin(mapping.GetBinNumber(), 0);
    for (size_t i = 0; i < w * h; i++)
        if (pixel_to_bin[i] < mapping.GetBinNumber())
            pixels_per_bin[pixel_to_bin[i]]++;

    // Plant 200 bright pixels spread over one well-populated ring - far more than a real spot, so a
    // plain mean would move visibly.
    const uint16_t ring = static_cast<uint16_t>(
        std::max_element(pixels_per_bin.begin(), pixels_per_bin.end()) - pixels_per_bin.begin());
    REQUIRE(pixels_per_bin[ring] > 10000);
    int planted = 0;
    for (size_t i = 0; i < w * h && planted < 200; i++)
        if (pixel_to_bin[i] == ring) {
            buffer[i] = 100000;
            planted++;
        }
    REQUIRE(planted == 200);

    finder.Detect(buffer, AdaptiveSettings());
    const std::vector<float> spiked = finder.GetRingBackground();

    REQUIRE(std::isfinite(clean[ring]));
    REQUIRE(std::isfinite(spiked[ring]));
    CHECK(spiked[ring] == Catch::Approx(clean[ring]).epsilon(0.01));

    for (size_t b = 0; b < clean.size(); b++)
        CHECK(std::isfinite(clean[b]) == (pixels_per_bin[b] >= 40));
}

// The footprint of a bright reflection must not grow with its brightness. A whole-ring threshold is
// an absolute contour, so the area a Gaussian puts above it grows as sigma^2 ln(peak/threshold) - the
// same spot detected a hundred times brighter comes back tens of pixels larger, and an upper bound on
// spot size becomes an upper bound on spot INTENSITY. Intersecting with the local-box test removes
// that: the spot inflates the box's own variance, the peak divides out of the acceptance test, and the
// contour lands at a fixed fraction of the spot's own height whatever that height is.
TEST_CASE("AdaptiveSpotFinderCPU_FootprintDoesNotGrowWithBrightness", "[AdaptiveSpotFinder]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080);
    x.QSpacingForAzimInt_recipA(0.05).QRangeForAzimInt_recipA(0.05, 5.0);
    x.GeometryTransformation(false);

    PixelMask pixel_mask(x);
    AzimuthalIntegrationMapping mapping(x, pixel_mask);
    const auto &pixel_to_bin = mapping.GetPixelToBin();

    const size_t w = x.GetXPixelsNum();
    const size_t h = x.GetYPixelsNum();

    size_t spot_row = 0, spot_col = 0;
    for (size_t row = 400; row < h - 400 && spot_row == 0; row++)
        for (size_t col = 400; col < w - 400; col++)
            if (pixel_to_bin[row * w + col] != UINT16_MAX) {
                spot_row = row;
                spot_col = col;
                break;
            }
    REQUIRE(spot_row > 0);

    std::vector<bool> res_mask(x.GetPixelsNum(), false);
    auto settings = AdaptiveSettings();
    settings.min_pix_per_spot = 2;
    settings.max_pix_per_spot = 100000;   // no bound, so the footprint itself is what is measured

    // One Gaussian of width 1.5 px on a flat background of 10, at three amplitudes a hundred apart.
    auto footprint = [&](double amplitude) {
        ImagePreprocessorBuffer buffer(x.GetPixelsNum());
        for (size_t i = 0; i < w * h; i++)
            buffer[i] = 10;
        constexpr double sigma = 1.5;
        for (int dr = -12; dr <= 12; dr++)
            for (int dc = -12; dc <= 12; dc++) {
                const double r2 = dr * dr + dc * dc;
                buffer[(spot_row + dr) * w + spot_col + dc] =
                    10 + static_cast<int32_t>(amplitude * std::exp(-r2 / (2 * sigma * sigma)));
            }
        AdaptiveSpotFinderCPU finder(mapping);
        finder.SetResolutionMask(res_mask);
        const auto spots = finder.Run(buffer, settings);
        REQUIRE(spots.size() == 1);
        return spots[0].PixelCount();
    };

    const int64_t small = footprint(300.0);
    const int64_t large = footprint(30000.0);
    CHECK(small > 0);
    // A hundredfold in peak is 4.6 e-folds. An absolute contour would add sigma^2 ln(100) ~ 10 pixels
    // per e-fold of AREA here, tens of pixels in all; a peak-relative one adds nothing.
    CHECK(large - small <= 4);
}