From c55e91ccfa801311e2d4b9e9ab4faf133d3faaaf Mon Sep 17 00:00:00 2001
From: leonarski_f <filip.leonarski@psi.ch>
Date: Wed, 13 May 2026 07:48:17 +0200
Subject: [PATCH] Merge: Simplify (but not yet gaining performance)

---
 image_analysis/scale_merge/Merge.cpp | 63 ++++++++++------------------
 1 file changed, 22 insertions(+), 41 deletions(-)

diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp
index 9f049362..e5ef98b1 100644
--- a/image_analysis/scale_merge/Merge.cpp
+++ b/image_analysis/scale_merge/Merge.cpp
@@ -71,46 +71,15 @@ namespace {
         return out;
     }
 
-    MergeResult InitResult(const std::vector<HKLKey> &slot_to_hkl,
-                                const std::vector<Obs> &obs) {
-        MergeResult out;
-        out.merged.resize(slot_to_hkl.size());
-
-        for (int i = 0; i < static_cast<int>(slot_to_hkl.size()); ++i) {
-            out.merged[i].h = slot_to_hkl[i].h;
-            out.merged[i].k = slot_to_hkl[i].k;
-            out.merged[i].l = slot_to_hkl[i].l;
-            out.merged[i].I = 0.0;
-            out.merged[i].sigma = 0.0;
-            out.merged[i].d = 0.0;
-        }
-
-        std::vector<std::vector<double>> d_values(slot_to_hkl.size());
-        for (const auto &o: obs) {
-            if (std::isfinite(o.r->d) && o.r->d > 0.0f)
-                d_values[o.hkl].push_back(o.r->d);
-        }
-
-        for (int h = 0; h < static_cast<int>(d_values.size()); ++h) {
-            auto &v = d_values[h];
-            if (v.empty())
-                continue;
-
-            std::nth_element(v.begin(), v.begin() + static_cast<long>(v.size() / 2), v.end());
-            out.merged[h].d = v[v.size() / 2];
-        }
-
-        return out;
-    }
-
-    void Merge(size_t nhkl, MergeResult &out, const std::vector<Obs> &obs) {
+    MergeResult Merge(const std::vector<HKLKey> &slot_to_hkl, const std::vector<Obs> &obs) {
         struct Accum {
             double sum_wI = 0.0;
             double sum_w = 0.0;
             double sum_wsigma2 = 0.0;
+            std::vector<double> d_values;
         };
 
-        std::vector<Accum> acc(nhkl);
+        std::vector<Accum> acc(slot_to_hkl.size());
 
         for (const auto &o: obs) {
             const double I_corr = static_cast<double>(o.r->I) * o.r->scaling_correction;
@@ -126,16 +95,31 @@ namespace {
             a.sum_wI += w * I_corr;
             a.sum_w += w;
             a.sum_wsigma2 += w * w * sigma_corr * sigma_corr;
+            if (std::isfinite(o.r->d) && o.r->d > 0.0f)
+                a.d_values.push_back(o.r->d);
+
         }
 
-        for (int h = 0; h < static_cast<int>(nhkl); ++h) {
-            const auto &a = acc[h];
+        MergeResult out;
+        out.merged.resize(slot_to_hkl.size());
+
+        for (int h = 0; h < static_cast<int>(slot_to_hkl.size()); ++h) {
+            auto &a = acc[h];
             if (a.sum_w <= 0.0)
                 continue;
-
+            out.merged[h].h = slot_to_hkl[h].h;
+            out.merged[h].k = slot_to_hkl[h].k;
+            out.merged[h].l = slot_to_hkl[h].l;
             out.merged[h].I = a.sum_wI / a.sum_w;
             out.merged[h].sigma = std::sqrt(a.sum_wsigma2) / a.sum_w;
+            if (a.d_values.empty())
+                continue;
+
+            std::ranges::nth_element(a.d_values,
+                a.d_values.begin() + static_cast<long>(a.d_values.size() / 2));
+            out.merged[h].d = a.d_values[a.d_values.size() / 2];
         }
+        return out;
     }
 
     void Stats(const DiffractionExperiment &x, MergeResult &out, const std::vector<Obs> &obs) {
@@ -281,10 +265,7 @@ MergeResult MergeReflections(const std::vector<std::vector<Reflection> > &observ
                              const DiffractionExperiment &x) {
     std::vector<HKLKey> slot_to_hkl;
     auto obs = BuildObservations(observations, x, slot_to_hkl);
-
-    auto out = InitResult(slot_to_hkl, obs);
-
-    Merge(slot_to_hkl.size(), out, obs);
+    auto out = Merge(slot_to_hkl, obs);
     Stats(x, out, obs);
 
     return out;