From c00970db0a2461a6a59a41c758188ea55100d81b Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Sat, 4 Oct 2025 15:10:02 +0200
Subject: [PATCH] ImageAnalysisCPU: Optimize azimuthal integration calculations
 to increase performance

---
 common/AzimuthalIntegrationSettings.cpp | 16 ++++++++++++++--
 common/AzimuthalIntegrationSettings.h   |  6 ++++++
 image_analysis/ImageAnalysisCPU.cpp     | 15 ++++++++-------
 image_analysis/ImageAnalysisCPU.h       |  2 ++
 4 files changed, 30 insertions(+), 9 deletions(-)

diff --git a/common/AzimuthalIntegrationSettings.cpp b/common/AzimuthalIntegrationSettings.cpp
index 717a70a6..ebdc4275 100644
--- a/common/AzimuthalIntegrationSettings.cpp
+++ b/common/AzimuthalIntegrationSettings.cpp
@@ -10,6 +10,10 @@
 #define check_min(param, val, min) if ((val) < (min)) throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin, param)
 #define check_finite(param, val) if (!std::isfinite(val)) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, param)
 
+AzimuthalIntegrationSettings::AzimuthalIntegrationSettings() {
+    UpdateBinCount();
+}
+
 AzimuthalIntegrationSettings &AzimuthalIntegrationSettings::SolidAngleCorrection(bool input) {
     solid_angle_correction = input;
     return *this;
@@ -26,6 +30,7 @@ AzimuthalIntegrationSettings &AzimuthalIntegrationSettings::QRange_recipA(float
 
     high_q_recipA = high;
     low_q_recipA = low;
+    UpdateBinCount();
     return *this;
 }
 
@@ -33,6 +38,7 @@ AzimuthalIntegrationSettings &AzimuthalIntegrationSettings::QSpacing_recipA(floa
     check_finite("Q spacing for azimuthal integration", input);
     check_min("Q spacing for azimuthal integration", input, 0.00999);
     q_spacing = input;
+    UpdateBinCount();
     return *this;
 }
 
@@ -88,11 +94,12 @@ AzimuthalIntegrationSettings &AzimuthalIntegrationSettings::AzimuthalBinCount(in
     check_min("Azimuthal bin count", input, 1);
     check_max("Azimuthal bin count", input, 512);
     azim_bins = input;
+    UpdateBinCount();
     return *this;
 }
 
 int32_t AzimuthalIntegrationSettings::GetQBinCount() const {
-    return std::ceil((high_q_recipA - low_q_recipA) / q_spacing);
+    return q_bins;
 }
 
 int32_t AzimuthalIntegrationSettings::GetAzimuthalBinCount() const {
@@ -100,7 +107,7 @@ int32_t AzimuthalIntegrationSettings::GetAzimuthalBinCount() const {
 }
 
 int32_t AzimuthalIntegrationSettings::GetBinCount() const {
-    return GetQBinCount() * GetAzimuthalBinCount();
+    return total_bins;
 }
 
 uint16_t AzimuthalIntegrationSettings::QToBin(float q) const {
@@ -117,3 +124,8 @@ uint16_t AzimuthalIntegrationSettings::GetBin(float q, float phi_deg) const {
     int16_t phi_bin = std::floor(phi_deg / 360.0f * azim_bins);
     return q_bin + phi_bin * GetQBinCount();
 }
+
+void AzimuthalIntegrationSettings::UpdateBinCount() {
+    q_bins = std::ceil((high_q_recipA - low_q_recipA) / q_spacing);
+    total_bins = q_bins * azim_bins;
+}
diff --git a/common/AzimuthalIntegrationSettings.h b/common/AzimuthalIntegrationSettings.h
index d19b96e0..02ddea55 100644
--- a/common/AzimuthalIntegrationSettings.h
+++ b/common/AzimuthalIntegrationSettings.h
@@ -16,7 +16,13 @@ class AzimuthalIntegrationSettings {
     float bkg_estimate_low_q_recipA = 2.0f * M_PI / 5.0;
     float q_spacing = 0.05;
     int32_t azim_bins = 1;
+
+    int32_t q_bins= 0;
+    int32_t total_bins = 0;
+
+    void UpdateBinCount();
 public:
+    AzimuthalIntegrationSettings();
     AzimuthalIntegrationSettings& SolidAngleCorrection(bool input);
     AzimuthalIntegrationSettings& PolarizationCorrection(bool input);
     AzimuthalIntegrationSettings& QRange_recipA(float low, float high);
diff --git a/image_analysis/ImageAnalysisCPU.cpp b/image_analysis/ImageAnalysisCPU.cpp
index 5021662e..79d1f164 100644
--- a/image_analysis/ImageAnalysisCPU.cpp
+++ b/image_analysis/ImageAnalysisCPU.cpp
@@ -22,7 +22,8 @@ ImageAnalysisCPU::ImageAnalysisCPU(const DiffractionExperiment &in_experiment,
       spotFinder(in_integration),
       indexer(in_indexer),
       saturation_limit(experiment.GetSaturationLimit()),
-      mask(in_mask) {
+      mask(in_mask),
+      azint_bins(in_integration.GetBinNumber()) {
 
     roi_map = experiment.ExportROIMap();
     roi_count = experiment.ROI().size();
@@ -77,9 +78,9 @@ void ImageAnalysisCPU::Analyze(DataMessage &output,
 
     std::vector<ROIMessage> roi(roi_count);
 
-    std::vector<float> azim_sum(integration.GetBinNumber(), 0.0f);
-    std::vector<float> azim_sum2(integration.GetBinNumber(), 0.0f);
-    std::vector<uint32_t> azim_count(integration.GetBinNumber(), 0);
+    std::vector<float> azim_sum(azint_bins, 0.0f);
+    //std::vector<float> azim_sum2(integration.GetBinNumber(), 0.0f);
+    std::vector<uint32_t> azim_count(azint_bins, 0);
 
     size_t err_pixels = 0;
     size_t masked_pixels = 0;
@@ -130,11 +131,11 @@ void ImageAnalysisCPU::Analyze(DataMessage &output,
                     }
                 }
             }
-            uint16_t bin = pixel_to_bin[i];
-            if (bin < integration.GetBinNumber()) {
+            const uint16_t bin = pixel_to_bin[i];
+            if (bin < azint_bins) {
                 float val = image[i] * corrections[i];
                 azim_sum[bin] += val;
-                azim_sum2[bin] += val * val;
+                //azim_sum2[bin] += val * val;
                 ++azim_count[bin];
             }
         }
diff --git a/image_analysis/ImageAnalysisCPU.h b/image_analysis/ImageAnalysisCPU.h
index 96f81d61..39a27f13 100644
--- a/image_analysis/ImageAnalysisCPU.h
+++ b/image_analysis/ImageAnalysisCPU.h
@@ -31,6 +31,8 @@ class ImageAnalysisCPU {
     ImageSpotFinder spotFinder;
     IndexerThreadPool *indexer;
 
+    uint16_t azint_bins;
+
     const int64_t saturation_limit;
 
     const PixelMask &mask;