diff --git a/common/AzimuthalIntegrationMapping.cpp b/common/AzimuthalIntegrationMapping.cpp index bb74eda4..03480451 100644 --- a/common/AzimuthalIntegrationMapping.cpp +++ b/common/AzimuthalIntegrationMapping.cpp @@ -2,15 +2,22 @@ // SPDX-License-Identifier: GPL-3.0-only #include +#include +#include #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(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 &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 &mask) { + const std::vector &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> 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 &mask) { + const std::vector &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(nthreads, experiment.GetModulesNum()); + std::vector> 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(); -} \ No newline at end of file +} + +size_t AzimuthalIntegrationMapping::GetNThreads() const { + return nthreads; +} diff --git a/common/AzimuthalIntegrationMapping.h b/common/AzimuthalIntegrationMapping.h index 9a84dd4d..996e4b32 100644 --- a/common/AzimuthalIntegrationMapping.h +++ b/common/AzimuthalIntegrationMapping.h @@ -25,22 +25,21 @@ protected: std::optional 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 &mask); - - void SetupConvGeom(const DiffractionGeometry &geom, - const std::vector &mask); + void SetupRawGeom(const DiffractionExperiment& experiment, const std::vector &mask, size_t nthreads); + void SetupConvGeomRows(const DiffractionGeometry &geom, const std::vector &mask, size_t row0, size_t row_end); + void SetupConvGeom(const DiffractionGeometry &geom, const std::vector &mask, size_t nthreads); void SetupPixel(const DiffractionGeometry &geom, const std::vector &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&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; }; diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp index 5e560d18..857d4de3 100644 --- a/receiver/JFJochReceiver.cpp +++ b/receiver/JFJochReceiver.cpp @@ -49,7 +49,7 @@ JFJochReceiver::JFJochReceiver(const DiffractionExperiment &in_experiment, az_int_mapping = std::make_unique(experiment, pixel_mask); auto end_time_point = std::chrono::steady_clock::now(); auto duration = std::chrono::duration(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); diff --git a/tests/AzimuthalIntegrationTest.cpp b/tests/AzimuthalIntegrationTest.cpp index 9dd030db..ed48c14d 100644 --- a/tests/AzimuthalIntegrationTest.cpp +++ b/tests/AzimuthalIntegrationTest.cpp @@ -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 +static void RequireVectorsEqual(const std::vector &ref, + const std::vector &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); +} \ No newline at end of file