diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index f9f16dad5..ce7262e5c 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -604,7 +604,18 @@ per-image channels call normal, and what is left has to carry the verdict again that cannot be taken without healthy frames holds more than one thing and is not removed at all. A rejected stretch is never narrower than one rocking event, because the partials of one event are combined into the same intensities and inside it no frame can be judged apart from its neighbours. Never more than -a quarter of the sweep is removed. +a quarter of the sweep is removed this way. + +A stretch the per-image channels flagged — a range of this table, with neighbours closer than one +rocking event joined — is then judged once more, **taken whole**, and here neither the quarter-sweep +limit nor the run's own batches apply: when the crystal dies half-way, or a second crystal or a slipped +orientation takes over the second half of the sweep, the defect *is* most of the run, and the run's own +batches are no reference. The stretch is removed when keeping it lowers the **plain** CC1/2 (every +reflection at weight 1, as the CC1/2 the report prints) by several standard errors, by at least 0.01, +and by more than the spread of the batches neither channel flags (the same robust floor as above, +measured on those batches alone; with fewer than three of them nothing is removed). The worst stretch +goes first and the rest are measured again against what is left. Weak stretches that agree with the +rest of the run — a crystal that leaves the beam and comes back — cost nothing on this measure and stay. **What ΔCC1/2 cannot do**, because the report must not imply otherwise: diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index e3a838a2e..42133bbbe 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -3274,6 +3274,9 @@ namespace { constexpr double DELTA_CC_MAX_REJECTED = 0.25; // ... and never more than this fraction of the sweep constexpr int DELTA_CC_MIN_REFLECTIONS = 50; // fewer than this and the pair of CC1/2 means nothing constexpr int DELTA_CC_MIN_NULL_BATCHES = 8; // fewer measured batches than this and there is no null + // A stretch the per-image channels flagged, taken whole (the end of MeasureBatchDeltaCCHalf): + constexpr double DELTA_CC_MIN_HARM = 0.01; // the least cost in CC1/2 that is worth removing + constexpr int DELTA_CC_MIN_HEALTHY_BATCHES = 3; // fewer unflagged batches than this and there is no null // Median, and the median absolute deviation scaled to a standard deviation. Both destroy the order // of `v`, which no caller needs. @@ -3515,6 +3518,8 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) { bool ok = false; double delta_cc = 0.0, delta_z = 0.0, se = 0.0, se_z = 0.0; }; + // Every reflection at weight 1 instead of the information it carries: the plain CC1/2 (see below). + bool unit_weight = false; // Add the frames [f0, f1] still in to a range's contribution. Frame by frame in frame order, so a // longer range that starts on the same frame is the shorter one carried on: the same sums, added // in the same order. @@ -3546,9 +3551,9 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) { if (t.n >= 2 && n_without >= 2) { double mean, eps; ReflectionMoments(t.n, t.s1, t.s2, mean, eps); - with.Add(mean, eps, t.w / t.n); + with.Add(mean, eps, unit_weight ? 1.0 : t.w / t.n); ReflectionMoments(n_without, t.s1 - d.s1, t.s2 - d.s2, mean, eps); - without.Add(mean, eps, (t.w - d.w) / n_without); + without.Add(mean, eps, unit_weight ? 1.0 : (t.w - d.w) / n_without); } } const double cc_with = with.CC(), cc_without = without.CC(); @@ -3814,6 +3819,75 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) { candidate_delta[open_index[i]] = open_delta[i]; } + // Then the stretches the per-image channels flagged (the ledger MeasureSweepQuality drew), each taken + // whole. The loop above cannot remove a stretch that is most of the sweep: its candidates stop at + // the cap, and its null is the run's own batches, which is meaningless when half of them are the + // defect - a crystal dead after 40% of a room-temperature sweep, a second crystal or a slipped + // orientation taking over the second half. Here the evidence that something is wrong is already + // there, from channels that owe nothing to the merge, and the stretch is known whole, so neither + // the cap nor the run-wide null is needed: + // * the cost is read on the plain CC1/2, every reflection at weight 1 - the statistic the run + // reports - and not on the information-weighted one, which by construction cannot see a weak + // stretch whose noise reaches the reported CC1/2; + // * the null is the batches nothing flags - every frame normal by both channels and still in - + // measured the same way, DELTA_CC_SIGMA robust standard deviations below their median; with + // fewer than DELTA_CC_MIN_HEALTHY_BATCHES of them there is no null and nothing is removed; + // * and the cost has to be at least DELTA_CC_MIN_HARM in CC: on a crystal merging at 0.999 a tail + // that costs 0.0002 is many standard errors, and harmless. + // The worst stretch goes first and the rest are measured again against what is left, as above. + unit_weight = true; + double healthy_floor = -std::numeric_limits::infinity(); + { + std::vector> healthy; + for (int b = 0; b < n_batch; ++b) { + bool normal = true; + for (int f = batch_first(b); f <= batch_last(b); ++f) + normal = normal && !frame_rejected[f] && normal_frame(f); + if (normal) healthy.emplace_back(batch_first(b), batch_last(b)); + } + std::vector v; + for (const RangeDelta &h : measure_all(healthy)) + if (h.ok) v.push_back(h.delta_cc); + if (static_cast(v.size()) >= DELTA_CC_MIN_HEALTHY_BATCHES) { + const double median = MedianOf(v); + for (double &d : v) d = std::fabs(d - median); + healthy_floor = std::min(median - DELTA_CC_SIGMA * 1.4826 * MedianOf(v), -DELTA_CC_MIN_HARM); + } + } + // Ledger ranges closer together than one rocking event are one stretch. + std::vector> stretches; + for (const auto &r : sweep_quality.ranges) { + if (!stretches.empty() && r.first_image - stretches.back().second - 1 < min_frames) + stretches.back().second = std::max(stretches.back().second, r.last_image); + else + stretches.emplace_back(r.first_image, r.last_image); + } + while (true) { + // A stretch already wholly out touches no reflection and reads as not measured. + const std::vector stretch_delta = measure_all(stretches); + int worst = -1; + for (size_t i = 0; i < stretches.size(); ++i) { + const RangeDelta &r = stretch_delta[i]; + if (r.ok && r.delta_z < -DELTA_CC_SIGMA * r.se_z && r.delta_cc < healthy_floor + && (worst < 0 || r.delta_z < stretch_delta[worst].delta_z)) + worst = static_cast(i); + } + if (worst < 0) + break; + const auto [f0, f1] = stretches[worst]; + for (int f = f0; f <= f1; ++f) { + n_rejected += frame_rejected[f] ? 0 : 1; + frame_rejected[f] = 1; + } + rejected_ranges.emplace_back(f0, f1); + range_cost.emplace_back(stretch_delta[worst].delta_cc, stretch_delta[worst].se); + logger->Info("delta-CC1/2: frames {}-{}, flagged by the per-image channels, removed whole - keeping " + "them lowered the plain CC1/2 by {:.4f} +/- {:.4f}", f0, f1, -stretch_delta[worst].delta_cc, + stretch_delta[worst].se); + build_totals(); + } + unit_weight = false; + // The ledger: one range per contiguous stretch of one disposition, with what keeping it costs. SplitSweepRanges(); FillSweepRanges(/*assign_reasons=*/false);