mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests 1 Actions Packages Projects Releases 64 Wiki Activity

AzimuthalIntegrationMapping: Parallel calculation of azimuthal integration mapping

Browse Source
This commit is contained in:
leonarski_f
2026-04-30 11:15:39 +02:00
parent 37495c0902
commit c7f2cad73f
4 changed files with 183 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files
common/AzimuthalIntegrationMapping.cpp
+78 -24
View File
@@ -2,15 +2,22 @@
// SPDX-License-Identifier: GPL-3.0-only
#include <cmath>
#include <thread>
#include <future>
#include "AzimuthalIntegrationMapping.h"
#include "JFJochException.h"
#include "DiffractionGeometry.h"
#include "RawToConvertedGeometry.h"
AzimuthalIntegrationMapping::AzimuthalIntegrationMapping(const AzimuthalIntegrationSettings& in_settings, size_t width, size_t height,
float wavelength)
: settings(in_settings), width(width), height(height), wavelength(wavelength) {
AzimuthalIntegrationMapping::AzimuthalIntegrationMapping(const DiffractionExperiment &experiment,
const PixelMask& mask,
size_t in_nthreads)
: settings(experiment.GetAzimuthalIntegrationSettings()),
wavelength(experiment.GetWavelength_A()),
width(experiment.GetXPixelsNumConv()),
height(experiment.GetYPixelsNumConv()) {
if (width <= 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Detector width must be above 0");
@@ -22,40 +29,60 @@ AzimuthalIntegrationMapping::AzimuthalIntegrationMapping(const AzimuthalIntegrat
if (settings.GetBinCount() >= UINT16_MAX)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Cannot handle more than 65534 az. int. bins");
}
AzimuthalIntegrationMapping::AzimuthalIntegrationMapping(const DiffractionExperiment &experiment,
const PixelMask& mask)
: AzimuthalIntegrationMapping(experiment.GetAzimuthalIntegrationSettings(),
experiment.GetXPixelsNum(),
experiment.GetYPixelsNum(),
experiment.GetWavelength_A()) {
polarization_factor = experiment.GetPolarizationFactor();
if (!experiment.IsGeometryTransformed())
SetupRawGeom(experiment, mask.GetMaskRaw(experiment));
if (in_nthreads == 0)
nthreads = std::thread::hardware_concurrency();
else
SetupConvGeom(experiment.GetDiffractionGeometry(),mask.GetMask());
nthreads = in_nthreads;
nthreads = std::clamp<size_t>(nthreads, 1, 64);
if (!experiment.IsGeometryTransformed())
SetupRawGeom(experiment, mask.GetMaskRaw(experiment), nthreads);
else
SetupConvGeom(experiment.GetDiffractionGeometry(),mask.GetMask(), nthreads);
UpdateMaxBinNumber();
}
void AzimuthalIntegrationMapping::SetupConvGeomRows(const DiffractionGeometry &geom, const std::vector<uint32_t> &mask,
size_t row0, size_t row_end) {
for (size_t row = row0; row < row_end && row < height; row++) {
for (size_t col = 0; col < width; col++)
SetupPixel(geom, mask, row * width + col, col, row);
}
}
void AzimuthalIntegrationMapping::SetupConvGeom(const DiffractionGeometry &geom,
const std::vector<uint32_t> &mask) {
const std::vector<uint32_t> &mask,
size_t nthreads) {
pixel_to_bin.resize(width * height, UINT16_MAX);
pixel_resolution.resize(width * height, 0);
corrections.resize(width * height, 0);
if (mask.size() != width * height)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Mask size invalid");
for (int row = 0; row < height; row++)
for (int col = 0; col < width; col++)
SetupPixel(geom, mask, row * width + col, col, row);
if (nthreads <= 1) {
SetupConvGeomRows(geom, mask, 0, height);
} else {
nthreads = std::min(nthreads, height);
std::vector<std::future<void>> futures;
for (size_t t = 0; t < nthreads; ++t)
futures.emplace_back(std::async(std::launch::async,
&AzimuthalIntegrationMapping::SetupConvGeomRows,
this, std::cref(geom), std::cref(mask), t * height / nthreads, (t + 1) * height / nthreads));
for (auto &f: futures)
f.get();
}
}
void AzimuthalIntegrationMapping::SetupRawGeom(const DiffractionExperiment &experiment,
const std::vector<uint32_t> &mask) {
const std::vector<uint32_t> &mask, size_t nthreads) {
if (mask.size() != RAW_MODULE_SIZE * experiment.GetModulesNum())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Mask size invalid");
@@ -65,11 +92,34 @@ void AzimuthalIntegrationMapping::SetupRawGeom(const DiffractionExperiment &expe
auto geom = experiment.GetDiffractionGeometry();
for (int m = 0; m < experiment.GetModulesNum(); m++) {
for (int pxl = 0; pxl < RAW_MODULE_SIZE; pxl++) {
auto [x,y] = RawToConvertedCoordinate(experiment, m, pxl);
SetupPixel(geom, mask, m * RAW_MODULE_SIZE + pxl, x, y);
if (nthreads <= 1) {
for (int m = 0; m < experiment.GetModulesNum(); m++) {
for (int pxl = 0; pxl < RAW_MODULE_SIZE; pxl++) {
auto [x,y] = RawToConvertedCoordinate(experiment, m, pxl);
SetupPixel(geom, mask, m * RAW_MODULE_SIZE + pxl, x, y);
}
}
} else {
nthreads = std::min<size_t>(nthreads, experiment.GetModulesNum());
std::vector<std::future<void>> futures;
futures.reserve(nthreads);
for (size_t t = 0; t < nthreads; ++t) {
const size_t module_begin = t * experiment.GetModulesNum() / nthreads;
const size_t module_end = (t + 1) * experiment.GetModulesNum() / nthreads;
futures.emplace_back(std::async(std::launch::async, [&, module_begin, module_end] {
for (size_t m = module_begin; m < module_end; ++m) {
for (int pxl = 0; pxl < RAW_MODULE_SIZE; ++pxl) {
auto [x, y] = RawToConvertedCoordinate(experiment, m, pxl);
SetupPixel(geom, mask, m * RAW_MODULE_SIZE + pxl, x, y);
}
}
}));
}
for (auto &f: futures)
f.get();
}
}
@@ -170,4 +220,8 @@ int32_t AzimuthalIntegrationMapping::GetAzimuthalBinCount() const {
int32_t AzimuthalIntegrationMapping::GetQBinCount() const {
return settings.GetQBinCount();
}
}
size_t AzimuthalIntegrationMapping::GetNThreads() const {
return nthreads;
}
common/AzimuthalIntegrationMapping.h
+8 -8
View File
@@ -25,22 +25,21 @@ protected:
std::optional<float> polarization_factor;
size_t nthreads;
void UpdateMaxBinNumber();
AzimuthalIntegrationMapping(const AzimuthalIntegrationSettings& settings,
size_t width, size_t height, float wavelength);
void SetupRawGeom(const DiffractionExperiment& experiment,
const std::vector<uint32_t> &mask);
void SetupConvGeom(const DiffractionGeometry &geom,
const std::vector<uint32_t> &mask);
void SetupRawGeom(const DiffractionExperiment& experiment, const std::vector<uint32_t> &mask, size_t nthreads);
void SetupConvGeomRows(const DiffractionGeometry &geom, const std::vector<uint32_t> &mask, size_t row0, size_t row_end);
void SetupConvGeom(const DiffractionGeometry &geom, const std::vector<uint32_t> &mask, size_t nthreads);
void SetupPixel(const DiffractionGeometry &geom, const std::vector<uint32_t> &mask,
uint32_t pxl, uint32_t col, uint32_t row);
public:
AzimuthalIntegrationMapping(const DiffractionExperiment& experiment,
const PixelMask& mask);
const PixelMask& mask,
size_t nthreads = 0);
[[nodiscard]] uint16_t GetBinNumber() const;
[[nodiscard]] const std::vector<uint16_t>&GetPixelToBin() const;
@@ -56,6 +55,7 @@ public:
[[nodiscard]] const AzimuthalIntegrationSettings& Settings() const;
[[nodiscard]] int32_t GetAzimuthalBinCount() const;
[[nodiscard]] int32_t GetQBinCount() const;
[[nodiscard]] size_t GetNThreads() const;
};
receiver/JFJochReceiver.cpp
+1 -1
View File
@@ -49,7 +49,7 @@ JFJochReceiver::JFJochReceiver(const DiffractionExperiment &in_experiment,
az_int_mapping = std::make_unique<AzimuthalIntegrationMapping>(experiment, pixel_mask);
auto end_time_point = std::chrono::steady_clock::now();
auto duration = std::chrono::duration<float>(end_time_point - start_time_point);
logger.Info("Azimuthal integration mapping done in {:.5f} s", duration.count());
logger.Info("Azimuthal integration mapping done in {:.5f} s with {} threads", duration.count(), az_int_mapping->GetNThreads());
plots.Setup(experiment, *az_int_mapping);
tests/AzimuthalIntegrationTest.cpp
+96
View File
@@ -318,3 +318,99 @@ TEST_CASE("AzimuthalIntegrationProfile_GetResult1D","[AzimuthalIntegration]") {
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);
}