Merge stills-satspot into rc178: per-spot saturation on stills (stuck pixels masked, lower bound only with pixels at the limit), run-level gate removed

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

common/Reflection.h

@@ -62,9 +62,6 @@ struct Reflection {
     // intensity at all - its brightest pixels are missing - and a rotation merge drops the whole
     // rocking event it belongs to, as XDS does with an overloaded reflection.
     bool overloaded = false;
-    // A still's overloaded reflection kept as a lower bound instead (BraggIntegrationEngine::Finalize);
-    // the run may still hand it back to the overload drop (RugnuxScaleMerge).
-    bool overload_kept = false;
 };
 
 // One full reflection of a rotation sweep as the merge used it - its partials summed into one

image_analysis/bragg_integration/BraggIntegrationEngine.cpp

@@ -284,11 +284,14 @@ std::vector<Reflection> BraggIntegrationEngine::Finalize(const std::vector<Refle
         // So a still keeps the observation as a lower bound: each saturated pixel counted at the
         // saturation limit. (A fit to the unsaturated wings was tried and read lower still: the
         // profile learned from unsaturated spots does not describe the wings of a clipped one.)
-        if (fr.overloaded && keep_overloads) {
+        // That takes a pixel AT the limit: an observation overloaded only by pixels unreadable on
+        // this frame has nothing to count at the limit, and counting those pixels as zero is no
+        // lower bound of anything - it stays overloaded and is dropped. (A pixel stuck at the limit
+        // is no clipped spot either; the pre-scan masks it, HotPixelFinder.)
+        if (fr.overloaded && keep_overloads && fr.saturated > 0) {
             refl.I = fr.I + static_cast<float>(fr.saturated * (saturation_limit - fr.bkg));
             refl.sigma = std::sqrt(fr.var_bkg + std::max(refl.I, 1.0f));
             refl.overloaded = false;
-            refl.overload_kept = true;
             ++counts.overloads_kept;
         }
         refl.observed = true;

image_analysis/hot_pixels/HotPixels.cpp

@@ -203,7 +203,8 @@ void HotPixelFinder::AddImage(const int32_t *image, std::vector<int32_t> &scratc
         for (size_t i = begin; i < end; i++) {
             const int32_t k = key[i], v = image[i];
             if (k < 0) continue;
-            if (v == INT32_MIN) {
+            // The error value and the saturation limit alike: no count (see the class comment).
+            if (v == INT32_MIN || v == INT32_MAX) {
                 n_error[i]++;
                 if (ring_ok[k / SECTORS]) {   // counted in its ring-sector's frames, but not valid here
                     n_error_ring_ok[i]++;
@@ -213,7 +214,7 @@ void HotPixelFinder::AddImage(const int32_t *image, std::vector<int32_t> &scratc
             }
             if (!ring_ok[k / SECTORS]) continue;
             sum_value[i] += v;
-            if (v == INT32_MAX || static_cast<float>(v) > threshold[k])
+            if (static_cast<float>(v) > threshold[k])
                 n_lit[i]++;
         }
     }

image_analysis/hot_pixels/HotPixels.h

@@ -31,8 +31,10 @@
 // is above the Poisson bound: a weaker one cannot make an outlier. Alone or as part of a patch - a
 // patch lit through the whole sweep is stationary in the lab, so it is not a reflection of the
 // rotating crystal, and a reflection crossing it reads the patch. A pixel holding the
-// detector's error value on most frames is masked with it: every frame already treats it as invalid,
-// but the static mask did not know.
+// detector's error value or its saturation limit on most frames is masked with it: every frame already
+// treats it as invalid, but the static mask did not know. A pixel stuck at the limit is no overload -
+// a CCD can have thousands, at 65535 on every frame - and unmasked, every reflection whose disk it
+// falls in is flagged as overloaded.
 //
 // A DEAD pixel is the mirror of a hot one: it reads on average at most half of its own sector's level,
 // and its mean deficit is above the bound a hot pixel's excess has to clear. Nothing but a defect keeps
@@ -85,7 +87,7 @@ public:
     struct Result {
         std::vector<uint32_t> mask;   // non-zero = masked (1 hot, 2 error value, 3 dead), converted geometry
         size_t hot = 0;               // persistent above the ring
-        size_t error = 0;             // holding the error value on most frames
+        size_t error = 0;             // holding the error value or the saturation limit on most frames
         size_t dead = 0;              // reading at most half of its sector's level
         uint32_t frames = 0;
     };

image_analysis/hot_pixels/HotPixelsGPU.cu

@@ -132,7 +132,7 @@ __global__ void accumulate_kernel(const int32_t *__restrict__ image, const int32
          i += static_cast<size_t>(blockDim.x) * gridDim.x) {
         const int32_t k = key[i], v = image[i];
         if (k < 0) continue;
-        if (v == INT32_MIN) {
+        if (v == INT32_MIN || v == INT32_MAX) {   // the error value or the saturation limit: no count
             n_error[i]++;
             if (ring_ok[k / sectors]) {
                 n_error_ring_ok[i]++;
@@ -142,7 +142,7 @@ __global__ void accumulate_kernel(const int32_t *__restrict__ image, const int32
         }
         if (!ring_ok[k / sectors]) continue;
         sum_value[i] += v;
-        if (v == INT32_MAX || static_cast<float>(v) > threshold[k])
+        if (static_cast<float>(v) > threshold[k])
             n_lit[i]++;
     }
 }

rugnux/RugnuxPreScan.cpp

@@ -1218,7 +1218,7 @@ void Rugnux::MaskDefectivePixels(int start_image, const std::vector<int> &sample
         pixel_mask_.LoadHotPixelMask(experiment_, mask);
         const size_t masked = defects->hot + defects->error + defects->dead;
         logger.Info("Defective pixels: {} masked from {} pre-scan frames - {} lit above their ring on more "
-                    "frames than one reflection or chance explains, {} holding the detector's error value, {} dead",
+                    "frames than one reflection or chance explains, {} holding the detector's error value or saturation limit, {} dead",
                     masked, defects->frames, defects->hot, defects->error, defects->dead);
         Note(fmt::format("{} defective pixel{} masked", masked, masked == 1 ? "" : "s"));
     };

rugnux/RugnuxScaleMerge.cpp

@@ -549,45 +549,6 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
     size_t ice_flagged = 0, ice_total = 0;
     {
         auto &outcomes = indexer->GetIntegrationOutcome();
-        // Saturation clips the strongest spots, and those are at low resolution: kept overloads of a
-        // still are a measurement only where they concentrate there. Where they do not, the saturated
-        // pixels are something else (hot or flagged pixels across the pattern), and an observation
-        // with such a pixel set to the saturation limit is not a lower bound but an outlier - so they
-        // go back to the overload drop. A pragmatic gate until saturation is told from bad pixels per
-        // spot.
-        {
-            std::atomic<int64_t> k_all{0}, k_low{0}, n_all{0}, n_low{0};
-            ParallelFor(static_cast<int>(outcomes.size()), config_.nthreads, [&](int o) {
-                int64_t all = 0, low = 0, kept = 0, kept_low = 0;
-                for (const auto &r : outcomes[o].reflections) {
-                    const bool low_res = r.d >= OVERLOAD_LOW_RES_A;
-                    ++all;
-                    low += low_res;
-                    if (r.overload_kept) {
-                        ++kept;
-                        kept_low += low_res;
-                    }
-                }
-                n_all += all; n_low += low; k_all += kept; k_low += kept_low;
-            });
-            if (k_all > 0 && n_low > 0) {
-                const double frac_all = static_cast<double>(k_all) / static_cast<double>(n_all);
-                const double frac_low = static_cast<double>(k_low) / static_cast<double>(n_low);
-                if (frac_low < 1.5 * frac_all) {
-                    for (auto &oc : outcomes)
-                        for (auto &r : oc.reflections)
-                            if (r.overload_kept) {
-                                r.overloaded = true;
-                                r.overload_kept = false;
-                            }
-                    result.overloads_kept = 0;
-                    logger.Info("Overloads: the {} kept as a lower bound are not concentrated at low resolution "
-                                "({:.3f}% of the reflections coarser than {:.0f} A, {:.3f}% of all) - not "
-                                "saturated spots, so they are dropped as overloads",
-                                k_all.load(), 100.0 * frac_low, OVERLOAD_LOW_RES_A, 100.0 * frac_all);
-                }
-            }
-        }
         std::atomic<int64_t> n_all{0}, n_low{0}, n_over{0}, n_over_low{0}, n_ice{0};
         ParallelFor(static_cast<int>(outcomes.size()), config_.nthreads, [&](int o) {
             int64_t all = 0, low = 0, over = 0, over_low = 0, on_ice = 0;

tests/BraggIntegrationEngineCPUTest.cpp

@@ -307,3 +307,47 @@ TEST_CASE("BraggIntegrationEngineCPU_StillKeepsOverloadAsLowerBound", "[Integrat
     CHECK(flagged[0].overloaded);
     CHECK(rotation.Counts().overloads_kept == 0);
 }
+
+// A still keeps an overloaded observation as a lower bound only where it has pixels at the saturation
+// limit. One overloaded by a pixel unreadable on this frame alone has nothing to count at the limit:
+// it stays flagged and is dropped.
+TEST_CASE("BraggIntegrationEngineCPU_StillKeepsOnlySaturatedOverloads", "[Integration][portable]") {
+    DiffractionExperiment experiment(DetJF(2));
+    experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f);
+    experiment.ImportBraggIntegrationSettings(BraggIntegrationSettings());
+    const int32_t limit = 500, bkg = 12;
+    experiment.Detector().SaturationLimit(limit);
+    const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum();
+
+    ImagePreprocessorBuffer image(npixel);
+    for (size_t i = 0; i < npixel; ++i)
+        image[i] = bkg;
+    auto spot = [&](int cx, int cy, double amplitude) {
+        for (int dy = -6; dy <= 6; ++dy)
+            for (int dx = -6; dx <= 6; ++dx) {
+                const int32_t v = bkg + static_cast<int32_t>(std::lround(amplitude * std::exp(-(dx * dx + dy * dy) / (2 * 1.3 * 1.3))));
+                image[(cy + dy) * width + cx + dx] = v >= limit ? INT32_MAX : v;
+            }
+    };
+    spot(600, 300, 2000.0);                  // clipped: a plateau of saturated pixels in a bright spot
+    spot(300, 300, 40.0);
+    image[300 * width + 300] = INT32_MIN;    // its peak unreadable on this frame, not in the mask
+
+    std::vector<Reflection> predicted;
+    for (const auto &[x, y] : {std::pair{600, 300}, std::pair{300, 300}}) {
+        Reflection r{};
+        r.h = x; r.k = y; r.l = 1;
+        r.predicted_x = x; r.predicted_y = y;
+        r.d = 2.0f;
+        r.prescaling_corr = 1.0f;
+        r.partiality = 1.0f;
+        predicted.push_back(r);
+    }
+
+    BraggIntegrationEngineCPU still(experiment, PixelMask(experiment));
+    const auto out = still.Run(image, predicted, predicted.size(), 0);
+    REQUIRE(out.size() == 2);
+    CHECK_FALSE(out[0].overloaded);
+    CHECK(out[1].overloaded);
+    CHECK(still.Counts().overloads_kept == 1);
+}

tests/BraggIntegrationEngineGPUTest.cpp

@@ -364,6 +364,58 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") {
 // next one, so a reused engine has to give the same answer as a fresh one. A first frame whose spots
 // are somewhere else entirely is what would show a leftover mark: a stale mask pixel is read as a
 // neighbour's signal and dropped from the background ring, a stale owner steals a pixel outright.
+// The still's overload decision on the device: a clipped spot is kept as a lower bound, an overload by
+// an unreadable pixel is not - as on the CPU (BraggIntegrationEngineCPU_StillKeepsOnlySaturatedOverloads).
+TEST_CASE("BraggIntegrationEngineGPU_StillOverloadsMatchCPU", "[Integration]") {
+    if (get_gpu_count() == 0) {
+        WARN("No CUDA GPU present. Skipping BraggIntegrationEngineGPU_StillOverloadsMatchCPU");
+        return;
+    }
+    DiffractionExperiment experiment(DetJF(2));
+    experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f);
+    experiment.ImportBraggIntegrationSettings(BraggIntegrationSettings());
+    const int32_t limit = 500, bkg = 12;
+    experiment.Detector().SaturationLimit(limit);
+    const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum();
+
+    std::vector<int32_t> image(npixel, bkg);
+    auto spot = [&](int cx, int cy, double amplitude) {
+        for (int dy = -6; dy <= 6; ++dy)
+            for (int dx = -6; dx <= 6; ++dx) {
+                const int32_t v = bkg + static_cast<int32_t>(std::lround(amplitude * std::exp(-(dx * dx + dy * dy) / (2 * 1.3 * 1.3))));
+                image[(cy + dy) * width + cx + dx] = v >= limit ? INT32_MAX : v;
+            }
+    };
+    spot(600, 300, 2000.0);
+    spot(300, 300, 40.0);
+    image[300 * width + 300] = INT32_MIN;
+    std::vector<Reflection> predicted;
+    for (const auto &[x, y] : {std::pair{600, 300}, std::pair{300, 300}})
+        predicted.push_back(MakeReflection(static_cast<float>(x), static_cast<float>(y), 2.0f, x));
+
+    ImagePreprocessorBuffer cpu_image(npixel);
+    for (size_t i = 0; i < npixel; ++i)
+        cpu_image[i] = image[i];
+    BraggIntegrationEngineCPU cpu(experiment, PixelMask(experiment));
+    const auto out_cpu = cpu.Run(cpu_image, predicted, predicted.size(), 0);
+
+    auto stream = std::make_shared<CudaStream>();
+    ImagePreprocessorBufferGPU gpu_image(npixel);
+    for (size_t i = 0; i < npixel; ++i)
+        gpu_image[i] = image[i];
+    REQUIRE(cudaMemcpyAsync(gpu_image.getGPUBuffer(), gpu_image.getBuffer().data(),
+                            npixel * sizeof(int32_t), cudaMemcpyHostToDevice, *stream) == cudaSuccess);
+    BraggIntegrationEngineGPU gpu(experiment, stream, PixelMask(experiment));
+    const auto out_gpu = gpu.Run(gpu_image, predicted, predicted.size(), 0);
+
+    REQUIRE(out_cpu.size() == 2);
+    REQUIRE(out_gpu.size() == 2);
+    CHECK_FALSE(out_gpu[0].overloaded);
+    CHECK(out_gpu[1].overloaded);
+    CHECK(gpu.Counts().overloads_kept == 1);
+    CHECK(out_gpu[0].I == Catch::Approx(out_cpu[0].I).epsilon(0.001));
+}
+
 TEST_CASE("BraggIntegrationEngineGPU_ReusedEngineMatchesFresh") {
     if (get_gpu_count() == 0) {
         WARN("No CUDA GPU present. Skipping BraggIntegrationEngineGPU_ReusedEngineMatchesFresh");

tests/HotPixelFinderTest.cpp

@@ -25,9 +25,10 @@ namespace {
 
 // A Poisson background with a powder ring, and on it four kinds of pixel: one that reads 60 counts
 // high on every frame, one only 20 high - persistent, but too weak to make an outlier - one that holds
-// the error value on every frame, and one crossed by a genuine reflection. The last sits close to the rotation axis, where one reflection stays on a pixel for a
+// the error value on every frame, one stuck at the saturation limit, and one crossed by a genuine
+// reflection, saturated on a few frames. The last sits close to the rotation axis, where one reflection stays on a pixel for a
 // long stretch of rotation - here nine consecutive sampled frames, more than the chance bound allows
-// but fewer than the one-reflection bound at that zeta. Only the first and the third are masked - and
+// but fewer than the one-reflection bound at that zeta. Only the first, third and fourth are masked - and
 // a 3x3 patch reading high on every frame, which is stationary in the lab and so no reflection either.
 // Three dark pixels besides: one dead on the powder ring, which is masked, one on the ring reading 70%
 // of it, and one dead on the 5-count background, whose deficit is too small to make an outlier.
@@ -40,7 +41,7 @@ TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") {
     HotPixelFinder finder(x, pixel_mask, 4);
 
     constexpr int HOT_X = 200, HOT_Y = 200, WARM_X = 190, WARM_Y = 60, ERR_X = 60, ERR_Y = 190;
-    constexpr int BRAGG_X = 40, BRAGG_Y = 115;
+    constexpr int BRAGG_X = 40, BRAGG_Y = 115, STUCK_X = 100, STUCK_Y = 30;
     constexpr int PATCH_X = 150, PATCH_Y = 225;
     constexpr int DEAD_X = C + 82, DEAD_Y = C, DIM_X = C, DIM_Y = C + 82, DEAD_BKG_X = 20, DEAD_BKG_Y = 20;
     std::mt19937 rng(1);
@@ -57,11 +58,12 @@ TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") {
         frame[I(HOT_X, HOT_Y)] += 60;
         frame[I(WARM_X, WARM_Y)] += 20;
         frame[I(ERR_X, ERR_Y)] = INT32_MIN;
+        frame[I(STUCK_X, STUCK_Y)] = INT32_MAX;
         for (int dy = -1; dy <= 1; dy++)
             for (int dx = -1; dx <= 1; dx++)
                 frame[I(PATCH_X + dx, PATCH_Y + dy)] += 60;
         if (f >= 20 && f < 29)
-            frame[I(BRAGG_X, BRAGG_Y)] += 500;
+            frame[I(BRAGG_X, BRAGG_Y)] = f < 23 ? INT32_MAX : frame[I(BRAGG_X, BRAGG_Y)] + 500;
         frame[I(DEAD_X, DEAD_Y)] = 0;
         frame[I(DIM_X, DIM_Y)] = static_cast<int32_t>(0.7 * frame[I(DIM_X, DIM_Y)]);
         frame[I(DEAD_BKG_X, DEAD_BKG_Y)] = 0;
@@ -72,6 +74,7 @@ TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") {
     CHECK(result.frames == NFRAMES);
     CHECK(result.mask[I(HOT_X, HOT_Y)] == 1);   // 1 = hot, 2 = error value, 3 = dead
     CHECK(result.mask[I(ERR_X, ERR_Y)] == 2);
+    CHECK(result.mask[I(STUCK_X, STUCK_Y)] == 2);
     CHECK(result.mask[I(WARM_X, WARM_Y)] == 0);
     CHECK(result.mask[I(BRAGG_X, BRAGG_Y)] == 0);
     for (int dy = -1; dy <= 1; dy++)
@@ -81,7 +84,7 @@ TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") {
     CHECK(result.mask[I(DIM_X, DIM_Y)] == 0);
     CHECK(result.mask[I(DEAD_BKG_X, DEAD_BKG_Y)] == 0);
     CHECK(result.hot == 10);
-    CHECK(result.error == 1);
+    CHECK(result.error == 2);
     CHECK(result.dead == 1);
 }