From 8c490db59c611f356a077dff65dd17fd64b91aa5 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 13 Nov 2026 00:55:31 +0100 Subject: [PATCH] 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 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK --- image_analysis/hot_pixels/HotPixels.cpp | 18 +++++--- image_analysis/hot_pixels/HotPixels.h | 4 +- image_analysis/hot_pixels/HotPixelsGPU.cu | 10 ++--- image_analysis/hot_pixels/HotPixelsGPU.h | 8 ++-- rugnux/RugnuxPreScan.cpp | 9 ++-- tests/HotPixelFinderTest.cpp | 55 ++++++++++++++++++++++- 6 files changed, 80 insertions(+), 24 deletions(-) diff --git a/image_analysis/hot_pixels/HotPixels.cpp b/image_analysis/hot_pixels/HotPixels.cpp index 2f56213f1..998056e6c 100644 --- a/image_analysis/hot_pixels/HotPixels.cpp +++ b/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 persistent(width * height, 0); ParallelChunks(static_cast(height), nthreads, [&](int y0, int y1) { ForEachPixelInRows(modules, static_cast(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(px), static_cast(py)).Normalize(); - const double zeta = std::max(1e-4, static_cast(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(px), static_cast(py)).Normalize(); + const double zeta = std::max(1e-4, static_cast(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(std::max(k_bragg, k_chance[key[i] / SECTORS]), n_valid(i)); persistent[i] = n_lit[i] >= k; diff --git a/image_analysis/hot_pixels/HotPixels.h b/image_analysis/hot_pixels/HotPixels.h index 270a9c911..4f0265731 100644 --- a/image_analysis/hot_pixels/HotPixels.h +++ b/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. diff --git a/image_analysis/hot_pixels/HotPixelsGPU.cu b/image_analysis/hot_pixels/HotPixelsGPU.cu index 672196f10..5d36f5cdc 100644 --- a/image_analysis/hot_pixels/HotPixelsGPU.cu +++ b/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(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 &lit, std::vector &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<<>>(npixels, sectors, rules, frames, min_valid, spacing_ok, key, key_frames, + candidate_kernel<<>>(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 g_lit(n, ALLOC), g_error(n, ALLOC), g_error_ring_ok(n, ALLOC), g_neg(n, ALLOC); CudaDevicePtr g_sum(n, ALLOC), g_error_level(n, ALLOC); cuda_err(cudaMemsetAsync(d_count, 0, sizeof(uint32_t), stream)); - candidate_kernel<<>>(npixels, sectors, rules, frames, min_valid, spacing_ok, key, + candidate_kernel<<>>(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()); diff --git a/image_analysis/hot_pixels/HotPixelsGPU.h b/image_analysis/hot_pixels/HotPixelsGPU.h index 15d0d2e4b..76b450b1d 100644 --- a/image_analysis/hot_pixels/HotPixelsGPU.h +++ b/image_analysis/hot_pixels/HotPixelsGPU.h @@ -96,9 +96,9 @@ public: void ChanceCounts(std::vector &lit, std::vector &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 key_frames; std::vector key_level_sum, key_local_level_sum; }; - Candidates GetCandidates(uint32_t frames, int min_valid, bool spacing_ok, const std::vector &k_chance); + Candidates GetCandidates(uint32_t frames, int min_valid, const std::vector &k_chance); }; diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 00ce896dd..1cd0d3031 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/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 &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 diff --git a/tests/HotPixelFinderTest.cpp b/tests/HotPixelFinderTest.cpp index 44e64a065..1d30dea29 100644 --- a/tests/HotPixelFinderTest.cpp +++ b/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(C)).BeamY_pxl(static_cast(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 frame(static_cast(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(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 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) {