Parallelise the incident-flux divide, drop a redundant sync

DivideOutIncidentFlux was still the last fully serial pass in Ingest: a sweep
over every observation to take each frame's mean background, and another to
divide every rlp by its frame's flux. Ten gigabytes of traffic on one thread.

The per-frame means go a frame at a time rather than an observation at a time,
so each frame's running sum stays in one thread and in the order it had -
splitting by observation would cut a frame across two threads and the partial
sums would have to be recombined, which is a different sequence of roundings.
The divide is per-element and splits anywhere.

The adaptive spot finder synchronised after flagging strong pixels. The
extractor that reads those pixels runs on the same stream, so the ordering
already guaranteed the flagging had finished; the wait only idled the host,
once per image.

Measured on a crystal with 66 million partial observations: Ingest 8.5 s and
7.7 s -> 7.1 s and 6.6 s, whole crystal 1m24s -> 1m17s. Merged statistics
unchanged.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu

@@ -767,6 +767,8 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) {
     h_var_bkg.resize(new_capacity);
     h_obs_x.resize(new_capacity); h_obs_y.resize(new_capacity);
     h_ok.resize(new_capacity); h_has_obs.resize(new_capacity);
+
+
     capacity = new_capacity;
 }
 

image_analysis/bragg_integration/BraggIntegrationEngineGPU.h

@@ -60,6 +60,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine {
     std::vector<float>   h_I, h_sigma, h_bkg, h_var_bkg, h_obs_x, h_obs_y;
     std::vector<uint8_t> h_ok, h_has_obs;
 
+
     void EnsureCapacity(size_t n);
 
 public:

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -203,16 +203,24 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment,
 void RotationScaleMerge::DivideOutIncidentFlux() {
     frame_flux.assign(n_frames, 1.0);
 
-    std::vector<double> mean_bkg(n_frames, NAN), finite;
-    for (int f = 0; f < n_frames; ++f) {
+    // A frame at a time, not an observation at a time: each frame's mean is a running sum over its
+    // own contiguous range, so giving a thread whole frames keeps every sum in one thread and in the
+    // order it had. Splitting by observation would cut a frame across two threads and the sum would
+    // have to be recombined, which is a different sequence of roundings.
+    std::vector<double> mean_bkg(n_frames, NAN);
+    ParallelFor(n_frames, nthreads, [&](int f) {
         double sum = 0.0;
         int n = 0;
         for (int i = frame_start[f]; i < frame_start[f] + frame_count[f]; ++i)
             if (std::isfinite(partials[i].bkg) && partials[i].bkg > 0.0f) { sum += partials[i].bkg; ++n; }
-        if (n == 0) continue;
-        mean_bkg[f] = sum / n;
-        finite.push_back(mean_bkg[f]);
-    }
+        if (n > 0) mean_bkg[f] = sum / n;
+    });
+    // Collected in frame order, serially: it is one value per frame, and the median only depends on
+    // the set, but keeping the order removes the question.
+    std::vector<double> finite;
+    finite.reserve(n_frames);
+    for (int f = 0; f < n_frames; ++f)
+        if (std::isfinite(mean_bkg[f])) finite.push_back(mean_bkg[f]);
     if (finite.empty())
         return;
 
@@ -224,8 +232,12 @@ void RotationScaleMerge::DivideOutIncidentFlux() {
     for (int f = 0; f < n_frames; ++f)
         if (std::isfinite(mean_bkg[f])) frame_flux[f] = mean_bkg[f] / median;
 
-    for (auto &o : partials)
-        o.rlp = static_cast<float>(o.rlp / frame_flux[o.frame]);
+    ParallelChunks(static_cast<int>(partials.size()), nthreads, [&](int lo, int hi) {
+        for (int i = lo; i < hi; ++i) {
+            auto &o = partials[i];
+            o.rlp = static_cast<float>(o.rlp / frame_flux[o.frame]);
+        }
+    });
 }
 
 void RotationScaleMerge::Ingest() {

image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu

@@ -441,8 +441,8 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image,
     cuda_err(cudaMemsetAsync(gpu_strong, 0, OutputByteSize(), *stream));
     flag_strong<<<flag_blocks, flag_threads, 0, *stream>>>(
         image.getGPUBuffer(), gpu_pixel_to_bin->get(), gpu_thr, gpu_strong, npix, nbins);
-    // The bit buffer stays on the device - ExtractComponents reads it there.
-    cuda_err(cudaStreamSynchronize(*stream));
+    // The bit buffer stays on the device and ExtractComponents reads it there, on this same stream,
+    // so the ordering already guarantees flag_strong has finished. Waiting here only idled the host.
 }
 
 void AdaptiveSpotFinderGPU::SetResolutionMask(const std::vector<bool> &mask) {