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/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/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"); 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); }