diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index c08cbf3b..b37f3237 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1099,6 +1099,43 @@ void RotationScaleMerge::SmoothG(std::vector &obs, std::vector &g, g[f] = g_smooth[f]; } +bool RotationScaleMerge::RejectCollapsedFullScales() { + // The Unity model leaves corr = 1/G, and every full of a frame carries that frame's G. + std::vector g_frame(n_frames, NAN); + for (const auto &o : fulls) + if (o.frame >= 0 && o.frame < n_frames && std::isfinite(o.corr) && o.corr > 0.0f) + g_frame[o.frame] = 1.0 / static_cast(o.corr); + + std::vector fitted; + fitted.reserve(g_frame.size()); + for (const double gf : g_frame) + if (std::isfinite(gf) && gf > 0.0) + fitted.push_back(gf); + if (fitted.empty()) + return false; + const size_t mid = fitted.size() / 2; + std::nth_element(fitted.begin(), fitted.begin() + mid, fitted.end()); + const double g_floor = fitted[mid] * MIN_CREDIBLE_SCALE_RATIO; + + int n_rejected = 0; + double worst = 1.0; + for (auto &o : fulls) { + if (o.frame < 0 || o.frame >= n_frames) + continue; + const double gf = g_frame[o.frame]; + if (!std::isfinite(gf) || gf >= g_floor) + continue; + worst = std::max(worst, 1.0 / gf); + o.corr = 1.0f; // as if the frame's scale had never been fitted + ++n_rejected; + } + if (n_rejected > 0) + logger.Warning("Rejected the fitted scale of {} full(s) whose frame scaled to less than {:.0f}x " + "below the run median (worst {:.0f}x amplification); those frames merge unscaled", + n_rejected, 1.0 / MIN_CREDIBLE_SCALE_RATIO, worst / fitted[mid]); + return n_rejected > 0; +} + void RotationScaleMerge::Combine() { fulls.clear(); g_full.assign(n_frames, 1.0); @@ -2028,6 +2065,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, } logger.Info("Scaled fulls (XDS order, Unity model)"); } + const bool rejected_full_scales = scale_fulls && RejectCollapsedFullScales(); // --- 4b. Optional correction surfaces (decay = resolution x time; absorption = goniometer-frame // diffracted-beam direction), each an alternating multiplicative fit of the fulls' corr against @@ -2051,7 +2089,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, #ifdef JFJOCH_USE_CUDA // The corrections mutate the host fulls' corr; when the merge runs on the resident (GPU) fulls, push // the corrected corr back to the device so the merge reads it. - if (corrections && combined_on_gpu && scaled_fulls_on_gpu) { + if ((corrections || rejected_full_scales) && combined_on_gpu && scaled_fulls_on_gpu) { std::vector fcorr(fulls.size()); for (size_t i = 0; i < fulls.size(); ++i) fcorr[i] = fulls[i].corr; gpu_->SetFullsCorr(fcorr.data()); diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 1ca5c5a5..2dd3be38 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -198,6 +198,14 @@ private: void Combine(); // partials -> fulls (CPU) + // Undo the scale on fulls whose frame's scale collapsed toward zero. The fulls are scaled with the + // Unity model, so their corr IS 1/G and a collapsed G multiplies every intensity on that frame + // without bound. This stage refits G from scratch with no smoothing to fall back on, so a rejected + // frame simply keeps the unscaled corr the combine gave it. Reads the host fulls, so it covers the + // CPU and GPU scaling paths alike. Returns true if anything was rejected (the caller then has to + // push the corrected corr back to the device). + bool RejectCollapsedFullScales(); + // Post-scale-fulls correction surfaces, each an alternating multiplicative fit of the host fulls' corr // against the merged reference (cheap host loops; the corrected corr is re-uploaded to the resident // fulls afterwards). Each is cross-validated (fit even frames, keep only if held-out odd equivalents