Jungfraujoch/tests/RugnuxTest.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/DiffractionExperiment.h"
#include "../common/ScanResultGenerator.h"
#include "../writer/FileWriter.h"
#include "../reader/JFJochHDF5Reader.h"
#include "../rugnux/Rugnux.h"
#include "../rugnux/RugnuxCommandLine.h"
#include "../rugnux/SpotWidth.h"

namespace {
    // Write a small VDS dataset of `n` flat images and return nothing (prefix_master.h5 +
    // prefix_data_000001.h5 land in the test working directory).
    void WriteTestDataset(const std::string &prefix, int n) {
        RegisterHDF5Filter();

        DiffractionExperiment x(DetJF(1));
        x.FilePrefix(prefix).ImagesPerTrigger(n).OverwriteExistingFiles(true);
        x.BitDepthImage(16).ImagesPerFile(n).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true);
        x.Compression(CompressionAlgorithm::NO_COMPRESSION);
        x.BeamX_pxl(512).BeamY_pxl(256).DetectorDistance_mm(150).IncidentEnergy_keV(WVL_1A_IN_KEV)
                .FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));

        std::vector<int16_t> image(x.GetPixelsNum(), 5);
        StartMessage start_message;
        x.FillMessage(start_message);
        FileWriter file_set(start_message);
        ScanResultGenerator generator(x);
        for (int i = 0; i < n; i++) {
            DataMessage message{};
            message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
            message.number = i;
            REQUIRE_NOTHROW(file_set.WriteHDF5(message));
            generator.Add(message);
        }
        EndMessage end_message;
        end_message.max_image_number = n;
        generator.FillEndMessage(end_message);
        file_set.WriteHDF5(end_message);
        file_set.Finalize();
    }
}

TEST_CASE("Rugnux_AzInt", "[HDF5][Full]") {
    WriteTestDataset("process_azint_in", 8);

    JFJochHDF5Reader reader;
    REQUIRE_NOTHROW(reader.ReadFile("process_azint_in_master.h5"));
    auto dataset = reader.GetDataset();
    REQUIRE(dataset);

    ProcessConfig config;
    config.mode = ProcessMode::AzimuthalIntegration;
    config.nthreads = 2;
    config.output_prefix = "process_azint_out";

    Rugnux process(reader, dataset->experiment, *dataset->pixel_mask, config);
    ProcessResult result;
    REQUIRE_NOTHROW(result = process.Run());

    CHECK_FALSE(result.cancelled);
    CHECK(result.images_processed == 8);
    REQUIRE(result.written_master_path.has_value());

    {
        // The _process.h5 links back to the source images and carries an azimuthal profile per image.
        JFJochHDF5Reader out;
        REQUIRE_NOTHROW(out.ReadFile("process_azint_out_process.h5"));
        CHECK(out.GetNumberOfImages() == 8);
        std::shared_ptr<JFJochReaderImage> img;
        REQUIRE_NOTHROW(img = out.LoadImage(0));
        REQUIRE(img);
        CHECK_FALSE(img->ImageData().az_int_profile.empty());
    }

    reader.Close();
    remove("process_azint_in_master.h5");
    remove("process_azint_in_data_000001.h5");
    remove("process_azint_out_process.h5");
    REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}

TEST_CASE("Rugnux_NoOutput", "[HDF5][Full]") {
    WriteTestDataset("process_noout_in", 6);

    JFJochHDF5Reader reader;
    REQUIRE_NOTHROW(reader.ReadFile("process_noout_in_master.h5"));
    auto dataset = reader.GetDataset();

    // Empty output prefix => process without writing any file.
    ProcessConfig config;
    config.mode = ProcessMode::AzimuthalIntegration;
    config.nthreads = 3;

    Rugnux process(reader, dataset->experiment, *dataset->pixel_mask, config);
    auto result = process.Run();

    CHECK_FALSE(result.cancelled);
    CHECK(result.images_processed == 6);
    CHECK_FALSE(result.written_master_path.has_value());

    reader.Close();
    remove("process_noout_in_master.h5");
    remove("process_noout_in_data_000001.h5");
    REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}

