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

#include <catch2/catch_all.hpp>

#include "../common/AzimuthalIntegrationProfile.h"
#include "../common/AzimuthalIntegrationMapping.h"
#include "../common/JFJochMath.h"

TEST_CASE("AzimuthalIntegrationMapping_Constructor","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());

    REQUIRE(x.GetPixelsNum() == 2164*2068);

    std::unique_ptr<AzimuthalIntegrationMapping> radial;

    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 5);

    PixelMask pixel_mask(x);
    REQUIRE_NOTHROW(radial = std::make_unique<AzimuthalIntegrationMapping>(x, pixel_mask));
}

TEST_CASE("AzimuthalIntegrationMapping_GetBinNumber","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 4);

    PixelMask pixel_mask(x);
    AzimuthalIntegrationMapping mapping(x, pixel_mask);
    REQUIRE(mapping.GetBinNumber() == 39);
}

TEST_CASE("AzimuthalIntegrationMapping_GetBinNumber_mask","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 9);

    std::vector<uint32_t> pixel_mask(x.GetPixelsNum(), 0);
    auto geom = x.GetDiffractionGeometry();
    for (int row = 0; row < x.GetYPixelsNum(); row++) {
        for (int col = 0; col < x.GetXPixelsNum(); col++) {
            float d = geom.PxlToRes(col, row);
            float q = 2 * PI / d;
            if (q >= 3.1)
                pixel_mask[row * x.GetXPixelsNum() + col] = 1;
        }
    }

    PixelMask pixel_mask_obj(x);
    pixel_mask_obj.LoadUserMask(x, pixel_mask);
    AzimuthalIntegrationMapping mapping(x, pixel_mask_obj);
    REQUIRE(mapping.GetBinNumber() == 89);
}

TEST_CASE("AzimuthalIntegrationMapping_GetBinNumber_DetectorLimit","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);

    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 9.9);
    PixelMask pixel_mask(x);
    AzimuthalIntegrationMapping mapping(x, pixel_mask);

    REQUIRE(mapping.GetBinNumber() == 98);
}

TEST_CASE("AzimuthalIntegrationMapping_GetBinToQ","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 4);

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

    auto bin_to_q = mapping.GetBinToQ();

    REQUIRE(bin_to_q.size() == 39);
    CHECK(bin_to_q[0] == Catch::Approx(0.15));
    CHECK(bin_to_q[1] == Catch::Approx(0.25));
    CHECK(bin_to_q[15] == Catch::Approx(1.65));
    CHECK(bin_to_q[38] == Catch::Approx(3.95));
}

TEST_CASE("AzimuthalIntegrationMapping_GetBinToPhi","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);

    AzimuthalIntegrationSettings settings;

    settings.QSpacing_recipA(0.1).QRange_recipA(0.1, 4.1).AzimuthalBinCount(4);
    x.ImportAzimuthalIntegrationSettings(settings);

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

    auto &bin_to_q = mapping.GetBinToQ();

    REQUIRE(bin_to_q.size() == 40 * 4);
    CHECK(bin_to_q[0] == Catch::Approx(0.15));
    CHECK(bin_to_q[1] == Catch::Approx(0.25));
    CHECK(bin_to_q[15] == Catch::Approx(1.65));
    CHECK(bin_to_q[38] == Catch::Approx(3.95));

    auto &bin_to_phi = mapping.GetBinToPhi();

    REQUIRE(bin_to_phi.size() == 40 * 4);
    CHECK(bin_to_phi[0] == Catch::Approx(0));
    CHECK(bin_to_phi[1] == Catch::Approx(0));
    CHECK(bin_to_phi[38] == Catch::Approx(0));
    CHECK(bin_to_phi[40] == Catch::Approx(90));
}

TEST_CASE("AzimuthalIntegrationMapping_GetBin","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);

    AzimuthalIntegrationSettings settings;

    settings.QSpacing_recipA(0.1).QRange_recipA(0.1, 4.1).AzimuthalBinCount(4);

    CHECK(settings.GetBin(0.11f, 0.0f) == 0);
    CHECK(settings.GetBin(0.11f, 95.0f) == 40);
    CHECK(settings.GetBin(0.11f, 185.0f) == 80);
    CHECK(settings.GetBin(4.02f, 280.0f) == 120 + 39);
}

#include "../preview/JFJochTIFF.h"

TEST_CASE("AzimuthalIntegrationMapping_GetMapping","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(100).BeamX_pxl(1000).BeamY_pxl(1000);

    AzimuthalIntegrationSettings settings;

    settings.QSpacing_recipA(0.1).QRange_recipA(0.1, 4.1).AzimuthalBinCount(4);
    x.ImportAzimuthalIntegrationSettings(settings);

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

    auto map = mapping.GetPixelToBin();

    CompressedImage image(map, x.GetXPixelsNum(), x.GetYPixelsNum());
    WriteTIFFToFile("test.tiff", image);

    CHECK(map[x.GetXPixelsNum() * 500 + 1500] / 40 == 3);
    CHECK(map[x.GetXPixelsNum() * 500 + 500] / 40 == 2);
    CHECK(map[x.GetXPixelsNum() * 1500 + 500] / 40 == 1);
    CHECK(map[x.GetXPixelsNum() * 1500 + 1500] / 40 == 0);
}

