diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp
index 37518848c..1c50c3626 100644
--- a/image_analysis/scale_merge/RotationScaleMerge.cpp
+++ b/image_analysis/scale_merge/RotationScaleMerge.cpp
@@ -5300,17 +5300,20 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
     // not rejected by either outlier test, inside the cut (judged on the group's d, as the merged
     // reflections are) - each with the sigma the merge weighted it by.
     if (export_scaled_fulls && full_stats && !for_search) {
-        result.scaled_fulls.clear();
-        for (int i = 0; i < n_full; ++i) {
+        // Which fulls, in order, and then each one's record on its own - on all threads.
+        const std::vector<int32_t> kept = KeptIndices(n_full, nthreads, [&](int i) {
             const int g = mf.group[i];
-            if (g < 0 || rejected_obs[i] || !std::isfinite(merged_I[g]))
-                continue;
-            if (effective_d_min && acc[g].d < *effective_d_min)
-                continue;
-            const Obs &o = fulls[i];
-            const float I_corr = o.I * o.corr;
-            result.scaled_fulls.push_back({o.h, o.k, o.l, I_corr, corrected_sigma(o, I_corr, o.sigma * o.corr)});
-        }
+            return g >= 0 && !rejected_obs[i] && std::isfinite(merged_I[g])
+                   && !(effective_d_min && acc[g].d < *effective_d_min);
+        });
+        result.scaled_fulls.resize(kept.size());
+        ParallelChunks(static_cast<int>(kept.size()), ThreadsForWork(kept.size(), nthreads), [&](int lo, int hi) {
+            for (int j = lo; j < hi; ++j) {
+                const Obs &o = fulls[kept[j]];
+                const float I_corr = o.I * o.corr;
+                result.scaled_fulls[j] = {o.h, o.k, o.l, I_corr, corrected_sigma(o, I_corr, o.sigma * o.corr)};
+            }
+        });
     }
 
     // Asymptotic I/sigma. ISa is by definition the I -> infinity limit of the signal-to-noise, i.e. the