mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 92 Wiki Activity
Files
rc184
Jungfraujoch/tests/AzimuthalIntegrationTest.cpp
T
leonarski_fandClaude Opus 5.5 8d1157f711 Azimuthal integration: review follow-ups for the Q binning 
The low-Q headroom check and the detector-edge high Q are floored by the linear
spacing only where one was given, and by the smallest allowed Q otherwise, instead
of by a spacing log10 bins do not use. GetQSpacing_recipA says that under log10 it
means nothing about the bins. GetQBinEdges_recipA, which only the tests used, is
dropped: the tests state the edges from the documented formula. A test pins the
far-edge clamp of linear-by-spacing bins on a two-row setting. The broker docs note
that a GET carries q_bins and no q_spacing once bins are configured.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK
2026-11-10 16:55:48 +01:00

869 lines
36 KiB
C++
Raw Permalink Blame History

// 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));
}
Reference in new issue View Git Blame Copy Permalink