TEST_CASE("AzimuthalIntegrationMapping_QToBin","[AzimuthalIntegration][portable]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 4);

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

    REQUIRE(mapping.QToBin(0.0) == 0);
    REQUIRE(std::floor(mapping.QToBin(0.200001)) == 1);
    REQUIRE(mapping.QToBin(0.6) == Catch::Approx(5));
    REQUIRE(mapping.QToBin(50.0) == Catch::Approx(38));
}

TEST_CASE("AzimuthalIntegrationProfile","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 4);

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

    AzimuthalIntegrationProfile profile(mapping);

    std::vector<float> sum(mapping.GetBinNumber()), sum2(mapping.GetBinNumber());
    std::vector<uint32_t> count(mapping.GetBinNumber());

    for (int i = 0; i < mapping.GetBinNumber(); i++) {
        sum[i] = i * i * 4;
        sum2[i] = i * i * i * i * 4;
        count[i] = i;
    }
    REQUIRE_NOTHROW(profile.Add(sum, sum2, count));
    REQUIRE_NOTHROW(profile.Add(sum, sum2, count));

    std::vector<float> sum_wr(mapping.GetBinNumber() - 1);
    REQUIRE_THROWS(profile.Add(sum_wr, sum2, count));

    auto plot = profile.GetPlot();
    REQUIRE(plot.GetPlots().size() == 1);

    REQUIRE(plot.GetPlots()[0].x.size() == mapping.GetBinNumber());
    REQUIRE(plot.GetPlots()[0].y.size() == mapping.GetBinNumber());
    REQUIRE(plot.GetPlots()[0].x[0] == Catch::Approx(mapping.GetBinToQ()[0]));
    REQUIRE(std::isnan(plot.GetPlots()[0].y[0]));

    for (int i = 1; i < mapping.GetBinNumber(); i++) {
        REQUIRE(plot.GetPlots()[0].x[i] == Catch::Approx(mapping.GetBinToQ()[i]));
        REQUIRE(plot.GetPlots()[0].y[i] == Catch::Approx(i * 4));
    }
}

TEST_CASE("AzimuthalIntegrationProfile_GetStd","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 4);

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

    AzimuthalIntegrationProfile profile(mapping);

    REQUIRE(mapping.GetBinNumber() >= 4);

    std::vector<float> sum(mapping.GetBinNumber()), sum2(mapping.GetBinNumber());
    std::vector<uint32_t> count(mapping.GetBinNumber());

    sum[0] = 2 + 3 + 4 + 5;
    sum2[0] = 4 + 9 + 16 + 25;
    count[0] = 4;

    sum[1] = 1 + 1;
    sum2[1] = 1 + 1;
    count[1] = 2;

    sum[2] = 1;
    sum2[2] = 1;
    count[2] = 1;

    sum[3] = 0;
    sum2[3] = 0;
    count[3] = 0;

    REQUIRE_NOTHROW(profile.Add(sum, sum2, count));

    auto ret_mean = profile.GetResult();
    auto ret_stddev = profile.GetStd();
    auto ret_count = profile.GetPixelCount();
    REQUIRE(ret_mean.size() == mapping.GetBinNumber());
    REQUIRE(ret_stddev.size() == mapping.GetBinNumber());
    REQUIRE(ret_count.size() == mapping.GetBinNumber());

    CHECK(ret_mean[0] == Catch::Approx(3.5));
    CHECK(ret_stddev[0] == Catch::Approx(std::sqrt(5.0/ 3.0)));
    CHECK(ret_count[0] == 4);
    
    CHECK(ret_mean[1] == Catch::Approx(1.0));
    CHECK(ret_stddev[1] == 0.0f);
    CHECK(ret_count[1] == 2);
    
    CHECK(ret_mean[2] == Catch::Approx(1.0));
    CHECK(std::isnan(ret_stddev[2]));
    CHECK(ret_count[2] == 1);
    
    CHECK(std::isnan(ret_mean[3]));
    CHECK(std::isnan(ret_stddev[3]));
    CHECK(ret_count[3] == 0);
}

TEST_CASE("AzimuthalIntegrationMapping_DimensionsMatchPixelToBin","[AzimuthalIntegration]") {
    // The reported dimensions have to describe the image pixel_to_bin was built for, in both
    // geometry modes - the adaptive spot finders walk the image with them and index pixel_to_bin.
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080);
    x.QSpacingForAzimInt_recipA(0.05).QRangeForAzimInt_recipA(0.05, 5.0);

    PixelMask converted_mask(x);
    AzimuthalIntegrationMapping converted(x, converted_mask);
    CHECK(converted.GetWidth() * converted.GetHeight() == converted.GetPixelToBin().size());

    x.GeometryTransformation(false);
    PixelMask raw_mask(x);
    AzimuthalIntegrationMapping raw(x, raw_mask);
    CHECK(raw.GetWidth() * raw.GetHeight() == raw.GetPixelToBin().size());
}