TEST_CASE("Rugnux_Cancel", "[HDF5][Full]") {
    WriteTestDataset("process_cancel_in", 8);

    JFJochHDF5Reader reader;
    REQUIRE_NOTHROW(reader.ReadFile("process_cancel_in_master.h5"));
    auto dataset = reader.GetDataset();

    ProcessConfig config;
    config.mode = ProcessMode::AzimuthalIntegration;
    config.nthreads = 2;

    Rugnux process(reader, dataset->experiment, *dataset->pixel_mask, config);
    process.Cancel();   // cancel before running: the worker loop stops immediately
    auto result = process.Run();

    CHECK(result.cancelled);
    CHECK(result.images_processed == 0);
    CHECK_FALSE(result.written_master_path.has_value());

    reader.Close();
    remove("process_cancel_in_master.h5");
    remove("process_cancel_in_data_000001.h5");
    REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}

TEST_CASE("RugnuxCommandLine_Full", "[process]") {
    DiffractionExperiment x(DetJF(1));
    IndexingSettings idx;
    idx.Algorithm(IndexingAlgorithmEnum::FFT);
    idx.GeomRefinementAlgorithm(GeomRefinementAlgorithmEnum::BeamCenter);
    x.ImportIndexingSettings(idx);
    x.SpaceGroupNumber(96);

    ProcessConfig config;
    config.mode = ProcessMode::FullAnalysis;
    config.nthreads = 8;
    config.output_prefix = "run1";
    config.end_image = 500;
    config.rotation_indexing = true;
    config.two_pass_rotation = true;
    config.rotation_indexing_image_count = 30;
    config.spot_finding = DiffractionExperiment::DefaultDataProcessingSettings();

    const std::string cmd = RugnuxCommandLine(config, x, "/data/test_master.h5");
    CHECK(cmd.rfind("rugnux", 0) == 0);
    CHECK(cmd.find("-N 8") != std::string::npos);
    CHECK(cmd.find("-e 500") != std::string::npos);
    CHECK(cmd.find("-o run1") != std::string::npos);
    CHECK(cmd.find("-X fft") != std::string::npos);
    CHECK(cmd.find("-S 96") != std::string::npos);
    // -R takes an optional argument, so its value must be attached (-R30); a separate "-R 30" token
    // would not re-parse (getopt would leave 30 as a positional and drop the count).
    CHECK(cmd.find("-R30") != std::string::npos);
    CHECK(cmd.find("-R 30") == std::string::npos);
    CHECK(cmd.find("/data/test_master.h5") != std::string::npos);
}

TEST_CASE("RugnuxCommandLine_AzInt", "[process]") {
    DiffractionExperiment x(DetJF(1));
    AzimuthalIntegrationSettings a;
    a.AzimuthalBinCount(4);
    x.ImportAzimuthalIntegrationSettings(a);

    ProcessConfig config;
    config.mode = ProcessMode::AzimuthalIntegration;
    config.nthreads = 2;
    config.output_prefix = "az";

    const std::string cmd = RugnuxCommandLine(config, x, "in.h5");
    CHECK(cmd.rfind("rugnux", 0) == 0);
    CHECK(cmd.find("--mode azint") != std::string::npos);
    CHECK(cmd.find("--azim-phi-bins 4") != std::string::npos);
    CHECK(cmd.find("--azim-min-q") != std::string::npos);
    CHECK(cmd.find("in.h5") != std::string::npos);
}

namespace {
    // A field of identical round Gaussian spots on three rings, so that the width estimator sees
    // several resolution bands with the same true width and its 1/d fit has to come back flat.
    void PaintGaussianSpots(ImagePreprocessorBuffer &image, int w, double sigma, double total_counts,
                            std::vector<DiffractionSpot> &spots) {
        constexpr int BKG = 3;
        for (size_t i = 0; i < image.size(); i++) image[i] = BKG;
        const double amp = total_counts / (2.0 * M_PI * sigma * sigma);
        for (int radius : {150, 350, 550})
            for (int k = 0; k < 20; k++) {
                const double phi = 2.0 * M_PI * k / 20.0 + 0.1 * radius;
                const int cx = static_cast<int>(std::lround(600 + radius * std::cos(phi)));
                const int cy = static_cast<int>(std::lround(600 + radius * std::sin(phi)));
                for (int dy = -14; dy <= 14; dy++)
                    for (int dx = -14; dx <= 14; dx++)
                        image[static_cast<size_t>(cy + dy) * w + (cx + dx)] +=
                                static_cast<int32_t>(std::lround(
                                        amp * std::exp(-(dx * dx + dy * dy) / (2.0 * sigma * sigma))));
                spots.emplace_back(static_cast<uint32_t>(cx), static_cast<uint32_t>(cy),
                                   static_cast<int64_t>(total_counts));
            }
    }
}

