leonarski_fandClaude Opus 5
0162afd324
The mapping took its width and height from the CONVERTED geometry unconditionally, while pixel_to_bin is sized per mode: converted when the geometry is transformed, raw module layout when it is not. In raw mode the two disagreed - 2068x2162 reported against a 1024x4096 map on a JF4M. Only the adaptive spot finders read those dimensions, and they read them for exactly the thing that breaks: the CPU finder derives npix = w*h and then indexes the image, pixel_to_bin and the resolution mask with it, so it walked ~277k pixels past the end of all three; the GPU finder stays in bounds but decodes the strong-pixel bit index with the wrong row stride and reports spots at wrong coordinates. Nothing combines raw geometry with adaptive detection today, so this was latent rather than live. Take them from GetXPixelsNum()/GetYPixelsNum(), which already follow the geometry mode. The converted path is unchanged - it is the same number there - and every internal use is inside SetupConvGeom, which only runs when the geometry is transformed. Covered by two tests: the mapping's dimensions must match pixel_to_bin in both modes, and the CPU adaptive finder must return a spot planted on a raw-geometry image at that raw pixel. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
548 lines
21 KiB
C++
548 lines
21 KiB
C++
// 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"
|
|
|
|
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 * M_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]") {
|
|
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]") {
|
|
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);
|
|
} |