diff --git a/common/JFJochReceiverPlots.cpp b/common/JFJochReceiverPlots.cpp index 17116054..23960661 100644 --- a/common/JFJochReceiverPlots.cpp +++ b/common/JFJochReceiverPlots.cpp @@ -472,6 +472,16 @@ std::optional JFJochReceiverPlots::GetBkgEstimate() const { return {}; } +std::optional JFJochReceiverPlots::GetResolutionEstimate() const { + std::vector v = resolution_estimate.ExportArray(); + std::erase_if(v, [](float x) { return !std::isfinite(x); }); + if (v.empty()) + return {}; + const size_t mid = v.size() / 2; + std::ranges::nth_element(v, v.begin() + mid); + return v[mid]; +} + std::optional JFJochReceiverPlots::GetIceRingSpotRatio() const { // A ratio of MEANS, not a mean of ratios: one image holds a handful of control spots, so a // per-image ratio is dominated by its own denominator. Pooling over the run is the measurement. diff --git a/common/JFJochReceiverPlots.h b/common/JFJochReceiverPlots.h index 3b711a40..a7ebdddf 100644 --- a/common/JFJochReceiverPlots.h +++ b/common/JFJochReceiverPlots.h @@ -126,6 +126,10 @@ public: // Pooled over the run: spots on the hexagonal rings over the same q width of ice-free control // flanks. 1 = spots spread evenly, > 1 = they pile up on the rings (textured ice). [[nodiscard]] std::optional GetIceRingSpotRatio() const; + // The run's spot-finding resolution estimate: the MEDIAN of the per-image estimates. A median and + // not a mean - a few blank or badly-shot frames pull a mean a long way and leave the median where + // the run actually is. + [[nodiscard]] std::optional GetResolutionEstimate() const; std::vector GetIceRingScoreArray() const; std::vector GetAzIntProfile() const; diff --git a/docs/CBOR.md b/docs/CBOR.md index b961e8ff..c887b6f3 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -213,7 +213,7 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | strong_pixel_count | uint64 | Strong pixel count (first stage of spot finding) | | | | min_viable_pixel_value | int64 | Minimal pixel value, excl. error and saturation | | | | max_viable_pixel_value | int64 | Maximal pixel value, excl. error and saturation | | | -| resolution_estimate | float | Diffraction resolution estimation \[Angstrom\] | | X | +| resolution_estimate | float | Resolution the merged data are predicted to reach, from this image's spots alone \[Angstrom\] | | X | | data_collection_efficiency | float | Image collection efficiency \[\] | | | | packets_expected | uint64 | Number of packets expected per image (in units of 2 kB) | | | | packets_received | uint64 | Number of packets received per image (in units of 2 kB) | | | @@ -268,7 +268,7 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | - strong_pixel_count | uint64 | Strong pixel count (first stage of spot finding) | | | | - data_collection_efficiency | float | Image collection efficiency \[\] | | | | - bkg_estimate | float | Mean value for pixels in resolution range from 3.0 to 5.0 A \[photons\] (with solid angle/polarization corrections, if applied) | | X | -| - resolution_estimate | float | Diffraction resolution estimation | | X | +| - resolution_estimate | float | Predicted merged resolution, from spots alone | | X | | - adu_histogram | Array(uint64) | ADU histogram | | X | | - roi_integrals | object | Results of ROI calculation | | X | | - - sum | int64 | Sum of pixels in ROI area \[photons\] | | | @@ -316,7 +316,7 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | profile_radius | Array(float) | Per-image profile radius \[Angstrom^-1\] | | | mosaicity | Array(float) | Per-image mosaicity \[degree\] | | | bFactor | Array(float) | Per-image estimated B-factor \[Angstrom^2\] | | -| resolution_estimate | Array(float) | Per-image diffraction resolution estimate \[Angstrom\] | | +| resolution_estimate | Array(float) | Per-image predicted merged resolution, from spots alone \[Angstrom\] | | | min_viable_pixel_value | Array(int64) | Per-image minimum valid pixel value, excluding error/saturated pixels | | | max_viable_pixel_value | Array(int64) | Per-image maximum valid pixel value, excluding error/saturated pixels | | | saturated_pixel_count | Array(int32) | Per-image saturated pixel count | | diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 1e0d37c3..c5bf88fe 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -18,6 +18,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * HDF5 and image stream: `mirror_y` records whether the assembled image is mirrored in Y relative to the detector's raw readout. * rugnux: an image integrated in pyFAI through the `.poni` file written by `--mode calibration` now comes out with the correct azimuth. Radial integration is unchanged. * rugnux: the `.poni` file declares pyFAI's `orientation`, which needs pyFAI 2024.01 or newer. +* The per-image resolution estimate now predicts the resolution the merged data reach, rather than reporting the highest-resolution spot found; rugnux reports the run's value as `SPOT_RESOLUTION_ESTIMATE` in its report. * rugnux: fixing the space group with `-S` no longer prevents the lattice from being found; the group is applied to scaling and merging rather than to the indexing search. * rugnux: the detector geometry is also logged in XDS's convention (`ORGX`/`ORGY`, detector axis vectors, rotation axis), so it can be compared with an XDS refinement. * rugnux: `_process.h5` describes the pixel format of the images it links to. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 9391fe29..4bb910a8 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -318,6 +318,10 @@ Both channels are used offline as a **gate** on ice handling: unless the run rea A further optional safeguard removes isolated high-resolution “spur” spots by detecting large gaps in $1/d$ (or $q$) space and discarding spots beyond the gap. This is intended for macromolecular diffraction where edge-of-detector backgrounds can be extremely low. +A per-image **resolution estimate** is read off the same spot list, and it predicts the resolution the *merged* data will reach rather than the highest-resolution spot on this image. Each non-ice spot is weighted by $\sqrt{I}$ — the intensity is a summed photon count, so $\sqrt{I}$ is its Poisson significance — and the quantile is taken beyond which a fraction $f=0.30$ of that weight lies; the estimate is that resolution taken $2.25\times$ further in $1/d$, clamped so it can never beat the corner of the detector. The dataset value is the median over images. + +Two things make this work where reading off the highest-resolution spot does not. A quantile near the middle of the distribution measures the *shape* of the fall-off, which is the crystal's own decay $\exp(-B/2d^{2})$; the extreme end measures where detection stops, which moves with the exposure and with how many reflections the unit cell puts on a frame, and is a poor guide to the crystal. And merging averages many observations of each reflection, so intensities go on being measurable a fixed factor past the point where one image's spot finder still detects them — that factor is what the $2.25$ carries. Both constants are calibrated on rotation data against the resolution at which per-shell CC1/2 falls through 0.30. The estimate is a prediction, good to about 0.2 Å; it is not a measurement of what a run achieved, and nothing downstream is cut on it. + ### 3.4 Connected-component labeling (CCL) Strong pixels are grouped into connected components (adjacent strong pixels) using a CCL algorithm. Each component yields a candidate spot with: diff --git a/docs/HDF5.md b/docs/HDF5.md index 88a6e336..116b2e2d 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -280,7 +280,7 @@ In legacy/VDS mode these live in the data files and are linked/virtual-stacked i | `profileRadius` | Å⁻¹ | crystal profile radius | | `mosaicity` | deg | mosaicity estimate | | `bFactor` | Ų | per-image B-factor estimate | -| `resolutionEstimate` | Å | diffraction resolution estimate | +| `resolutionEstimate` | Å | resolution the merged data are predicted to reach, from this image's spots alone | | `integratedReflections` | | number of integrated reflections | | `bkgEstimate` | photons | mean background in the 3–5 Å resolution band | | `iceRingScore` | ratio | strongest hexagonal-ice ring intensity over the smooth radial background (1 = no ice) | diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index de5063db..a1b910c4 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -152,6 +152,11 @@ Sections, in order: `1. DATA SET`, `2. INDEXING`, `3. GEOMETRY POST-REFINEMENT` `4. SPACE GROUP DETERMINATION`, `5. SCALING AND MERGING`, `6. TWINNING`, `7. RADIATION DAMAGE`, `8. SWEEP QUALITY`, `9. WARNINGS`. +**`SPOT_RESOLUTION_ESTIMATE=`** in section 1 is how far the merged data are expected to reach, read +off the found spots alone — no lattice, no integration, no merge — so it is there on a run that never +merges, and on a run that does it can be read against `INCLUDE_RESOLUTION_RANGE` in section 5. It is a +prediction, good to about 0.2 Å on rotation data; nothing is cut on it. + **Which pass.** A rotation run integrates twice — once at the geometry in the input file, then again at the post-refined geometry — and can integrate a third time if a guard rejects the second pass. There is **one** report, for the pass that became the canonical output, and `PASS=` / diff --git a/image_analysis/spot_finding/SpotUtils.cpp b/image_analysis/spot_finding/SpotUtils.cpp index 641ebd26..04fe8966 100644 --- a/image_analysis/spot_finding/SpotUtils.cpp +++ b/image_analysis/spot_finding/SpotUtils.cpp @@ -114,22 +114,57 @@ void FilterSpuriousHighResolutionSpots(std::vector &spots, float thr } } -std::optional GetResolution(const std::vector &spots) { - std::vector resolutions; - resolutions.reserve(spots.size()); +namespace { + // Fraction of the image's weighted spot signal that is allowed to lie beyond the quantile read + // off below. A quantile near the middle of the distribution measures the shape of the fall-off, + // which is the crystal's own; the extreme end of it measures the detection threshold and how many + // reflections the unit cell puts on the frame, which are not. + constexpr float SPOT_RESOLUTION_TAIL_FRACTION = 0.30f; + // How much further in 1/d the merged data reach than that quantile. Merging averages many + // observations of each reflection, so intensities go on being measurable well past the point where + // one image's spot finder still detects them. Calibrated on rotation data against the resolution at + // which per-shell CC1/2 falls through 0.30. + constexpr float SPOT_RESOLUTION_MERGE_REACH = 2.25f; + // Fewer spots than this and the quantile is not a fall-off, it is a handful of points. + constexpr size_t SPOT_RESOLUTION_MIN_SPOTS = 4; +} +std::optional GetResolution(const std::vector &spots, float detector_d_min_A) { + // Each spot enters weighted by its own signal-to-noise. The intensity is a summed photon count, so + // it is Poisson and its significance is sqrt(I): that keeps a marginal high-resolution detection + // from counting for as much as a real reflection, without letting the handful of very strong + // low-resolution reflections - which say nothing about how far the crystal diffracts - decide the + // answer, as weighting by intensity itself would. + std::vector> spot_1_over_d2_weight; // (1/d^2, sqrt(intensity)) + spot_1_over_d2_weight.reserve(spots.size()); + float total_weight = 0.0f; for (const auto &spot: spots) { - if (!spot.ice_ring) - resolutions.push_back(spot.d_A); + if (spot.ice_ring || !(spot.d_A > 0.0f) || !(spot.intensity > 0.0f)) + continue; + const float weight = std::sqrt(spot.intensity); + spot_1_over_d2_weight.emplace_back(1.0f / (spot.d_A * spot.d_A), weight); + total_weight += weight; } - std::ranges::sort(resolutions); - - if (resolutions.size() < 4) + if (spot_1_over_d2_weight.size() < SPOT_RESOLUTION_MIN_SPOTS || !(total_weight > 0.0f)) return std::nullopt; - if (resolutions.size() < 20) - return resolutions[2]; - return resolutions[static_cast(resolutions.size() * 0.05)]; + + // Walk in from the highest-resolution spot until the tail fraction of the weight is behind us. + std::ranges::sort(spot_1_over_d2_weight, std::ranges::greater{}, + [](const std::pair &s) { return s.first; }); + float walked = 0.0f; + float one_over_d2 = spot_1_over_d2_weight.front().first; + for (const auto &[s, weight]: spot_1_over_d2_weight) { + walked += weight; + one_over_d2 = s; + if (walked >= SPOT_RESOLUTION_TAIL_FRACTION * total_weight) + break; + } + + const float d_A = 1.0f / (SPOT_RESOLUTION_MERGE_REACH * std::sqrt(one_over_d2)); + + // However far the crystal diffracts, no merge reaches past the corner of the detector. + return detector_d_min_A > 0.0f ? std::max(d_A, detector_d_min_A) : d_A; } void GenerateSpotPlot(DataMessage &msg, const std::vector &spots, float d_min_A) { @@ -197,7 +232,7 @@ void SpotAnalyze(const DiffractionExperiment &experiment, GenerateSpotPlot(output, spots_out, spot_d_min.value_or(0.0f) > 0 ? *spot_d_min : experiment.GetDetectorMaxResolution_A()); - output.resolution_estimate = GetResolution(spots_out); + output.resolution_estimate = GetResolution(spots_out, experiment.GetDetectorMaxResolution_A()); // One decision drives both: if indexing is to use the ice-band spots, the spot budget must not // throw them away before it gets the chance. diff --git a/image_analysis/spot_finding/SpotUtils.h b/image_analysis/spot_finding/SpotUtils.h index a6b79e42..f0b1aa3d 100644 --- a/image_analysis/spot_finding/SpotUtils.h +++ b/image_analysis/spot_finding/SpotUtils.h @@ -26,11 +26,20 @@ void MarkIceRings(std::vector &spots, float tolerance_q_recipA); void FilterSpotsByCount(std::vector &input, int64_t count, bool deprioritise_ice); void FilterSpuriousHighResolutionSpots(std::vector &spots, float threshold); -// Ignore high res. spots if there is a gap in (1/d) between two spots of dist_threshold (default: 0.25 A^-1) -// For less than 4 spots - don't provide resolution estimation -// For less than 20 spots - take the third best high resolution -// For 20 spots or more - take the 95% percentile best high resolution spot -std::optional GetResolution(const std::vector &spots); + +// How far the SCALED AND MERGED data from a run of images like this one are expected to reach, in A, +// from the spots of this one image alone. Nothing else - no lattice, no integration, no merge. +// +// The spots' intensity-weighted resolution distribution falls off at the crystal's own rate, and a +// quantile taken near the middle of that fall-off measures that rate. The far end of the distribution +// does not: the highest-resolution spot found says where DETECTION stops, which moves with the +// exposure and with how many reflections the unit cell puts on the frame. Merging then reaches a fixed +// factor further in 1/d than the quantile, because averaging many observations goes on measuring +// intensities that one image cannot detect. +// +// detector_d_min_A is the corner of the detector, which the answer is never allowed to beat; pass 0 to +// leave it unclamped. Returns nothing when the image has too few spots to have a fall-off at all. +std::optional GetResolution(const std::vector &spots, float detector_d_min_A = 0.0f); void SpotAnalyze(const DiffractionExperiment &experiment, const SpotFindingSettings &settings, diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 8b6fd1d4..6157dfce 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -78,6 +78,16 @@ std::string RenderResultReport(const std::string &output_prefix, << " The distance and beam centre above are the ones this result was integrated at, which on\n" << " a rotation run is the post-refined geometry rather than the values in the input file.\n"; + if (result.spot_resolution_estimate_A.has_value()) { + os << "\n"; + Key(os, "SPOT_RESOLUTION_ESTIMATE", fmt::format("{:.2f}", *result.spot_resolution_estimate_A)); + os << "\n" + << " How far the merged data are expected to reach, read off the found spots alone - no\n" + << " lattice, no integration, no merge. It is a prediction, good to about 0.2 A on the\n" + << " rotation data it was calibrated on, and it is not what the run achieved: compare it\n" + << " with INCLUDE_RESOLUTION_RANGE in section 5.\n"; + } + if (result.pass_count > 1) { os << "\n"; Key(os, "PASS", fmt::format("{} of {}", result.pass_number, result.pass_count)); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 50d89b81..d0cff4c9 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1925,6 +1925,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b end_msg.run_number = experiment_.GetRunNumber(); end_msg.run_name = experiment_.GetRunName(); end_msg.bkg_estimate = plots.GetBkgEstimate(); + result.spot_resolution_estimate_A = plots.GetResolutionEstimate(); end_msg.ice_ring_score = plots.GetIceRingScoreArray(); end_msg.ice_ring_score_mean = plots.GetIceRingScore(); end_msg.az_int_result["dataset"] = plots.GetAzIntProfile(); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index ed4a8307..e8fee939 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -140,6 +140,10 @@ struct ProcessResult { double frame_rate_hz = 0.0; double throughput_MBs = 0.0; std::optional indexing_rate; + // How far the merged data are expected to reach, predicted from spot finding alone (the median of + // the per-image estimates, see SpotUtils.h GetResolution). Independent of the merge, so it can be + // read against INCLUDE_RESOLUTION_RANGE - or on a run that never merges. + std::optional spot_resolution_estimate_A; std::optional consensus_cell; bool rotation_lattice_found = false; MeanProcessingTime mean_processing_time{}; diff --git a/tests/SpotUtilsTest.cpp b/tests/SpotUtilsTest.cpp index 54470895..17f70b70 100644 --- a/tests/SpotUtilsTest.cpp +++ b/tests/SpotUtilsTest.cpp @@ -27,3 +27,32 @@ TEST_CASE("FilterSpuriousHighResolutionSpots") { CHECK(spots[3].d_A == Catch::Approx(6.0)); } + +TEST_CASE("GetResolution") { + // Eleven equally strong spots at 1/d^2 = 0.1, 0.2, ... 1.1. Walking in from the highest-resolution + // one, four of the eleven are the first to carry 30% of the weight, so the quantile is the fourth + // spot in, 1/d^2 = 0.8. The estimate is that resolution taken 2.25x further in 1/d. + std::vector spots; + for (int i = 1; i <= 11; i++) + spots.push_back(SpotToSave{.intensity = 100.0f, .d_A = 1.0f / std::sqrt(0.1f * static_cast(i))}); + const auto d = GetResolution(spots); + REQUIRE(d.has_value()); + CHECK(*d == Catch::Approx(1.0 / (2.25 * std::sqrt(0.8))).epsilon(1e-4)); + + // The merged data cannot beat the corner of the detector. + CHECK(*GetResolution(spots, 2.0f) == Catch::Approx(2.0)); + + // Ice-flagged spots take no part, however strong they are. + std::vector with_ice = spots; + with_ice.push_back(SpotToSave{.intensity = 1e6f, .d_A = 0.5f, .ice_ring = true}); + CHECK(*GetResolution(with_ice) == Catch::Approx(*d)); + + // A weak high-resolution spot moves the answer far less than a strong one, which is the point of + // weighting by sqrt(I) rather than counting: the old order statistic would follow it entirely. + std::vector with_spur = spots; + with_spur.push_back(SpotToSave{.intensity = 1.0f, .d_A = 0.5f}); + CHECK(*GetResolution(with_spur) == Catch::Approx(*d).epsilon(0.02)); + + // Too few spots to have a fall-off at all. + CHECK_FALSE(GetResolution(std::vector(3)).has_value()); +}