From b7e4ea23a3e8d5e81eff54772ed0c113491c80e1 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 22 Oct 2026 17:27:43 +0200 Subject: [PATCH 1/2] Pre-scan masks a pixel at the saturation limit on most frames, as one holding the error value A pixel stuck at the saturation limit was never masked: the defective-pixel scan counted only the error value as "invalid on most frames", and a saturated pixel merely as lit, which the hot rule does not judge on stills. A Rayonix stills set has 3,627 pixels at 65535 on every frame, at high resolution; each flagged every reflection whose disk it fell in as overloaded, and the spot finder found them on every frame. The saturation limit now counts with the error value (CPU and CUDA), so such a pixel is masked like one holding the error value. On that set the indexed fraction rises from 10% to 33% (and 7% to 30% on the second lattice of the same run), R_split 0.120 -> 0.106 and the unrefined R-free against the deposition -0.001 / -0.024. Rotation merges checked bit-identical (three in-house sets with overloads). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- image_analysis/hot_pixels/HotPixels.cpp | 5 +++-- image_analysis/hot_pixels/HotPixels.h | 8 +++++--- image_analysis/hot_pixels/HotPixelsGPU.cu | 4 ++-- rugnux/RugnuxPreScan.cpp | 2 +- tests/HotPixelFinderTest.cpp | 13 ++++++++----- 5 files changed, 19 insertions(+), 13 deletions(-) diff --git a/image_analysis/hot_pixels/HotPixels.cpp b/image_analysis/hot_pixels/HotPixels.cpp index 29937caec..79736c80b 100644 --- a/image_analysis/hot_pixels/HotPixels.cpp +++ b/image_analysis/hot_pixels/HotPixels.cpp @@ -203,7 +203,8 @@ void HotPixelFinder::AddImage(const int32_t *image, std::vector &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 &scratc } if (!ring_ok[k / SECTORS]) continue; sum_value[i] += v; - if (v == INT32_MAX || static_cast(v) > threshold[k]) + if (static_cast(v) > threshold[k]) n_lit[i]++; } } diff --git a/image_analysis/hot_pixels/HotPixels.h b/image_analysis/hot_pixels/HotPixels.h index e364871e0..7320979bc 100644 --- a/image_analysis/hot_pixels/HotPixels.h +++ b/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 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; }; diff --git a/image_analysis/hot_pixels/HotPixelsGPU.cu b/image_analysis/hot_pixels/HotPixelsGPU.cu index a5725533b..61dbc3474 100644 --- a/image_analysis/hot_pixels/HotPixelsGPU.cu +++ b/image_analysis/hot_pixels/HotPixelsGPU.cu @@ -132,7 +132,7 @@ __global__ void accumulate_kernel(const int32_t *__restrict__ image, const int32 i += static_cast(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(v) > threshold[k]) + if (static_cast(v) > threshold[k]) n_lit[i]++; } } diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 9603581df..6786be5be 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -1218,7 +1218,7 @@ void Rugnux::MaskDefectivePixels(int start_image, const std::vector &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")); }; diff --git a/tests/HotPixelFinderTest.cpp b/tests/HotPixelFinderTest.cpp index a4015e8b8..da00f00cd 100644 --- a/tests/HotPixelFinderTest.cpp +++ b/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(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); } From 5905d34cba276788f810b2cfaa02e465f7446e58 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 22 Oct 2026 17:27:44 +0200 Subject: [PATCH 2/2] Stills keep an overload as a lower bound only where it has pixels at the limit; run-level gate removed The lower bound counts each saturated pixel at the saturation limit, so it needs one. An observation overloaded only by pixels unreadable on this frame (a per-frame bad-pixel flag, data value 0) was kept as a box sum with those pixels counted as zero, which bounds nothing; it now stays overloaded and is dropped, per observation. With the stuck pixels masked by the pre-scan, the run-level gate (kept overloads must be enriched coarser than 4 A) has nothing left to catch: without it every still set gives the gated result (clipped CCD and Rayonix sets identical to the fourth decimal, photon-counting sets unchanged), so it is removed with Reflection::overload_kept. A neighbour-brightness test per saturated pixel was tried and not kept: after the mask it only dropped 9-20% of the Rayonix overloads and worsened R_split on one set (0.245 -> 0.278). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- common/Reflection.h | 3 -- .../BraggIntegrationEngine.cpp | 7 ++- rugnux/RugnuxScaleMerge.cpp | 39 -------------- tests/BraggIntegrationEngineCPUTest.cpp | 44 ++++++++++++++++ tests/BraggIntegrationEngineGPUTest.cpp | 52 +++++++++++++++++++ 5 files changed, 101 insertions(+), 44 deletions(-) diff --git a/common/Reflection.h b/common/Reflection.h index 1ab59fc19..1b8b5caec 100644 --- a/common/Reflection.h +++ b/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 diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index 8ed7b99a6..25769d19f 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -284,11 +284,14 @@ std::vector BraggIntegrationEngine::Finalize(const std::vector 0) { refl.I = fr.I + static_cast(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; diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index b3e876321..815b0edec 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/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 k_all{0}, k_low{0}, n_all{0}, n_low{0}; - ParallelFor(static_cast(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(k_all) / static_cast(n_all); - const double frac_low = static_cast(k_low) / static_cast(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 n_all{0}, n_low{0}, n_over{0}, n_over_low{0}, n_ice{0}; ParallelFor(static_cast(outcomes.size()), config_.nthreads, [&](int o) { int64_t all = 0, low = 0, over = 0, over_low = 0, on_ice = 0; diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp index 1d1068987..e4789cc3e 100644 --- a/tests/BraggIntegrationEngineCPUTest.cpp +++ b/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(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 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); +} diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index 5fcf3d48e..233449cc5 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/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 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(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 predicted; + for (const auto &[x, y] : {std::pair{600, 300}, std::pair{300, 300}}) + predicted.push_back(MakeReflection(static_cast(x), static_cast(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(); + 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");