TEST_CASE("AzimuthalIntegrationProfile_GetStd_AfterClear","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 4);

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

    AzimuthalIntegrationProfile profile(mapping);

    std::vector<float> sum(mapping.GetBinNumber()), sum2(mapping.GetBinNumber());
    std::vector<uint32_t> count(mapping.GetBinNumber());

    sum[0] = 2 + 3 + 4 + 5;
    sum2[0] = 4 + 9 + 16 + 25;
    count[0] = 4;

    // The same frame twice, with a Clear() in between: a profile is reused for every image of a
    // dataset, so the second image has to give exactly the first one's standard deviation.
    profile.Add(sum, sum2, count);
    const auto first = profile.GetStd();

    profile.Clear(mapping);
    profile.Add(sum, sum2, count);
    const auto second = profile.GetStd();

    CHECK(first[0] == Catch::Approx(std::sqrt(5.0 / 3.0)));
    CHECK(second[0] == Catch::Approx(first[0]));
}

TEST_CASE("AzimuthalIntegrationProfile_ClearToLargerMapping","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 4);

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

    // A wider q range - the viewer re-uses one profile across datasets, so every vector Clear()
    // touches has to end up the size the new mapping asks for.
    x.QRangeForAzimInt_recipA(0.1, 9);
    AzimuthalIntegrationMapping large(x, pixel_mask);
    REQUIRE(large.GetBinNumber() > small.GetBinNumber());

    AzimuthalIntegrationProfile profile(small);
    profile.Clear(large);

    std::vector<float> sum(large.GetBinNumber()), sum2(large.GetBinNumber());
    std::vector<uint32_t> count(large.GetBinNumber());

    // A bin that exists only in the wider mapping.
    const auto bin = small.GetBinNumber();
    sum[bin] = 2 + 3 + 4 + 5;
    sum2[bin] = 4 + 9 + 16 + 25;
    count[bin] = 4;

    profile.Add(sum, sum2, count);

    const auto stddev = profile.GetStd();
    REQUIRE(stddev.size() == large.GetBinNumber());
    CHECK(stddev[bin] == Catch::Approx(std::sqrt(5.0 / 3.0)));
}

TEST_CASE("AzimuthalIntegrationProfile_operatorAdd","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 4);

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

    AzimuthalIntegrationProfile profile0(mapping), profile1(mapping);

    std::vector<float> sum(mapping.GetBinNumber()), sum2(mapping.GetBinNumber());
    std::vector<uint32_t> count(mapping.GetBinNumber());

    for (int i = 0; i < mapping.GetBinNumber(); i++) {
        sum[i] = (i + 1) * i * 4;
        sum2[i] = (i+ 1) * i * 5;
        count[i] = i + 1;
    }
    REQUIRE_NOTHROW(profile0.Add(sum, sum2, count));
    REQUIRE_NOTHROW(profile1 += profile0);

    auto plot = profile1.GetPlot();

    REQUIRE(plot.GetPlots().size() == 1);
    REQUIRE(plot.GetPlots()[0].x.size() == mapping.GetBinNumber());
    REQUIRE(plot.GetPlots()[0].y.size() == mapping.GetBinNumber());
    for (int i = 0; i < mapping.GetBinNumber(); i++) {
        REQUIRE(plot.GetPlots()[0].x[i] == Catch::Approx(mapping.GetBinToQ()[i]));
        REQUIRE(plot.GetPlots()[0].y[i] == Catch::Approx(i * 4));
    }
}

TEST_CASE("AzimuthalIntegrationProfile_GetMeanValueOfBins","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 4);

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

    AzimuthalIntegrationProfile profile(mapping);

    std::vector<float> sum(mapping.GetBinNumber());
    std::vector<float> sum2(mapping.GetBinNumber());
    std::vector<uint32_t> count(mapping.GetBinNumber());

    for (int i = 0; i < mapping.GetBinNumber(); i++) {
        sum[i] = i * i * 4;
        sum2[i] = i * i * i * i * 4;
        count[i] = i;
    }
    REQUIRE_NOTHROW(profile.Add(sum, sum2, count));

    REQUIRE(profile.GetMeanValueOfBins(0,2) == Catch::Approx((sum[0] + sum[1] + sum[2])  / double(count[0] + count[1] + count[2])));
    REQUIRE(profile.GetMeanValueOfBins(5,7) == Catch::Approx((sum[5] + sum[6] + sum[7])  / double (count[5] + count[6] + count[7])));

    x.BkgEstimateQRange_recipA(0.7, 0.8);

    REQUIRE(profile.GetBkgEstimate(x.GetAzimuthalIntegrationSettings()) == Catch::Approx((sum[5] + sum[6] + sum[7])  / double (count[5] + count[6] + count[7])));

    x.BkgEstimateQRange_recipA(0.01, 0.345);
    REQUIRE(profile.GetBkgEstimate(x.GetAzimuthalIntegrationSettings()) == Catch::Approx((sum[0] + sum[1] + sum[2])  / double(count[0] + count[1] + count[2])));
}

