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.