From a609f6ba1cb075526ee478f65d3b7fe3cc9befc2 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 11 Oct 2026 09:54:47 +0200 Subject: [PATCH] rugnux two-pass: the integration radius decided on the pre-pass frames, before the canonical pass Where the measured spot width widened the integration radius, the canonical pass integrated the whole sweep at the widened radius, and only then did the starvation guard read how many background rings the neighbouring predictions crowd; over the bound (1.13%) the pass was thrown away and run again at the fixed radius. The geometry pre-pass now counts that same quantity at the widened radii without integrating there (BraggIntegrationEngineCPU::CountNeighbourCrowding: the neighbour mask and the ring pixels the run's mask leaves readable, no image) on 20 frames spread over the pre-pass, with the standard error of the ratio over those frames. Over the bound, the canonical pass integrates at the fixed radius from the outset. A count within NEAR_TIE_SIGMA of the bound is a near tie in the part-B sense: the canonical pass counts the widened radii again on 20 frames over the whole sweep, and where that puts them under the bound the run is made again deciding late. A run deciding late counts nothing early and decides as before. Under the bound the canonical pass integrates at the widened radius and the measured guard stays as it was, so a pre-pass count that is too low (its prediction holds fewer tails than the canonical pass's, or it stood on another lattice class) costs only the pass it cost before. No new threshold. Battery against 46be3f647 (the 14 open sets that re-ran at the fixed radius, the 9 open/in-house sets closest under the bound, the ci tier, 2 private sets): - p.hkl and MTZ data identical on all 51 open/in-house sets and on one private set. - The 14 re-running sets: rugnux wall 534 -> 447 s (7qij -27, 4nwv -12, 7yzx -10, 9fhc -9, 5src -6 s, on a shared machine). 8sqq and 8xtf count under the bound on the pre-pass and re-run as before. 4nwv (1.69 +- 0.28%) and one other set were near ties, re-checked on the whole sweep (1.72%, 1.61%), ending where they did. - One private, densely crowded set: same space group and d_min, R_meas 28.6 -> 39.3%. Its old re-run re-indexed starting from the spot budget the thrown-away pass had measured (152 of 1000). With that budget given (--max-spots 152), the new run gives 30.6%. The difference comes from that path, not from the radius. - The count costs ~13 CPU-s on the heaviest set (65k predictions a frame). Over the 22 ci sets that count, pre-pass integration took 45.2 -> 45.8 s. - A build forcing the near tie and its overturn on 4nwv re-checks on the whole sweep (1.80% against the engine's 1.81%) and makes the run again deciding late. It ends with the same p.hkl. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- image_analysis/MXAnalysisWithoutFPGA.cpp | 12 +++ image_analysis/MXAnalysisWithoutFPGA.h | 12 +++ .../BraggIntegrationEngineCPU.cpp | 86 +++++++++++++------ .../BraggIntegrationEngineCPU.h | 9 ++ rugnux/Rugnux.cpp | 2 + rugnux/Rugnux.h | 18 ++++ rugnux/RugnuxImagePass.cpp | 51 +++++++++++ rugnux/RugnuxPasses.cpp | 41 ++++++++- tests/BraggIntegrationEngineCPUTest.cpp | 35 ++++++++ 9 files changed, 240 insertions(+), 26 deletions(-) diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 566782e62..9f6354a7c 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -188,6 +188,8 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, : FixedThresholdFinder(); const auto integrate_fn = [this](const std::vector &predicted, size_t npredicted, int64_t image_number) { + if (count_crowding) + crowding_engine->CountNeighbourCrowding(predicted, npredicted); return bragg_engine->Run(*preprocessor_buffer, predicted, npredicted, image_number); }; // The supercell probe integrates on the same engine, but what the engine counts is the run's @@ -446,6 +448,16 @@ BraggIntegrationCounts MXAnalysisWithoutFPGA::BraggCounts() const { return c; } +void MXAnalysisWithoutFPGA::CountCrowdingAt(const BraggIntegrationSettings &settings) { + DiffractionExperiment x = experiment; + x.ImportBraggIntegrationSettings(settings); + crowding_engine = std::make_unique(x, mask); +} + +BraggIntegrationCounts MXAnalysisWithoutFPGA::CrowdingCounts() const { + return crowding_engine ? crowding_engine->Counts() : BraggIntegrationCounts{}; +} + void MXAnalysisWithoutFPGA::UpdateMaskResolution(const SpotFindingSettings &settings) { mask_low_res = settings.low_resolution_limit; mask_high_res = settings.high_resolution_limit; diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index 670a8b1e5..b12d1c38d 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -12,6 +12,7 @@ #include "../common/AzimuthalIntegrationProfile.h" #include "bragg_prediction/BraggPrediction.h" #include "bragg_integration/BraggIntegrationEngine.h" +#include "bragg_integration/BraggIntegrationEngineCPU.h" #include "spot_finding/ImageSpotFinder.h" #include "spot_finding/AdaptiveSpotFinderCPU.h" #include "indexing/IndexerThreadPool.h" @@ -64,6 +65,10 @@ class MXAnalysisWithoutFPGA { std::unique_ptr bragg_engine; // What the supercell probe's integrations added to bragg_engine's counts (see BraggCounts). BraggIntegrationCounts probe_counts; + // Counts the neighbour crowding of the integrated reflections at another stencil (CountCrowdingAt), + // on the images for which count_crowding is on. + std::unique_ptr crowding_engine; + bool count_crowding = false; std::unique_ptr preprocessor_buffer; const PixelMask &mask; @@ -111,6 +116,13 @@ public: // What this worker's Bragg integrator counted (BraggIntegrationCounts). Each worker builds its own // analysis, so a caller that wants the run's totals sums this over the workers it started. [[nodiscard]] BraggIntegrationCounts BraggCounts() const; + + // Also count, on the images Analyze is given while CountCrowding(true), how crowded the background + // rings of the integrated reflections would be at `settings` rather than the experiment's own + // (BraggIntegrationEngineCPU::CountNeighbourCrowding), read back with CrowdingCounts. + void CountCrowdingAt(const BraggIntegrationSettings &settings); + void CountCrowding(bool on) { count_crowding = on; } + [[nodiscard]] BraggIntegrationCounts CrowdingCounts() const; }; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 3bb37b74d..146e850e6 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -69,6 +69,66 @@ BraggIntegrationEngineCPU::BraggIntegrationEngineCPU(const DiffractionExperiment refl_mask(static_cast(xpixel), static_cast(ypixel), 0), owner(static_cast(xpixel), static_cast(ypixel), BRAGG_OWNER_NONE) {} +void BraggIntegrationEngineCPU::MarkNeighbours(const std::vector &predicted, size_t npredicted) { + const int W = static_cast(xpixel), H = static_cast(ypixel); + // --- Reflection mask: mark the r2 signal region of every predicted reflection, so a neighbour's + // signal is excluded from this reflection's r2..r3 background ring. MASK_FLUX where the + // neighbour puts flux on this frame (NEIGHBOUR_MASK_MIN_PARTIALITY) - only those are excluded + // - and MASK_TAIL where only the tail of a rocking curve does, which is counted for the density + // guard but read as background. The region is the INNER stencil ellipse, taken in each + // neighbour's OWN frame - an elongated neighbour whose streak is still marked as a disk would + // leak its tails into this reflection's ring. --- + refl_mask.Clear(); + for (size_t i = 0; i < npredicted; ++i) { + const auto &r = predicted[i]; + const uint8_t level = r.partiality >= NEIGHBOUR_MASK_MIN_PARTIALITY ? MASK_FLUX : MASK_TAIL; + const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil); + const int x0 = std::max(0, static_cast(std::floor(r.predicted_x - st.ex_in - 1.0f))); + const int x1 = std::min(W - 1, static_cast(std::ceil(r.predicted_x + st.ex_in + 1.0f))); + const int y0 = std::max(0, static_cast(std::floor(r.predicted_y - st.ey_in - 1.0f))); + const int y1 = std::min(H - 1, static_cast(std::ceil(r.predicted_y + st.ey_in + 1.0f))); + for (int y = y0; y <= y1; ++y) + for (int x = x0; x <= x1; ++x) { + const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y); + if (d.inner < r2_sq) { + uint8_t &m = refl_mask.At(x, y); + m = std::max(m, level); + } + } + } +} + +void BraggIntegrationEngineCPU::CountNeighbourCrowding(const std::vector &predicted, + size_t npredicted) { + const int W = static_cast(xpixel), H = static_cast(ypixel); + counts.predicted += npredicted; + MarkNeighbours(predicted, npredicted); + for (size_t i = 0; i < npredicted; ++i) { + const auto &r = predicted[i]; + const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil); + const int x0 = std::max(0, static_cast(std::floor(r.predicted_x - st.ex_out - 1.0f))); + const int x1 = std::min(W - 1, static_cast(std::ceil(r.predicted_x + st.ex_out + 1.0f))); + const int y0 = std::max(0, static_cast(std::floor(r.predicted_y - st.ey_out - 1.0f))); + const int y1 = std::min(H - 1, static_cast(std::ceil(r.predicted_y + st.ey_out + 1.0f))); + int n_ring = 0, n_clean = 0; + for (int y = y0; y <= y1; ++y) + for (int x = x0; x <= x1; ++x) { + const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y); + if (BraggStencilInSignal(st, x - r.predicted_x, y - r.predicted_y, d.signal, r1_sq) + || d.inner < r2_sq || d.outer >= r3_sq) + continue; + const size_t idx = static_cast(y) * W + x; + if (idx / 32 < static_mask.size() && ((static_mask[idx / 32] >> (idx % 32)) & 1U)) + continue; + ++n_ring; + if (refl_mask.Get(x, y) == 0) + ++n_clean; + } + if (n_clean <= 5 && n_ring > 5) + ++counts.bkg_starved_by_neighbour; + } +} + template std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const std::vector &predicted, @@ -90,32 +150,8 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, auto grid_idx = [this](int dx, int dy) { return (dy + R) * G + (dx + R); }; - // --- Reflection mask: mark the r2 signal region of every predicted reflection, so a neighbour's - // signal is excluded from this reflection's r2..r3 background ring. MASK_FLUX where the - // neighbour puts flux on this frame (NEIGHBOUR_MASK_MIN_PARTIALITY) - only those are excluded - // - and MASK_TAIL where only the tail of a rocking curve does, which is counted for the density - // guard but read as background. The region is the INNER stencil ellipse, taken in each - // neighbour's OWN frame - an elongated neighbour whose streak is still marked as a disk would - // leak its tails into this reflection's ring. --- - refl_mask.Clear(); + MarkNeighbours(predicted, npredicted); owner.Clear(); - for (size_t i = 0; i < npredicted; ++i) { - const auto &r = predicted[i]; - const uint8_t level = r.partiality >= NEIGHBOUR_MASK_MIN_PARTIALITY ? MASK_FLUX : MASK_TAIL; - const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil); - const int x0 = std::max(0, static_cast(std::floor(r.predicted_x - st.ex_in - 1.0f))); - const int x1 = std::min(W - 1, static_cast(std::ceil(r.predicted_x + st.ex_in + 1.0f))); - const int y0 = std::max(0, static_cast(std::floor(r.predicted_y - st.ey_in - 1.0f))); - const int y1 = std::min(H - 1, static_cast(std::ceil(r.predicted_y + st.ey_in + 1.0f))); - for (int y = y0; y <= y1; ++y) - for (int x = x0; x <= x1; ++x) { - const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y); - if (d.inner < r2_sq) { - uint8_t &m = refl_mask.At(x, y); - m = std::max(m, level); - } - } - } // --- Signal-region ownership: a pixel inside two reflections' signal regions belongs to the // NEARER predicted centre. The union mask above cannot answer that - it also marks a diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h index 4060c2f06..95b0a943a 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h @@ -57,6 +57,9 @@ class BraggIntegrationEngineCPU : public BraggIntegrationEngine { // it does not explain was unreadable on this frame only - an overload (see Reflection::overloaded). std::vector static_mask; + // Marks every predicted reflection's r2 region in refl_mask (MASK_FLUX / MASK_TAIL). + void MarkNeighbours(const std::vector &predicted, size_t npredicted); + template std::vector RunImpl(const Sampler &img, const std::vector &predicted, size_t npredicted, int64_t image_number); @@ -70,6 +73,12 @@ public: const std::vector &predicted, size_t npredicted, int64_t image_number) override; + // Only the neighbour count of these predicted reflections at this engine's stencil + // (BraggIntegrationCounts::predicted and bkg_starved_by_neighbour), reading no image: what a pass + // integrating at one radius measures of another. A ring pixel counts as readable unless the run's + // mask excludes it, so an overload on one frame is read as readable here. + void CountNeighbourCrowding(const std::vector &predicted, size_t npredicted); + // FPGA workflow: integrate straight off the assembled detector image, reading only the pixels // inside each reflection disk (no whole-image conversion - the FPGA host cannot afford one at its // frame rate). Masked pixels carry the type minimum and saturated the type maximum. diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 4d5e408e2..72bf72174 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -461,6 +461,8 @@ bool Rugnux::NearTieOverturned() { logger.Info("Near ties re-checked on the rest of the sweep: {} first pass(es), {} ending elsewhere", near_tie_rechecks_->pending.size(), elsewhere); near_tie_rechecks_->pending.clear(); + if (near_tie_rechecks_->radius_elsewhere) + elsewhere++; // The metric-symmetry checks decided early, where no first pass has already sent the run again. for (auto &[demoted, promoted, lost] : near_tie_rechecks_->arms) { if (elsewhere > 0) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index a476554ac..efe8986cc 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -899,6 +899,21 @@ class Rugnux { // reflections left their background ring below six clean pixels. Measured, not predicted - it is // the only thing that tells a radius the pattern can take from one it cannot (see RunAllPasses). std::optional bkg_starved_fraction_; + // The same neighbour crowding at the widened radii, counted without integrating there + // (BraggIntegrationEngineCPU::CountNeighbourCrowding) on CROWDING_FRAMES frames spread over a + // pass: over the geometry pre-pass, where the canonical pass's radius is decided on it before that + // pass runs (RefinedPasses), and over the whole sweep in a canonical pass that a near tie left at + // the fixed radius (radius_recheck_). se is the standard error of the fraction over those frames. + struct RingCrowding { + double fraction = NAN; + double se = NAN; + size_t images = 0; + }; + static constexpr int CROWDING_FRAMES = 20; + std::optional widened_crowding_; + // The widened radii a near tie left out of the canonical pass: that pass counts their crowding, + // and where all the frames put it under the bound the run is made again deciding late. + std::optional radius_recheck_; // The header-geometry pass, for as long as the refined passes that are judged against it run (see // the quality guard in RunAllPasses). Held here so a refined pass can ask the guard's own question // of itself as soon as its search merge is in, which is where it learns that its report and its @@ -991,6 +1006,9 @@ class Rugnux { std::atomic carried = false; // a first pass decided a near tie early std::vector pending; std::vector arms; + // The integration radius decided on the pre-pass frames within its noise ends elsewhere on + // all of them (see RefinedPasses). + std::atomic radius_elsewhere = false; }; std::shared_ptr near_tie_rechecks_ = std::make_shared(); // Leaves the first pass the pass just run made to be re-checked; `run` is the copy of this run taken diff --git a/rugnux/RugnuxImagePass.cpp b/rugnux/RugnuxImagePass.cpp index 1498143b2..4c95d0593 100644 --- a/rugnux/RugnuxImagePass.cpp +++ b/rugnux/RugnuxImagePass.cpp @@ -144,6 +144,18 @@ void Rugnux::ImagePass(PipelineLocals &p) { std::mutex harmonic_m; HarmonicEvidence harmonic_evidence; + // The crowding of the background rings at the widened radii, counted on frames spread over the + // pass (widened_crowding_): on the geometry pre-pass, which integrates at the fixed radius while the + // widened one is held back (bragg_adaptive_), and on a canonical pass a near tie left at the fixed + // radius (radius_recheck_). A run deciding late decides the radius on the canonical pass itself. + const std::optional crowding_at = + geometry_prepass ? (decide_late_ ? std::nullopt : bragg_adaptive_) + : prepass_frames_only_ ? std::nullopt : radius_recheck_; + const int crowding_stride = std::max(1, (geometry_prepass ? prepass_images : images_to_process) + / CROWDING_FRAMES); + std::mutex crowding_m; + std::vector> crowding_frames; // predicted, crowded + // Where the indexed spots sit against their predicted reflections, measured on the geometry // pre-pass for the canonical pass's footprint (SpotFootprint.h). const bool measure_offsets = geometry_prepass && prescan_footprint_.has_value(); @@ -228,6 +240,9 @@ void Rugnux::ImagePass(PipelineLocals &p) { AzimuthalIntegrationProfile profile(mapping); HarmonicEvidence worker_harmonic; std::vector worker_offsets; + std::vector> worker_crowding; + if (crowding_at) + analysis.CountCrowdingAt(*crowding_at); while (!cancelled_) { const int ordinal = next_ordinal.fetch_add(1); @@ -261,6 +276,9 @@ void Rugnux::ImagePass(PipelineLocals &p) { msg.image_collection_efficiency = dataset->efficiency[image_idx]; total_uncompressed_bytes += msg.image.GetUncompressedSize(); + const bool count_crowding = crowding_at && ordinal % crowding_stride == 0; + analysis.CountCrowding(count_crowding); + const BraggIntegrationCounts crowding_before = analysis.CrowdingCounts(); const auto t0 = std::chrono::steady_clock::now(); try { analysis.Analyze(msg, profile, config_.spot_finding); @@ -283,6 +301,11 @@ void Rugnux::ImagePass(PipelineLocals &p) { calibration_spot_list.push_back(msg.spots[i]); } + if (count_crowding) { + const BraggIntegrationCounts c = analysis.CrowdingCounts(); + worker_crowding.emplace_back(c.predicted - crowding_before.predicted, + c.bkg_starved_by_neighbour - crowding_before.bkg_starved_by_neighbour); + } AddHarmonicEvidence(msg.spots, experiment_.GetDiffractionGeometry(), worker_harmonic); if (measure_offsets) MeasureFootprintOffsets(msg.spots, msg.reflections, experiment_.GetBeamX_pxl(), @@ -304,6 +327,10 @@ void Rugnux::ImagePass(PipelineLocals &p) { std::lock_guard lock(footprint_m); footprint_offsets.insert(footprint_offsets.end(), worker_offsets.begin(), worker_offsets.end()); } + { + std::lock_guard lock(crowding_m); + crowding_frames.insert(crowding_frames.end(), worker_crowding.begin(), worker_crowding.end()); + } const auto bragg = analysis.BraggCounts(); bragg_predicted += bragg.predicted; @@ -540,6 +567,30 @@ void Rugnux::ImagePass(PipelineLocals &p) { 100.0 * experiment_.GetBraggIntegrationSettings().GetOverlapMinPeak()); } + // The fraction is a ratio of sums over the counted frames, and its standard error is that ratio's, + // from the scatter of the frames about it. + if (crowding_at) { + widened_crowding_.reset(); + uint64_t predicted = 0, crowded = 0; + for (const auto &[p_i, c_i] : crowding_frames) { + predicted += p_i; + crowded += c_i; + } + if (predicted > 0 && !cancelled_) { + const double f = static_cast(crowded) / static_cast(predicted); + const double n = static_cast(crowding_frames.size()); + double ss = 0.0; + for (const auto &[p_i, c_i] : crowding_frames) + ss += std::pow(static_cast(c_i) - f * static_cast(p_i), 2); + const double se = n > 1 ? std::sqrt(ss * n / (n - 1.0)) / static_cast(predicted) : NAN; + widened_crowding_ = RingCrowding{f, se, crowding_frames.size()}; + logger.Info("Counted at r1={:.1f} r2={:.2f} r3={:.2f} without integrating there, on {} images: the " + "neighbouring predictions would crowd {:.3f}% +- {:.3f}% of the {} rings", + crowding_at->GetR1(), crowding_at->GetR2(), crowding_at->GetR3(), crowding_frames.size(), + 100.0 * f, 100.0 * se, predicted); + } + } + // Every image failing is a total failure, not a run that produced nothing: it used to be reported // only as per-image log lines while the process still exited 0 with no output file. if (!cancelled_ && images_to_process > 0 && result.images_processed == 0) diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index 93bcbb6bb..3c1af7b85 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -12,6 +12,7 @@ #include "../image_analysis/hot_pixels/HotPixels.h" #include "../image_analysis/SensorAbsorption.h" #include "DiagnosticOutput.h" +#include "RugnuxFirstPass.h" #include #include @@ -429,6 +430,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { const auto file_bragg_adaptive = bragg_adaptive_; const auto file_bragg_before_adaptive = bragg_before_adaptive_; const auto file_bkg_starved_fraction = bkg_starved_fraction_; + const auto file_widened_crowding = widened_crowding_; const auto file_last_rotation_result = last_rotation_result_; const bool file_supercell_settled = supercell_settled_; @@ -515,6 +517,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { bragg_adaptive_ = file_bragg_adaptive; bragg_before_adaptive_ = file_bragg_before_adaptive; bkg_starved_fraction_ = file_bkg_starved_fraction; + widened_crowding_ = file_widened_crowding; last_rotation_result_ = file_last_rotation_result; supercell_settled_ = file_supercell_settled; } @@ -820,7 +823,31 @@ std::optional Rugnux::RefinedPasses(RugnuxObserver *observer, Pro // bragg_adaptive_), so everything pass 2 consumes from it - the post-refined distance and beam, // the space group, the smoothed mosaicity - is what the fixed radius produces, and the widened // disk changes only the intensities of the canonical pass. - if (bragg_adaptive_) { + // + // Unless the pre-pass already counted that its neighbours would crowd too many of the widened + // background rings (widened_crowding_, the count the starvation guard below measures): the + // canonical pass then integrates at the radius the run started with from the outset, rather than + // integrating at the widened one only to learn that and run again. A count within NEAR_TIE_SIGMA of + // the bound is a near tie: the canonical pass counts the widened radii again, over the whole sweep, + // and where that puts them under the bound the run is made again deciding late (radius_recheck_). + if (bragg_adaptive_ && widened_crowding_ && widened_crowding_->fraction > spot_width::BKG_STARVED_MAX_FRACTION) { + const RingCrowding c = *widened_crowding_; + const bool near_tie = c.fraction - spot_width::BKG_STARVED_MAX_FRACTION + <= FirstPassRun::NEAR_TIE_SIGMA * c.se; + logger.Info("Two-pass: at the widened integration radius r1={:.1f} the neighbouring reflections would leave " + "{:.2f}% +- {:.2f}% of the predicted reflections on {} pre-pass images without a background " + "ring (bound {:.2f}%), so the canonical pass integrates at r1={:.1f} r2={:.1f} r3={:.2f}{}", + bragg_adaptive_->GetR1(), 100.0 * c.fraction, 100.0 * c.se, c.images, + 100.0 * spot_width::BKG_STARVED_MAX_FRACTION, experiment_.GetBraggIntegrationSettings().GetR1(), + experiment_.GetBraggIntegrationSettings().GetR2(), experiment_.GetBraggIntegrationSettings().GetR3(), + near_tie ? " - a near tie, re-checked on the whole sweep" : ""); + note(fmt::format("{:.1f}% of the background rings would starve at r1 {:.1f} px: the canonical pass " + "integrates at r1 {:.1f} px", 100.0 * c.fraction, bragg_adaptive_->GetR1(), + experiment_.GetBraggIntegrationSettings().GetR1())); + if (near_tie) + radius_recheck_ = bragg_adaptive_; + bragg_adaptive_.reset(); + } else if (bragg_adaptive_) { experiment_.ImportBraggIntegrationSettings(*bragg_adaptive_); logger.Info("Two-pass: the canonical pass integrates at the measured spot width, " "r1={:.1f} r2={:.1f} r3={:.2f}", bragg_adaptive_->GetR1(), @@ -994,6 +1021,18 @@ std::optional Rugnux::RefinedPasses(RugnuxObserver *observer, Pro rerun_if_starved_ = true; const DiffractionExperiment before_pass2 = experiment_; auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + if (radius_recheck_) { + // The radius decided on the pre-pass frames within its noise, re-checked on the whole sweep. + const bool elsewhere = widened_crowding_ + && widened_crowding_->fraction <= spot_width::BKG_STARVED_MAX_FRACTION; + logger.Info("Near-tie re-check: on {} images over the whole sweep the widened radius r1={:.1f} would leave " + "{:.2f}% of the predicted reflections without a background ring (bound {:.2f}%) - the " + "integration radius {}", widened_crowding_ ? widened_crowding_->images : 0, + radius_recheck_->GetR1(), widened_crowding_ ? 100.0 * widened_crowding_->fraction : NAN, + 100.0 * spot_width::BKG_STARVED_MAX_FRACTION, elsewhere ? "ends elsewhere" : "ends where it did"); + near_tie_rechecks_->radius_elsewhere = elsewhere; + radius_recheck_.reset(); + } if (lattice_arm_check_) { // The metric-symmetry check above, decided inside the pass just run (the promoted arm). const LatticeArmCheck arm = *lattice_arm_check_; diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp index 8c68bb5e6..067723b69 100644 --- a/tests/BraggIntegrationEngineCPUTest.cpp +++ b/tests/BraggIntegrationEngineCPUTest.cpp @@ -77,6 +77,41 @@ TEST_CASE("BraggIntegrationEngineCPU_NeighbourMaskFollowsPartiality", "[Integrat CHECK(on_frame.bkg_starved_by_neighbour == tails.bkg_starved_by_neighbour); } +// The neighbour count taken without an image (CountNeighbourCrowding) is the one the integrator +// takes while integrating, on a frame whose pixels are all readable, whatever the partialities. +TEST_CASE("BraggIntegrationEngineCPU_CrowdingCountedWithoutImage", "[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().R1(5.0f).R2(7.0f).R3(13.0f)); + const size_t npixel = experiment.GetPixelsNum(); + ImagePreprocessorBuffer image(npixel); + for (size_t i = 0; i < npixel; ++i) + image[i] = 12; + + for (float spacing : {7.0f, 11.0f, 16.0f}) { + std::vector predicted; + for (int gy = -6; gy <= 6; ++gy) + for (int gx = -6; gx <= 6; ++gx) { + Reflection r{}; + r.h = gx; r.k = gy; r.l = 1; + r.predicted_x = 600.3f + spacing * gx + 0.37f * gy; + r.predicted_y = 300.2f + spacing * gy; + r.d = 2.0f; + r.prescaling_corr = 1.0f; + r.partiality = (gx + gy) % 3 == 0 ? 0.01f : 0.8f; + predicted.push_back(r); + } + BraggIntegrationEngineCPU integrated(experiment, PixelMask(experiment)); + integrated.Run(image, predicted, predicted.size(), 0); + BraggIntegrationEngineCPU counted(experiment, PixelMask(experiment)); + counted.CountNeighbourCrowding(predicted, predicted.size()); + CHECK(counted.Counts().predicted == integrated.Counts().predicted); + CHECK(counted.Counts().bkg_starved_by_neighbour == integrated.Counts().bkg_starved_by_neighbour); + if (spacing == 7.0f) + CHECK(counted.Counts().bkg_starved_by_neighbour > 0); + } +} + // A saturated pixel (INT32_MAX, as the preprocessor marks it) at the peak of a spot: the reflection is // not measured but kept, flagged overloaded, so that a rotation merge can drop its whole rocking event. // The same spot without the saturated pixel is an ordinary measurement.