TEST_CASE("AzimuthalIntegrationProfile_GetResult1D","[AzimuthalIntegration][portable]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);

    AzimuthalIntegrationSettings settings;
    settings.QSpacing_recipA(0.1f).QRange_recipA(0.1f, 0.4f).AzimuthalBinCount(3);
    x.ImportAzimuthalIntegrationSettings(settings);

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

    REQUIRE(mapping.GetQBinCount() == 3);
    REQUIRE(mapping.GetAzimuthalBinCount() == 3);
    REQUIRE(mapping.GetBinNumber() == 9);

    std::vector<float> sum(mapping.GetBinNumber(), 0.0f);
    std::vector<float> sum2(mapping.GetBinNumber(), 20.0f);
    std::vector<uint32_t> count(mapping.GetBinNumber(), 0);

    // Layout is [azimuth][q], flattened:
    // az0: q0 q1 q2
    // az1: q0 q1 q2
    // az2: q0 q1 q2
    //
    // Choose values so the correct collapsed result is easy to verify:
    // q0 -> (10 + 20 + 30) / 3 = 20
    // q1 -> (11 + 21 + 31) / 3 = 21
    // q2 -> (12 + 22 + 32) / 3 = 22
    sum[0] = 10; count[0] = 1; // az0 q0
    sum[1] = 11; count[1] = 1; // az0 q1
    sum[2] = 12; count[2] = 1; // az0 q2

    sum[3] = 20; count[3] = 1; // az1 q0
    sum[4] = 21; count[4] = 1; // az1 q1
    sum[5] = 22; count[5] = 1; // az1 q2

    sum[6] = 30; count[6] = 1; // az2 q0
    sum[7] = 31; count[7] = 1; // az2 q1
    sum[8] = 32; count[8] = 1; // az2 q2

    REQUIRE_NOTHROW(profile.Add(sum, sum2, count));

    auto result_1d = profile.GetResult1D();

    REQUIRE(result_1d.size() == 3);
    CHECK(result_1d[0] == Catch::Approx(20.0f));
    CHECK(result_1d[1] == Catch::Approx(21.0f));
    CHECK(result_1d[2] == Catch::Approx(22.0f));
}

template <class T>
static void RequireVectorsEqual(const std::vector<T> &ref,
                                const std::vector<T> &other,
                                const std::string &name,
                                int nthreads) {
    INFO(name << ", threads=" << nthreads);
    REQUIRE(ref.size() == other.size());
    CHECK(memcmp(ref.data(), other.data(), sizeof(T) * ref.size()) == 0);
}

static void CheckAzimuthalIntegrationMappingThreadingExact(const DiffractionExperiment &experiment) {
    PixelMask pixel_mask(experiment);

    AzimuthalIntegrationMapping mapping_1(experiment, pixel_mask, 1);
    AzimuthalIntegrationMapping mapping_2(experiment, pixel_mask, 2);
    AzimuthalIntegrationMapping mapping_16(experiment, pixel_mask, 16);
    AzimuthalIntegrationMapping mapping_0(experiment, pixel_mask, 0);

    REQUIRE(mapping_1.GetBinNumber() == mapping_2.GetBinNumber());
    REQUIRE(mapping_1.GetBinNumber() == mapping_16.GetBinNumber());
    REQUIRE(mapping_1.GetBinNumber() == mapping_0.GetBinNumber());

    REQUIRE(mapping_1.GetQBinCount() == mapping_2.GetQBinCount());
    REQUIRE(mapping_1.GetQBinCount() == mapping_16.GetQBinCount());
    REQUIRE(mapping_1.GetQBinCount() == mapping_0.GetQBinCount());

    REQUIRE(mapping_1.GetAzimuthalBinCount() == mapping_2.GetAzimuthalBinCount());
    REQUIRE(mapping_1.GetAzimuthalBinCount() == mapping_16.GetAzimuthalBinCount());
    REQUIRE(mapping_1.GetAzimuthalBinCount() == mapping_0.GetAzimuthalBinCount());

    RequireVectorsEqual(mapping_1.GetPixelToBin(), mapping_2.GetPixelToBin(), "pixel_to_bin", 2);
    RequireVectorsEqual(mapping_1.GetPixelToBin(), mapping_16.GetPixelToBin(), "pixel_to_bin", 16);
    RequireVectorsEqual(mapping_1.GetPixelToBin(), mapping_0.GetPixelToBin(), "pixel_to_bin", 0);

    RequireVectorsEqual(mapping_1.Resolution(), mapping_2.Resolution(), "resolution", 2);
    RequireVectorsEqual(mapping_1.Resolution(), mapping_16.Resolution(), "resolution", 16);
    RequireVectorsEqual(mapping_1.Resolution(), mapping_0.Resolution(), "resolution", 0);

    RequireVectorsEqual(mapping_1.Corrections(), mapping_2.Corrections(), "corrections", 2);
    RequireVectorsEqual(mapping_1.Corrections(), mapping_16.Corrections(), "corrections", 16);
    RequireVectorsEqual(mapping_1.Corrections(), mapping_0.Corrections(), "corrections", 0);

    RequireVectorsEqual(mapping_1.GetBinToQ(), mapping_2.GetBinToQ(), "bin_to_q", 2);
    RequireVectorsEqual(mapping_1.GetBinToQ(), mapping_16.GetBinToQ(), "bin_to_q", 16);
    RequireVectorsEqual(mapping_1.GetBinToQ(), mapping_0.GetBinToQ(), "bin_to_q", 0);

    RequireVectorsEqual(mapping_1.GetBinToD(), mapping_2.GetBinToD(), "bin_to_d", 2);
    RequireVectorsEqual(mapping_1.GetBinToD(), mapping_16.GetBinToD(), "bin_to_d", 16);
    RequireVectorsEqual(mapping_1.GetBinToD(), mapping_0.GetBinToD(), "bin_to_d", 0);

    RequireVectorsEqual(mapping_1.GetBinToTwoTheta(), mapping_2.GetBinToTwoTheta(), "bin_to_2theta", 2);
    RequireVectorsEqual(mapping_1.GetBinToTwoTheta(), mapping_16.GetBinToTwoTheta(), "bin_to_2theta", 16);
    RequireVectorsEqual(mapping_1.GetBinToTwoTheta(), mapping_0.GetBinToTwoTheta(), "bin_to_2theta", 0);

    RequireVectorsEqual(mapping_1.GetBinToPhi(), mapping_2.GetBinToPhi(), "bin_to_phi", 2);
    RequireVectorsEqual(mapping_1.GetBinToPhi(), mapping_16.GetBinToPhi(), "bin_to_phi", 16);
    RequireVectorsEqual(mapping_1.GetBinToPhi(), mapping_0.GetBinToPhi(), "bin_to_phi", 0);
}

