From 2eac43b92588c8efd01851623227d4582d05a58b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 8 Aug 2023 10:20:26 +0200 Subject: [PATCH] RadialIntegration: Enable pixel splitting with CPU based routine --- image_analysis/RadialIntegration.cpp | 67 +++++++++++++++------ image_analysis/RadialIntegration.h | 22 ++++--- image_analysis/RadialIntegrationMapping.cpp | 11 +++- image_analysis/RadialIntegrationMapping.h | 1 + tests/RadialIntegrationTest.cpp | 28 +++++++++ tools/DataAnalysisPerfTest.cpp | 38 ++++++++---- 6 files changed, 126 insertions(+), 41 deletions(-) diff --git a/image_analysis/RadialIntegration.cpp b/image_analysis/RadialIntegration.cpp index 7118af15..73895449 100644 --- a/image_analysis/RadialIntegration.cpp +++ b/image_analysis/RadialIntegration.cpp @@ -4,12 +4,24 @@ #include "RadialIntegration.h" #include "../common/JFJochException.h" -RadialIntegration::RadialIntegration(const std::vector& in_mapping, uint16_t in_nbins) : - pixel_to_bin(in_mapping), nbins(in_nbins), sum(in_nbins, 0), count(in_nbins, 0) -{} +RadialIntegration::RadialIntegration(const std::vector& in_mapping, uint32_t in_nbins, uint32_t in_pixel_split) : + pixel_to_bin(in_mapping), nbins(in_nbins), sum(in_nbins, 0), count(in_nbins, 0), coeff(in_mapping.size(), 1.0f), + pixel_split(in_pixel_split) +{ + if (pixel_split == 1) + pixel_spit_weight = 1.0f; + else if (pixel_split == 4) { + pixel_spit_weight = 0.25f; + if (pixel_to_bin.size() % 4 != 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "With pixel split of 4 input array must be of size multiple of 4"); + } else + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Only pixel split of 1 and 4 allowed at the moment for radial integration"); +} -RadialIntegration::RadialIntegration(const RadialIntegrationMapping &mapping) : - RadialIntegration(mapping.GetPixelToBinMapping(), mapping.GetBinNumber()) +RadialIntegration::RadialIntegration(const RadialIntegrationMapping &mapping, uint32_t in_pixel_split) : + RadialIntegration(mapping.GetPixelToBinMapping(), mapping.GetBinNumber(), in_pixel_split) {} void RadialIntegration::Clear() { @@ -20,17 +32,34 @@ void RadialIntegration::Clear() { i = 0; } -void RadialIntegration::Process(const int16_t *data, size_t npixel) { - if (npixel != pixel_to_bin.size()) +void RadialIntegration::Process(const int16_t *__restrict data, size_t npixel) { + if (npixel != pixel_to_bin.size() / pixel_split) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Mismatch in size of pixel-to-bin mapping and image"); - for (int i = 0; i < npixel; i++) { - auto bin = pixel_to_bin[i]; - auto value = data[i]; - if ((value > INT16_MIN + 4) && (value < INT16_MAX - 4) && (bin < nbins)) { - sum[bin] += value; - count[bin] += 1; + if (pixel_split == 1) { + for (int i = 0; i < npixel; i++) { + auto value = data[i]; + if ((value > INT16_MIN + 4) && (value < INT16_MAX - 4)) { + auto bin = pixel_to_bin[i]; + if (bin < nbins) { + sum[bin] += coeff[i] * value; + count[bin] += 1.0f; + } + } + } + } else if (pixel_split == 4) { + for (int i = 0; i < npixel; i++) { + auto value = data[i]; + if ((value > INT16_MIN + 4) && (value < INT16_MAX - 4)) { + for (int p = 0; p < 4; p++) { + auto bin = pixel_to_bin[i * pixel_split + p]; + if (bin < nbins) { + sum[bin] += coeff[i] * value * 0.25f; + count[bin] += 0.25f; + } + } + } } } } @@ -51,19 +80,19 @@ void RadialIntegration::GetResult(std::vector &result) const { } } -const std::vector &RadialIntegration::GetSum() const { +const std::vector &RadialIntegration::GetSum() const { return sum; } -const std::vector &RadialIntegration::GetCount() const { +const std::vector &RadialIntegration::GetCount() const { return count; } -float RadialIntegration::GetRangeValue(uint16_t min_bin, uint16_t max_bin) { - int64_t ret_sum = 0; - int64_t ret_count = 0; +float RadialIntegration::GetRangeValue(uint32_t min_bin, uint32_t max_bin) { + float ret_sum = 0; + float ret_count = 0; - for (int i = std::min(nbins,min_bin); i <= std::min((uint16_t)(nbins-1),max_bin); i++) { + for (int i = std::min(nbins,min_bin); i <= std::min(nbins-1,max_bin); i++) { ret_sum += sum[i]; ret_count += count[i]; } diff --git a/image_analysis/RadialIntegration.h b/image_analysis/RadialIntegration.h index 80f18bbb..d77ddb88 100644 --- a/image_analysis/RadialIntegration.h +++ b/image_analysis/RadialIntegration.h @@ -12,20 +12,24 @@ #include "RadialIntegrationMapping.h" class RadialIntegration { - const std::vector& pixel_to_bin; - const uint16_t nbins; - std::vector sum; - std::vector count; + const std::vector pixel_to_bin; + const std::vector coeff; + const uint32_t nbins; + const uint32_t pixel_split; + float pixel_spit_weight; + float one_over_pixel_spit_weight; + std::vector sum; + std::vector count; public: - RadialIntegration(const RadialIntegrationMapping& mapping); - RadialIntegration(const std::vector& mapping, uint16_t nbins); + RadialIntegration(const RadialIntegrationMapping& mapping, uint32_t pixel_split = 1); + RadialIntegration(const std::vector& mapping, uint32_t nbins, uint32_t pixel_split = 1); void Clear(); void Process(const int16_t *data, size_t npixel); void ProcessOneImage(const int16_t *data, size_t npixel); // Process + Clear void GetResult(std::vector &result) const; - [[nodiscard]] float GetRangeValue(uint16_t min_bin, uint16_t max_bin); - [[nodiscard]] const std::vector& GetSum() const; - [[nodiscard]] const std::vector& GetCount() const; + [[nodiscard]] float GetRangeValue(uint32_t min_bin, uint32_t max_bin); + [[nodiscard]] const std::vector& GetSum() const; + [[nodiscard]] const std::vector& GetCount() const; }; diff --git a/image_analysis/RadialIntegrationMapping.cpp b/image_analysis/RadialIntegrationMapping.cpp index d7c2f883..2c65c887 100644 --- a/image_analysis/RadialIntegrationMapping.cpp +++ b/image_analysis/RadialIntegrationMapping.cpp @@ -79,4 +79,13 @@ const std::vector &RadialIntegrationMapping::GetSolidAngleCorr() const { double RadialIntegrationMapping::QToBin(double q) const { return std::min(static_cast(max_bin_number), std::max(0.0, (q - low_q) / q_spacing)); -} \ No newline at end of file +} + +std::vector RadialIntegrationMapping::GetPixelToBinMappingSplitTo4() const { + std::vector ret(pixel_to_bin.size() * 4); + for (int i = 0; i < pixel_to_bin.size(); i++) { + for (int j = 0; j < 4; j++) + ret[i * 4 + j] = pixel_to_bin[i]; + } + return ret; +} diff --git a/image_analysis/RadialIntegrationMapping.h b/image_analysis/RadialIntegrationMapping.h index 4c185b19..fb948591 100644 --- a/image_analysis/RadialIntegrationMapping.h +++ b/image_analysis/RadialIntegrationMapping.h @@ -17,6 +17,7 @@ public: RadialIntegrationMapping(const DiffractionExperiment& experiment, const uint8_t *one_byte_mask = nullptr); [[nodiscard]] uint16_t GetBinNumber() const; [[nodiscard]] const std::vector &GetPixelToBinMapping() const; + [[nodiscard]] std::vector GetPixelToBinMappingSplitTo4() const; [[nodiscard]] const std::vector &GetBinToQ() const; [[nodiscard]] const std::vector &GetSolidAngleCorr() const; [[nodiscard]] double QToBin(double q) const; diff --git a/tests/RadialIntegrationTest.cpp b/tests/RadialIntegrationTest.cpp index 1a816904..02e8e91b 100644 --- a/tests/RadialIntegrationTest.cpp +++ b/tests/RadialIntegrationTest.cpp @@ -133,6 +133,34 @@ TEST_CASE("RadialIntegration_Process","[RadialIntegration]") { REQUIRE(radial.GetSum()[1] == 6+2); } +TEST_CASE("RadialIntegration_Process_PxlSplit4","[RadialIntegration]") { + std::vector pixel_to_bin = {0, 1, 2, 3, 2, 2, 1, 4, 0, 0, 0, 0, 5, 5, 5, 5}; + std::vector test_image = {5, 7, INT16_MIN, 123}; + std::vector result; + + RadialIntegration radial(pixel_to_bin, 5, 4); + radial.ProcessOneImage(test_image.data(), 4); + + REQUIRE(radial.GetCount().size() == 5); + REQUIRE(radial.GetSum().size() == 5); + + float count = 0; + for (const auto &i: radial.GetCount()) + count += i; + REQUIRE(count == Approx(2.0f)); + + REQUIRE(radial.GetCount()[0] == Approx(0.25f)); + REQUIRE(radial.GetCount()[1] == Approx(0.5f)); + REQUIRE(radial.GetCount()[2] == Approx(0.75f)); + REQUIRE(radial.GetCount()[3] == Approx(0.25f)); + REQUIRE(radial.GetCount()[4] == Approx(0.25f)); + + REQUIRE(radial.GetSum()[0] == Approx(5 * 0.25f)); + REQUIRE(radial.GetSum()[1] == Approx(5 * 0.25f + 7 * 0.25f)); + REQUIRE(radial.GetSum()[2] == Approx(5 * 0.25f + 7 * 0.5f)); + REQUIRE(radial.GetSum()[4] == Approx(7 * 0.25f)); +} + TEST_CASE("RadialIntegration_GetResult","[RadialIntegration]") { std::vector pixel_to_bin = {0,1,2,4,3,1,2,3}; std::vector test_image = {7,6,5,4,3,2,1,0}; diff --git a/tools/DataAnalysisPerfTest.cpp b/tools/DataAnalysisPerfTest.cpp index d3c2ac3e..499bde6d 100644 --- a/tools/DataAnalysisPerfTest.cpp +++ b/tools/DataAnalysisPerfTest.cpp @@ -45,7 +45,7 @@ auto TestAll(const DiffractionExperiment &experiment, const JFJochProtoBuf::Data auto elapsed = std::chrono::duration_cast(end_time - start_time); std::ostringstream strstream; - logger.Info("{:20s} {:8.1f} ms/image", "Full", + logger.Info("{:30s} {:8.1f} ms/image", "Full", elapsed.count() / (1000.0 * (double) nimages)); return strstream.str(); @@ -85,7 +85,7 @@ auto TestAllWithROI(const DiffractionExperiment &experiment, const JFJochProtoBu auto elapsed = std::chrono::duration_cast(end_time - start_time); std::ostringstream strstream; - logger.Info("{:20s} {:8.1f} ms/image", "Full+ROI", + logger.Info("{:30s} {:8.1f} ms/image", "Full+ROI", elapsed.count() / (1000.0 * (double) nimages)); return strstream.str(); @@ -127,7 +127,7 @@ void TestIndexing() { auto end_time = std::chrono::system_clock::now(); auto elapsed = std::chrono::duration_cast(end_time - start_time); - logger.Info("Fast feedback index. {:8.1f} ms/image", elapsed.count() / (1000.0 * nexec)); + logger.Info("{:30s} {:8.1f} ms/image", "Fast feedback index.", elapsed.count() / (1000.0 * nexec)); } } @@ -149,7 +149,7 @@ auto TestSpotFinder(const DiffractionExperiment &experiment, const JFJochProtoBu auto elapsed = std::chrono::duration_cast(end_time - start_time); std::ostringstream strstream; - logger.Info("{:20s} {:8.1f} ms/image {:5d} spots", "Spot finding", + logger.Info("{:30s} {:8.1f} ms/image {:5d} spots", "Spot finding", elapsed.count() / (1000.0 * (double) nimages), spots.size()); return strstream.str(); @@ -173,7 +173,7 @@ auto TestSpotFinderWithoutCopyToDevice(const DiffractionExperiment &experiment, auto elapsed = std::chrono::duration_cast(end_time - start_time); std::ostringstream strstream; - logger.Info("{:20s} {:8.1f} ms/image {:5d} spots", "Spot finding", + logger.Info("{:30s} {:8.1f} ms/image {:5d} spots", "Spot finding", elapsed.count() / (1000.0 * (double) nimages), spots.size()); return strstream.str(); @@ -237,18 +237,29 @@ auto TestRadialIntegrationGPUWithoutCopyToDevice(const DiffractionExperiment &x, } -auto TestRadialIntegration(const DiffractionExperiment &experiment, const JFJochProtoBuf::DataProcessingSettings &settings, - int16_t* image, size_t nimages) { +auto TestRadialIntegration(const DiffractionExperiment &experiment, + const JFJochProtoBuf::DataProcessingSettings &settings, + int16_t* image, size_t nimages, + uint32_t pixel_split = 1) { uint32_t nredo = 20; RadialIntegrationMapping mapping(experiment); - RadialIntegration integration(mapping); + std::unique_ptr integration; + + if (pixel_split == 1) { + integration = std::make_unique(mapping); + } else { + integration = std::make_unique(mapping.GetPixelToBinMappingSplitTo4(), + mapping.GetBinNumber(), + 4); + } + std::vector result; auto start_time = std::chrono::system_clock::now(); for (int redo = 0; redo < nredo; redo++) { for (int i = 0; i < nimages; i++) { - integration.ProcessOneImage(image + i * experiment.GetPixelsNum(), experiment.GetPixelsNum()); + integration->ProcessOneImage(image + i * experiment.GetPixelsNum(), experiment.GetPixelsNum()); //integration.GetResult(result); } } @@ -332,14 +343,17 @@ int main(int argc, char **argv) { TestSpotFinder(x, settings, local_peakfinder_gpu,image_conv.data(), nimages); } - logger.Info("{:20s} {:8.1f} ms/image", "Radial int. (GPU)", TestRadialIntegrationGPU(x, settings, + logger.Info("{:30s} {:8.1f} ms/image", "Radial int. (GPU)", TestRadialIntegrationGPU(x, settings, image_conv.data(), nimages)); - logger.Info("{:20s} {:8.1f} ms/image", "Radial int. (GPU/nocopy)", TestRadialIntegrationGPUWithoutCopyToDevice(x, settings, + logger.Info("{:30s} {:8.1f} ms/image", "Radial int. (GPU/nocpy)", TestRadialIntegrationGPUWithoutCopyToDevice(x, settings, image_conv.data(), nimages)); - logger.Info("{:20s} {:8.1f} ms/image", "Radial int. (CPU)", TestRadialIntegration(x, settings, + logger.Info("{:30s} {:8.1f} ms/image", "Radial int. (CPU)", TestRadialIntegration(x, settings, image_conv.data(), nimages)); + logger.Info("{:30s} {:8.1f} ms/image", "Radial int. pxlspl 2 (CPU)", TestRadialIntegration(x, settings, + image_conv.data(), nimages, 4)); + TestIndexing(); logger.Info("Full package");