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