diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 50e64e775..9b198fa66 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -2,6 +2,7 @@ ## 1.0.0 ### 1.0.0-rc.166 * The de-novo space-group search ignores the frames whose fitted per-frame scale came out below a tenth of the run median. Its merge is in P1, where a reflection has too few observations for anything to catch an intensity that scale amplified; the production merge still keeps every frame. +* rugnux measures its own spot budget on rotation data instead of keeping the strongest 1000 spots of every image: the first pass reads how deep into an image's spot list the spots still lie on the lattice it found, and the run keeps that many. `--max-spots` still pins the budget, and stills are unchanged. * Self-calibrating spot detection intersects its per-resolution-ring threshold with the local signal-to-noise test again, instead of replacing it. The ring threshold takes the place of the fixed photon floor and nothing else; standing alone it followed a bright reflection's skirt outwards, so on a strongly diffracting crystal the brightest reflections were detected as 100-500 pixel blobs and then discarded for being too large. * The self-calibrating threshold no longer steps where it switches from the exact Poisson tail to a normal approximation. * A spot may hold up to 200 connected pixels rather than 50. Under the self-calibrating threshold a reflection's footprint grows with its brightness, so the old bound acted as an intensity ceiling and discarded the strongest reflections of a strongly diffracting crystal - on one such set, every one of the ten brightest on an image. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 68ed963a0..481c85a71 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -434,6 +434,20 @@ A spot is indexed if $\delta^2 < \tau^2$, where $\tau$ is the configured toleran For indexed spots, the reciprocal lattice point $\mathbf{p} = h\mathbf{a}^*+k\mathbf{b}^*+l\mathbf{c}^*$ is used to compute $\Delta_\mathrm{Ewald}(\mathbf{p})$ (stored as a diagnostic and later used in profile-radius estimation). +A frame is taken to be this crystal's when at least a fraction $g = 0.20$ of its in-resolution, non-ice spots index. On rotation data that decision is what admits the frame to integration, so its denominator matters: every spot handed to it that is not a reflection of this crystal argues against the frame. + +### 4.2 The spot budget + +Only the strongest `--max-spots` spots of an image are kept (`FilterSpotsByCount`), and that budget therefore sets the denominator above. Detections are not all reflections — background structure, unlisted ice and detector artefacts are found too — so a budget deeper than an image's reflections makes the test above a measurement of the background rather than of the crystal, and a *larger* budget can integrate *fewer* images. + +rugnux measures the budget instead of fixing it. With the sweep's lattice in hand, the first pass tallies the spots of a sample of frames by their rank in the intensity-ordered list: how many images carried a spot at that rank, and on how many of them it indexed. Weighting each indexed spot by $1-g$ and each unindexed one by $-g$ — the same weighing the frame test applies to the list as a whole — the running sum over ranks + +$$ E(N) = n_\mathrm{indexed}(N) - g\, n_\mathrm{counted}(N) $$ + +rises exactly while the spots at that depth lie on the lattice more often than $g$, and falls after. The budget is $\arg\max_N E(N)$. Its meaning is "as deep into the list as the image is still showing reflections of this crystal": deeper spots cannot help the frame test and can only push a frame towards rejection. On crystals whose spot lists are reflections all the way down the maximum is at the end of the list and the budget is unchanged. + +The peak has to be one. Under the null — the spots lie on the lattice at the same rate at every depth — $E$ is a driftless random walk in the counted spots, with per-spot variance $g(1-g)$, and the maximum of such a walk is positive whatever the data; an $\arg\max$ taken on its own would shorten every dataset, including one with nothing to shorten. What the budget acts on is the fall from the peak to the end of the list, $E(N^*) - E(L)$, which is the maximum of the same walk read backwards; by the reflection principle its null law is $P(\mathrm{fall} > z\sqrt{g(1-g)T}) = 2(1-\Phi(z))$ over $T$ counted spots in all, so the search over ranks is already accounted for and no further multiple-comparison correction applies. The budget is taken only where the fall clears that bar at $z = 3.29$, one false shortening in a thousand measurements; otherwise the whole list is kept. + --- ## 5. FFT indexing (unknown unit cell) diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 77508371e..e91a44de2 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -957,7 +957,7 @@ Spot finding: | `--spot-high-resolution ` | High-resolution limit for spot finding, Å. Omitted (or 0): no resolution clipping — spot finding extends as far as the detector reaches, for rotation data as well as stills | | `--spot-low-resolution ` | Low-resolution limit for spot finding, Å (default: 50; lower it, e.g. 24, to exclude the direct-beam halo on weak serial data; 0 removes the limit) | | `--min-pix-per-spot ` | Minimum connected strong pixels per spot. **If omitted, min-pix is chosen per image** (stills indexing): the frame is indexed at min-pix 3/2/1 and the one maximising indexed-spot count × indexed fraction is kept. Give an explicit value to force a fixed min-pix instead. | -| `--max-spots ` | Maximum spots kept per image (the strongest ones) and handed to indexing (default: 1000) | +| `--max-spots ` | Maximum spots kept per image (the strongest ones) and handed to indexing. **If omitted, the budget is measured** on rotation data: the first pass reads how deep into an image's spot list its spots still lie on the lattice it found, and the run keeps that many (never more than 1000). Give a value to pin it. Stills always use the fixed 1000. | | `--detect-ice-rings[=on\|off]` | Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling. Default: the master file's `detect_ice_rings`, or — where the file carries no such key — **on for rotation and off for stills** | Azimuthal integration (the radial profile behind the per-image ice-ring score): diff --git a/image_analysis/indexing/AnalyzeIndexing.cpp b/image_analysis/indexing/AnalyzeIndexing.cpp index 056b8fe95..a7adc108b 100644 --- a/image_analysis/indexing/AnalyzeIndexing.cpp +++ b/image_analysis/indexing/AnalyzeIndexing.cpp @@ -420,14 +420,13 @@ bool AnalyzeIndexing(DataMessage &message, int64_t indexing_lattice_count = 0; bool outcome = false; - // Minimum fraction of the in-resolution spots a candidate lattice must index to be accepted. - // Lowering it admits weaker/sparser crystals (more real ones on flooded XFEL frames, but also more - // spurious lattices that a downstream merge-consistency gate must remove). The gate is a stills - // notion (CrystFEL, White et al., J. Appl. Cryst. 45, 335-341 (2012)): XDS, MOSFLM and DIALS index - // once over the sweep and then integrate every frame from that lattice - none of them re-decides - // per frame whether a frame may be integrated. - constexpr float min_frac = 0.20f; - const bool lattice_fits = nspots_indexed >= std::lround(min_frac * nspots_ref); + // LATTICE_MIN_INDEXED_FRACTION is the minimum fraction of the in-resolution spots a candidate + // lattice must index to be accepted. Lowering it admits weaker/sparser crystals (more real ones on + // flooded XFEL frames, but also more spurious lattices that a downstream merge-consistency gate must + // remove). The gate is a stills notion (CrystFEL, White et al., J. Appl. Cryst. 45, 335-341 (2012)): + // XDS, MOSFLM and DIALS index once over the sweep and then integrate every frame from that lattice - + // none of them re-decides per frame whether a frame may be integrated. + const bool lattice_fits = nspots_indexed >= std::lround(LATTICE_MIN_INDEXED_FRACTION * nspots_ref); // Two different questions. "Does this frame index?" - reported as the indexing rate, and what the // rotation first pass scores candidate lattices on - needs the absolute floor too, because a // handful of spots sit on almost any lattice by chance. "Is this frame worth integrating?" needs @@ -522,3 +521,38 @@ bool AnalyzeIndexing(DataMessage &message, message.indexing_result = outcome && frame_indexes; return outcome; } + +void AddSpotBudgetEvidence(const std::vector &spots, bool index_ice_rings, + std::vector &indexed, std::vector &counted) { + const size_t n = std::min({spots.size(), indexed.size(), counted.size()}); + for (size_t i = 0; i < n; i++) { + if (!index_ice_rings && spots[i].ice_ring) + continue; + counted[i]++; + if (spots[i].indexed) + indexed[i]++; + } +} + +int64_t SpotBudgetFromEvidence(const std::vector &indexed, const std::vector &counted) { + int64_t indexed_sum = 0; + int64_t counted_sum = 0; + int64_t budget = 0; + float best = 0.0f; + float at_end = 0.0f; + for (size_t i = 0; i < std::min(indexed.size(), counted.size()); i++) { + indexed_sum += indexed[i]; + counted_sum += counted[i]; + at_end = static_cast(indexed_sum) + - LATTICE_MIN_INDEXED_FRACTION * static_cast(counted_sum); + if (at_end > best) { + best = at_end; + budget = static_cast(i) + 1; + } + } + + const float g = LATTICE_MIN_INDEXED_FRACTION; + const float noise = SPOT_BUDGET_SIGNIFICANCE_Z + * std::sqrt(g * (1.0f - g) * static_cast(counted_sum)); + return (best - at_end > noise) ? budget : 0; +} diff --git a/image_analysis/indexing/AnalyzeIndexing.h b/image_analysis/indexing/AnalyzeIndexing.h index 1e612343d..3fb2d1f1b 100644 --- a/image_analysis/indexing/AnalyzeIndexing.h +++ b/image_analysis/indexing/AnalyzeIndexing.h @@ -7,6 +7,45 @@ #include "../../common/DiffractionExperiment.h" #include "../../common/JFJochMessages.h" +// Minimum fraction of a frame's in-resolution spots that must lie on a candidate lattice for the +// frame to be that crystal's. See the frame gate in AnalyzeIndexing, which is where it is applied. +constexpr float LATTICE_MIN_INDEXED_FRACTION = 0.20f; + +// Tally one image's spots by their rank in its intensity-ordered spot list: how many images had a spot +// at that rank at all (`counted`) and on how many of them it lay on the lattice (`indexed`). Ice spots +// are skipped, as they are in the frame gate. Both are added to, and their length bounds the ranks +// considered. Counts rather than weights so that the tally is exact whatever order the images are +// summed in, which is what makes the budget below independent of the thread schedule. +void AddSpotBudgetEvidence(const std::vector &spots, bool index_ice_rings, + std::vector &indexed, std::vector &counted); + +// How far the fall from the peak must exceed the counting noise of the spots for the peak to be one. +// +// Under the null - the spots lie on the lattice at the same rate at every depth - the running sum +// below is a driftless random walk in the counted spots: each is worth 1 - g with probability g and +// -g otherwise, so its step has mean zero and variance g(1-g). The MAXIMUM of such a walk is positive +// whatever the data, so an argmax taken on its own cuts every dataset, including one with nothing to +// cut. What is acted on is the FALL from the peak to the end of the list, which is the maximum of the +// same walk read backwards from the end, and the reflection principle gives that maximum's null law +// exactly: P(fall > z sqrt(g(1-g)T)) = 2(1 - Phi(z)) over T counted spots in all. The search over the +// ranks is therefore already paid for and no further multiple-comparison correction is due. z is set +// for one false cut in a thousand measurements, which over a corpus the size of a rotation test set +// (tens of crystals, a measurement per pass) expects none at all. +constexpr float SPOT_BUDGET_SIGNIFICANCE_Z = 3.29f; // 2(1 - Phi(z)) = 0.001 + +// The spot budget those tallies support: the rank at which indexed - LATTICE_MIN_INDEXED_FRACTION * +// counted, summed over the ranks down to it, peaks. Each spot that lies on the lattice is worth +// 1 - LATTICE_MIN_INDEXED_FRACTION and each one that does not costs LATTICE_MIN_INDEXED_FRACTION - the +// same weighing the frame gate applies to a spot list as a whole - so the sum rises exactly while the +// spots at that depth are on the lattice more often than the gate's floor. Deeper than the peak they +// are not: they are no longer this crystal's reflections, and they can only push a frame towards +// rejection while adding nothing the lattice recognises. +// +// Zero - keep the whole list - when the fall from that peak to the end of the list is no larger than +// the counting noise above, which is the case whenever the spots go on lying on the lattice at the +// same rate all the way down, and the case a bare argmax gets wrong. +int64_t SpotBudgetFromEvidence(const std::vector &indexed, const std::vector &counted); + bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experiment, const CrystalLattice &latt, diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index ed1e51b63..db0f96f86 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -31,6 +31,7 @@ #include "../image_analysis/MXAnalysisWithoutFPGA.h" #include "../image_analysis/beam_stop/ShadowFinder.h" #include "../image_analysis/IndexAndRefine.h" +#include "../image_analysis/indexing/AnalyzeIndexing.h" #include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" #include "../image_analysis/geom_refinement/BeamCenterFromSpots.h" #include "../image_analysis/geom_refinement/GeometryRefiner.h" @@ -996,6 +997,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { Logger logger("Rugnux"); const std::string base_prefix = config_.output_prefix; const auto gonio_snapshot = experiment_.GetGoniometer(); + const int64_t max_spot_count_snapshot = experiment_.GetDatasetSettings().GetMaxSpotCount(); prepass_detector_geometry_.reset(); prepass_rotation_scale_.reset(); prepass_result_.reset(); @@ -1014,6 +1016,10 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { pass1.pass_count = 2; if (cancelled_) { config_.output_prefix = base_prefix; return pass1; } if (gonio_snapshot) experiment_.Goniometer(*gonio_snapshot); // undo the pre-pass goniometer shift + // A measured spot budget is the pass's own, not the run's: give the second pass the same list to + // measure from, so the two passes cannot ratchet each other down. The canonical pass then reads + // the budget off the refined geometry, where the spots that do lie on the lattice actually do. + experiment_.MaxSpotCount(max_spot_count_snapshot); // Apply the post-refined detector geometry for the second pass, keeping the header geometry so // the run can go back to it if the refined pass turns out worse (see the quality guard below). @@ -1662,6 +1668,32 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b return count.load(); }; + // How deep into an image's intensity-ordered spot list this lattice is still being seen - the + // measured spot budget. Same frames, same per-image path as count_indexed above, but scoring the + // spots rather than the frames: each frame's spots are tallied by rank and the run keeps the + // ranks the tally supports (AddSpotBudgetEvidence / SpotBudgetFromEvidence). A budget deeper than + // that admits detections that are not this crystal's reflections, and the frame gate then counts + // them against the frame - which is how a larger budget can integrate FEWER images. + auto measure_budget = [&](IndexAndRefine &idx, const RotationIndexerResult &r) -> int64_t { + idx.ForceRotationIndexerResult(r); + prefetch_spots(validation); + std::mutex evidence_mutex; + std::vector indexed(experiment_.GetMaxSpotCount(), 0); + std::vector counted(indexed.size(), 0); + const bool index_ice_rings = experiment_.GetIndexingSettings().GetIndexIceRings(); + ParallelFor(static_cast(validation.size()), + std::min(validation.size(), config_.nthreads), [&](int i) { + DataMessage m{}; + m.number = validation[i]; + m.spots = spot_cache.at(validation[i]); + if (!idx.IndexFrameOnly(m, validation_settings)) + return; + std::unique_lock ul(evidence_mutex); + AddSpotBudgetEvidence(m.spots, index_ice_rings, indexed, counted); + }); + return SpotBudgetFromEvidence(indexed, counted); + }; + // Feed one first-pass scheme (a set of image ordinals) into its own rotation indexer, ready to // be indexed. Spots are pulled from the cache here; the FFT + refinement runs separately // (RunIndexing) so the two schemes' indexing can overlap. @@ -1922,6 +1954,19 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b logger.Info("Two-pass rotation indexing found lattice (scheme '{}': {}/{} validation frames)", best.name, best.score, static_cast(validation.size())); + // The spots in hand were already cut to the budget in force, so the measurement can only + // shorten a budget, never lengthen one - which is why RunAllPasses hands the second pass the + // list the first one started from rather than the list the first one settled on. + if (config_.measure_spot_budget) { + const int64_t budget = measure_budget(*indexer, *best.result); + if (budget > 0 && budget < experiment_.GetMaxSpotCount()) { + logger.Info("Spot budget: keeping the strongest {} spots per image " + "(measured on {} frames; was {})", budget, + static_cast(validation.size()), experiment_.GetMaxSpotCount()); + experiment_.MaxSpotCount(budget); + } + } + // Second pass: compare the de-novo lattice with pass 1's HERE, before integrating every image // with it. A markedly larger cell is the bistable supercell collapse, and a centring the group // carried over from pass 1 cannot describe is the same disagreement seen from the symmetry side. diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 4b49e0a12..1e863dc5c 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -68,6 +68,13 @@ struct ProcessConfig { CalibrationMethod calibration_method = CalibrationMethod::Rings; std::vector calibrant_ring_q; + // Measure the spot budget instead of taking experiment.GetMaxSpotCount() as given (FullAnalysis, + // rotation only; the rugnux CLI sets it unless --max-spots pinned the budget). The first pass reads + // how deep into each image's spot list the spots still lie on the lattice it found, and the run + // keeps that many - see AddSpotBudgetEvidence. Off by default, so the online receiver and the viewer + // keep the fixed budget they are configured with. + bool measure_spot_budget = false; + // Rotation indexing (FullAnalysis) bool rotation_indexing = false; bool two_pass_rotation = true; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 790e7d03a..5019e0a05 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -116,7 +116,7 @@ void print_usage() { std::cout << " --spot-false-pixels Adaptive detection operating point: expected noise pixels tolerated per frame (default: 100; implies --adaptive-spots)" << std::endl; std::cout << " --spot-high-resolution High resolution limit for spot finding. If omitted (or 0), spot finding is not clipped in resolution and extends as far as the detector reaches" << std::endl; std::cout << " --spot-low-resolution Low resolution limit for spot finding, in A (default: 50; lower it, e.g. 24, to exclude the direct-beam halo on weakly-diffracting serial data; 0 removes the limit)" << std::endl; - std::cout << " --max-spots Max spot count per image, the strongest ones, handed to indexing (default: 1000)" << std::endl; + std::cout << " --max-spots Max spot count per image, the strongest ones, handed to indexing. Default: measured from the data on rotation (as deep into each image's spot list as its spots still lie on the lattice), capped at 1000" << std::endl; std::cout << " --detect-ice-rings[=on|off] Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling. Default: the master file's setting, or - where the file says nothing - on for rotation and off for stills" << std::endl; std::cout << std::endl; @@ -2312,6 +2312,8 @@ static int RunRugnux(int argc, char **argv) { config.fit_spindle = fit_spindle; config.adaptive_integration_radius = adaptive_integration_radius; config.rotation_postrefine_geometry = rotation_postrefine_geometry; + // Measure the spot budget from the data unless the user pinned it. + config.measure_spot_budget = !max_spot_count_override.has_value(); config.rotation_scale = rotation_scale; config.rotation_indexing_image_count = rotation_indexing_image_count; config.forced_rotation_lattice = forced_rotation_lattice; diff --git a/tests/SpotUtilsTest.cpp b/tests/SpotUtilsTest.cpp index 0a804c9f0..28abb93c5 100644 --- a/tests/SpotUtilsTest.cpp +++ b/tests/SpotUtilsTest.cpp @@ -4,6 +4,7 @@ #include #include "../image_analysis/spot_finding/SpotUtils.h" +#include "../image_analysis/indexing/AnalyzeIndexing.h" TEST_CASE("FilterSpuriousHighResolutionSpots") { std::vector spots; @@ -59,3 +60,62 @@ TEST_CASE("GetResolution") { // Too few spots to have a fall-off at all. CHECK_FALSE(GetResolution(std::vector(3)).has_value()); } + +TEST_CASE("SpotBudgetFromEvidence") { + // One image's worth of spots, repeated over 60 frames as the first pass does: the first 100 index + // and the next 100 do not. Every indexed spot adds 1 - 0.2 and every unindexed one takes 0.2 away, + // so the running tally rises to rank 100 and falls after it. + std::vector spots(200); + for (size_t i = 0; i < spots.size(); i++) + spots[i].indexed = i < 100; + + constexpr int frames = 60; + std::vector indexed(spots.size(), 0), counted(spots.size(), 0); + for (int f = 0; f < frames; f++) + AddSpotBudgetEvidence(spots, false, indexed, counted); + CHECK(SpotBudgetFromEvidence(indexed, counted) == 100); + + // Spots that go on indexing all the way down: the tally never falls, so there is nothing to cut. + for (auto &s: spots) + s.indexed = true; + std::vector all_hit(spots.size(), 0), all_seen(spots.size(), 0); + for (int f = 0; f < frames; f++) + AddSpotBudgetEvidence(spots, false, all_hit, all_seen); + CHECK(SpotBudgetFromEvidence(all_hit, all_seen) == 0); + + // A budget already cut to its peak has no fall left in it, so a second measurement takes nothing + // further off: the rule does not ratchet down on repetition. + CHECK(SpotBudgetFromEvidence({indexed.begin(), indexed.begin() + 100}, + {counted.begin(), counted.begin() + 100}) == 0); + + // Ice-flagged spots take no part, so a run of them neither ends the budget nor moves it: the peak + // stays at the last indexed non-ice rank before them. + std::vector with_ice(300); + for (size_t i = 0; i < with_ice.size(); i++) { + with_ice[i].ice_ring = (i >= 100 && i < 160); + with_ice[i].indexed = i < 100; + } + std::vector ice_indexed(with_ice.size(), 0), ice_counted(with_ice.size(), 0); + for (int f = 0; f < frames; f++) + AddSpotBudgetEvidence(with_ice, false, ice_indexed, ice_counted); + CHECK(SpotBudgetFromEvidence(ice_indexed, ice_counted) == 100); + + // Nothing indexes: no rank carries evidence and there is no budget to report. + for (auto &s: spots) + s.indexed = false; + std::vector none_indexed(spots.size(), 0), none_counted(spots.size(), 0); + for (int f = 0; f < frames; f++) + AddSpotBudgetEvidence(spots, false, none_indexed, none_counted); + CHECK(SpotBudgetFromEvidence(none_indexed, none_counted) == 0); + + // The case a bare argmax gets wrong: the spots index at exactly the gate's own fraction at every + // depth, so there is no depth at which the list stops being reflections. The tally still has a + // maximum - it always does - but the fall from it is inside the counting noise, and nothing is cut. + std::vector flat(1000); + for (size_t i = 0; i < flat.size(); i++) + flat[i].indexed = (i % 5 == 0); + std::vector flat_indexed(flat.size(), 0), flat_counted(flat.size(), 0); + for (int f = 0; f < frames; f++) + AddSpotBudgetEvidence(flat, false, flat_indexed, flat_counted); + CHECK(SpotBudgetFromEvidence(flat_indexed, flat_counted) == 0); +}