diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index a1dac33e0..ecbec9c6e 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,7 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* The spot-finding resolution estimate is no longer capped at the corner of the detector, so a crystal that diffracts past the edge is reported as reaching past it, and a run where the estimate is finer than what was merged is a detector-limited run. * Spot finding no longer throws away a whole image when it holds many strong pixels: the limit follows the detector (one pixel in 64) instead of standing at the 65535 that suited a 4-megapixel detector, which a strongly diffracting crystal on an 18-megapixel one passes on its best frames. * The connected-component search is linear in the strong pixels rather than quadratic in how many of them a detector line holds; a fully lit image is labelled in 0.16 s instead of 76 s, and the spots it finds are unchanged. * In `jfjoch_viewer`, Alt and the mouse wheel step through the dataset one image at a time. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index b53529b9e..29650cc68 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -389,7 +389,7 @@ is kept; the frame is then integrated once at that min-pix. The fraction factor ### 3.6 Predicting the resolution the merged data will reach -A per-image **resolution estimate** is read off the finished spot list. It predicts how far the *merged* data will reach, not how far the furthest spot on this image lies. Each non-ice spot is weighted by $\sqrt{I}$ — the intensity is a summed photon count, so $\sqrt{I}$ is its Poisson significance — the $1/d^2$ is found beyond which a fraction $f=0.30$ of that weight lies, and 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. +A per-image **resolution estimate** is read off the finished spot list. It predicts how far the *merged* data will reach, not how far the furthest spot on this image lies. Each non-ice spot is weighted by $\sqrt{I}$ — the intensity is a summed photon count, so $\sqrt{I}$ is its Poisson significance — the $1/d^2$ is found beyond which a fraction $f=0.30$ of that weight lies, and the estimate is that resolution taken $2.25\times$ further in $1/d$. It is deliberately **not** limited to what the detector records: the quantile sits in the middle of the fall-off, well inside the recorded range, so it goes on measuring the crystal where the detector stops before the diffraction does, and on such a run it reads finer than the detector corner. The dataset value is the median over images. Both constants carry a mechanism. A quantile from the middle of the distribution measures the *shape* of the fall-off, which is the crystal's own $\exp(-B/2d^{2})$, where the extreme end of it measures where detection stops — a threshold that moves with the exposure and with how many reflections the unit cell puts on a frame. And merging averages many observations of each reflection, so intensities go on being measurable a fixed factor in $1/d$ past the point at which one image's spot finder still detects them; that factor is the $2.25$. Both are calibrated on rotation data against the resolution at which per-shell CC1/2 falls through 0.30, and the estimate is good to about 0.2 Å there. It is a prediction and not a measurement of what a run achieved: nothing downstream is cut on it, and it is reported alone (rugnux `SPOT_RESOLUTION_ESTIMATE`, and per image in the stream, the plots and HDF5). diff --git a/image_analysis/spot_finding/SpotUtils.cpp b/image_analysis/spot_finding/SpotUtils.cpp index 04fe89661..433eb6610 100644 --- a/image_analysis/spot_finding/SpotUtils.cpp +++ b/image_analysis/spot_finding/SpotUtils.cpp @@ -129,7 +129,7 @@ namespace { constexpr size_t SPOT_RESOLUTION_MIN_SPOTS = 4; } -std::optional GetResolution(const std::vector &spots, float detector_d_min_A) { +std::optional GetResolution(const std::vector &spots) { // 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 @@ -161,10 +161,10 @@ std::optional GetResolution(const std::vector &spots, float d 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; + // Not clamped at the corner of the detector. The quantile is read from the middle of the + // fall-off, so it still measures the crystal where the detector cuts that fall-off short; + // clamping reported where the detector stops instead, which is the one thing this is not for. + return 1.0f / (SPOT_RESOLUTION_MERGE_REACH * std::sqrt(one_over_d2)); } void GenerateSpotPlot(DataMessage &msg, const std::vector &spots, float d_min_A) { @@ -232,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, experiment.GetDetectorMaxResolution_A()); + output.resolution_estimate = GetResolution(spots_out); // 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 f0b1aa3d3..1eaa5f03e 100644 --- a/image_analysis/spot_finding/SpotUtils.h +++ b/image_analysis/spot_finding/SpotUtils.h @@ -37,9 +37,11 @@ void FilterSpuriousHighResolutionSpots(std::vector &spots, float thr // 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); +// The answer is deliberately NOT limited to what this detector records. The quantile sits in the +// middle of the fall-off, well inside the recorded range, so it goes on measuring the crystal when +// the detector stops before the diffraction does - which is the case the number is most wanted for. +// Returns nothing when the image has too few spots to have a fall-off at all. +std::optional GetResolution(const std::vector &spots); void SpotAnalyze(const DiffractionExperiment &experiment, const SpotFindingSettings &settings, diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 8b6ef13fb..876121f8d 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -102,7 +102,9 @@ std::string RenderResultReport(const std::string &output_prefix, << " 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"; + << " with INCLUDE_RESOLUTION_RANGE in section 5. It is not limited to what this detector\n" + << " records: where it reads finer than the high-resolution end of that range, the crystal\n" + << " diffracts past the corner and the run is detector-limited.\n"; } if (result.pass_count > 1) { diff --git a/tests/SpotUtilsTest.cpp b/tests/SpotUtilsTest.cpp index 17f70b701..0a804c9f0 100644 --- a/tests/SpotUtilsTest.cpp +++ b/tests/SpotUtilsTest.cpp @@ -39,8 +39,11 @@ TEST_CASE("GetResolution") { 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)); + // The answer is not limited to what a detector records. Keeping only the five spots a detector + // reaching 1/d^2 = 0.5 would have recorded leaves the quantile at 0.4, and the estimate still + // extrapolates 2.25x past it instead of stopping at the cut. + const std::vector cut(spots.begin(), spots.begin() + 5); + CHECK(*GetResolution(cut) == Catch::Approx(1.0 / (2.25 * std::sqrt(0.4))).epsilon(1e-4)); // Ice-flagged spots take no part, however strong they are. std::vector with_ice = spots;