// The width the adaptive integration radius is set from. A round Gaussian of width sigma holds 80 %
// of its flux inside sqrt(2 ln 5) * sigma = 1.794 * sigma, and that is what the estimator has to
// return - over an aperture that owes nothing to the integrator's r1, which is the whole point of
// measuring it here rather than reading the integrator's own second moment.
TEST_CASE("SpotWidth_Gaussian", "[process]") {
    constexpr int W = 1200, H = 1200;
    DiffractionGeometry geometry;
    geometry.BeamX_pxl(600).BeamY_pxl(600).DetectorDistance_mm(200).PixelSize_mm(0.075)
            .Wavelength_A(1.0);

    for (double sigma : {1.0, 2.2}) {
        ImagePreprocessorBuffer image(static_cast<size_t>(W) * H);
        std::vector<DiffractionSpot> spots;
        PaintGaussianSpots(image, W, sigma, 20000.0, spots);

        std::vector<spot_width::FluxCurve> curves;
        MeasureSpotFluxCurves(image, W, H, geometry, spots, curves);
        REQUIRE(curves.size() >= 45);

        const auto r80 = spot_width::R80AtReference(curves);
        REQUIRE(r80.has_value());
        CHECK(*r80 == Catch::Approx(1.794 * sigma).margin(0.3));
    }

    // The rule the measurement drives: the shipped radius below the line, the capped one above it.
    CHECK(spot_width::R1ForWidth(1.0f) == 4.0f);
    CHECK(spot_width::R1ForWidth(1.794f) == 4.0f);
    CHECK(spot_width::R1ForWidth(2.4f) == 5.0f);
    CHECK(spot_width::R1ForWidth(3.947f) == 6.0f);
    CHECK(spot_width::R1ForWidth(9.0f) == 6.0f);
}

// The beam-centre arbitration between a first pass at the file's centre and one at the measured
// centre. Two arms on the same lattice are two geometries of one hypothesis: they are judged on the
// signal each measured, not on a CC1/2 that reads the same on both.
TEST_CASE("MeasuredCentreWins", "[process]") {
    const auto arm = [](const UnitCell &cell, gemmi::CrystalSystem system, double cc_half, int64_t strong) {
        ProcessResult r;
        r.has_merge_statistics = true;
        r.consensus_cell = cell;
        r.consensus_centering = 'P';
        r.rotation_lattice_type = LatticeMessage{'P', 0, system};
        r.search_merge_cc_half = cc_half;
        r.search_merge_strong_reflections = strong;
        r.search_merge_reflections = 40000;
        return r;
    };
    const UnitCell triclinic{40.0f, 41.0f, 100.0f, 86.0f, 84.0f, 72.0f};
    // The same lattice in the setting with a and b exchanged, as a noisy a ~ b can come out.
    const UnitCell swapped{41.02f, 39.98f, 100.1f, 84.0f, 86.0f, 72.0f};
    const auto tri = gemmi::CrystalSystem::Triclinic;

    const ProcessResult file = arm(triclinic, tri, 0.87, 9000);
    REQUIRE(ArmsHoldSameLattice(file, arm(swapped, tri, 0.87, 9000)));

    // Same lattice, same CC1/2: the arm that measured more signal wins, whichever centre it is...
    CHECK(!MeasuredCentreWins(file, arm(swapped, tri, 0.87, 9500)).empty());
    CHECK(MeasuredCentreWins(file, arm(swapped, tri, 0.87, 8500)).empty());
    // ...a tie keeps the file's centre, and a higher CC1/2 alone does not move it.
    CHECK(MeasuredCentreWins(file, arm(swapped, tri, 0.99, 9000)).empty());

    // Different lattices: the metric question, decided on the search merges.
    const UnitCell monoclinic{57.0f, 42.0f, 100.0f, 90.0f, 95.0f, 90.0f};
    const auto mono = gemmi::CrystalSystem::Monoclinic;
    REQUIRE(!ArmsHoldSameLattice(file, arm(monoclinic, mono, 0.87, 9000)));
    CHECK(!MeasuredCentreWins(file, arm(monoclinic, mono, 0.95, 9000)).empty());
    CHECK(MeasuredCentreWins(file, arm(monoclinic, mono, 0.88, 20000)).empty());
}