TEST_CASE("AzimuthalIntegrationMapping_Threading_FixedGeometry_2000x2000", "[AzimuthalIntegration]") {
    DiffractionExperiment x(DetDECTRIS(2000, 2000, "E16M", ""));
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 10);
    x.PolarizationFactor(0.99f);

    REQUIRE(x.IsGeometryTransformed());

    CheckAzimuthalIntegrationMappingThreadingExact(x);
}

TEST_CASE("AzimuthalIntegrationMapping_Threading_RawGeometry_18Modules", "[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF9M());
    x.Raw();
    x.DetectorDistance_mm(100).BeamX_pxl(1500).BeamY_pxl(1500);
    x.QSpacingForAzimInt_recipA(0.05).QRangeForAzimInt_recipA(0.1, 8.0);
    x.PolarizationFactor(0.99f);

    REQUIRE(!x.IsGeometryTransformed());
    REQUIRE(x.GetModulesNum() == 18);

    CheckAzimuthalIntegrationMappingThreadingExact(x);
}

TEST_CASE("AzimuthalIntegrationMapping_Threading_ConvertedGeometry_18Modules", "[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF9M());
    x.Conversion();
    x.DetectorDistance_mm(100).BeamX_pxl(1500).BeamY_pxl(1500);
    x.QSpacingForAzimInt_recipA(0.05).QRangeForAzimInt_recipA(0.1, 8.0);
    x.PolarizationFactor(0.99f);

    REQUIRE(x.IsGeometryTransformed());
    REQUIRE(x.GetModulesNum() == 18);

    CheckAzimuthalIntegrationMappingThreadingExact(x);
}

// A mapping that takes its geometry from the cache - built there under another mask, or at another
// centre before - is the mapping built from scratch, bit for bit.
TEST_CASE("AzimuthalIntegrationMapping_GeometryCacheExact", "[AzimuthalIntegration]") {
    DiffractionExperiment x(DetDECTRIS(2000, 2000, "E16M", ""));
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 10);
    x.PolarizationFactor(0.99f);

    PixelMask plain_mask(x);
    PixelMask user_mask(x);
    std::vector<uint32_t> blocked(x.GetPixelsNum(), 0);
    for (int row = 900; row < 1100; row++)
        for (int col = 0; col < 1100; col++)
            blocked[row * x.GetXPixelsNum() + col] = 1;
    user_mask.LoadUserMask(x, blocked);

    AzimuthalIntegrationGeometryCache cache;
    const auto check = [&](const DiffractionExperiment &experiment, const PixelMask &mask) {
        const AzimuthalIntegrationMapping ref(experiment, mask);
        const AzimuthalIntegrationMapping cached(experiment, mask, cache);
        RequireVectorsEqual(ref.GetPixelToBin(), cached.GetPixelToBin(), "pixel_to_bin", 0);
        RequireVectorsEqual(ref.Resolution(), cached.Resolution(), "resolution", 0);
        RequireVectorsEqual(ref.Corrections(), cached.Corrections(), "corrections", 0);
        CHECK(ref.GetPixelToBinChecksum() == cached.GetPixelToBinChecksum());
        CHECK(ref.GetCorrectionsChecksum() == cached.GetCorrectionsChecksum());
    };
    check(x, plain_mask);   // computed
    check(x, user_mask);    // taken from the cache, under another mask
    check(x, plain_mask);

    DiffractionExperiment moved(x);
    moved.BeamX_pxl(1003.5f);
    check(moved, user_mask);   // another centre: computed again
    check(x, user_mask);
}

// The ice score on an explicit profile: flat means no ice (1.0), and a bump planted on a hexagonal
// ring is reported at its own height over the background. Uses the static entry point, which is what
// feeds the score the peak-excluded per-ring background instead of the plain profile.
TEST_CASE("AzimuthalIntegrationProfile_IceRingScore","[AzimuthalIntegration]") {
    AzimuthalIntegrationSettings settings;
    settings.QSpacing_recipA(0.01f).QRange_recipA(0.1f, 4.5f);

    const int q_bins = settings.GetQBinCount();
    REQUIRE(q_bins > 400);

    auto bin_of = [&](float d_A) {
        const float q = 6.283185307f / d_A;
        return static_cast<int>(std::lround((q - settings.GetLowQ_recipA()) / settings.GetQSpacing_recipA() - 0.5f));
    };

    std::vector<float> profile(q_bins, 100.0f);
    CHECK(AzimuthalIntegrationProfile::IceRingScore(profile, q_bins, settings, 0.03f) == Catch::Approx(1.0f));

    // A bump at the 2.249 A hexagonal ring, on a background the running median still reads as 100.
    profile[bin_of(2.249f)] = 250.0f;
    CHECK(AzimuthalIntegrationProfile::IceRingScore(profile, q_bins, settings, 0.03f) == Catch::Approx(2.5f));

    // A bump of the same size well away from every ring is not ice and must not be reported.
    std::vector<float> off_ring(q_bins, 100.0f);
    off_ring[bin_of(2.500f)] = 250.0f;
    CHECK(AzimuthalIntegrationProfile::IceRingScore(off_ring, q_bins, settings, 0.03f) == Catch::Approx(1.0f));
}

