Jungfraujoch/tests/SpotFootprintTest.cpp

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

#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>

#include <cmath>
#include <vector>

#include "../image_analysis/bragg_integration/SpotFootprint.h"

// Spots drawn as Gaussians elongated along and across the radius are measured back at their widths,
// whatever their azimuth, and tabulated by distance from the beam.
TEST_CASE("SpotFootprint_MeasuresRadialAndTangentialWidths", "[Integration][portable]") {
    const int W = 1200, H = 1200;
    const float bx = 600.0f, by = 600.0f;
    std::vector<int32_t> img(static_cast<size_t>(W) * H, 10);   // flat background
    std::vector<float> xs, ys;
    // Width grows with distance: sigma_rad = 1 + r/200, sigma_tan = 1 + r/100.
    for (int k = 0; k < 400; ++k) {
        const float r = 60.0f + 480.0f * static_cast<float>(k % 20) / 20.0f;
        const float phi = 0.61f * static_cast<float>(k);
        const float x = bx + r * std::cos(phi), y = by + r * std::sin(phi);
        const float ux = std::cos(phi), uy = std::sin(phi);
        const float sr = 1.0f + r / 200.0f, st = 1.0f + r / 100.0f;
        for (int py = static_cast<int>(y) - 30; py <= static_cast<int>(y) + 30; ++py)
            for (int px = static_cast<int>(x) - 30; px <= static_cast<int>(x) + 30; ++px) {
                if (px < 0 || py < 0 || px >= W || py >= H) continue;
                const float dx = px - x, dy = py - y;
                const float rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux;
                img[static_cast<size_t>(py) * W + px] += static_cast<int32_t>(
                        std::lround(2000.0f * std::exp(-rad * rad / (2 * sr * sr) - tn * tn / (2 * st * st))));
            }
        xs.push_back(x);
        ys.push_back(y);
    }
    // Overlapping spots are not what this checks: keep those far from every other one.
    std::vector<float> kx, ky;
    for (size_t i = 0; i < xs.size(); ++i) {
        bool alone = true;
        for (size_t j = 0; j < xs.size(); ++j)
            if (i != j && std::hypot(xs[i] - xs[j], ys[i] - ys[j]) < 45.0f) alone = false;
        if (alone) { kx.push_back(xs[i]); ky.push_back(ys[i]); }
    }
    REQUIRE(kx.size() > 40);
    std::vector<FootprintSpot> spots;
    MeasureFootprintSpots(img.data(), W, H, bx, by, kx, ky, spots);
    REQUIRE(spots.size() > 30);
    for (const auto &s : spots) {
        CHECK_THAT(s.sigma_rad, Catch::Matchers::WithinRel(1.0f + s.r_px / 200.0f, 0.12f));
        CHECK_THAT(s.sigma_tan, Catch::Matchers::WithinRel(1.0f + s.r_px / 100.0f, 0.12f));
    }
}

TEST_CASE("SpotFootprint_TableFillsSparseBinsFromNeighbours", "[Integration][portable]") {
    std::vector<FootprintSpot> spots;
    for (int i = 0; i < FOOTPRINT_MIN_SPOTS_PER_BIN; ++i) {
        spots.push_back({50.0f, 1.0f, 1.5f});    // bin 0
        spots.push_back({1150.0f, 3.0f, 4.0f});  // last bin
    }
    const SpotFootprint fp = FootprintFromSpots(spots, 1200.0f);
    REQUIRE(fp.sigma_rad.size() == static_cast<size_t>(FOOTPRINT_BINS));
    REQUIRE(fp.sigma_rad.front() == 1.0f);
    REQUIRE(fp.sigma_tan.back() == 4.0f);
    REQUIRE(fp.sigma_rad[2] == 1.0f);                    // nearer the first bin
    REQUIRE(fp.sigma_rad[FOOTPRINT_BINS - 3] == 3.0f);   // nearer the last
    REQUIRE(FootprintFromSpots({}, 1200.0f).empty());
}

// A reflection recorded as a doublet: two indexed spots on either side of the prediction, along the
// radius. The offsets are found against the reflection of the same hkl, and their mean square is
// added to the widths of the bin they fall in.
TEST_CASE("SpotFootprint_OffsetsFromPredictionWidenTheTable", "[Integration][portable]") {
    const float bx = 600.0f, by = 600.0f;
    std::vector<Reflection> refl(1);
    refl[0].h = 1; refl[0].k = 2; refl[0].l = 3;
    refl[0].predicted_x = bx + 1000.0f;   // on +x, so radial = x
    refl[0].predicted_y = by;
    std::vector<SpotToSave> spots;
    spots.push_back({.x = bx + 1004.0f, .y = by, .lattice = 0, .h = 1, .k = 2, .l = 3, .indexed = true});
    spots.push_back({.x = bx + 996.0f, .y = by, .lattice = 0, .h = 1, .k = 2, .l = 3, .indexed = true});
    spots.push_back({.x = bx + 1000.0f, .y = by + 30.0f, .lattice = 0, .h = 3, .k = 2, .l = 1, .indexed = true});   // no such prediction
    spots.push_back({.x = bx + 1000.0f, .y = by + 30.0f, .lattice = -1, .h = 1, .k = 2, .l = 3, .indexed = false});  // not indexed
    std::vector<FootprintOffset> offsets;
    MeasureFootprintOffsets(spots, refl, bx, by, offsets);
    REQUIRE(offsets.size() == 2);
    REQUIRE_THAT(offsets[0].off_rad, Catch::Matchers::WithinAbs(4.0, 1e-4));
    REQUIRE_THAT(offsets[0].off_tan, Catch::Matchers::WithinAbs(0.0, 1e-4));

    std::vector<FootprintOffset> pool;
    for (int i = 0; i < FOOTPRINT_MIN_SPOTS_PER_BIN; ++i)
        pool.insert(pool.end(), offsets.begin(), offsets.end());
    SpotFootprint widths;
    widths.bin_px = 100.0f;
    widths.sigma_rad.assign(12, 1.0f);
    widths.sigma_tan.assign(12, 2.0f);
    const SpotFootprint fp = FootprintWithOffsets(widths, pool);
    for (int b = 0; b < 12; ++b) {   // one filled bin (10) serves them all
        REQUIRE_THAT(fp.sigma_rad[b], Catch::Matchers::WithinAbs(std::sqrt(17.0), 1e-4));
        REQUIRE_THAT(fp.sigma_tan[b], Catch::Matchers::WithinAbs(2.0, 1e-4));
    }
    REQUIRE(FootprintWithOffsets(widths, {}).sigma_rad[3] == 1.0f);
    REQUIRE(FootprintWithOffsets(SpotFootprint{}, pool).empty());
}