Jungfraujoch/tests/DetectorModulesTest.cpp

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

#include <catch2/catch_all.hpp>
#include <cmath>
#include <filesystem>

#include "../common/DetectorGeometryModular.h"
#include "../common/DetectorModules.h"
#include "../common/DiffractionExperiment.h"
#include "../common/JFJochMath.h"
#include "../common/PixelMask.h"
#include "../reader/NXmxModules.h"
#include "../writer/NXmxModules.h"

// Two JUNGFRAU modules side by side, 8 pixels apart: each a 1030 x 514 box of the assembled image, placed
// as imaged, so every pixel maps onto itself.
TEST_CASE("DetectorModules_JungfrauAsImaged", "[DetectorModules]") {
    const DetectorGeometryModular geom(2, 2, 8, 36, false);
    const DetectorModules modules = ModulesAsImaged(geom);
    REQUIRE(modules.size() == 2);
    CHECK(modules[0].x0 == 0);
    CHECK(modules[1].x0 == 1038);
    for (const auto &m: modules.Get()) {
        CHECK(m.width == 1030);
        CHECK(m.height == 514);
    }
    CHECK_FALSE(modules.Moved());

    CHECK(modules.ModuleOfImagePixel(1033, 100) == -1);   // in the gap
    CHECK(modules.ModuleOfImagePixel(1040, 100) == 1);
    const auto [px, py] = modules.ImageToPlane(1040.25, 100.5);
    CHECK(px == Catch::Approx(1040.25));
    CHECK(py == Catch::Approx(100.5));

    // Not modular: no modules.
    CHECK(ModulesAsImaged(DetectorGeometryModular(1)).empty());
}

// A module moved by a rigid step: the image -> plane map moves its pixels, and plane -> image undoes it.
TEST_CASE("DetectorModules_Apply", "[DetectorModules]") {
    DetectorModules modules({DetectorModule::AsImaged(0, 0, 100, 50), DetectorModule::AsImaged(108, 0, 100, 50)});
    const auto centres = modules.Centres();
    std::vector<ModuleShift> step = centres;
    step[1].dx = 1.5;
    step[1].dy = -0.25;
    step[1].rot_rad = 0.2 * PI / 180.0;
    modules.Apply(step);
    CHECK(modules.Moved());

    // Module 0 is where it was.
    const auto [ax, ay] = modules.ImageToPlane(10, 10);
    CHECK(ax == Catch::Approx(10));
    CHECK(ay == Catch::Approx(10));

    // Module 1's centre moved by the shift, a corner by the shift and the turn.
    const auto [cx, cy] = modules[1].PlaneCentre();
    CHECK(cx == Catch::Approx(centres[1].x0 + 1.5));
    CHECK(cy == Catch::Approx(centres[1].y0 - 0.25));
    const double c = std::cos(step[1].rot_rad), s = std::sin(step[1].rot_rad);
    const double ux = 108 - centres[1].x0, uy = 0 - centres[1].y0;
    const auto [qx, qy] = modules.ImageToPlane(108, 0);
    CHECK(qx == Catch::Approx(centres[1].x0 + c * ux - s * uy + 1.5));
    CHECK(qy == Catch::Approx(centres[1].y0 + s * ux + c * uy - 0.25));

    // And back: a plane point inside module 1's footprint goes to the image pixel that put it there.
    const auto [ix, iy] = modules.PlaneToImage(qx, qy);
    CHECK(ix == Catch::Approx(108).margin(1e-9));
    CHECK(iy == Catch::Approx(0).margin(1e-9));
    const auto [jx, jy] = modules.PlaneToImage(modules.ImageToPlane(150.3, 20.7).first,
                                               modules.ImageToPlane(150.3, 20.7).second);
    CHECK(jx == Catch::Approx(150.3));
    CHECK(jy == Catch::Approx(20.7));
}

// Written as NXdetector_module groups and read back, the modules come back as they were - here two
// JUNGFRAU modules, one of them moved off the image's grid, on a detector stored mirrored.
TEST_CASE("DetectorModules_NXmxRoundTrip", "[DetectorModules]") {
    DetectorModules modules = ModulesAsImaged(DetectorGeometryModular(2, 2, 8, 36, false));
    auto step = modules.Centres();
    step[0].dx = -0.7;
    step[0].dy = 0.35;
    step[0].rot_rad = -0.05 * PI / 180.0;
    modules.Apply(step);

    DiffractionGeometry geom;
    geom.BeamX_pxl(1100).BeamY_pxl(250).DetectorDistance_mm(120).PixelSize_mm(0.075)
        .Orientation(DetectorOrientation(true, 1));

    const std::string path = "detector_modules_roundtrip.h5";
    WriteModuleGeometryFile(path, modules, geom);
    DetectorModules back;
    {
        HDF5ReadOnlyFile file(path);
        back = ReadNXmxModules(file, 0.075e-3, geom.GetOrientation());
    }
    std::filesystem::remove(path);

    REQUIRE(back.size() == 2);
    for (size_t i = 0; i < 2; i++) {
        const auto &a = modules[i], &b = back[i];
        CHECK(b.x0 == a.x0);
        CHECK(b.y0 == a.y0);
        CHECK(b.width == a.width);
        CHECK(b.height == a.height);
        CHECK(b.ox == Catch::Approx(a.ox).margin(1e-3));
        CHECK(b.oy == Catch::Approx(a.oy).margin(1e-3));
        CHECK(b.fx == Catch::Approx(a.fx).margin(1e-6));
        CHECK(b.fy == Catch::Approx(a.fy).margin(1e-6));
        CHECK(b.sx == Catch::Approx(a.sx).margin(1e-6));
        CHECK(b.sy == Catch::Approx(a.sy).margin(1e-6));
    }
}

// A JUNGFRAU file that marks only the gaps between its modules in the mask: the modules come back from
// the mask, one per gap-bounded region; a region that is not one whole module gives none at all.
TEST_CASE("DetectorModules_FromGapMask", "[DetectorModules]") {
    const DetectorGeometryModular geom(4, 2, 8, 36, false);
    const DetectorModules expected = ModulesAsImaged(geom);
    const int64_t width = geom.GetWidth(true), height = geom.GetHeight(true);
    std::vector<uint32_t> mask(width * height, 1u << PixelMask::ModuleGapPixelBit);
    expected.ForEachPixel(width, height, [&](int64_t, size_t i, double, double) { mask[i] = 0; });
    mask[5] |= 1u << PixelMask::ErrorPixelBit;   // other bits do not matter

    const DetectorModules modules = ModulesFromGapMask(mask, width, height);
    REQUIRE(modules.size() == expected.size());
    for (size_t i = 0; i < modules.size(); i++) {
        CHECK(modules[i].x0 == expected[i].x0);
        CHECK(modules[i].y0 == expected[i].y0);
        CHECK(modules[i].width == expected[i].width);
        CHECK(modules[i].height == expected[i].height);
    }

    // A gap pixel inside a module: no longer whole modules.
    auto holed = mask;
    holed[100 * width + 100] |= 1u << PixelMask::ModuleGapPixelBit;
    CHECK(ModulesFromGapMask(holed, width, height).empty());
    // No gaps at all: one region the size of the image, not a module.
    CHECK(ModulesFromGapMask(std::vector<uint32_t>(width * height, 0), width, height).empty());
}