diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index a4063ba29..2ac99041c 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -3853,10 +3853,14 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) { // 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. + // that costs 0.0002 is many standard errors, and harmless - unless the information-weighted CC1/2 + // says the same, as significantly: then the merged intensities do not want the stretch either, and + // there is nothing its removal could cost. The floor protects a stretch the merge still learns + // from, not one whose frames carry nothing but noise (a sample that has left the beam). // 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(); + double healthy_null = -std::numeric_limits::infinity(); + double healthy_floor = healthy_null; { std::vector> healthy; for (int b = 0; b < n_batch; ++b) { @@ -3871,7 +3875,8 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) { 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); + healthy_null = median - DELTA_CC_SIGMA * 1.4826 * MedianOf(v); + healthy_floor = std::min(healthy_null, -DELTA_CC_MIN_HARM); } } // Ledger ranges closer together than one rocking event are one stretch. @@ -3885,10 +3890,22 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) { while (true) { // A stretch already wholly out touches no reflection and reads as not measured. const std::vector stretch_delta = measure_all(stretches); + unit_weight = false; + const std::vector weighted_delta = measure_all(stretches); + unit_weight = true; + for (size_t i = 0; i < stretches.size(); ++i) + if (stretch_delta[i].ok) + logger->Debug("delta-CC1/2 flagged stretch {}-{}: plain {:+.4f} +/- {:.4f} (z {:+.1f}), weighted {:+.4f} " + "(z {:+.1f}), healthy null {:+.4f}", stretches[i].first, stretches[i].second, + stretch_delta[i].delta_cc, stretch_delta[i].se, stretch_delta[i].delta_z / stretch_delta[i].se_z, + weighted_delta[i].delta_cc, weighted_delta[i].ok ? weighted_delta[i].delta_z / weighted_delta[i].se_z : 0.0, + healthy_null); 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 + const RangeDelta &r = stretch_delta[i], &w = weighted_delta[i]; + const bool weighted_agrees = w.ok && w.delta_z < -DELTA_CC_SIGMA * w.se_z; + if (r.ok && r.delta_z < -DELTA_CC_SIGMA * r.se_z + && (r.delta_cc < healthy_floor || (weighted_agrees && r.delta_cc < healthy_null)) && (worst < 0 || r.delta_z < stretch_delta[worst].delta_z)) worst = static_cast(i); }