diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index b1b9dd74d..92e197e0e 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -408,23 +408,38 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, // the scaling and the sweep ledger to judge as they judge any other. Only from the ends inward: a frame // inside the sweep without beam is news, and is left to the ledger too. A frame with nothing integrated // on it (missing from disk, not processed) is walked past and not counted. -std::pair TrimBlankSweepEnds(std::vector &outcomes, size_t nthreads) { +namespace { constexpr double BLANK_TAIL_SIGMA = 3.0; +} + +SweepEndTally TallySweepEnd(const IntegrationOutcome &outcome) { + SweepEndTally t; + for (const auto &r : outcome.reflections) { + if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f) { + t.upper += r.I > BLANK_TAIL_SIGMA * r.sigma; + t.lower += r.I < -BLANK_TAIL_SIGMA * r.sigma; + } + if (std::isfinite(r.bkg)) { t.bkg_sum += r.bkg; ++t.bkg_n; } + } + return t; +} + +std::pair TrimBlankSweepEnds(std::vector &outcomes, size_t nthreads) { + std::vector tally(outcomes.size()); + ParallelFor(static_cast(outcomes.size()), nthreads, [&](int f) { tally[f] = TallySweepEnd(outcomes[f]); }); + return TrimBlankSweepEnds(outcomes, tally); +} + +std::pair TrimBlankSweepEnds(std::vector &outcomes, + const std::vector &tally) { const int n = static_cast(outcomes.size()); std::vector upper(n, 0), lower(n, 0); std::vector bkg(n, 0.0); - ParallelFor(n, nthreads, [&](int f) { - double bkg_sum = 0.0; - int64_t bkg_n = 0; - for (const auto &r : outcomes[f].reflections) { - if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f) { - upper[f] += r.I > BLANK_TAIL_SIGMA * r.sigma; - lower[f] += r.I < -BLANK_TAIL_SIGMA * r.sigma; - } - if (std::isfinite(r.bkg)) { bkg_sum += r.bkg; ++bkg_n; } - } - if (bkg_n > 0) bkg[f] = bkg_sum / bkg_n; - }); + for (int f = 0; f < n; ++f) { + upper[f] = tally[f].upper; + lower[f] = tally[f].lower; + if (tally[f].bkg_n > 0) bkg[f] = tally[f].bkg_sum / tally[f].bkg_n; + } std::vector excess, tails, background; for (int f = 0; f < n; ++f) if (!outcomes[f].reflections.empty()) { diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index ced81eccb..ef4b730f5 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -67,7 +67,17 @@ template using FillVector = std::vector>; // Stills use the per-image ScaleOnTheFly (fixed partiality) instead. // // Before any of it: TrimBlankSweepEnds clears the reflections of the frames at either end of the sweep -// that the beam never reached, and returns how many it cleared at the start and at the end. +// that the beam never reached, and returns how many it cleared at the start and at the end. It decides +// on each frame's SweepEndTally, one walk over the frame's reflections in their order; a caller walking +// them anyway for something else hands the tallies over, the other form makes them on all threads. +struct SweepEndTally { + int64_t upper = 0, lower = 0; // reflections standing above / below zero by more than the tail sigma + double bkg_sum = 0.0; + int64_t bkg_n = 0; +}; +SweepEndTally TallySweepEnd(const IntegrationOutcome &outcome); +std::pair TrimBlankSweepEnds(std::vector &outcomes, + const std::vector &tally); std::pair TrimBlankSweepEnds(std::vector &outcomes, size_t nthreads); class RotationScaleMerge { diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 51d157475..af29b845f 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -507,8 +507,31 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // and what a model's low-resolution terms rest on. The ice flags are set in the same sweep. Integer // counts, one image per task: the totals do not depend on the split. TimingMark ends_mark(logger, "sweep ends, overloads, ice rings"); + // One walk over the reflections, per frame, for all three: what the sweep ends are decided on, the + // overload counts and the ice flags. The counts are added up below over the frames still in. + auto &outcomes = indexer->GetIntegrationOutcome(); + struct FrameCounts { int64_t all = 0, low = 0, over = 0, over_low = 0, on_ice = 0; }; + std::vector tally(outcomes.size()); + std::vector counts(outcomes.size()); + ParallelFor(static_cast(outcomes.size()), config_.nthreads, [&](int o) { + tally[o] = TallySweepEnd(outcomes[o]); + FrameCounts &c = counts[o]; + for (auto &r : outcomes[o].reflections) { + const bool low_res = r.d >= OVERLOAD_LOW_RES_A; + ++c.all; + c.low += low_res; + if (r.overloaded) { + ++c.over; + c.over_low += low_res; + } + if (flag_ice) { + r.on_ice_ring = IsOnIceRing(r.d, ice_width); + c.on_ice += r.on_ice_ring; + } + } + }); if (experiment_.IsRotationIndexing()) { - const auto [at_start, at_end] = TrimBlankSweepEnds(indexer->GetIntegrationOutcome(), config_.nthreads); + const auto [at_start, at_end] = TrimBlankSweepEnds(outcomes, tally); result.blank_start_frames = at_start; result.blank_end_frames = at_end; if (at_start + at_end > 0) @@ -518,7 +541,6 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // Stills given a cell: only the frames whose refined cell belongs to its population are merged // (see GivenCellPopulation) - a second crystal form, fitted to the same cell, is left out. if (!experiment_.IsRotationIndexing() && stills_user_cell_) { - auto &outcomes = indexer->GetIntegrationOutcome(); const double given = CrystalLattice(*stills_user_cell_).CalcVolume(); std::vector rel; for (const auto &o : outcomes) @@ -550,25 +572,13 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { } size_t ice_flagged = 0, ice_total = 0; { - auto &outcomes = indexer->GetIntegrationOutcome(); - 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; - for (auto &r : outcomes[o].reflections) { - const bool low_res = r.d >= OVERLOAD_LOW_RES_A; - ++all; - low += low_res; - if (r.overloaded) { - ++over; - over_low += low_res; - } - if (flag_ice) { - r.on_ice_ring = IsOnIceRing(r.d, ice_width); - on_ice += r.on_ice_ring; - } - } - n_all += all; n_low += low; n_over += over; n_over_low += over_low; n_ice += on_ice; - }); + int64_t n_all = 0, n_low = 0, n_over = 0, n_over_low = 0, n_ice = 0; + for (size_t o = 0; o < outcomes.size(); ++o) { + if (outcomes[o].reflections.empty()) // a frame left out above counts nothing + continue; + const FrameCounts &c = counts[o]; + n_all += c.all; n_low += c.low; n_over += c.over; n_over_low += c.over_low; n_ice += c.on_ice; + } result.integrated_reflections += n_all; result.integrated_reflections_low_res += n_low; result.overloaded_reflections += n_over;