Gain: the frames of the gain measurement are read and measured side by side

RugnuxMeasureGain read and measured its 16 frames one after another, 1.25-1.5 s
of step 1 on every CCD/SMV input. Each frame now runs on its own thread with its
own levels and histogram, put together in frame order: the same integer sums and
the same list of levels, so the same gain. Step 1 on marCCD/SMV sweeps 1.3 -> 0.5 s.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UqXYdL4qBGYRiSXSGXfCLu

rugnux/RugnuxDefaults.cpp

@@ -5,6 +5,9 @@
 
 #include <algorithm>
 #include <cstring>
+#include <optional>
+
+#include "../common/ParallelFor.h"
 
 namespace {
     // Smoothing range for the per-frame scale G on rotation data (XDS DELPHI-like), in degrees.
@@ -125,18 +128,20 @@ bool RugnuxDetectIceRings(std::optional<bool> file_detect_ice_rings, bool rotati
     return file_detect_ice_rings.value_or(rotation);
 }
 
-RugnuxMeasuredGain RugnuxMeasureGain(JFJochReader &reader) {
+RugnuxMeasuredGain RugnuxMeasureGain(JFJochReader &reader, size_t nthreads) {
     const auto dataset = reader.GetDataset();
     const size_t n = reader.GetNumberOfImages();
     const size_t frames = std::min<size_t>(n, background_gain::GAIN_MAX_FRAMES);
-    std::vector<background_gain::Levels> levels;
-    std::vector<uint64_t> histogram(background_gain::PEAK_HISTOGRAM_BINS, 0);
-    JFJochReaderRawImage raw;
-    std::vector<uint8_t> buffer;
-    std::vector<int32_t> pixels;
-    for (size_t i = 0; i < frames; i++) {
+    // The frames are read and measured side by side; each keeps its own levels and histogram, which are
+    // put together in frame order below.
+    std::vector<std::optional<background_gain::Levels>> frame_levels(frames);
+    std::vector<std::vector<uint64_t>> frame_histogram(frames);
+    ParallelFor(static_cast<int>(frames), std::min(frames, nthreads), [&](int i) {
+        JFJochReaderRawImage raw;
+        std::vector<uint8_t> buffer;
+        std::vector<int32_t> pixels;
         if (!reader.ReadRawImage(static_cast<int64_t>(i * n / frames), raw))
-            continue;
+            return;
         const size_t w = raw.image.GetWidth(), h = raw.image.GetHeight();
         pixels.resize(w * h);
         const uint8_t *ptr = raw.image.GetUncompressedPtr(buffer);
@@ -147,11 +152,20 @@ RugnuxMeasuredGain RugnuxMeasureGain(JFJochReader &reader) {
         else if (raw.image.GetMode() == CompressedImageMode::Uint16)
             std::copy_n(reinterpret_cast<const uint16_t *>(ptr), w * h, pixels.begin());
         else
-            continue;
-        levels.push_back(background_gain::MeasureFrame(pixels, w, h, dataset->pixel_mask->GetMask(),
-                                                       dataset->experiment.GetSaturationLimit()));
+            return;
+        frame_levels[i] = background_gain::MeasureFrame(pixels, w, h, dataset->pixel_mask->GetMask(),
+                                                        dataset->experiment.GetSaturationLimit());
+        frame_histogram[i].assign(background_gain::PEAK_HISTOGRAM_BINS, 0);
         background_gain::AddToHistogram(pixels, dataset->pixel_mask->GetMask(),
-                                        dataset->experiment.GetSaturationLimit(), histogram);
+                                        dataset->experiment.GetSaturationLimit(), frame_histogram[i]);
+    });
+    std::vector<background_gain::Levels> levels;
+    std::vector<uint64_t> histogram(background_gain::PEAK_HISTOGRAM_BINS, 0);
+    for (size_t i = 0; i < frames; i++) {
+        if (frame_levels[i])
+            levels.push_back(*frame_levels[i]);
+        for (size_t b = 0; b < frame_histogram[i].size(); b++)
+            histogram[b] += frame_histogram[i][b];
     }
     return {background_gain::PhotonPeak(histogram), background_gain::Gain(levels)};
 }

rugnux/RugnuxDefaults.h

@@ -84,7 +84,7 @@ struct RugnuxMeasuredGain {
     std::optional<background_gain::GainEstimate> peaks;
     std::optional<float> variance;
 };
-RugnuxMeasuredGain RugnuxMeasureGain(JFJochReader &reader);
+RugnuxMeasuredGain RugnuxMeasureGain(JFJochReader &reader, size_t nthreads);
 
 // The gain the frames are analysed at, into everything that takes a photon count from them: the
 // integration (BraggIntegrationSettings::pixel_gain, which the defective-pixel tests read too) and the

rugnux/rugnux_cli.cpp

@@ -2156,7 +2156,7 @@ int RugnuxInputRun::PrepareRugnuxInput() {
         const float unit = float_frames ? 1.0f / 16.0f : 1.0f;
         if (float_frames)
             xfel_reader.ApplyGain(unit);
-        const RugnuxMeasuredGain g = RugnuxMeasureGain(reader);
+        const RugnuxMeasuredGain g = RugnuxMeasureGain(reader, nthreads);
         if (float_frames)
             xfel_reader.ClearGain();
         if (g.peaks) {

viewer/JFJochProcessController.cpp

@@ -83,7 +83,7 @@ void JFJochProcessController::run_(QString file_path, DiffractionExperiment expe
         JFJochReader &reader = *reader_ptr;
         // A CCD's pixel values are ADU: analysed at the gain its background shows, as rugnux does.
         if (reader_ptr == &marccd_reader || reader_ptr == &smv_reader) {
-            const RugnuxMeasuredGain g = RugnuxMeasureGain(reader);
+            const RugnuxMeasuredGain g = RugnuxMeasureGain(reader, static_cast<size_t>(std::max(1, config.nthreads)));
             const std::optional<float> gain = g.peaks ? std::optional<float>(g.peaks->gain) : g.variance;
             ApplyRugnuxPixelGain(experiment, config.spot_finding, gain.value_or(1.0f));
         }