// A profile given per (q, azimuth) bin is averaged over azimuth first, so a ring seen in every sector
// scores the same as the equivalent 1-D profile.
TEST_CASE("AzimuthalIntegrationProfile_IceRingScore_Azimuthal","[AzimuthalIntegration]") {
    AzimuthalIntegrationSettings settings;
    settings.QSpacing_recipA(0.01f).QRange_recipA(0.1f, 4.5f).AzimuthalBinCount(4);

    const int q_bins = settings.GetQBinCount();
    const float q = 6.283185307f / 2.249f;
    const int ring = static_cast<int>(std::lround((q - settings.GetLowQ_recipA()) / settings.GetQSpacing_recipA() - 0.5f));

    std::vector<float> flat(q_bins, 100.0f);
    std::vector<float> sectors(static_cast<size_t>(q_bins) * 4, 100.0f);
    flat[ring] = 250.0f;
    for (int az = 0; az < 4; az++)
        sectors[static_cast<size_t>(az) * q_bins + ring] = 250.0f;

    CHECK(AzimuthalIntegrationProfile::IceRingScore(sectors, q_bins, settings, 0.03f)
          == Catch::Approx(AzimuthalIntegrationProfile::IceRingScore(flat, q_bins, settings, 0.03f)));
}

TEST_CASE("AzimuthalIntegrationSettings_QBinning_Validation","[AzimuthalIntegration]") {
    AzimuthalIntegrationSettings s;
    s.QRange_recipA(0.01f, 1.0f);

    // Linear takes either the spacing or the bin count, not both.
    REQUIRE_THROWS(s.QBinning(QBinSpacing::Linear, 0.01f, 100));
    // Log10 takes only the bin count, and needs it.
    REQUIRE_THROWS(s.QBinning(QBinSpacing::Log10, 0.01f, 100));
    REQUIRE_THROWS(s.QBinning(QBinSpacing::Log10, 0.01f, std::nullopt));
    REQUIRE_THROWS(s.QBinning(QBinSpacing::Log10, std::nullopt, std::nullopt));
    REQUIRE_THROWS(s.QBinning(QBinSpacing::Log10, std::nullopt, 0));
    REQUIRE_THROWS(s.QBinning(QBinSpacing::Log10, std::nullopt, 65535));
    REQUIRE_THROWS(s.QBinning(QBinSpacing::Linear, 0.0f, std::nullopt));
    // log10(q) needs q > 0: a low Q of zero or below is refused (minimum 1e-5 1/A, for every spacing).
    {
        AzimuthalIntegrationSettings log_s;
        log_s.QBinning(QBinSpacing::Log10, std::nullopt, 100);
        REQUIRE_THROWS(log_s.QRange_recipA(0.0f, 1.0f));
        REQUIRE_THROWS(log_s.QRange_recipA(-0.1f, 1.0f));
        REQUIRE_NOTHROW(log_s.QRange_recipA(1e-5f, 1.0f));
        CHECK(std::isfinite(log_s.GetQBinCentres_recipA().front()));
    }

    // A refused call leaves the settings as they were.
    CHECK(s.GetQBinSpacing() == QBinSpacing::Linear);
    CHECK(s.GetRequestedQSpacing_recipA() == 0.01f);
    CHECK(!s.GetRequestedQBinCount());

    REQUIRE_NOTHROW(s.QBinning(QBinSpacing::Log10, std::nullopt, 200));
    CHECK(s.GetQBinSpacing() == QBinSpacing::Log10);
    CHECK(s.GetQBinCount() == 200);
    CHECK(s.GetRequestedQBinCount() == 200);
    CHECK(!s.GetRequestedQSpacing_recipA());

    REQUIRE_NOTHROW(s.QBinning(QBinSpacing::Linear, std::nullopt, 99));
    CHECK(s.GetQBinCount() == 99);
    CHECK(s.GetQSpacing_recipA() == Catch::Approx(0.01));
    CHECK(!s.GetRequestedQSpacing_recipA());

    // Linear with neither keeps the spacing it had.
    REQUIRE_NOTHROW(s.QSpacing_recipA(0.02f));
    REQUIRE_NOTHROW(s.QBinning(QBinSpacing::Linear, std::nullopt, std::nullopt));
    CHECK(s.GetRequestedQSpacing_recipA() == 0.02f);
    CHECK(s.GetQBinCount() == 50);
}

