Pre-scan: the hot-pixel rule applies on stills too

A pixel lit on more frames than one reflection or chance explains is a defect on stills as on rotation:
a still's reflection is on a pixel for one frame, so two frames are already more than one explains, and
the chance bound measured on the ring covers a lattice that recurs. The rule was switched off without a
rotation (spacing 0), so a live or offline stills run whose input mask was incomplete analysed every
frame with its hot pixels: on a JUNGFRAU 8M stills set with no detector mask ~350 of them made ~350
fixed false spots on every frame, junk lattices on the empty frames and a merge far below one with the
mask. The error-value, dead and unstable rules are unchanged; the GPU candidate pass no longer takes
the spacing flag it only tested.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK

image_analysis/hot_pixels/HotPixels.cpp

@@ -367,7 +367,7 @@ HotPixelFinder::Result HotPixelFinder::GetMask(double oscillation_deg, double sp
     // tests below could pass (see GetCandidates) - into the host arrays, which are zero everywhere else,
     // so the tests below run on them unchanged.
     if (gpu) {
-        const auto c = gpu->GetCandidates(frames, min_valid, spacing_deg > 0.0, k_chance);
+        const auto c = gpu->GetCandidates(frames, min_valid, k_chance);
         std::fill(sums_known.begin(), sums_known.end(), 0);
         for (size_t j = 0; j < c.index.size(); j++) {
             const size_t i = c.index[j];
@@ -389,13 +389,17 @@ HotPixelFinder::Result HotPixelFinder::GetMask(double oscillation_deg, double sp
     std::vector<uint8_t> persistent(width * height, 0);
     ParallelChunks(static_cast<int>(height), nthreads, [&](int y0, int y1) {
         ForEachPixelInRows(modules, static_cast<int64_t>(width), y0, y1, [&](int64_t, size_t i, double px, double py) {
-            if (key[i] < 0 || n_valid(i) < min_valid || n_lit[i] == 0 || !(spacing_deg > 0.0)) return;
+            if (key[i] < 0 || n_valid(i) < min_valid || n_lit[i] == 0) return;
             // One reflection's stay on this pixel, in sampled frames: 1 / |zeta| turns the rocking
-            // width into rotation, and zeta = |axis . (s1 x s0)| with s0 along z.
-            const Coord s1 = geometry.LabCoord(static_cast<float>(px), static_cast<float>(py)).Normalize();
-            const double zeta = std::max(1e-4, static_cast<double>(std::fabs(axis * (s1 % Coord(0, 0, 1)))));
-            const double stay = (oscillation_deg + ROCKING_WIDTH_DEG) / zeta;
-            const double k_bragg = 1.0 + std::ceil(stay / spacing_deg);
+            // width into rotation, and zeta = |axis . (s1 x s0)| with s0 along z. On stills (no
+            // spacing) a reflection is on a pixel for one frame, so two frames are more than one explains.
+            double k_bragg = 2.0;
+            if (spacing_deg > 0.0) {
+                const Coord s1 = geometry.LabCoord(static_cast<float>(px), static_cast<float>(py)).Normalize();
+                const double zeta = std::max(1e-4, static_cast<double>(std::fabs(axis * (s1 % Coord(0, 0, 1)))));
+                const double stay = (oscillation_deg + ROCKING_WIDTH_DEG) / zeta;
+                k_bragg = 1.0 + std::ceil(stay / spacing_deg);
+            }
             const double k = std::min<double>(std::max<double>(k_bragg, k_chance[key[i] / SECTORS]),
                                               n_valid(i));
             persistent[i] = n_lit[i] >= k;

image_analysis/hot_pixels/HotPixels.h

@@ -22,7 +22,9 @@
 // either explanation that is not a defect allows:
 //
 //   * one Bragg reflection: it stays on a pixel for (oscillation + rocking width) / |zeta| of
-//     rotation, i.e. on at most 1 + ceil(that / frame spacing) of frames sampled that far apart;
+//     rotation, i.e. on at most 1 + ceil(that / frame spacing) of frames sampled that far apart; on
+//     stills (no spacing) a reflection is on a pixel for one frame, and a lattice at the same
+//     orientation on frame after frame is what the chance rate below measures;
 //   * chance: independent reflections and noise light a pixel of this ring on a fraction q of the
 //     frames (measured, on the ring's pixels that are not candidates), and the number of lit frames
 //     that fewer than 0.01 of all the detector's pixels would reach by chance is the binomial bound.

image_analysis/hot_pixels/HotPixelsGPU.cu

@@ -208,7 +208,7 @@ __global__ void chance_kernel(size_t npixels, int sectors, const int32_t *__rest
 }
 
 // Writes the candidates' indices from `out` on when `out` is given, and counts them either way.
-__global__ void candidate_kernel(size_t npixels, int sectors, HotPixelLevelRules rules, uint32_t frames, int min_valid, bool spacing_ok,
+__global__ void candidate_kernel(size_t npixels, int sectors, HotPixelLevelRules rules, uint32_t frames, int min_valid,
                                  const int32_t *__restrict__ key, const uint32_t *__restrict__ key_frames,
                                  const uint16_t *__restrict__ n_lit, const uint16_t *__restrict__ n_error,
                                  const uint16_t *__restrict__ n_error_ring_ok, const int64_t *__restrict__ sum_value,
@@ -219,7 +219,7 @@ __global__ void candidate_kernel(size_t npixels, int sectors, HotPixelLevelRules
     const int nv = valid_frames(key, key_frames, n_error_ring_ok, i);
     const bool error = 2 * static_cast<uint32_t>(n_error[i]) > frames;
     const int lit = n_lit[i];
-    const bool persistent = spacing_ok && nv >= min_valid && lit > 0
+    const bool persistent = nv >= min_valid && lit > 0
                             && lit >= min(max(2, k_chance[key[i] / sectors]), nv);
     // Can be dead (HotPixelFinder::GetMask): at most half of its sector's level, no error frames.
     const bool dark = n_error[i] == 0 && nv >= min_valid && 2 * sum_value[i] <= key_local_level_sum[key[i]];
@@ -385,7 +385,7 @@ void HotPixelFinderGPU::ChanceCounts(std::vector<int64_t> &lit, std::vector<int6
     cuda_err(cudaStreamSynchronize(stream));
 }
 
-HotPixelFinderGPU::Candidates HotPixelFinderGPU::GetCandidates(uint32_t frames, int min_valid, bool spacing_ok,
+HotPixelFinderGPU::Candidates HotPixelFinderGPU::GetCandidates(uint32_t frames, int min_valid,
                                                                const std::vector<int> &k_chance) {
     CudaStream stream;
     cuda_err(cudaStreamWaitEvent(stream, last_accumulate, 0));
@@ -396,7 +396,7 @@ HotPixelFinderGPU::Candidates HotPixelFinderGPU::GetCandidates(uint32_t frames,
         cuda_err(cudaMemcpyAsync(d_k_chance, k_chance.data(), k_chance.size() * sizeof(int), cudaMemcpyHostToDevice,
                                  stream));
     cuda_err(cudaMemsetAsync(d_count, 0, sizeof(uint32_t), stream));
-    candidate_kernel<<<blocks, THREADS, 0, stream>>>(npixels, sectors, rules, frames, min_valid, spacing_ok, key, key_frames,
+    candidate_kernel<<<blocks, THREADS, 0, stream>>>(npixels, sectors, rules, frames, min_valid, key, key_frames,
                                                      n_lit, n_error, n_error_ring_ok, sum_value, key_local_level_sum,
                                                      d_k_chance, n_neg, d_count, nullptr);
     cuda_err(cudaGetLastError());
@@ -421,7 +421,7 @@ HotPixelFinderGPU::Candidates HotPixelFinderGPU::GetCandidates(uint32_t frames,
         CudaDevicePtr<uint16_t> g_lit(n, ALLOC), g_error(n, ALLOC), g_error_ring_ok(n, ALLOC), g_neg(n, ALLOC);
         CudaDevicePtr<int64_t> g_sum(n, ALLOC), g_error_level(n, ALLOC);
         cuda_err(cudaMemsetAsync(d_count, 0, sizeof(uint32_t), stream));
-        candidate_kernel<<<blocks, THREADS, 0, stream>>>(npixels, sectors, rules, frames, min_valid, spacing_ok, key,
+        candidate_kernel<<<blocks, THREADS, 0, stream>>>(npixels, sectors, rules, frames, min_valid, key,
                                                          key_frames, n_lit, n_error, n_error_ring_ok, sum_value,
                                                          key_local_level_sum, d_k_chance, n_neg, d_count, index);
         cuda_err(cudaGetLastError());

image_analysis/hot_pixels/HotPixelsGPU.h

@@ -96,9 +96,9 @@ public:
     void ChanceCounts(std::vector<int64_t> &lit, std::vector<int64_t> &seen);
 
     // The pixels that can be masked: those holding the error value on more than half of `frames`;
-    // where spacing_ok, those lit on at least min(max(2, k_chance[ring]), valid frames) of at least
-    // min_valid valid frames - a persistent pixel is lit on at least that many, because one reflection
-    // explains at least two; and those with no error frame, at least min_valid valid frames and a sum
+    // those lit on at least min(max(2, k_chance[ring]), valid frames) of at least min_valid valid
+    // frames - a persistent pixel is lit on at least that many, because one reflection explains at
+    // least two; and those with no error frame, at least min_valid valid frames and a sum
     // at most half of their sector's level sum, which a dead pixel has. For each, its index and
     // per-pixel sums; and every key's sums.
     struct Candidates {
@@ -108,5 +108,5 @@ public:
         std::vector<uint32_t> key_frames;
         std::vector<int64_t> key_level_sum, key_local_level_sum;
     };
-    Candidates GetCandidates(uint32_t frames, int min_valid, bool spacing_ok, const std::vector<int> &k_chance);
+    Candidates GetCandidates(uint32_t frames, int min_valid, const std::vector<int> &k_chance);
 };

rugnux/RugnuxPreScan.cpp

@@ -240,9 +240,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
     // radius the first one measured rather than re-reading the frames.
     const bool want_width = config_.adaptive_integration_radius && !spot_width_measured_;
     // Defective pixels, for the same reason on the same pass: each pixel is compared against its
-    // resolution ring. Only under rotation, where a reflection moves off a pixel and a defect does not -
-    // except for the pixels holding the detector's error value and the dead ones, which need no
-    // rotation and are measured on stills too (MaskDefectivePixels).
+    // resolution ring, and one lit on more frames than a reflection or chance explains is masked
+    // (MaskDefectivePixels) - on stills as on rotation: a still's reflection is on a pixel for one frame.
     const auto goniometer = experiment_.GetGoniometer();
     const bool want_hot = !hot_pixels_measured_
                           && (!experiment_.IsRotationIndexing() || (goniometer && goniometer->IsScanning()));
@@ -1414,8 +1413,8 @@ void Rugnux::MaskDefectivePixels(int start_image, const std::vector<int> &sample
         }
     });
 
-    // On stills there is no rotation to count a reflection's stay on a pixel in, so spacing stays 0 and
-    // the hot-pixel rule is not applied (HotPixelFinder::GetMask); the error-value and dead pixels are.
+    // On stills there is no rotation to count a reflection's stay on a pixel in: spacing stays 0 and a
+    // pixel lit on two frames is already more than one reflection explains (HotPixelFinder::GetMask).
     const auto goniometer = experiment_.GetGoniometer();
     const double oscillation = goniometer ? std::abs(goniometer->GetIncrement_deg()) : 0.0;
     const double spacing = goniometer && sample.size() > 1

tests/HotPixelFinderTest.cpp

@@ -93,6 +93,57 @@ TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") {
     CHECK(result.unstable == 1);
 }
 
+// Stills: no rotation, so a reflection is on a pixel for one frame and nothing but a defect lights the
+// same pixel on frame after frame. The same background and ring; a pixel 60 high on every frame and a 3x3
+// patch like it are masked, a pixel 20 high on every frame is not (too weak for an outlier), and a pixel
+// that a strong reflection lands on on four of the sixty frames - a lattice that recurs - is not: four
+// frames are within what chance allows a ring lit as often as this one. The error-value and dead rules as
+// on rotation.
+TEST_CASE("HotPixelFinder_StillsPersistentPixel", "[HotPixelFinder]") {
+    DiffractionExperiment x(DetDECTRIS(W, H, "Test detector", ""));
+    x.IncidentEnergy_keV(WVL_1A_IN_KEV).DetectorDistance_mm(10.0f);
+    x.BeamX_pxl(static_cast<float>(C)).BeamY_pxl(static_cast<float>(C));
+    const PixelMask pixel_mask(x);
+    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, PATCH_X = 150, PATCH_Y = 225, DEAD_X = C + 82, DEAD_Y = C;
+    std::mt19937 rng(2);
+    std::vector<int32_t> frame(static_cast<size_t>(W) * H), scratch;
+    for (int f = 0; f < NFRAMES; f++) {
+        for (int y = 0; y < H; y++)
+            for (int x_ = 0; x_ < W; x_++) {
+                const double r = std::hypot(x_ - C, y - C);
+                const double mean = 5.0 + 45.0 * std::exp(-0.5 * (r - 82.0) * (r - 82.0) / 16.0);
+                frame[I(x_, y)] = std::poisson_distribution<int>(mean)(rng);
+            }
+        frame[I(HOT_X, HOT_Y)] += 60;
+        frame[I(WARM_X, WARM_Y)] += 20;
+        frame[I(ERR_X, ERR_Y)] = INT32_MIN;
+        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 % 15 == 7)
+            frame[I(BRAGG_X, BRAGG_Y)] += 500;
+        frame[I(DEAD_X, DEAD_Y)] = 0;
+        finder.AddImage(frame.data(), scratch);
+    }
+
+    const auto result = finder.GetMask(0.0, 0.0, 4);
+    CHECK(result.frames == NFRAMES);
+    CHECK(result.mask[I(HOT_X, HOT_Y)] == 1);
+    for (int dy = -1; dy <= 1; dy++)
+        for (int dx = -1; dx <= 1; dx++)
+            CHECK(result.mask[I(PATCH_X + dx, PATCH_Y + dy)] == 1);
+    CHECK(result.mask[I(WARM_X, WARM_Y)] == 0);
+    CHECK(result.mask[I(BRAGG_X, BRAGG_Y)] == 0);
+    CHECK(result.mask[I(ERR_X, ERR_Y)] == 2);
+    CHECK(result.mask[I(DEAD_X, DEAD_Y)] == 3);
+    CHECK(result.hot == 10);
+    CHECK(result.error == 1);
+    CHECK(result.dead == 1);
+}
+
 #ifdef JFJOCH_USE_CUDA
 // The device path against the host one, on frames that take every branch of both: a background bright
 // enough on one annulus that the host reads its rings' statistics off the exact selection rather than
@@ -227,8 +278,8 @@ TEST_CASE("HotPixelFinderGPU_FoldedCardsMatchOneCard", "[HotPixelFinder]") {
     CHECK(lit_folded == lit_one);
     CHECK(seen_folded == seen_one);
     const std::vector<int> k_chance(NRINGS, 5);
-    const auto expected = one.GetCandidates(NFRAMES, 10, true, k_chance);
-    const auto result = first.GetCandidates(NFRAMES, 10, true, k_chance);
+    const auto expected = one.GetCandidates(NFRAMES, 10, k_chance);
+    const auto result = first.GetCandidates(NFRAMES, 10, k_chance);
     CHECK(expected.index.size() >= 40);
     // The candidates are written by atomic slot, so their order is the device's; compare them by pixel.
     const auto by_pixel = [](const HotPixelFinderGPU::Candidates &c) {