From b3039245f97d8d8fb19ca1343b02df7ddfbe6ecb Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Mon, 19 Oct 2026 11:21:25 +0200
Subject: [PATCH] BeamCenterFFT: the 99.9th percentile and median of the
 projection by radix select instead of nth_element

PrepareImage copied the finite pixels and ran two nth_element over them (18 M pixels on a 16M
detector, ~200 ms single-threaded). Two counting passes over the float bits (upper, then lower
16 bits of an order-preserving key) give the same k-th values in ~55 ms without the copy. Values
are identical; only the sign of a zero picked as a quantile may differ, which the clipping maps
to the same prepared image.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4
---
 .../geom_refinement/BeamCenterFFT.cpp         | 58 +++++++++++++++----
 1 file changed, 46 insertions(+), 12 deletions(-)

diff --git a/image_analysis/geom_refinement/BeamCenterFFT.cpp b/image_analysis/geom_refinement/BeamCenterFFT.cpp
index ece2d6fb1..53c4468ac 100644
--- a/image_analysis/geom_refinement/BeamCenterFFT.cpp
+++ b/image_analysis/geom_refinement/BeamCenterFFT.cpp
@@ -4,7 +4,9 @@
 #include "BeamCenterFFT.h"
 
 #include <algorithm>
+#include <bit>
 #include <cmath>
+#include <cstdint>
 #include <limits>
 #include <stdexcept>
 
@@ -30,26 +32,58 @@ struct Prepared {
     double variance = 0.0;  // of a over the valid pixels, after the whole preparation
 };
 
+// Order-preserving: unsigned comparison of the keys orders the floats as < does (-0 below +0).
+uint32_t SortKey(float v) {
+    const uint32_t u = std::bit_cast<uint32_t>(v);
+    return (u & 0x80000000u) ? ~u : (u | 0x80000000u);
+}
+
+// For each k, the value nth_element would leave at position k of the finite values of x - without
+// copying or reordering the image. Counting the finite values by the upper 16 bits of their key
+// finds the bucket holding the k-th; counting that bucket by the lower 16 bits finds its key,
+// which is the value. Two passes over x for all k together.
+std::vector<float> KthFinite(const std::vector<float> &x, std::vector<size_t> k) {
+    std::vector<size_t> upper(1 << 16, 0);
+    for (const float v : x)
+        if (std::isfinite(v))
+            upper[SortKey(v) >> 16]++;
+    std::vector<uint32_t> bucket(k.size(), 0);
+    for (size_t i = 0; i < k.size(); i++)
+        while (k[i] >= upper[bucket[i]])
+            k[i] -= upper[bucket[i]++];
+
+    std::vector<std::vector<size_t>> lower(k.size(), std::vector<size_t>(1 << 16, 0));
+    for (const float v : x)
+        if (std::isfinite(v))
+            for (size_t i = 0; i < k.size(); i++)
+                if (SortKey(v) >> 16 == bucket[i])
+                    lower[i][SortKey(v) & 0xFFFF]++;
+    std::vector<float> out;
+    for (size_t i = 0; i < k.size(); i++) {
+        uint32_t b = 0;
+        while (k[i] >= lower[i][b])
+            k[i] -= lower[i][b++];
+        const uint32_t key = (bucket[i] << 16) | b;
+        out.push_back(std::bit_cast<float>((key & 0x80000000u) ? (key & 0x7FFFFFFFu) : ~key));
+    }
+    return out;
+}
+
 Prepared PrepareImage(const std::vector<float> &mean) {
     Prepared out;
     const size_t n = mean.size();
     out.a.assign(n, 0.0f);
     out.m.assign(n, 0.0f);
 
-    std::vector<float> valid;
-    valid.reserve(n);
-    for (float v : mean)
-        if (std::isfinite(v))
-            valid.push_back(v);
-    if (valid.size() < 2)
+    const size_t n_valid = std::count_if(mean.begin(), mean.end(), [](float v) { return std::isfinite(v); });
+    if (n_valid < 2)
         return out;
 
-    auto kth = [&valid](size_t k) {
-        std::nth_element(valid.begin(), valid.begin() + k, valid.end());
-        return valid[k];
-    };
-    const float clip_hi = kth(static_cast<size_t>(0.999 * (valid.size() - 1)));
-    const float median = kth(valid.size() / 2);
+    // Where these fall among equal values a -0 can come out for a +0, which the clipping below
+    // turns into the same image.
+    const std::vector<float> q = KthFinite(mean, {static_cast<size_t>(0.999 * (n_valid - 1)), n_valid / 2});
+    const float clip_hi = q[0];
+    const float median = q[1];
 
     double sum = 0.0;
     size_t count = 0;