diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 3a53277c2..a644b54a9 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -5439,16 +5439,29 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st } #endif if (!scaled_on_gpu) { + // RunScalingLoop rescales after every iteration before anything reads corr, so the iteration's + // UpdateCorr is folded into that pass over the partials: the G and the fitted frames as the + // iteration left them, then the rescale on top - the same two roundings, one pass instead of two. + std::vector g_fitted; + std::vector frame_fitted; partial_loop = RunScalingLoop("partials", g_partial, frame_count, smooth_window, /*release_restraint=*/false, [&] { ReducePartialGroupMeans(n_groups, partial_mean); FitPerFrameG(partials, frame_start, frame_count, partial_mean, /*unity=*/false, g_partial); - UpdateCorr(partials, g_partial, frame_scaled_scratch); + g_fitted = g_partial; + frame_fitted = frame_scaled_scratch; }, [&](const std::vector &apply, const std::vector &ratio) { ParallelChunks(static_cast(partials.size()), ThreadsForWork(partials.size(), nthreads), [&](int lo, int hi) { for (int i = lo; i < hi; ++i) { Obs &o = partials[i]; + if (frame_fitted[o.frame]) { // UpdateCorr + const double denom = static_cast(o.partiality) * g_fitted[o.frame]; + if (std::isfinite(o.prescaling_corr) && std::isfinite(denom) && denom > 0.0) + o.corr = static_cast(o.prescaling_corr / denom); + else + o.corr = NAN; + } if (apply[o.frame] && std::isfinite(o.corr)) o.corr = static_cast(o.corr * ratio[o.frame]); }