Jungfraujoch/tests/ImageSpotFinderCPUTest.cpp

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

#include <catch2/catch_all.hpp>
#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h"

TEST_CASE("ImageSpotFinderCPU_SignalToNoise") {
    size_t width = 100;
    size_t height = 100;

    ImageSpotFinderCPU s(width, height);

    auto &input = s.GetInputBuffer();
    for (int i = 0; i < width * height; i++)
        input[i] = (i % 2) * 5 + 5;

    input[width * 50 + 50] = 20;
    input[width * 25 + 26] = 16;
    input[width * 75 + 25] = 12;

    // Mean is 7.5, std. dev is 2.5
    SpotFindingSettings settings{
            .signal_to_noise_threshold = 3.0,
            .photon_count_threshold = 0,
            .min_pix_per_spot = 1,
            .max_pix_per_spot = 20,
            .high_resolution_limit = 0.5,
            .low_resolution_limit = 3.0,
    };

    std::vector<bool> mask_resolution(width * height, false);
    auto spots = s.Run(settings, mask_resolution);

    REQUIRE(spots.size() == 2);
    REQUIRE(spots[0].RawCoord().y == 25);
    REQUIRE(spots[1].RawCoord().y == 50);
}

TEST_CASE("ImageSpotFinderCPU_SignalToNoise_Resolution") {
    size_t width = 100;
    size_t height = 100;

    ImageSpotFinderCPU s(width, height);

    auto &input = s.GetInputBuffer();
    for (int i = 0; i < width * height; i++)
        input[i] = (i % 2) * 5 + 5;

    input[width * 50 + 50] = 20;
    input[width * 25 + 26] = 16;
    input[width * 75 + 25] = 12;

    // Mean is 7.5, std. dev is 2.5
    SpotFindingSettings settings{
            .signal_to_noise_threshold = 3.0,
            .photon_count_threshold = 0,
            .min_pix_per_spot = 1,
            .max_pix_per_spot = 20,
            .high_resolution_limit = 1.5,
            .low_resolution_limit = 3.0,
    };

    std::vector<bool> mask_resolution(width * height, false);
    mask_resolution[width * 50 + 50] = true;
    auto spots = s.Run(settings, mask_resolution);

    REQUIRE(spots.size() == 1);
    REQUIRE(spots[0].RawCoord().x == 26);
    REQUIRE(spots[0].RawCoord().y == 25);
}

TEST_CASE("ImageSpotFinderCPU_CountThreshold_Resolution") {
    size_t width = 100;
    size_t height = 100;

    ImageSpotFinderCPU s(width, height);

    auto &input = s.GetInputBuffer();
    for (int i = 0; i < width * height; i++)
        input[i] = (i % 2) * 5 + 5;

    input[width * 50 + 50] = 20;
    input[width * 25 + 26] = 16;
    input[width * 75 + 25] = 12;

    // Mean is 7.5, std. dev is 2.5
    SpotFindingSettings settings{
            .signal_to_noise_threshold = -1,
            .photon_count_threshold = 11,
            .min_pix_per_spot = 1,
            .max_pix_per_spot = 20,
            .high_resolution_limit = 1.5,
            .low_resolution_limit = 3.0,
    };

    std::vector<bool> mask_resolution(width * height, false);
    mask_resolution[width * 50 + 50] = true;
    auto spots = s.Run(settings, mask_resolution);

    REQUIRE(spots.size() == 2);
    REQUIRE(spots[0].RawCoord().y == 25);
    REQUIRE(spots[1].RawCoord().y == 75);
}

TEST_CASE("ImageSpotFinderCPU_CountThreshold_Mask") {
    size_t width = 100;
    size_t height = 100;

    ImageSpotFinderCPU s( width, height);

    auto &input = s.GetInputBuffer();
    for (int i = 0; i < width * height; i++)
        input[i] = (i % 2) * 5 + 5;

    input[width * 50 + 50] = INT32_MIN;
    input[width * 50 + 51] = INT32_MAX;
    input[width * 25 + 26] = 16;
    input[width * 75 + 25] = 12;

    // Mean is 7.5, std. dev is 2.5
    SpotFindingSettings settings{
        .signal_to_noise_threshold = -1,
        .photon_count_threshold = 11,
        .min_pix_per_spot = 1,
        .max_pix_per_spot = 20,
        .high_resolution_limit = 1.5,
        .low_resolution_limit = 3.0,
    };

    std::vector<bool> mask_resolution(width * height, false);
    auto spots = s.Run(settings, mask_resolution);

    REQUIRE(spots.size() == 3);
    REQUIRE(spots[0].RawCoord().x == 26);
    REQUIRE(spots[0].RawCoord().y == 25);
    REQUIRE(spots[1].RawCoord().x == 51);
    REQUIRE(spots[1].RawCoord().y == 50);
    REQUIRE(spots[2].RawCoord().x == 25);
    REQUIRE(spots[2].RawCoord().y == 75);
}

TEST_CASE("ImageSpotFinderCPU_SignalToNoise_Mask") {
    size_t width = 100;
    size_t height = 100;
    ImageSpotFinderCPU s(width, height);

    auto &input = s.GetInputBuffer();
    for (int i = 0; i < width * height; i++)
        input[i] = (i % 2) * 5 + 5;

    input[width * 25 + 25] = INT32_MIN;
    input[width * 25 + 26] = 16;
    input[width * 74 + 25] = INT32_MIN;
    input[width * 74 + 26] = INT32_MIN;
    input[width * 75 + 25] = 15;
    input[width * 76 + 27] = INT32_MAX;

    // Mean is 7.5, std. dev is 2.5
    SpotFindingSettings settings{
        .signal_to_noise_threshold = 3.0,
        .photon_count_threshold = 0,
        .min_pix_per_spot = 1,
        .max_pix_per_spot = 20,
        .high_resolution_limit = 1.5,
        .low_resolution_limit = 3.0,
};

    std::vector<bool> mask_resolution(width * height, false);
    auto spots = s.Run(settings, mask_resolution);

    REQUIRE(spots.size() == 3);
    REQUIRE(spots[0].RawCoord().x == 26);
    REQUIRE(spots[0].RawCoord().y == 25);
    REQUIRE(spots[1].RawCoord().x == 25);
    REQUIRE(spots[1].RawCoord().y == 75);
    REQUIRE(spots[2].RawCoord().x == 27);
    REQUIRE(spots[2].RawCoord().y == 76);
}