diff --git a/image_analysis/UpdateReflectionResolution.cpp b/image_analysis/UpdateReflectionResolution.cpp
index 456d56f3..6cf1bf13 100644
--- a/image_analysis/UpdateReflectionResolution.cpp
+++ b/image_analysis/UpdateReflectionResolution.cpp
@@ -2,9 +2,11 @@
 // SPDX-License-Identifier: GPL-3.0-only
 
 #include "UpdateReflectionResolution.h"
+
+#include "IntegrationOutcome.h"
 #include "../common/CrystalLattice.h"
 
-ResolutionStats UpdateReflectionResolution(const UnitCell &cell, std::vector<std::vector<Reflection> > &reflections) {
+ResolutionStats UpdateReflectionResolution(const UnitCell &cell, std::vector<IntegrationOutcome > &reflections) {
     ResolutionStats ret;
 
     CrystalLattice latt(cell);
@@ -13,12 +15,12 @@ ResolutionStats UpdateReflectionResolution(const UnitCell &cell, std::vector<std
     const auto cstar = latt.Cstar();
 
     for (auto &image: reflections) {
-        if (!image.empty()) {
+        if (!image.reflections.empty()) {
             ret.n_images++;
-            ret.n_reflections += image.size();
+            ret.n_reflections += image.reflections.size();
         }
 
-        for (auto &r: image) {
+        for (auto &r: image.reflections) {
             Coord q = r.h * astar + r.k * bstar + r.l * cstar;
             auto qlen = q.Length();
             if (qlen > 1e-6) {
diff --git a/image_analysis/UpdateReflectionResolution.h b/image_analysis/UpdateReflectionResolution.h
index 53c26e04..46f909e1 100644
--- a/image_analysis/UpdateReflectionResolution.h
+++ b/image_analysis/UpdateReflectionResolution.h
@@ -6,6 +6,7 @@
 #include <vector>
 #include "../common/Reflection.h"
 #include "../common/UnitCell.h"
+#include "IntegrationOutcome.h"
 
 struct ResolutionStats {
     float d_low = std::numeric_limits<float>::max();
@@ -14,5 +15,4 @@ struct ResolutionStats {
     int n_images = 0;
 };
 
-// Returns the highest (best) resolution in Angstrom observed by reflection
-ResolutionStats UpdateReflectionResolution(const UnitCell &cell, std::vector<std::vector<Reflection> > &reflections);
+ResolutionStats UpdateReflectionResolution(const UnitCell &cell, std::vector<IntegrationOutcome> &reflections);
diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp
index 418383e4..d455129e 100644
--- a/image_analysis/scale_merge/Merge.cpp
+++ b/image_analysis/scale_merge/Merge.cpp
@@ -17,60 +17,14 @@
 #include "../../common/ResolutionShells.h"
 #include "HKLKey.h"
 
-size_t CalcMergeMaskCC(const DiffractionExperiment &x,
-                   const std::vector<float> &scale_cc,
-                   std::vector<uint8_t> &result_mask) {
-
-    if (scale_cc.size() != result_mask.size())
-        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Scale CC size mismatch");
-
-    size_t n_rejected = 0;
-    auto min_cc = x.GetScalingSettings().GetMinCCForImage();
-    if (min_cc > 0.0) {
-        for (int i = 0; i < result_mask.size(); ++i) {
-            if (result_mask[i] && (scale_cc[i] < min_cc)) {
-                result_mask[i] = 0;
-                ++n_rejected;
-            }
-        }
-    }
-    return n_rejected;
-}
-
-size_t CalcMergeMaskUnitCell(const DiffractionExperiment &x,
-                             const UnitCell &ref_cell,
-                             const std::vector<std::optional<UnitCell> > &unit_cells,
-                             std::vector<uint8_t> &mask) {
-    if (mask.size() != unit_cells.size())
-        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Mismatch in vector size");
-
-    size_t n_rejected = 0;
-
-        const auto dtol = x.GetIndexingSettings().GetUnitCellDistTolerance();
-        const auto atol = x.GetIndexingSettings().GetUnitCellAngleTolerance_deg();
-    for (size_t i = 0; i < unit_cells.size(); ++i) {
-        if (mask[i]) {
-            // Counter should not include trivial case (non-indexed image or weird things)
-            // It should only count where indexing did happen, but cell was too far from reference
-            if (!unit_cells[i] || !unit_cells[i]->is_finite())
-                mask[i] = 0;
-            else if (!unit_cells[i]->is_close(ref_cell, dtol, atol)) {
-                mask[i] = 0;
-                ++n_rejected;
-            }
-        }
-    }
-
-    return n_rejected;
-}
-
 MergeOnTheFly::MergeOnTheFly(const DiffractionExperiment &x)
-    : scaling_settings(x.GetScalingSettings()),
+    : space_group_number(x.GetSpaceGroupNumber().value_or(1)),
+      scaling_settings(x.GetScalingSettings()),
       indexing_settings(x.GetIndexingSettings()),
       high_resolution_limit(scaling_settings.GetHighResolutionLimit_A()),
       image_cc_limit(scaling_settings.GetMinCCForImage()),
       min_partiality(scaling_settings.GetMinPartiality()),
-      generator(scaling_settings.GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)) {
+      generator(scaling_settings.GetMergeFriedel(), space_group_number) {
 }
 
 MergeOnTheFly &MergeOnTheFly::ReferenceCell(const std::optional<UnitCell> &cell) {
@@ -78,8 +32,14 @@ MergeOnTheFly &MergeOnTheFly::ReferenceCell(const std::optional<UnitCell> &cell)
     return *this;
 }
 
-void MergeOnTheFly::ProcessImage_i(const std::vector<Reflection> &v, int half) {
-    for (const auto &r: v) {
+void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, bool cc_mask) {
+    std::unique_lock ul(merged_mutex);
+    const int half = half_dist(rng);
+
+    if (Mask(outcome, cc_mask))
+        return;
+
+    for (const auto &r: outcome.reflections) {
         if (generator.IsSystematicallyAbsent(r))
             continue;
 
@@ -121,52 +81,22 @@ void MergeOnTheFly::ProcessImage_i(const std::vector<Reflection> &v, int half) {
     }
 }
 
-void MergeOnTheFly::AddImage(const std::vector<Reflection> &v, double image_cc, const UnitCell &cell) {
-    std::unique_lock ul(merged_mutex);
-    const int half = half_dist(rng);
-
-    if (image_cc_limit && image_cc < image_cc_limit)
-        return;
-
-    if (reference_cell
-        && !cell.is_close(*reference_cell,
+bool MergeOnTheFly::Mask(const IntegrationOutcome &outcome, bool cc_mask) {
+    if (reference_cell) {
+        auto cell = outcome.latt.GetUnitCell();
+        if (!cell.is_close(*reference_cell,
                           indexing_settings.GetUnitCellDistTolerance(),
                           indexing_settings.GetUnitCellAngleTolerance_deg()))
-        return;
-
-    ProcessImage_i(v, half);
-}
-
-void MergeOnTheFly::AddAll(const std::vector<std::vector<Reflection> > &reflections,
-                           const std::vector<uint8_t> &merge_mask) {
-    if (!merge_mask.empty() && merge_mask.size() != reflections.size())
-        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Merge mask size mismatch");
-
-    std::unique_lock ul(merged_mutex);
-    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;
-        ProcessImage_i(reflections[i], half);
+            return true;
     }
-}
 
-void MergeOnTheFly::AddAll(const std::vector<IntegrationOutcome> &integration_outcome,
-                           const std::vector<uint8_t> &merge_mask) {
-    if (!merge_mask.empty() && merge_mask.size() != integration_outcome.size())
-        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Merge mask size mismatch");
-
-    std::unique_lock ul(merged_mutex);
-    for (int i = 0; i < integration_outcome.size(); ++i) {
-        const int half = half_dist(rng);
-        if (!merge_mask.empty() && merge_mask[i] == 0)
-            continue;
-        if (integration_outcome[i].reflections.empty())
-            continue;
-        ProcessImage_i(integration_outcome[i].reflections, half);
+    if (cc_mask && image_cc_limit) {
+        if (!outcome.image_scale_cc
+            || std::isnan(outcome.image_scale_cc.value())
+            || outcome.image_scale_cc.value() < image_cc_limit.value())
+            return true;
     }
+    return false;
 }
 
 std::vector<MergedReflection> MergeOnTheFly::ExportReflections() {
@@ -242,18 +172,11 @@ std::vector<MergedReflection> MergeOnTheFly::ExportReflections() {
 }
 
 std::vector<MergedReflection> MergeAll(const DiffractionExperiment &x,
-                                       const std::vector<std::vector<Reflection> > &reflections,
-                                       const std::vector<uint8_t> &merge_mask) {
+                                       const std::vector<IntegrationOutcome> &integration_outcome,
+                                       bool mask) {
     MergeOnTheFly merge(x);
-    merge.AddAll(reflections, merge_mask);
-    return merge.ExportReflections();
-}
-
-std::vector<MergedReflection> MergeAll(const DiffractionExperiment &x,
-                                       const std::vector<IntegrationOutcome > &integration_outcome,
-                                       const std::vector<uint8_t> &merge_mask) {
-    MergeOnTheFly merge(x);
-    merge.AddAll(integration_outcome, merge_mask);
+    for (const auto &outcome: integration_outcome)
+        merge.AddImage(outcome, mask);
     return merge.ExportReflections();
 }
 
@@ -269,14 +192,17 @@ struct ShellAccum {
     CorrelationCoefficient cc_ref;
 };
 
-void CalcPossibleReflections(const DiffractionExperiment &x,
+void CalcPossibleReflections(int space_group_number ,
                              const UnitCell &cell,
                              double d_min,
                              double d_max,
                              const ResolutionShells &shells,
                              std::vector<ShellAccum> &acc) {
     gemmi::UnitCell gemmi_cell = cell;
-    const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(x.GetSpaceGroupNumber().value_or(1));
+    const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number);
+    if (sg == nullptr)
+        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+            "Invalid space group number " + std::to_string(space_group_number));
 
     // Generate unique reflections
     std::vector<gemmi::Miller> possible_hkls = gemmi::make_miller_vector(gemmi_cell, sg, d_min, d_max, true);
@@ -301,21 +227,13 @@ void CalcPossibleReflections(const DiffractionExperiment &x,
 }
 
 
-MergeStatistics MergeStats(const DiffractionExperiment &x,
-                           const std::vector<MergedReflection> &merged,
+MergeStatistics MergeOnTheFly::MergeStats(const std::vector<MergedReflection> &merged,
                            const std::vector<IntegrationOutcome > &integration_outcome,
-                           const UnitCell &cell,
-                           const std::vector<uint8_t> &merge_mask,
                            const std::vector<MergedReflection> &reference) {
-    if (!merge_mask.empty() && merge_mask.size() != integration_outcome.size())
-        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Merge mask size mismatch");
 
     constexpr int n_shells = 10;
 
-    auto min_partiality = x.GetScalingSettings().GetMinPartiality();
-    auto d_min_limit_A = x.GetScalingSettings().GetHighResolutionLimit_A();
-    auto scaling_settings = x.GetScalingSettings();
-    HKLKeyGenerator key_generator(scaling_settings.GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1));
+    auto d_min_limit_A = scaling_settings.GetHighResolutionLimit_A();
 
     std::unordered_map<uint64_t, float> reference_intensities;
     if (!reference.empty()) {
@@ -324,7 +242,7 @@ MergeStatistics MergeStats(const DiffractionExperiment &x,
             if (!std::isfinite(r.I))
                 continue;
 
-            const auto hkl = key_generator(r);
+            const auto hkl = generator(r);
             reference_intensities[hkl.pack()] = r.I;
         }
     }
@@ -355,7 +273,9 @@ MergeStatistics MergeStats(const DiffractionExperiment &x,
 
     std::vector<ShellAccum> acc(n_shells);
 
-    CalcPossibleReflections(x, cell, d_min_pad, d_max_pad, shells, acc);
+    if (reference_cell.has_value())
+        CalcPossibleReflections(space_group_number, reference_cell.value(),
+            d_min_pad, d_max_pad, shells, acc);
 
     CorrelationCoefficient cc_half_overall;
     CorrelationCoefficient cc_ref_overall;
@@ -373,7 +293,7 @@ MergeStatistics MergeStats(const DiffractionExperiment &x,
                 ++acc[s].n_i_over_sigma;
 
                 if (!reference_intensities.empty()) {
-                    const auto hkl = key_generator(m);
+                    const auto hkl = generator(m);
                     const auto ref_it = reference_intensities.find(hkl.pack());
                     if (ref_it != reference_intensities.end() && std::isfinite(ref_it->second)) {
                         acc[s].cc_ref.Add(m.I, ref_it->second);
@@ -391,11 +311,11 @@ MergeStatistics MergeStats(const DiffractionExperiment &x,
     }
 
     for (int i = 0; i < integration_outcome.size(); ++i) {
-        if (!merge_mask.empty() && merge_mask[i] == 0)
+        if (Mask(integration_outcome[i], true))
             continue;
 
         for (const auto &r: integration_outcome[i].reflections) {
-            if (key_generator.IsSystematicallyAbsent(r))
+            if (generator.IsSystematicallyAbsent(r))
                 continue;
             if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr))
                 continue;
@@ -467,173 +387,6 @@ MergeStatistics MergeStats(const DiffractionExperiment &x,
     return out;
 }
 
-MergeStatistics MergeStats(const DiffractionExperiment &x,
-                           const std::vector<MergedReflection> &merged,
-                           const std::vector<std::vector<Reflection> > &reflections,
-                           const UnitCell &cell,
-                           const std::vector<uint8_t> &merge_mask,
-                           const std::vector<MergedReflection> &reference) {
-    if (!merge_mask.empty() && merge_mask.size() != reflections.size())
-        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Merge mask size mismatch");
-
-    constexpr int n_shells = 10;
-
-    auto min_partiality = x.GetScalingSettings().GetMinPartiality();
-    auto d_min_limit_A = x.GetScalingSettings().GetHighResolutionLimit_A();
-    auto scaling_settings = x.GetScalingSettings();
-    HKLKeyGenerator key_generator(scaling_settings.GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1));
-
-    std::unordered_map<uint64_t, float> reference_intensities;
-    if (!reference.empty()) {
-        reference_intensities.reserve(reference.size());
-        for (const auto &r: reference) {
-            if (!std::isfinite(r.I))
-                continue;
-
-            const auto hkl = key_generator(r);
-            reference_intensities[hkl.pack()] = r.I;
-        }
-    }
-
-    float d_min = std::numeric_limits<float>::max();
-    float d_max = 0.0f;
-
-    for (const auto &m: merged) {
-        if (!std::isfinite(m.d) || m.d <= 0.0f)
-            continue;
-        if (d_min_limit_A && m.d < d_min_limit_A)
-            continue;
-
-        d_min = std::min(d_min, m.d);
-        d_max = std::max(d_max, m.d);
-    }
-
-    if (!(d_min < d_max && d_min > 0.0f))
-        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
-            "MergeStats: Error in resolution calculation");
-
-    const float d_min_pad = d_min * 0.999f;
-    const float d_max_pad = d_max * 1.001f;
-
-    ResolutionShells shells(d_min_pad, d_max_pad, n_shells);
-    const auto shell_mean_1_d2 = shells.GetShellMeanOneOverResSq();
-    const auto shell_min_res = shells.GetShellMinRes();
-
-    std::vector<ShellAccum> acc(n_shells);
-
-    CalcPossibleReflections(x, cell, d_min_pad, d_max_pad, shells, acc);
-
-    CorrelationCoefficient cc_half_overall;
-    CorrelationCoefficient cc_ref_overall;
-
-    for (const auto &m: merged) {
-        const auto shell = shells.GetShell(m.d);
-        if (!shell.has_value())
-            continue;
-
-        const int s = *shell;
-        if (s >= 0 && s < n_shells) {
-            if (std::isfinite(m.I) && std::isfinite(m.sigma) && m.sigma > 0.0) {
-                acc[s].unique++;
-                acc[s].sum_i_over_sigma += m.I / m.sigma;
-                ++acc[s].n_i_over_sigma;
-
-                if (!reference_intensities.empty()) {
-                    const auto hkl = key_generator(m);
-                    const auto ref_it = reference_intensities.find(hkl.pack());
-                    if (ref_it != reference_intensities.end() && std::isfinite(ref_it->second)) {
-                        acc[s].cc_ref.Add(m.I, ref_it->second);
-                        cc_ref_overall.Add(m.I, ref_it->second);
-                    }
-                }
-
-                if (std::isfinite(m.I_half[0]) && std::isfinite(m.I_half[1])) {
-                    acc[s].cc_half.Add(m.I_half[0], m.I_half[1]);
-                    cc_half_overall.Add(m.I_half[0], m.I_half[1]);
-                }
-
-            }
-        }
-    }
-
-    for (int i = 0; i < reflections.size(); ++i) {
-        if (!merge_mask.empty() && merge_mask[i] == 0)
-            continue;
-
-        for (const auto &r: reflections[i]) {
-            if (key_generator.IsSystematicallyAbsent(r))
-                continue;
-            if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr))
-                continue;
-            if (!AcceptReflection(r, d_min_limit_A))
-                continue;
-            if (r.partiality < min_partiality)
-                continue;
-
-            const float I_corr = r.I * r.image_scale_corr;
-            const float sigma_corr = r.sigma * r.image_scale_corr;
-            if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f)
-                continue;
-
-            const auto shell = shells.GetShell(r.d);
-            if (!shell.has_value())
-                continue;
-            const int s = *shell;
-            if (s >= 0 && s < n_shells)
-                acc[s].total_obs++;
-        }
-    }
-
-    MergeStatistics out;
-    out.shells.resize(n_shells);
-
-    for (int s = 0; s < n_shells; ++s) {
-        const auto &sa = acc[s];
-        auto &ss = out.shells[s];
-
-        ss.mean_one_over_d2 = shell_mean_1_d2[s];
-        ss.d_min = shell_min_res[s];
-        ss.d_max = s == 0 ? d_max_pad : shell_min_res[s - 1];
-        ss.total_observations = sa.total_obs;
-        ss.unique_reflections = sa.unique;
-        ss.possible_unique_reflections = sa.possible;
-        ss.mean_i_over_sigma = sa.n_i_over_sigma > 0
-                                   ? sa.sum_i_over_sigma / sa.n_i_over_sigma
-                                   : 0.0;
-
-        ss.cc_half = sa.cc_half.GetCC();
-        ss.cc_ref = sa.cc_ref.GetCC();
-    }
-
-    auto &overall = out.overall;
-    overall.d_min = d_min;
-    overall.d_max = d_max;
-
-    int all_possible = 0;
-    int all_unique = 0;
-    double sum_i_over_sigma = 0.0;
-    int n_i_over_sigma = 0;
-
-
-    for (const auto &sa: acc) {
-        overall.total_observations += sa.total_obs;
-        all_unique += sa.unique;
-        all_possible += sa.possible;
-        sum_i_over_sigma += sa.sum_i_over_sigma;
-        n_i_over_sigma += sa.n_i_over_sigma;
-
-    }
-
-    overall.possible_unique_reflections = all_possible;
-    overall.unique_reflections = all_unique;
-    overall.mean_i_over_sigma = n_i_over_sigma > 0 ? sum_i_over_sigma / n_i_over_sigma : 0.0;
-    overall.cc_half = cc_half_overall.GetCC();
-    overall.cc_ref = cc_ref_overall.GetCC();
-
-    return out;
-}
-
-
 std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in) {
     double completeness = in.possible_unique_reflections > 0
         ? static_cast<double>(in.unique_reflections) / in.possible_unique_reflections * 100.0 : 0.0;
diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h
index 0297767a..5f6ca6e8 100644
--- a/image_analysis/scale_merge/Merge.h
+++ b/image_analysis/scale_merge/Merge.h
@@ -53,6 +53,8 @@ struct MergeAccum {
 class MergeOnTheFly {
     mutable std::mutex merged_mutex;
 
+    const int space_group_number = 1;
+
     ScalingSettings scaling_settings;
     IndexingSettings indexing_settings;
 
@@ -70,50 +72,20 @@ class MergeOnTheFly {
 
     std::map<uint64_t, MergeAccum> accumulator;
 
-    void ProcessImage_i(const std::vector<Reflection> &v, int half);
+    bool Mask(const IntegrationOutcome &outcome, bool cc_mask);
 public:
     MergeOnTheFly(const DiffractionExperiment &x);
     MergeOnTheFly& ReferenceCell(const std::optional<UnitCell> &cell);
 
-    void AddAll(const std::vector<std::vector<Reflection> > &reflections,
-                const std::vector<uint8_t> &merge_mask = {});
+    void AddImage(const IntegrationOutcome& outcome, bool cc_mask = false);
 
-    void AddAll(const std::vector<IntegrationOutcome> &integration_outcome,
-        const std::vector<uint8_t> &merge_mask = {});
-    void AddImage(const std::vector<Reflection> &v, double image_cc, const UnitCell &cell);
+    MergeStatistics MergeStats(const std::vector<MergedReflection> &merged,
+                               const std::vector<IntegrationOutcome> &reflections,
+                               const std::vector<MergedReflection> &reference = {});
 
     std::vector<MergedReflection> ExportReflections();
 };
 
-size_t CalcMergeMaskCC(const DiffractionExperiment &x,
-                       const std::vector<float> &scale_cc,
-                       std::vector<uint8_t> &result_mask);
-
-size_t CalcMergeMaskUnitCell(const DiffractionExperiment &x,
-                             const UnitCell &reference_cell,
-                             const std::vector<std::optional<UnitCell> > &unit_cells,
-                             std::vector<uint8_t> &mask);
-
-std::vector<MergedReflection> MergeAll(const DiffractionExperiment &x,
-                                       const std::vector<std::vector<Reflection> > &reflections,
-                                       const std::vector<uint8_t> &merge_mask = {});
-
-
 std::vector<MergedReflection> MergeAll(const DiffractionExperiment &x,
                                        const std::vector<IntegrationOutcome> &reflections,
-                                       const std::vector<uint8_t> &merge_mask = {});
-
-MergeStatistics MergeStats(const DiffractionExperiment &x,
-                           const std::vector<MergedReflection> &merged,
-                           const std::vector<std::vector<Reflection> > &reflections,
-                           const UnitCell &cell,
-                           const std::vector<uint8_t> &merge_mask = {},
-                           const std::vector<MergedReflection> &reference = {});
-
-
-MergeStatistics MergeStats(const DiffractionExperiment &x,
-                           const std::vector<MergedReflection> &merged,
-                           const std::vector<IntegrationOutcome> &reflections,
-                           const UnitCell &cell,
-                           const std::vector<uint8_t> &merge_mask = {},
-                           const std::vector<MergedReflection> &reference = {});
+                                       bool mask = false);
diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp
index a70c9647..94ff070b 100644
--- a/image_analysis/scale_merge/ScaleOnTheFly.cpp
+++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp
@@ -504,159 +504,7 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const {
     }
 }
 
-ScaleOnTheFlyResult ScaleOnTheFly::Scale(std::vector<Reflection> &reflections,
-    std::optional<float> mosaicity_deg) const {
-    auto start = std::chrono::steady_clock::now();
-
-    ceres::Problem problem;
-
-    ScaleOnTheFlyResult result{
-        .B = 0.0,
-        .G = 1.0
-    };
-
-    if (model == PartialityModel::Rotation) {
-        if (mosaicity_deg && std::isfinite(*mosaicity_deg) && *mosaicity_deg > 0.0)
-            result.mos = *mosaicity_deg;
-        else
-            result.mos = s.GetDefaultMosaicity();
-        result.wedge = rot_wedge_deg.value_or(0.0);
-    } else {
-        result.mos = NAN;
-        result.wedge = NAN;
-    }
-
-    size_t n_reflections = 0;
-    for (const auto &r: reflections) {
-        if (!Accept(r))
-            continue;
-
-        const HKLKey key = hkl_key_generator(r);
-
-        if (!reference_data.contains(key))
-            continue;
-
-        ++n_reflections;
-
-        const double Itrue = reference_data.at(key);
-        const double sigma = r.sigma;
-
-        switch (model) {
-            case PartialityModel::Fixed: {
-                auto *cost = new ceres::AutoDiffCostFunction<IntensityFixedResidual, 1, 1, 1>(
-                    new IntensityFixedResidual(r, Itrue, sigma, r.partiality));
-                problem.AddResidualBlock(cost, nullptr, &result.G, &result.B);
-            }
-            break;
-            case PartialityModel::Unity: {
-                auto *cost = new ceres::AutoDiffCostFunction<IntensityFixedResidual, 1, 1, 1>(
-                    new IntensityFixedResidual(r, Itrue, sigma, 1.0));
-                problem.AddResidualBlock(cost, nullptr, &result.G, &result.B);
-            }
-            break;
-            case PartialityModel::Rotation: {
-                auto *cost = new ceres::AutoDiffCostFunction<ScalingRotationResidual, 1, 1, 1, 1, 1>(
-                    new ScalingRotationResidual(r, Itrue, sigma));
-                problem.AddResidualBlock(cost, nullptr, &result.G, &result.B, &result.mos,
-                    &result.wedge);
-            }
-                break;
-            default:
-                throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
-                                      "Not supported partiality model");
-        }
-    }
-
-    if (n_reflections < MIN_REFLECTIONS) {
-        result.succesful = false;
-        return result;
-    }
-    result.succesful = true;
-
-    constexpr double SIGMA_SCALE_FACTOR = 2; // We assume scale factor is +/- 2.0
-    const double scale_factor_regularization_weight = std::sqrt(static_cast<double>(n_reflections) / SIGMA_SCALE_FACTOR);
-
-    if (s.GetScalingRegularize()) {
-        auto *scale_reg_cost = new ceres::AutoDiffCostFunction<RegularizationResidual, 1, 1>(
-            new RegularizationResidual(scale_factor_regularization_weight, 1.0));
-        problem.AddResidualBlock(scale_reg_cost, nullptr, &result.G);
-    }
-
-    problem.SetParameterLowerBound(&result.G, 0, 0.0);
-
-    if (s.GetRefineB()) {
-        constexpr double SIGMA_B = 10.0; // in A^2
-        const double b_regularization_weight = std::sqrt(static_cast<double>(n_reflections) / SIGMA_B);
-
-        if (s.GetScalingRegularize()) {
-            auto *b_reg_cost = new ceres::AutoDiffCostFunction<RegularizationResidual, 1, 1>(
-                new RegularizationResidual(b_regularization_weight, 0.0));
-            problem.AddResidualBlock(b_reg_cost, nullptr, &result.B);
-        }
-
-        problem.SetParameterLowerBound(&result.B, 0, s.GetMinB());
-        problem.SetParameterUpperBound(&result.B, 0, s.GetMaxB());
-    } else {
-        problem.SetParameterBlockConstant(&result.B);
-    }
-
-    if (model == PartialityModel::Rotation) {
-        if (refine_rot_wedge) {
-            problem.SetParameterLowerBound(&result.wedge, 0, s.GetMinWedge());
-            problem.SetParameterUpperBound(&result.wedge, 0, s.GetMaxWedge());
-        } else {
-            problem.SetParameterBlockConstant(&result.wedge);
-        }
-        problem.SetParameterLowerBound(&result.mos, 0, s.GetMinMosaicity());
-        problem.SetParameterUpperBound(&result.mos, 0, s.GetMaxMosaicity());
-    }
-
-    ceres::Solver::Options options;
-    options.linear_solver_type = ceres::DENSE_QR;
-    options.minimizer_progress_to_stdout = false;
-    options.num_threads = 1;
-
-    ceres::Solver::Summary summary;
-    ceres::Solve(options, &problem, &summary);
-
-    for (auto &r: reflections) {
-        const double B_term = exp(result.B * -SafeInv(4.0 * r.d * r.d, 0.0));
-
-        switch (model) {
-            case PartialityModel::Unity:
-                r.partiality = 1.0;
-                break;
-            case PartialityModel::Rotation: {
-                if (std::isfinite(r.delta_phi_deg) && std::isfinite(r.zeta) && result.mos > 1e-6)
-                    r.partiality = RotationPartiality(r.delta_phi_deg, r.zeta, result.mos, result.wedge);
-                break;
-            }
-            default:
-                // For fixed partiality there is no need to change anything
-                break;
-        }
-        const double denom = B_term * r.partiality * result.G;
-        if (std::isfinite(r.rlp) &&
-            std::isfinite(denom) &&
-            denom > 0.0) {
-                r.image_scale_corr = static_cast<float>(r.rlp / denom);
-            } else {
-                r.image_scale_corr = NAN;
-            }
-    }
-
-    const auto [cc, cc_n] = CalculateGlobalCC(reflections);
-    result.cc = cc;
-    result.cc_n = cc_n;
-
-    auto end = std::chrono::steady_clock::now();
-    result.time_s = std::chrono::duration<float>(end - start).count();
-    return result;
-}
-
 void ScaleOnTheFly::Scale(std::vector<IntegrationOutcome> &integration, size_t nthreads) const {
-    ScalingResult result(integration.size());
-
     if (nthreads == 0)
         nthreads = std::thread::hardware_concurrency();
 
@@ -682,63 +530,3 @@ void ScaleOnTheFly::Scale(std::vector<IntegrationOutcome> &integration, size_t n
             f.get();
     }
 }
-
-ScalingResult ScaleOnTheFly::Scale(std::vector<std::vector<Reflection> > &reflections,
-                        const std::vector<float> &mosaicity,
-                        size_t nthreads) {
-    ScalingResult result(reflections.size());
-
-    if (nthreads == 0)
-        nthreads = std::thread::hardware_concurrency();
-
-    if (nthreads <= 1) {
-        for (int i = 0; i < reflections.size(); i++) {
-            std::optional<float> mos_val;
-            if (model == PartialityModel::Rotation
-                && mosaicity.size() > i
-                && std::isfinite(mosaicity[i])
-                && mosaicity[i] > 0.0)
-                mos_val = mosaicity[i];
-
-            auto local_result = Scale(reflections[i], mos_val);
-            result.mosaicity_deg[i] = local_result.mos;
-            result.image_bfactor_Ang2[i] = local_result.B;
-            result.image_scale_g[i] = local_result.G;
-            result.rotation_wedge_deg[i] = local_result.wedge;
-            result.image_cc[i] = local_result.cc;
-            result.image_cc_n[i] = local_result.cc_n;
-        }
-    } else {
-        auto local_nthreads = std::min(nthreads, reflections.size());
-        std::vector<std::future<void>> futures;
-        futures.reserve(local_nthreads);
-        std::atomic<size_t> curr_image = 0;
-
-        for (size_t t = 0; t < local_nthreads; ++t)
-            futures.emplace_back(std::async(std::launch::async, [&] {
-                size_t i = curr_image.fetch_add(1);
-                while (i < reflections.size()) {
-                    std::optional<float> mos_val;
-                    if (model == PartialityModel::Rotation
-                        && mosaicity.size() > i
-                        && std::isfinite(mosaicity[i])
-                        && mosaicity[i] > 0.0)
-                        mos_val = mosaicity[i];
-
-                    auto local_result = Scale(reflections[i], mos_val);
-                    result.mosaicity_deg[i] = local_result.mos;
-                    result.image_bfactor_Ang2[i] = local_result.B;
-                    result.image_scale_g[i] = local_result.G;
-                    result.rotation_wedge_deg[i] = local_result.wedge;
-                    result.image_cc[i] = local_result.cc;
-                    result.image_cc_n[i] = local_result.cc_n;
-                    i = curr_image.fetch_add(1);
-                }
-            }));
-
-        for (auto &f: futures)
-            f.get();
-    }
-
-    return result;
-}
diff --git a/image_analysis/scale_merge/ScaleOnTheFly.h b/image_analysis/scale_merge/ScaleOnTheFly.h
index 11f685ae..6c2a28d2 100644
--- a/image_analysis/scale_merge/ScaleOnTheFly.h
+++ b/image_analysis/scale_merge/ScaleOnTheFly.h
@@ -6,7 +6,6 @@
 #include "HKLKey.h"
 #include "Merge.h"
 #include "../../common/DiffractionExperiment.h"
-#include "ScalingResult.h"
 #include "../IntegrationOutcome.h"
 
 #include <map>
@@ -40,12 +39,6 @@ public:
     ScaleOnTheFly(const DiffractionExperiment &x, const std::vector<MergedReflection> &ref);
 
     void Scale(IntegrationOutcome &integration_outcome) const;
-    ScaleOnTheFlyResult Scale(std::vector<Reflection> &r, std::optional<float> mosaicity_deg) const;
-
-    ScalingResult Scale(std::vector<std::vector<Reflection> > &reflections,
-                        const std::vector<float> &mosaicity,
-                        size_t nthreads = 0);
-
     void Scale(std::vector<IntegrationOutcome> &integration_outcome, size_t nthreads = 0) const;
 };
 
diff --git a/reader/JFJochHDF5Reader.cpp b/reader/JFJochHDF5Reader.cpp
index 62a78c31..6a5f4f05 100644
--- a/reader/JFJochHDF5Reader.cpp
+++ b/reader/JFJochHDF5Reader.cpp
@@ -613,6 +613,7 @@ void JFJochHDF5Reader::ReadFile(const std::string &filename) {
             dataset->pixel_mask = PixelMask(mask_tmp);
         }
         dataset->experiment.ImagesPerTrigger(number_of_images);
+        cached_geom = dataset->experiment.GetDiffractionGeometry();
         SetStartMessage(dataset);
     } catch (const std::exception &e) {
         master_file = {};
@@ -969,6 +970,7 @@ void JFJochHDF5Reader::Close() {
     vds_data_mappings.clear();
     master_file_directory.clear();
     master_filename.clear();
+    cached_geom = DiffractionGeometry{};
     SetStartMessage({});
 }
 
@@ -1162,38 +1164,91 @@ std::vector<HDF5DataSourceMessage> JFJochHDF5Reader::GetHDF5DataSource(
 }
 
 
-void JFJochHDF5Reader::ReadReflections(std::vector<std::vector<Reflection> > &output, size_t start_image, std::optional<size_t> end_image) const {
+std::vector<IntegrationOutcome> JFJochHDF5Reader::ReadReflections(size_t start_image, std::optional<size_t> end_image) const {
     std::unique_lock ul(hdf5_mutex);
+
     if (start_image >= number_of_images)
         throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                               "start_image must be less than number_of_images");
 
-    if (end_image.has_value() && end_image.value() < start_image)
-        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
-                              "end_image must be greater or equal than start_image if provided");
+    const size_t end_image_val = end_image.value_or(number_of_images - 1);
 
-    int end_image_val = end_image.value_or(number_of_images - 1);
-
-    if (end_image_val > number_of_images - 1)
+    if (end_image_val < start_image)
         throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
-            "end_image_val must be less than or equal to number_of_images");
+                              "end_image must be greater or equal to start_image if provided");
+
+    if (end_image_val >= number_of_images)
+        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+                              "end_image must be less than number_of_images");
 
     if (!master_file)
         throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
-                              "Cannot read all reflections if file not loaded");
+                              "Cannot read reflections if file not loaded");
 
-    size_t output_start = output.size();
+    // Cache of data files opened during this call (keyed by path) to avoid
+    // re-opening the same file for every image in a multi-image batch.
+    std::map<std::string, std::shared_ptr<HDF5ReadOnlyFile>> open_files;
 
-    output.resize(output.size() + end_image_val - start_image + 1);
+    auto get_data_file = [&](const std::string &path) -> std::shared_ptr<HDF5ReadOnlyFile> {
+        auto it = open_files.find(path);
+        if (it != open_files.end())
+            return it->second;
+        auto f = std::make_shared<HDF5ReadOnlyFile>(path);
+        open_files[path] = f;
+        return f;
+    };
 
-    for (int i = 0; i < end_image_val - start_image + 1; ++i)
-        ReadReflectionsFromGroup(*master_file,
-            fmt::format("/entry/reflections/image_{:06d}", start_image + i),
-            output[output_start + i]);
-}
+    std::vector<IntegrationOutcome> ret;
+    ret.reserve(end_image_val - start_image + 1);
+
+    for (size_t img = start_image; img <= end_image_val; img++) {
+        IntegrationOutcome outcome;
+
+        // Generic (non-image-specific) detector geometry from experiment setup.
+        outcome.geom = cached_geom;
+
+        // Determine which HDF5 file holds the reflections and MX metadata for
+        // this image, and what local image index to use within that file.
+        HDF5Object *meta_file = master_file.get();
+        size_t meta_image_id = img;
+
+        std::shared_ptr<HDF5ReadOnlyFile> owned_source; // keeps the file alive
+
+        if (format == FileWriterFormat::NXmxLegacy) {
+            // Reflections and per-image MX data are in the individual data files,
+            // addressed by the source-local image index.
+            const size_t file_id = img / images_per_file;
+            meta_image_id = img % images_per_file;
+            owned_source = get_data_file(legacy_format_files.at(file_id));
+            meta_file = owned_source.get();
+        }
+        // NXmxVDS with contiguous layout (integrated): everything is in master_file
+        // with global indices — the defaults set above are correct.
+
+        // ── reflections ──────────────────────────────────────────────────────
+        const std::string refl_group = fmt::format("/entry/reflections/image_{:06d}", meta_image_id);
+        ReadReflectionsFromGroup(*meta_file, refl_group, outcome.reflections);
+
+        // ── per-image mosaicity ───────────────────────────────────────────────
+        if (meta_file->Exists("/entry/MX/mosaicity")) {
+            try {
+                outcome.mosaicity_deg =
+                    meta_file->ReadElement<float>("/entry/MX/mosaicity", meta_image_id);
+            } catch (...) {}
+        }
+
+        // ── indexed lattice (stored as 9-element row-major matrix) ────────────
+        if (meta_file->Exists("/entry/MX/latticeIndexed")) {
+            try {
+                auto lattice_vec = meta_file->ReadOptVector<float>(
+                    "/entry/MX/latticeIndexed",{meta_image_id, 0}, {1, 9});
+                if (lattice_vec.size() == 9)
+                    outcome.latt = CrystalLattice(lattice_vec);
+            } catch (...) {}
+        }
+
+        ret.push_back(std::move(outcome));
+    }
 
-std::vector<std::vector<Reflection> > JFJochHDF5Reader::ReadReflections(size_t start_image, std::optional<size_t> end_image) const {
-    std::vector<std::vector<Reflection> > ret;
-    ReadReflections(ret, start_image, end_image);
     return ret;
-}
+}
\ No newline at end of file
diff --git a/reader/JFJochHDF5Reader.h b/reader/JFJochHDF5Reader.h
index 9531fd08..f952462a 100644
--- a/reader/JFJochHDF5Reader.h
+++ b/reader/JFJochHDF5Reader.h
@@ -6,6 +6,7 @@
 
 #include "JFJochReader.h"
 #include "../writer/HDF5Objects.h"
+#include "../image_analysis/IntegrationOutcome.h"
 
 class JFJochHDF5Reader : public JFJochReader {
     FileWriterFormat format = FileWriterFormat::NoFile;
@@ -13,6 +14,8 @@ class JFJochHDF5Reader : public JFJochReader {
 
     std::shared_ptr<HDF5ReadOnlyFile> master_file;
 
+    DiffractionGeometry cached_geom;
+
     std::vector<std::string> legacy_format_files;
     std::vector<HDF5VirtualDatasetMapping> vds_data_mappings;
     std::string master_file_directory;
@@ -52,8 +55,7 @@ public:
         std::optional<uint64_t> image_count = {}
     ) const;
 
-    std::vector<std::vector<Reflection>> ReadReflections(size_t start_image = 0, std::optional<size_t> end_image = {}) const;
-    void ReadReflections(std::vector<std::vector<Reflection>> &output, size_t start_image = 0, std::optional<size_t> end_image = {}) const;
+    std::vector<IntegrationOutcome> ReadReflections(size_t start_image = 0, std::optional<size_t> end_image = {}) const;
 
     CompressedImage ReadCalibration(std::vector<uint8_t> &tmp, const std::string &name) const;
 
@@ -61,4 +63,4 @@ public:
 };
 
 
-#endif //JFJOCHHDF5IMAGEREADER_H
+#endif //JFJOCHHDF5IMAGEREADER_H
\ No newline at end of file
diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp
index ebd7799f..332e4bf9 100644
--- a/tests/JFJochReaderTest.cpp
+++ b/tests/JFJochReaderTest.cpp
@@ -1970,6 +1970,9 @@ TEST_CASE("JFJochReader_ReadReflections_VDS", "[HDF5][Full]") {
                         .zeta = 0.02f
                     }
                 };
+
+                message.mosaicity_deg = i*0.15f;
+                message.indexing_lattice = CrystalLattice({100,0,0}, {0,50,0}, {0,0,30});
             }
 
             REQUIRE_NOTHROW(writer.WriteHDF5(message));
@@ -1992,23 +1995,28 @@ TEST_CASE("JFJochReader_ReadReflections_VDS", "[HDF5][Full]") {
 
         REQUIRE(reflections.size() == 4);
 
-        CHECK(reflections[0].empty());
+        CHECK(reflections[0].reflections.empty());
 
-        REQUIRE(reflections[1].size() == 2);
-        CHECK(reflections[1][0].h == 11);
-        CHECK(reflections[1][0].k == 20);
-        CHECK(reflections[1][0].l == 30);
-        CHECK(reflections[1][0].I == Catch::Approx(1001.0f));
-        CHECK(reflections[1][0].predicted_x == Catch::Approx(101.0f));
-        CHECK(reflections[1][0].predicted_y == Catch::Approx(201.0f));
+        REQUIRE(reflections[1].reflections.size() == 2);
+        CHECK(reflections[1].reflections[0].h == 11);
+        CHECK(reflections[1].reflections[0].k == 20);
+        CHECK(reflections[1].reflections[0].l == 30);
+        CHECK(reflections[1].reflections[0].I == Catch::Approx(1001.0f));
+        CHECK(reflections[1].reflections[0].predicted_x == Catch::Approx(101.0f));
+        CHECK(reflections[1].reflections[0].predicted_y == Catch::Approx(201.0f));
+        CHECK(reflections[1].mosaicity_deg == Catch::Approx(0.15f));
+        CHECK(reflections[1].latt.CalcVolume() == Catch::Approx(100*50*30));
 
-        CHECK(reflections[2].empty());
+        CHECK(reflections[2].reflections.empty());
 
-        REQUIRE(reflections[3].size() == 2);
-        CHECK(reflections[3][0].h == 13);
-        CHECK(reflections[3][0].I == Catch::Approx(1003.0f));
-        CHECK(reflections[3][1].h == 43);
-        CHECK(reflections[3][1].I == Catch::Approx(2003.0f));
+        REQUIRE(reflections[3].reflections.size() == 2);
+        CHECK(reflections[3].reflections[0].h == 13);
+        CHECK(reflections[3].reflections[0].I == Catch::Approx(1003.0f));
+        CHECK(reflections[3].reflections[1].h == 43);
+        CHECK(reflections[3].reflections[1].I == Catch::Approx(2003.0f));
+        CHECK(reflections[3].mosaicity_deg == Catch::Approx(0.45f));
+        CHECK(reflections[3].latt.Vec0().x == Catch::Approx(100.0f));
+        CHECK(reflections[3].latt.Vec1().y == Catch::Approx(50.0f));
     }
 
     {
@@ -2019,13 +2027,13 @@ TEST_CASE("JFJochReader_ReadReflections_VDS", "[HDF5][Full]") {
 
         REQUIRE(reflections.size() == 3);
 
-        REQUIRE(reflections[0].size() == 2); // original image 1
-        CHECK(reflections[0][0].h == 11);
+        REQUIRE(reflections[0].reflections.size() == 2); // original image 1
+        CHECK(reflections[0].reflections[0].h == 11);
 
-        CHECK(reflections[1].empty());       // original image 2
+        CHECK(reflections[1].reflections.empty());       // original image 2
 
-        REQUIRE(reflections[2].size() == 2); // original image 3
-        CHECK(reflections[2][0].h == 13);
+        REQUIRE(reflections[2].reflections.size() == 2); // original image 3
+        CHECK(reflections[2].reflections[0].h == 13);
     }
 
     remove("read_reflections_vds_master.h5");
diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp
index 0a8d72b9..53edfc10 100644
--- a/tools/jfjoch_process.cpp
+++ b/tools/jfjoch_process.cpp
@@ -696,12 +696,6 @@ int main(int argc, char **argv) {
         logger.Info("UC:         a={:.2f} b={:.2f} c={:.2f} alpha={:.2f} beta={:.2f} gamma={:.2f}",
                         consensus_cell->a, consensus_cell->b, consensus_cell->c,
                         consensus_cell->alpha, consensus_cell->beta, consensus_cell->gamma);
-
-        auto rejected_uc = CalcMergeMaskUnitCell(experiment, *consensus_cell, indexer.GetUnitCells(), merging_mask_uc);
-        if (rejected_uc > 0)
-            logger.Info("Rejected {} images for merging due to unit cell being too far from consensus", rejected_uc);
-        // UpdateReflectionResolution(consensus_cell.value(), indexer.GetReflections());
-        logger.Info("Reflection resolution updated based on consensus unit cell");
     } else
         logger.Info("Consensus unit cell not found - calculation tool {:.2f} ms", consensus_duration * 1e3);
     end_msg.unit_cell = consensus_cell;
@@ -709,8 +703,6 @@ int main(int argc, char **argv) {
     if (run_scaling || !reference_data.empty()) {
         logger.Info("Running scaling (mosaicity refinement) ...");
 
-        std::vector<float> scale_cc;
-
         if (reference_data.empty()) {
             // If reference data are given, there is live scaling (no need to go again)
             ScalingResult scale_result(0);
@@ -718,7 +710,7 @@ int main(int argc, char **argv) {
             auto scale_start = std::chrono::steady_clock::now();
             for (int i = 0; i < scaling_iter; i++) {
                 auto iter_start = std::chrono::steady_clock::now();
-                auto merge_result = MergeAll(experiment, indexer.GetIntegrationOutcome(), merging_mask_uc);
+                auto merge_result = MergeAll(experiment, indexer.GetIntegrationOutcome(), false);
                 scale_result = indexer.ScaleAllImages(merge_result);
                 scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat");
                 auto iter_end = std::chrono::steady_clock::now();
@@ -726,7 +718,6 @@ int main(int argc, char **argv) {
                 logger.Info("Scaling iteration {} took {:.3f} seconds", i, iter_time);
             }
 
-            scale_cc = scale_result.image_cc;
             end_msg.image_scale_factor = scale_result.image_scale_g;
             end_msg.image_scale_cc = scale_result.image_cc;
             end_msg.image_scale_mosaicity = scale_result.mosaicity_deg;
@@ -735,18 +726,18 @@ int main(int argc, char **argv) {
             auto scale_end = std::chrono::steady_clock::now();
             double scale_time = std::chrono::duration<double>(scale_end - scale_start).count();
             logger.Info("Scaling completed in {:.2f} s", scale_time);
-        } else
-            scale_cc = indexer.GetImageCC();
+        }
 
         auto merge_start = std::chrono::steady_clock::now();
 
-        auto rejected_cc = CalcMergeMaskCC(experiment, scale_cc, merging_mask_uc);
-        if (rejected_cc > 0)
-            logger.Info("Rejected {} images for merging due to low CC with reference", rejected_cc);
+        MergeOnTheFly merge_engine(experiment);
+        if (consensus_cell.has_value())
+            merge_engine.ReferenceCell(*consensus_cell);
+        for (auto &i : indexer.GetIntegrationOutcome())
+            merge_engine.AddImage(i);
 
-        auto merged_reflections = MergeAll(experiment, indexer.GetIntegrationOutcome(), merging_mask_uc);
-        auto merged_statistics = MergeStats(experiment, merged_reflections, indexer.GetIntegrationOutcome(), *consensus_cell, merging_mask_uc,
-            reference_data);
+        auto merged_reflections = merge_engine.ExportReflections();
+        auto merged_statistics = merge_engine.MergeStats(merged_reflections, indexer.GetIntegrationOutcome(), reference_data);
 
         auto merge_end = std::chrono::steady_clock::now();
         double merge_time = std::chrono::duration<double>(merge_end - merge_start).count();
diff --git a/tools/jfjoch_scale.cpp b/tools/jfjoch_scale.cpp
index cbfda64c..e6da83f2 100644
--- a/tools/jfjoch_scale.cpp
+++ b/tools/jfjoch_scale.cpp
@@ -208,7 +208,7 @@ int main(int argc, char **argv) {
         logger.Info("Using space group {} (number {})", space_group->hm, space_group_number.value());
     }
 
-    std::vector<std::vector<Reflection>> reflections;
+    std::vector<IntegrationOutcome> reflections;
     std::vector<float> mosaicity;
 
     if (optind >= argc) {
@@ -240,7 +240,7 @@ int main(int argc, char **argv) {
             for (int i = 0; i < end_image - start_image + 1; i++)
                 mosaicity[mosaicity_size + i] = dataset->mosaicity_deg[start_image + i];
 
-            reader.ReadReflections(reflections, start_image, end_image);
+            reflections = reader.ReadReflections(start_image, end_image);
             logger.Info("Loaded dataset from {}", input_file);
         } catch (const std::exception &e) {
             logger.Error("Error reading input file: {}", e.what());
@@ -250,7 +250,6 @@ int main(int argc, char **argv) {
         }
     }
 
-
     std::vector<MergedReflection> reference_data;
     if (!ref_mtz.empty()) {
         reference_data = LoadFCalcFromMtz(ref_mtz);
@@ -297,8 +296,6 @@ int main(int argc, char **argv) {
 
     experiment.ImagesPerTrigger(refl_stats.n_images);
 
-    ScalingResult scale_result(0);
-
     auto scale_start = std::chrono::steady_clock::now();
     for (int i = 0; i < scaling_iter; i++) {
         auto iter_start = std::chrono::steady_clock::now();
@@ -306,13 +303,13 @@ int main(int argc, char **argv) {
         if (reference_data.empty()) {
             auto merged = MergeAll(experiment, reflections);
             ScaleOnTheFly scaling(experiment, merged);
-            scale_result = scaling.Scale(reflections, mosaicity, nthreads);
+            scaling.Scale(reflections, nthreads);
         } else {
             ScaleOnTheFly scaling(experiment, reference_data);
-            scale_result = scaling.Scale(reflections, mosaicity, nthreads);
+            scaling.Scale(reflections, nthreads);
         }
 
-        scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat");
+        // scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat");
         auto iter_end = std::chrono::steady_clock::now();
         double iter_time = std::chrono::duration<double>(iter_end - iter_start).count();
         logger.Info("Scaling iteration {} took {:.3f} seconds", i, iter_time);
@@ -324,20 +321,21 @@ int main(int argc, char **argv) {
 
     auto merge_start = std::chrono::steady_clock::now();
 
-    std::vector<uint8_t> merge_mask(reflections.size(), 1);
-    auto rejected = CalcMergeMaskCC(experiment, scale_result.image_cc, merge_mask);
-    if (rejected > 0)
-        logger.Info("Rejected {} images due to low CC with reference", rejected);
+    MergeOnTheFly merge_engine(experiment);
+        merge_engine.ReferenceCell(experiment.GetUnitCell());
+    for (auto &i : reflections)
+        merge_engine.AddImage(i);
 
-    auto merge_result = MergeAll(experiment, reflections, merge_mask);
-    auto merge_stats = MergeStats(experiment, merge_result, reflections, experiment.GetUnitCell().value(), merge_mask,
-                                  reference_data);
+    auto merged_reflections = merge_engine.ExportReflections();
+    auto merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data);
 
     auto merge_end = std::chrono::steady_clock::now();
     double merge_time = std::chrono::duration<double>(merge_end - merge_start).count();
 
-    logger.Info("Merge completed in {:.2f} s ({} unique reflections)", merge_time,
-                merge_result.size());
+    logger.Info("Merge completed in {:.2f} s ({} unique reflections)", merge_time, merged_reflections.size());
+
+    // Print resolution-shell statistics table
+    std::cout << merged_statistics;
 
     const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value();
 
@@ -355,13 +353,10 @@ int main(int argc, char **argv) {
         sg_opts.min_total_compared = 100;
         sg_opts.test_systematic_absences = true;
 
-        const auto sg_search = SearchSpaceGroup(merge_result, sg_opts);
+        const auto sg_search = SearchSpaceGroup(merged_reflections, sg_opts);
         std::cout << std::endl << SearchSpaceGroupResultToText(sg_search) << std::endl;
     }
 
-    // Print resolution-shell statistics table
-    std::cout << merge_stats;
-
     if (!output_prefix.empty())
-        WriteReflections(merge_result, *experiment.GetUnitCell(), experiment, output_prefix);
+        WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, output_prefix);
 }
diff --git a/writer/HDF5Objects.h b/writer/HDF5Objects.h
index d011b51e..bd0b170d 100644
--- a/writer/HDF5Objects.h
+++ b/writer/HDF5Objects.h
@@ -171,6 +171,7 @@ public:
     std::string GetString(const std::string& name, const std::string& def = "") const;
 
     template <class T> std::optional<T> ReadElement(const std::string& name, size_t n) const;
+
     template <class T> std::vector<T> ReadVector(const std::string &name);
     template <class T> std::vector<T> ReadOptVector(const std::string &name);
     template <class T> std::vector<T> ReadVector(const std::string &name,