diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index bc1e74032..5d70637ce 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -734,8 +734,9 @@ against the ice-free flanks beside them (ice in large crystallites); 1 is no ice **`ICE_RINGS_DETECTED=`** is `TRUE` where either passes its gate (1.5 and 2.0), and **`ICE_REFLECTIONS_ON_RINGS_PCT=`** is then the share of the integrated reflections set aside from scaling. **`ICE_RING_BAND_nn=`** (rotation only) is one line per band of a ring - its core, within -0.01 Å$^{-1}$ in $q$ of the ring, and its low- and high-$q$ flanks out to the flag width: the ring's -$d$ (Å), the band, what the merge did with its reflections (`MERGE`, `DOWN_WEIGHT`, `REJECT`), how +0.01 Å$^{-1}$ in $q$ of the ring, and on each side an inner flank out to 0.02 and an outer one out to the +flag width (a ring's shoulder reaches into the inner flank, so it is judged apart from the outer): the ring's +$d$ (Å), the band (`CORE`, `LOW_FLANK`, `LOW_OUTER_FLANK`, ...), what the merge did with its reflections (`MERGE`, `DOWN_WEIGHT`, `REJECT`), how many merged reflections it holds, their bias against the off-ring reflections around the ring in standard errors, and the ratio of their variance to those reflections' (both on intensities normalised by a local Wilson fit). A band is rejected where it is biased beyond 3 standard errors, diff --git a/image_analysis/scale_merge/IceRingVerdict.cpp b/image_analysis/scale_merge/IceRingVerdict.cpp index 91e878ce3..240bf43ce 100644 --- a/image_analysis/scale_merge/IceRingVerdict.cpp +++ b/image_analysis/scale_merge/IceRingVerdict.cpp @@ -29,6 +29,16 @@ namespace { double OneOverD2(float d) { return 1.0 / (static_cast(d) * d); } } +const char *IceRingZoneName(int zone) { + switch (zone) { + case -2: return "low-q outer flank"; + case -1: return "low-q flank"; + case 1: return "high-q flank"; + case 2: return "high-q outer flank"; + default: return "core"; + } +} + const char *IceRingActionName(IceRingAction a) { switch (a) { case IceRingAction::DownWeight: return "down-weight"; @@ -98,20 +108,27 @@ IceRingDecision DecideIceRings(const std::vector &merged) { // Each band's reflections are counted whether or not there is a baseline to judge them // against, so a band merged for want of one is still reported as such. - Moments zone[3]; - int in_zone[3] = {0, 0, 0}; + // Five bands: the core and, on each side, an inner and an outer flank. A ring's shoulder reaches + // past the core into the inner flank; judged together with the outer flank it was diluted into a + // merge (a 5x shoulder was kept on one crystal), so each flank is judged on its own. + auto zone_of = [](double dq) { + return dq < -ICE_VERDICT_FLANK_Q ? 0 : dq < -ICE_VERDICT_CORE_Q ? 1 + : dq <= ICE_VERDICT_CORE_Q ? 2 : dq <= ICE_VERDICT_FLANK_Q ? 3 : 4; + }; + Moments zone[5]; + int in_zone[5] = {0, 0, 0, 0, 0}; for (size_t i = 0; i < merged.size(); ++i) { if (ring_of[i] != k) continue; - const int z = dq_of[i] < -ICE_VERDICT_CORE_Q ? 0 : (dq_of[i] > ICE_VERDICT_CORE_Q ? 2 : 1); + const int z = zone_of(dq_of[i]); ++in_zone[z]; if (nb.n > 0) zone[z].Add(t_of(merged[i])); } - IceRingZoneVerdict v[3]; - for (int z = 0; z < 3; ++z) { + IceRingZoneVerdict v[5]; + for (int z = 0; z < 5; ++z) { v[z].ring_d = ICE_RING_RES_A[k]; - v[z].zone = z - 1; + v[z].zone = z - 2; v[z].n = in_zone[z]; v[z].n_neighbours = nb.n; if (zone[z].n < ICE_VERDICT_MIN_ZONE || nb.n < ICE_VERDICT_MIN_NEIGHBOURS) @@ -128,7 +145,7 @@ IceRingDecision DecideIceRings(const std::vector &merged) { for (size_t i = 0; i < merged.size(); ++i) { if (ring_of[i] != k) continue; - const IceRingZoneVerdict &zv = v[dq_of[i] < -ICE_VERDICT_CORE_Q ? 0 : (dq_of[i] > ICE_VERDICT_CORE_Q ? 2 : 1)]; + const IceRingZoneVerdict &zv = v[zone_of(dq_of[i])]; if (zv.action == IceRingAction::Reject) out.weight_factor[i] = 0.0; else if (zv.action == IceRingAction::DownWeight) { @@ -136,7 +153,7 @@ IceRingDecision DecideIceRings(const std::vector &merged) { const auto &r = merged[i]; const double I_exp = std::exp(intercept + slope * OneOverD2(r.d)); const double var = static_cast(r.sigma) * r.sigma; - const double excess = (zone[zv.zone + 1].Var() - nb.Var()) * (var + I_exp * I_exp); + const double excess = (zone[zv.zone + 2].Var() - nb.Var()) * (var + I_exp * I_exp); out.weight_factor[i] = var / (var + excess); } } diff --git a/image_analysis/scale_merge/IceRingVerdict.h b/image_analysis/scale_merge/IceRingVerdict.h index 2b1ccb7a3..d132bf649 100644 --- a/image_analysis/scale_merge/IceRingVerdict.h +++ b/image_analysis/scale_merge/IceRingVerdict.h @@ -17,7 +17,8 @@ // from a log-linear Wilson fit (WilsonLineThrough, no I/sigma cut on its shells) over the off-ring // reflections within ICE_VERDICT_WINDOW_S2 of the ring, flat at their mean where that cannot be fitted, and each // reflection is read as t = (I - I_exp) / sqrt(sigma^2 + I_exp^2) - unit scatter whether the reflection -// is signal- or noise-dominated. Per band (the core of the ring and each flank): +// is signal- or noise-dominated. Per band (the core of the ring, and an inner and an outer flank on each +// side - a ring's shoulder reaches into the inner flank, and judged with the outer one it was diluted): // * the band's mean t differs from its neighbours' by more than ICE_VERDICT_SIGMA standard errors: // biased - REJECT, since a larger sigma would not remove a bias; // * otherwise its variance of t exceeds the neighbours' by more than ICE_VERDICT_SIGMA standard errors @@ -29,7 +30,7 @@ enum class IceRingAction { Merge, DownWeight, Reject }; struct IceRingZoneVerdict { float ring_d = 0.0f; // the ice ring (A) - int zone = 0; // -1 low-q flank, 0 core, +1 high-q flank + int zone = 0; // -2 low-q outer flank, -1 low-q flank, 0 core, +1 high-q flank, +2 high-q outer flank int n = 0; // merged reflections in the band int n_neighbours = 0; // off-ring reflections they are compared with double bias_z = 0.0; // (mean t band - mean t neighbours) / its standard error @@ -47,6 +48,7 @@ struct IceRingDecision { constexpr double ICE_VERDICT_SIGMA = 3.0; constexpr double ICE_VERDICT_CORE_Q = 0.01; // |q - q_ring| below this is the core (A^-1) +constexpr double ICE_VERDICT_FLANK_Q = 0.02; // ... below this the inner flank, beyond it the outer constexpr double ICE_VERDICT_WINDOW_S2 = 0.02; // neighbours within this of the ring in 1/d^2 (A^-2) constexpr int ICE_VERDICT_MIN_ZONE = 20; constexpr int ICE_VERDICT_MIN_NEIGHBOURS = 50; @@ -55,3 +57,4 @@ constexpr int ICE_VERDICT_MIN_NEIGHBOURS = 50; IceRingDecision DecideIceRings(const std::vector &merged); const char *IceRingActionName(IceRingAction a); +const char *IceRingZoneName(int zone); diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index e5e0e8d24..48d8c027a 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -6272,7 +6272,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool for (const auto &z : ice.zones) logger->Info("Ice ring {:.3f} A {}: {} reflections against {} off the ring, bias {:+.1f} sigma, " "variance x{:.2f} ({:+.1f} sigma) - {}", z.ring_d, - z.zone < 0 ? "low-q flank" : (z.zone > 0 ? "high-q flank" : "core"), z.n, + IceRingZoneName(z.zone), z.n, z.n_neighbours, z.bias_z, z.var_ratio, z.var_ratio_z, IceRingActionName(z.action)); result.ice_ring_verdicts = ice.zones; } diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 3ef220169..7deff7b35 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -2269,7 +2269,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, const auto &z = result.ice_ring_verdicts[i]; Add(s, KeyText(fmt::format("ICE_RING_BAND_{:02}", i + 1).c_str(), fmt::format("{:.3f} {} {} {} {:+.1f} {:.2f}", z.ring_d, - z.zone < 0 ? "LOW_FLANK" : (z.zone > 0 ? "HIGH_FLANK" : "CORE"), + z.zone == -2 ? "LOW_OUTER_FLANK" : z.zone == -1 ? "LOW_FLANK" : z.zone == 0 ? "CORE" + : z.zone == 1 ? "HIGH_FLANK" : "HIGH_OUTER_FLANK", z.action == IceRingAction::Reject ? "REJECT" : (z.action == IceRingAction::DownWeight ? "DOWN_WEIGHT" : "MERGE"), z.n, z.bias_z, z.var_ratio))); @@ -2325,8 +2326,9 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, ? " The ring handling above has treated what it could." : " The images' radial profile and spots showed no ring, so nothing was set\n" " aside for it: check the ice-ring shells.") - : fmt::format(" Ice in the intensities: none. The reflections on the ice rings agree with those beside\n" - " them ({:+.1f} sigma over {} rings).", ii.contrast_z, ii.n_rings))); + : fmt::format(" Ice in the intensities: no ring system. The ring cores stand {:+.1f} sigma from the\n" + " reflections beside them over {} rings, inside the 5 sigma a ring system would show.", + ii.contrast_z, ii.n_rings))); Add(s, Prose(" ICE_INTENSITY_RING_nn: ring d (A), merged reflections in its core, their mean intensity\n" " over a Wilson line fitted to the off-ring reflections within 0.15 A^-1 of the ring, that\n" " excess in standard errors, and the same mean over the ring's flanks. ICE_INTENSITY_Z fits\n" @@ -3653,8 +3655,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, : ii.Detected() ? fmt::format("; the merged intensities on the rings stand {:.1f} sigma {} their " "neighbours", std::fabs(ii.contrast_z), ii.contrast_z > 0 ? "above" : "below") - : fmt::format("; the merged intensities on the rings agree with their neighbours " - "({:+.1f} sigma)", ii.contrast_z); + : fmt::format("; no ring system in the merged intensities (core-flank contrast " + "{:+.1f} sigma)", ii.contrast_z); checked_row("Ice", {PathologyCode::ICE_RINGS}, (result.ice_rings_detected ? fmt::format("rings DETECTED: the {:.1f}% of the reflections on them were left out of " "scaling{}", 100.0 * result.ice_flagged_fraction, diff --git a/tests/IceRingVerdictTest.cpp b/tests/IceRingVerdictTest.cpp index c4564b2bc..d0d1de3d2 100644 --- a/tests/IceRingVerdictTest.cpp +++ b/tests/IceRingVerdictTest.cpp @@ -60,8 +60,8 @@ TEST_CASE("IceRingVerdict_biased_core_is_rejected") { }); const auto d = DecideIceRings(merged); CHECK(ActionOf(d, 2.667f, 0) == IceRingAction::Reject); - CHECK(ActionOf(d, 2.667f, -1) == IceRingAction::Reject); - CHECK(ActionOf(d, 2.667f, 1) == IceRingAction::Reject); + for (int zone : {-2, -1, 1, 2}) + CHECK(ActionOf(d, 2.667f, zone) == IceRingAction::Reject); CHECK(ActionOf(d, 2.249f, 0) == IceRingAction::Merge); size_t rejected = 0; for (double f : d.weight_factor)