// The linear spacing bins as it always has: the same count, the same bin for every q and the same
// centres, compared bit for bit with the formulas it used before log10 binning was added. The one
// exception, which this uniform sample of q cannot reach, is a q that rounds onto the far edge of the
// last bin: it now stays in that bin (AzimuthalIntegrationSettings_Linear_FarEdge).
TEST_CASE("AzimuthalIntegrationSettings_Linear_Unchanged","[AzimuthalIntegration]") {
    const float low = 0.1f, high = 4.5f, spacing = 0.01f;
    AzimuthalIntegrationSettings s;
    s.QSpacing_recipA(spacing).QRange_recipA(low, high).AzimuthalBinCount(3);

    const int32_t q_bins = std::ceil((high - low) / spacing);
    REQUIRE(s.GetQBinCount() == q_bins);
    REQUIRE(s.GetBinCount() == 3 * q_bins);

    const auto centres = s.GetQBinCentres_recipA();
    REQUIRE(centres.size() == q_bins);
    for (int i = 0; i < q_bins; i++)
        CHECK(static_cast<float>(centres[i]) == static_cast<float>(spacing * (i + 0.5) + low));

    for (int k = 0; k < 100000; k++) {
        const float q = low + (high - low) * static_cast<float>(k) / 100000.0f;
        const float phi = 359.9f * static_cast<float>(k % 97) / 97.0f;
        const int16_t q_bin = std::floor((q - low) / spacing);
        const int16_t phi_bin = std::floor(phi / 360.0f * 3);
        REQUIRE(s.GetBin(q, phi) == static_cast<uint16_t>(q_bin + phi_bin * q_bins));
        REQUIRE(s.QToBin(q) == std::min<uint16_t>(s.GetBinCount() - 1, std::floor(std::max(0.0f, (q - low) / spacing))));
    }
}

// A q just below high Q whose (q - low) / spacing rounds up to the bin count used to spill into the
// first bin of the next azimuthal row.
TEST_CASE("AzimuthalIntegrationSettings_Linear_FarEdge","[AzimuthalIntegration]") {
    const float low = 0.1f, high = 3.7f, spacing = 0.1f;
    AzimuthalIntegrationSettings s;
    s.QSpacing_recipA(spacing).QRange_recipA(low, high).AzimuthalBinCount(2);
    const float q = std::nextafter(high, 0.0f);
    REQUIRE(static_cast<int32_t>(std::floor((q - low) / spacing)) == s.GetQBinCount());
    CHECK(s.GetBin(q, 350.0f) == 2 * s.GetQBinCount() - 1);
}

TEST_CASE("AzimuthalIntegrationSettings_Linear_BinCount","[AzimuthalIntegration]") {
    AzimuthalIntegrationSettings s;
    s.QRange_recipA(0.1f, 4.1f).QBinning(QBinSpacing::Linear, std::nullopt, 40);
    CHECK(s.GetQSpacing_recipA() == Catch::Approx(0.1));

    const auto centres = s.GetQBinCentres_recipA();
    CHECK(centres[0] == Catch::Approx(0.15));
    CHECK(centres[39] == Catch::Approx(4.05));

    CHECK(s.GetBin(0.1f, 0.0f) == 0);
    CHECK(s.GetBin(std::nextafter(4.1f, 0.0f), 0.0f) == 39);
    CHECK(s.GetBin(4.1f, 0.0f) == UINT16_MAX);
    CHECK(s.GetBin(0.0999f, 0.0f) == UINT16_MAX);
}

TEST_CASE("AzimuthalIntegrationSettings_Log10","[AzimuthalIntegration]") {
    const float low = 0.001f, high = 1.0f;
    const int32_t n = 30;
    AzimuthalIntegrationSettings s;
    s.QRange_recipA(low, high).QBinning(QBinSpacing::Log10, std::nullopt, n).AzimuthalBinCount(2);
    REQUIRE(s.GetQBinCount() == n);
    REQUIRE(s.GetBinCount() == 2 * n);

    // The edges as the docs state them: q_i = low (high / low)^(i / n).
    std::vector<double> edges(n + 1);
    for (int i = 0; i <= n; i++)
        edges[i] = low * std::pow(static_cast<double>(high) / low, static_cast<double>(i) / n);
    const auto centres = s.GetQBinCentres_recipA();
    REQUIRE(centres.size() == n);
    // Equal steps in log10(q): three decades in 30 bins, every edge 10^0.1 times the one before.
    for (int i = 0; i < n; i++) {
        CHECK(edges[i + 1] > edges[i]);
        CHECK(edges[i + 1] / edges[i] == Catch::Approx(std::pow(10.0, 0.1)));
        // The centre is the mid-point on the log scale: the geometric mean of the edges.
        CHECK(centres[i] == Catch::Approx(std::sqrt(edges[i] * edges[i + 1])));
        // So the width of a bin grows with q, in proportion to its left edge.
        CHECK((edges[i + 1] - edges[i]) / edges[i] == Catch::Approx(std::pow(10.0, 0.1) - 1.0));
        if (i > 0)
            CHECK(edges[i + 1] - edges[i] > edges[i] - edges[i - 1]);
    }

    // Every q in range lands in the one bin whose edges hold it.
    for (int k = 0; k < 100000; k++) {
        const float q = low * std::pow(high / low, static_cast<float>(k) / 100000.0f);
        if (q >= high)
            continue;
        const uint16_t bin = s.GetBin(q, 0.0f);
        REQUIRE(bin < n);
        REQUIRE(q >= edges[bin] * (1 - 1e-5));
        REQUIRE(q < edges[bin + 1] * (1 + 1e-5));
        REQUIRE(s.QToBin(q) == bin);
        // The second azimuthal sector is the same q bins one row further.
        REQUIRE(s.GetBin(q, 270.0f) == bin + n);
    }
    CHECK(s.GetBin(std::nextafter(high, 0.0f), 0.0f) == n - 1);
    CHECK(s.GetBin(high, 0.0f) == UINT16_MAX);
    CHECK(s.GetBin(0.99f * low, 0.0f) == UINT16_MAX);
    CHECK(s.QToBin(0.0f) == 0);
}

