diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp
index 51702499..ac2d201a 100644
--- a/image_analysis/scale_merge/Merge.cpp
+++ b/image_analysis/scale_merge/Merge.cpp
@@ -6,6 +6,7 @@
 #include <algorithm>
 #include <cmath>
 #include <limits>
+#include <random>
 #include <unordered_map>
 
 #include "../../common/ResolutionShells.h"
@@ -25,7 +26,13 @@ std::vector<uint8_t> CalcMergeMask(const DiffractionExperiment &x,
 std::vector<MergedReflection> MergeAll(const DiffractionExperiment &x,
                                        const std::vector<std::vector<Reflection> > &reflections,
                                        const std::vector<uint8_t> &merge_mask) {
-    if (!merge_mask.empty() && merge_mask.size() < reflections.size())
+
+    // To select images for half-datasets to calculate CC1/2, I use a random number generator with a fixed seed.
+    // This makes sure that images are selected randomly, but in a fully reproducible manner (at least for the same binary)
+    std::mt19937 rng(123456789u);
+    std::bernoulli_distribution half_dist(0.5);
+
+    if (!merge_mask.empty() && merge_mask.size() != reflections.size())
         throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Merge mask size mismatch");
 
     auto scaling_settings = x.GetScalingSettings();
@@ -48,8 +55,11 @@ std::vector<MergedReflection> MergeAll(const DiffractionExperiment &x,
     std::unordered_map<uint64_t, Accum> acc;
 
     for (int i = 0; i < reflections.size(); ++i) {
+        const int half = half_dist(rng);
         if (!merge_mask.empty() && merge_mask[i] == 0)
             continue;
+        if (reflections[i].empty())
+            continue;
         for (const auto &r: reflections[i]) {
             if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr))
                 continue;
@@ -80,7 +90,6 @@ std::vector<MergedReflection> MergeAll(const DiffractionExperiment &x,
             it->second.sum_wI += wI;
             it->second.sum_w += w;
 
-            const int half = ((it->second.n_half[0] + it->second.n_half[1]) % 2 == 0) ? 0 : 1;
             it->second.sum_wI_half[half] += wI;
             it->second.sum_w_half[half] += w;
             it->second.n_half[half]++;
@@ -125,7 +134,7 @@ MergeStatistics MergeStats(const DiffractionExperiment &x,
                            const std::vector<std::vector<Reflection> > &reflections,
                            const std::vector<uint8_t> &merge_mask) {
 
-    if (!merge_mask.empty() && merge_mask.size() < reflections.size())
+    if (!merge_mask.empty() && merge_mask.size() != reflections.size())
         throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Merge mask size mismatch");
 
     constexpr int n_shells = 10;