// Every pixel of the detector whose q is within the range is assigned to exactly the bin whose edges
// hold its q, and the mapping's Q axis is the bins' true centres.
TEST_CASE("AzimuthalIntegrationMapping_Log10","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(200).BeamX_pxl(1000).BeamY_pxl(1000);
    AzimuthalIntegrationSettings settings;
    settings.QRange_recipA(0.01f, 2.0f).QBinning(QBinSpacing::Log10, std::nullopt, 100).AzimuthalBinCount(4);
    x.ImportAzimuthalIntegrationSettings(settings);

    PixelMask pixel_mask(x);
    AzimuthalIntegrationMapping mapping(x, pixel_mask);
    const auto &s = mapping.Settings();
    std::vector<double> edges(101);
    for (int i = 0; i <= 100; i++)
        edges[i] = 0.01 * std::pow(2.0 / 0.01, i / 100.0);
    const auto centres = s.GetQBinCentres_recipA();
    REQUIRE(mapping.GetQBinCount() == 100);

    const auto &bin_to_q = mapping.GetBinToQ();
    REQUIRE(bin_to_q.size() == 400);
    for (int az = 0; az < 4; az++)
        for (int i = 0; i < 100; i++) {
            CHECK(bin_to_q[az * 100 + i] == static_cast<float>(centres[i]));
            CHECK(mapping.GetBinToD()[az * 100 + i] == Catch::Approx(2.0 * PI / centres[i]));
        }

    const auto &pixel_to_bin = mapping.GetPixelToBin();
    const auto &resolution = mapping.Resolution();
    std::vector<int64_t> per_bin(mapping.GetBinNumber(), 0);
    int64_t in_range = 0;
    for (size_t i = 0; i < pixel_to_bin.size(); i++) {
        if (pixel_mask.GetMask()[i] != 0 || !(resolution[i] > 0)) {
            REQUIRE(pixel_to_bin[i] == UINT16_MAX);
            continue;
        }
        const float q = 2.0f * static_cast<float>(PI) / resolution[i];
        if (q < 0.01f || q >= 2.0f) {
            REQUIRE(pixel_to_bin[i] == UINT16_MAX);
            continue;
        }
        in_range++;
        REQUIRE(pixel_to_bin[i] < mapping.GetBinNumber());
        per_bin[pixel_to_bin[i]]++;
        const int q_bin = pixel_to_bin[i] % 100;
        REQUIRE(q >= edges[q_bin] * (1 - 1e-5));
        REQUIRE(q < edges[q_bin + 1] * (1 + 1e-5));
    }
    REQUIRE(in_range > 0);
    int64_t sum = 0;
    for (auto c: per_bin)
        sum += c;
    CHECK(sum == in_range);
}

// The linear mapping assigns every pixel to the bin the old formula gave it.
TEST_CASE("AzimuthalIntegrationMapping_Linear_Unchanged","[AzimuthalIntegration]") {
    DiffractionExperiment x(DetJF4M());
    x.DetectorDistance_mm(75).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForAzimInt_recipA(0.01).QRangeForAzimInt_recipA(0.05, 4.0);

    PixelMask pixel_mask(x);
    AzimuthalIntegrationMapping mapping(x, pixel_mask);
    const auto &pixel_to_bin = mapping.GetPixelToBin();
    const auto &resolution = mapping.Resolution();
    for (size_t i = 0; i < pixel_to_bin.size(); i++) {
        if (!(resolution[i] > 0))
            continue;
        const float q = 2.0f * static_cast<float>(PI) / resolution[i];
        uint16_t expected = UINT16_MAX;
        if (q >= 0.05f && q < 4.0f)
            expected = static_cast<int16_t>(std::floor((q - 0.05f) / 0.01f));
        REQUIRE(pixel_to_bin[i] == expected);
    }
    const auto &bin_to_q = mapping.GetBinToQ();
    for (size_t i = 0; i < bin_to_q.size(); i++)
        REQUIRE(bin_to_q[i] == static_cast<float>(0.01f * (i + 0.5) + 0.05f));
}

TEST_CASE("AzimuthalIntegrationProfile_IceRingScore_Log10","[AzimuthalIntegration]") {
    AzimuthalIntegrationSettings settings;
    settings.QRange_recipA(0.1f, 4.5f).QBinning(QBinSpacing::Log10, std::nullopt, 800);
    const int q_bins = settings.GetQBinCount();

    std::vector<float> profile(q_bins, 100.0f);
    CHECK(AzimuthalIntegrationProfile::IceRingScore(profile, q_bins, settings, 0.03f) == Catch::Approx(1.0f));

    profile[settings.GetBin(6.283185307f / 2.249f, 0.0f)] = 250.0f;
    CHECK(AzimuthalIntegrationProfile::IceRingScore(profile, q_bins, settings, 0.03f) == Catch::Approx(2.5f));
}