rugnux: automatic CC1/2-logistic high-resolution cutoff

Merge/report/write to the CC1/2 fall-off instead of the detector edge.
After the final merge, fit a logistic to a fine-binned CC1/2 profile and
cut one shell past cc_target (default 0.30) - generous, DIALS-style. Only
the written MTZ/CIF and the reported shell table are trimmed; scaling, the
3D combine and the error model still run over the full range, and the
per-image _process.h5 keeps full data for re-scaling. The de-novo P1 merge
that feeds the space-group search is never cut (gated on the final merge).

- ScalingSettings: resolution_cutoff {Off,CCHalfLogistic}, resolution_cc_target,
  report_shell_count (report shell count is now configurable).
- New ResolutionCutoff.{h,cpp}: ComputeCCHalfLogisticCutoff (fine-bin CC1/2 +
  grid/refine logistic fit + one-shell extension + degeneracy/low-res guards).
- Wired into RotationScaleMerge (rotation) and MergeOnTheFly/Rugnux (stills),
  plus the rugnux --scale path.
- CLI: --resolution-cutoff cc-logistic|off, --resolution-cc-target,
  --resolution-shells. --scaling-high-resolution stays the manual override
  (wins and disables the auto-cut).

Validated on the 24-crystal battery: no crashes (incl. Benas/cytC_10/Thau_9);
--resolution-cutoff off is bit-identical to baseline; cutoffs reproduce the
design table within ~0.05 A (lysoC 1.30, Ins_I 1.43, cytC_2 1.68, EcwtCQ244-E3
1.54); good-to-edge datasets (lyso_ref, EP_cs_02-10) correctly keep the full
range.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

common/ScalingSettings.cpp

@@ -209,4 +209,37 @@ ScalingSettings &ScalingSettings::ScalingRegularize(bool input) {
 
 bool ScalingSettings::GetScalingRegularize() const {
     return scaling_regularize;
-}
+}
+
+ScalingSettings &ScalingSettings::ResolutionCutoff(ResolutionCutoffMethod input) {
+    resolution_cutoff = input;
+    return *this;
+}
+
+ResolutionCutoffMethod ScalingSettings::GetResolutionCutoff() const {
+    return resolution_cutoff;
+}
+
+ScalingSettings &ScalingSettings::ResolutionCCTarget(double input) {
+    if (input <= 0.0 || input >= 1.0)
+        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+                              "Resolution CC target must be between 0 and 1");
+    resolution_cc_target = input;
+    return *this;
+}
+
+double ScalingSettings::GetResolutionCCTarget() const {
+    return resolution_cc_target;
+}
+
+ScalingSettings &ScalingSettings::ReportShellCount(int input) {
+    if (input < 1)
+        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+                              "Number of report shells must be at least 1");
+    report_shell_count = input;
+    return *this;
+}
+
+int ScalingSettings::GetReportShellCount() const {
+    return report_shell_count;
+}

common/ScalingSettings.h

@@ -8,6 +8,11 @@
 
 enum class IntensityFormat { Text, mmCIF, MTZ};
 
+// How the high-resolution cutoff for the written reflections and the reported shell table is chosen
+// when no explicit --scaling-high-resolution is given. Off = keep the full (detector-edge) range;
+// CCHalfLogistic = fit the CC1/2 fall-off and cut one shell past cc_target (DIALS-style, generous).
+enum class ResolutionCutoffMethod { Off, CCHalfLogistic };
+
 class ScalingSettings {
     bool refine_b = false;
     double max_b = 200.0;
@@ -45,6 +50,13 @@ class ScalingSettings {
     double rfree_fraction = 0.05;
     IntensityFormat intensity_format = IntensityFormat::mmCIF;
 
+    // Automatic high-resolution cutoff for the written reflections + reported shells (not the scaling
+    // or the error model, and not the per-image _process.h5). Applied only when no explicit
+    // high_resolution_limit_A is set - that manual limit always wins and disables the auto-cut.
+    ResolutionCutoffMethod resolution_cutoff = ResolutionCutoffMethod::CCHalfLogistic;
+    double resolution_cc_target = 0.30;   // CC1/2 value defining the fall-off limit before the +1 shell
+    int report_shell_count = 10;          // number of resolution shells in the reported statistics table
+
     bool scaling_regularize = false;
 public:
     ScalingSettings& RefineB(bool input);
@@ -66,6 +78,10 @@ public:
     ScalingSettings& FileFormat(IntensityFormat input);
     ScalingSettings& ScalingRegularize(bool input);
 
+    ScalingSettings& ResolutionCutoff(ResolutionCutoffMethod input);
+    ScalingSettings& ResolutionCCTarget(double input);
+    ScalingSettings& ReportShellCount(int input);
+
     [[nodiscard]] bool GetRefineB() const;
     [[nodiscard]] bool GetRefineWedge() const;
 
@@ -96,4 +112,8 @@ public:
     [[nodiscard]] double GetRfreeFraction() const;
     [[nodiscard]] IntensityFormat GetFileFormat() const;
     [[nodiscard]] bool GetScalingRegularize() const;
+
+    [[nodiscard]] ResolutionCutoffMethod GetResolutionCutoff() const;
+    [[nodiscard]] double GetResolutionCCTarget() const;
+    [[nodiscard]] int GetReportShellCount() const;
 };

image_analysis/scale_merge/CMakeLists.txt

@@ -9,6 +9,8 @@ ADD_LIBRARY(JFJochScaleMerge
         ScaleOnTheFly.h
         RotationScaleMerge.cpp
         RotationScaleMerge.h
+        ResolutionCutoff.cpp
+        ResolutionCutoff.h
         HKLKey.cpp
         HKLKey.h
         ScalingResult.h

image_analysis/scale_merge/Merge.cpp

@@ -572,11 +572,13 @@ void CalcPossibleReflections(int space_group_number ,
 
 MergeStatistics MergeOnTheFly::MergeStats(const std::vector<MergedReflection> &merged,
                            const std::vector<IntegrationOutcome > &integration_outcome,
-                           const std::vector<MergedReflection> &reference) {
+                           const std::vector<MergedReflection> &reference,
+                           std::optional<double> d_min_override) {
 
-    constexpr int n_shells = 10;
+    const int n_shells = scaling_settings.GetReportShellCount();
 
-    auto d_min_limit_A = scaling_settings.GetHighResolutionLimit_A();
+    auto d_min_limit_A = d_min_override.has_value()
+        ? d_min_override : scaling_settings.GetHighResolutionLimit_A();
 
     std::unordered_map<uint64_t, float> reference_intensities;
     if (!reference.empty()) {

image_analysis/scale_merge/Merge.h

@@ -153,9 +153,13 @@ public:
     [[nodiscard]] std::vector<char> DeltaCChalfReject(const std::vector<IntegrationOutcome> &outcomes,
                                                       double nsigma) const;
 
+    // d_min_override, when set, is the effective high-resolution limit for the shell table (used for
+    // the automatic resolution cutoff computed by the caller); otherwise the manual
+    // ScalingSettings high-resolution limit stands. The number of shells is ScalingSettings::ReportShellCount.
     MergeStatistics MergeStats(const std::vector<MergedReflection> &merged,
                                const std::vector<IntegrationOutcome> &reflections,
-                               const std::vector<MergedReflection> &reference = {});
+                               const std::vector<MergedReflection> &reference = {},
+                               std::optional<double> d_min_override = std::nullopt);
 
     std::vector<MergedReflection> ExportReflections();
 };

image_analysis/scale_merge/ResolutionCutoff.cpp

@@ -0,0 +1,174 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "ResolutionCutoff.h"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
+#include "../../common/CorrelationCoefficient.h"
+
+namespace {
+    // Fine CC1/2 bins for the fit (finer than the 10 reported shells, per the design). A bin needs
+    // this many merged reflections for its CC1/2 to be trusted.
+    constexpr int   N_FIT_BINS = 25;
+    constexpr int   MIN_BIN_COUNT = 10;
+    constexpr int   MIN_FIT_BINS = 5;          // need at least this many usable bins for a fit
+    constexpr int   EXTEND_BINS_PAST_FALLOFF = 2;   // bins kept beyond the first sub-target bin
+    constexpr double SHELLS_FOR_EXTENSION = 10.0;   // "+1 shell" = one 10-shell width in s (report-independent)
+
+    double Logistic(double s, double k, double s0) {
+        return 1.0 / (1.0 + std::exp(k * (s - s0)));
+    }
+
+    // Weighted (equal-weight) sum of squared residuals of the logistic against the binned CC1/2.
+    double FitSSE(const std::vector<double> &s, const std::vector<double> &cc, double k, double s0) {
+        double sse = 0.0;
+        for (size_t i = 0; i < s.size(); ++i) {
+            const double r = cc[i] - Logistic(s[i], k, s0);
+            sse += r * r;
+        }
+        return sse;
+    }
+}
+
+ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector<MergedReflection> &merged,
+                                                   double cc_target, Logger &logger) {
+    ResolutionCutoffResult result;
+    if (!(cc_target > 0.0 && cc_target < 1.0)) {
+        result.note = "invalid CC target";
+        return result;
+    }
+
+    // s = 1/d^2 range over the merged reflections that carry a half-set pair.
+    double s_lo = std::numeric_limits<double>::max(), s_hi = 0.0;
+    double d_data_min = std::numeric_limits<double>::max();
+    for (const auto &m : merged) {
+        if (!(m.d > 0.0f) || !std::isfinite(m.I_half[0]) || !std::isfinite(m.I_half[1]))
+            continue;
+        const double s = 1.0 / (static_cast<double>(m.d) * m.d);
+        s_lo = std::min(s_lo, s);
+        s_hi = std::max(s_hi, s);
+        d_data_min = std::min(d_data_min, static_cast<double>(m.d));
+    }
+    if (!(s_lo < s_hi)) {
+        result.note = "no half-set data for CC1/2 fit";
+        return result;
+    }
+
+    // Bin CC1/2 against s (equal width in s, matching the reporting shells which are equal in 1/d^2).
+    const double bin_w = (s_hi - s_lo) / N_FIT_BINS;
+    std::vector<CorrelationCoefficient> bin_cc(N_FIT_BINS);
+    std::vector<int> bin_n(N_FIT_BINS, 0);
+    for (const auto &m : merged) {
+        if (!(m.d > 0.0f) || !std::isfinite(m.I_half[0]) || !std::isfinite(m.I_half[1]))
+            continue;
+        const double s = 1.0 / (static_cast<double>(m.d) * m.d);
+        int b = static_cast<int>((s - s_lo) / bin_w);
+        b = std::clamp(b, 0, N_FIT_BINS - 1);
+        bin_cc[b].Add(m.I_half[0], m.I_half[1]);
+        ++bin_n[b];
+    }
+
+    // Usable bins (enough counts), in ascending-s order, with their bin-centre s.
+    std::vector<double> s_bin, cc_bin;
+    for (int b = 0; b < N_FIT_BINS; ++b) {
+        if (bin_n[b] < MIN_BIN_COUNT) continue;
+        const double cc = bin_cc[b].GetCC();
+        if (!std::isfinite(cc)) continue;
+        s_bin.push_back(s_lo + (b + 0.5) * bin_w);
+        cc_bin.push_back(cc);
+    }
+    if (static_cast<int>(s_bin.size()) < MIN_FIT_BINS) {
+        result.note = "too few usable CC1/2 bins";
+        return result;
+    }
+
+    // Restrict to the contiguous fall-off from low res: keep bins up to a couple past the first one
+    // that drops below cc_target, so a high-res noise blip cannot pull the fit back up. If the lowest
+    // bin is already below cc_target there is no low-res plateau to anchor on - bail out.
+    if (cc_bin.front() < cc_target) {
+        result.note = "no low-resolution CC1/2 plateau";
+        return result;
+    }
+    size_t keep = cc_bin.size();
+    for (size_t i = 0; i < cc_bin.size(); ++i) {
+        if (cc_bin[i] < cc_target) {
+            keep = std::min(cc_bin.size(), i + 1 + EXTEND_BINS_PAST_FALLOFF);
+            break;
+        }
+    }
+    s_bin.resize(keep);
+    cc_bin.resize(keep);
+    if (static_cast<int>(s_bin.size()) < MIN_FIT_BINS) {
+        result.note = "too few CC1/2 bins in the fall-off region";
+        return result;
+    }
+
+    // Fit the logistic by a grid search over (k>0, s0) then a local coordinate-descent refine
+    // (dependency-free; the fall-off is smooth and the grid lands close). s0 spans the s range; k
+    // spans transitions from very gradual to very sharp relative to that range.
+    const double s_range = s_hi - s_lo;
+    double best_k = 0.0, best_s0 = 0.0, best_sse = std::numeric_limits<double>::max();
+    constexpr int N_S0 = 60, N_K = 40;
+    const double k_min = 2.0 / s_range, k_max = 200.0 / s_range;
+    for (int ik = 0; ik < N_K; ++ik) {
+        const double k = k_min * std::pow(k_max / k_min, static_cast<double>(ik) / (N_K - 1));
+        for (int is = 0; is < N_S0; ++is) {
+            const double s0 = s_lo + s_range * static_cast<double>(is) / (N_S0 - 1);
+            const double sse = FitSSE(s_bin, cc_bin, k, s0);
+            if (sse < best_sse) { best_sse = sse; best_k = k; best_s0 = s0; }
+        }
+    }
+
+    double k = best_k, s0 = best_s0;
+    double step_s0 = s_range / N_S0, step_k = best_k * 0.5;
+    for (int iter = 0; iter < 200; ++iter) {
+        bool improved = false;
+        for (const double ds : {step_s0, -step_s0}) {
+            const double sse = FitSSE(s_bin, cc_bin, k, s0 + ds);
+            if (sse < best_sse) { best_sse = sse; s0 += ds; improved = true; }
+        }
+        for (const double dk : {step_k, -step_k}) {
+            const double kt = k + dk;
+            if (kt <= 0.0) continue;
+            const double sse = FitSSE(s_bin, cc_bin, kt, s0);
+            if (sse < best_sse) { best_sse = sse; k = kt; improved = true; }
+        }
+        if (!improved) { step_s0 *= 0.5; step_k *= 0.5; }
+        if (step_s0 < 1e-6 * s_range && step_k < 1e-6 * best_k) break;
+    }
+
+    // s where the fitted CC1/2 crosses cc_target, then "one shell too far".
+    const double s_cross = s0 + std::log(1.0 / cc_target - 1.0) / k;
+    const double delta_s = s_range / SHELLS_FOR_EXTENSION;
+    const double s_final = s_cross + delta_s;
+
+    // No cut if the fall-off is beyond the measured edge (CC1/2 still healthy at the highest s).
+    if (s_final >= s_hi) {
+        result.note = "CC1/2 does not fall off within the measured range";
+        return result;
+    }
+    // Low-resolution floor: never cut into good low-res data. A fit that puts the cutoff within two
+    // shells of the lowest-res data is not a real fall-off - keep the full range and warn.
+    if (s_final <= s_lo + 2.0 * delta_s) {
+        logger.Warning("Resolution cutoff fit landed at low resolution (degenerate CC1/2 fall-off); "
+                       "keeping the full resolution range");
+        result.note = "degenerate low-resolution fit";
+        return result;
+    }
+
+    double d_cut = 1.0 / std::sqrt(s_final);
+    d_cut = std::max(d_cut, d_data_min);   // cannot cut beyond the highest-resolution reflection
+
+    // A cut that is not meaningfully coarser than the data edge is a no-op.
+    if (d_cut <= d_data_min * 1.001) {
+        result.note = "CC1/2 healthy to the detector edge";
+        return result;
+    }
+
+    result.d_cut = d_cut;
+    result.note = "CC1/2 logistic fall-off, +1 shell";
+    return result;
+}

image_analysis/scale_merge/ResolutionCutoff.h

@@ -0,0 +1,29 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <optional>
+#include <string>
+#include <vector>
+
+#include "../../common/Logger.h"
+#include "../../common/Reflection.h"   // MergedReflection
+
+// Automatic high-resolution cutoff from the CC1/2 fall-off of the merged half-sets (DIALS-style).
+// The merge itself, the error model and the per-image _process.h5 are left untouched - only the
+// written reflections and the reported shell table should be trimmed to the returned d_cut.
+//
+// Method (see docs/rugnux_resolution_cutoff_design.md): bin CC1/2 against s = 1/d^2 in fine bins,
+// fit a logistic CC1/2(s) = 1/(1+exp(k*(s-s0))) to the contiguous-from-low-res fall-off, take the s
+// where the fit crosses cc_target, then extend by one mean (10-shell) shell width in s ("one shell
+// too far", generous). d_cut is nullopt when the fit is degenerate (too few bins, flat/non-monotone)
+// or CC1/2 never falls below cc_target inside the measured range - the caller then keeps the full
+// range. cc_target in (0,1); merged must carry finite I_half[0]/I_half[1] to contribute.
+struct ResolutionCutoffResult {
+    std::optional<double> d_cut;   // high-resolution limit (A); nullopt => keep the full range
+    std::string note;              // human-readable description of the decision (for logging)
+};
+
+ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector<MergedReflection> &merged,
+                                                   double cc_target, Logger &logger);

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -16,6 +16,7 @@
 #include <gemmi/symmetry.hpp>
 
 #include "HKLKey.h"
+#include "ResolutionCutoff.h"
 #include "../../common/CorrelationCoefficient.h"
 #include "../../common/CrystalLattice.h"
 #include "../../common/Definitions.h"
@@ -159,6 +160,9 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment,
     rfree_fraction = s.GetRfreeFraction();
     scale_fulls = s.GetScaleFulls();
     scaling_iter = std::max(1, scaling_iterations);
+    resolution_cutoff_method = s.GetResolutionCutoff();
+    resolution_cc_target = s.GetResolutionCCTarget();
+    report_shell_count = s.GetReportShellCount();
     if (const auto forced = s.GetForcedMosaicity(); forced.has_value())
         mosaicity_deg = *forced;
     else
@@ -1033,6 +1037,34 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
         result.merged.push_back(mr);
     }
 
+    // Automatic high-resolution cutoff (post-merge): trim the written reflections + reported shells to
+    // the CC1/2 fall-off. The scaling, combine and error model above already ran over the full range,
+    // and the per-image _process.h5 is written elsewhere from the partials, so no data is lost. A manual
+    // --scaling-high-resolution (d_min_limit) wins; the P1 search merge (for_search) is never cut, so
+    // the space-group search still sees the full range.
+    std::optional<double> effective_d_min = d_min_limit;
+    if (!effective_d_min && !for_search && resolution_cutoff_method == ResolutionCutoffMethod::CCHalfLogistic) {
+        const auto rc = ComputeCCHalfLogisticCutoff(result.merged, resolution_cc_target, logger);
+        if (rc.d_cut) {
+            effective_d_min = rc.d_cut;
+            logger.Info("Auto resolution cutoff: {:.2f} A ({}; override with --scaling-high-resolution)",
+                        *rc.d_cut, rc.note);
+        } else {
+            logger.Info("Auto resolution cutoff: none ({}); keeping the full resolution range", rc.note);
+        }
+    }
+    if (effective_d_min) {
+        result.merged.erase(std::remove_if(result.merged.begin(), result.merged.end(),
+            [&](const MergedReflection &m) { return std::isfinite(m.d) && m.d < *effective_d_min; }),
+            result.merged.end());
+        // Recompute the merged resolution span for the R-free binning below over the trimmed set.
+        d_min = std::numeric_limits<float>::max(); d_max = 0.0f;
+        for (const auto &m : result.merged) {
+            if (!std::isfinite(m.d) || m.d <= 0.0f) continue;
+            d_min = std::min(d_min, m.d); d_max = std::max(d_max, m.d);
+        }
+    }
+
     if (rfree_fraction > 0.0 && !result.merged.empty() && d_min < d_max && d_min > 0.0f) {
         constexpr int n_shells = 20;
         ResolutionShells shells(d_min * 0.999f, d_max * 1.001f, n_shells);
@@ -1049,12 +1081,12 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
 
     if (reject_count > 0)
         logger.Info("Merge outlier rejection: dropped {} observations", reject_count);
-    // ---- Statistics (10 shells): completeness, multiplicity, <I/sigma>, R_meas, CC1/2. ----
-    constexpr int n_shells = 10;
+    // ---- Statistics (report_shell_count shells): completeness, multiplicity, <I/sigma>, R_meas, CC1/2. ----
+    const int n_shells = report_shell_count;
     float sd_min = std::numeric_limits<float>::max(), sd_max = 0.0f;
     for (const auto &m : result.merged) {
         if (!std::isfinite(m.d) || m.d <= 0.0f) continue;
-        if (d_min_limit && m.d < *d_min_limit) continue;
+        if (effective_d_min && m.d < *effective_d_min) continue;
         sd_min = std::min(sd_min, m.d); sd_max = std::max(sd_max, m.d);
     }
     if (!(sd_min < sd_max && sd_min > 0.0f))

image_analysis/scale_merge/RotationScaleMerge.h

@@ -98,6 +98,11 @@ private:
     bool scale_fulls = true;
     double mosaicity_deg = 0.1;
     float ice_half_width_q = 0.0f;
+    // Automatic high-resolution cutoff for the written reflections + reported shells (post-merge; the
+    // scaling, combine and error model always run over the full range). Manual d_min_limit wins.
+    ResolutionCutoffMethod resolution_cutoff_method = ResolutionCutoffMethod::Off;
+    double resolution_cc_target = 0.30;
+    int report_shell_count = 10;
 
     // Flat buffers, allocated once by Ingest() and reused across Run() calls.
     std::vector<Obs> partials;                    // all per-frame partials, grouped by frame

rugnux/Rugnux.cpp

@@ -28,6 +28,7 @@
 #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h"
 #include "../image_analysis/scale_merge/Merge.h"
 #include "../image_analysis/scale_merge/RotationScaleMerge.h"
+#include "../image_analysis/scale_merge/ResolutionCutoff.h"
 #include "../image_analysis/scale_merge/ScalingResult.h"
 #include "../image_analysis/scale_merge/SearchSpaceGroup.h"
 #include "../image_analysis/scale_merge/TwinningAnalysis.h"
@@ -499,8 +500,32 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) {
 
             ScaleMergeResult out;
             out.merged = merge_engine.ExportReflections();
+
+            // Automatic high-resolution cutoff (post-merge): trim the written reflections + reported
+            // shells to the CC1/2 fall-off. The merge, scaling and error model above ran over the full
+            // range, and the _process.h5 is written from the per-image outcomes, so no data is lost. A
+            // manual --scaling-high-resolution wins; the P1 search merge (for_search) is never cut so the
+            // space-group search still sees the full range. (Rotation is cut inside RotationScaleMerge.)
+            const auto &cut_ss = experiment_.GetScalingSettings();
+            std::optional<double> effective_d_min = cut_ss.GetHighResolutionLimit_A();
+            if (!effective_d_min && !for_search
+                && cut_ss.GetResolutionCutoff() == ResolutionCutoffMethod::CCHalfLogistic) {
+                const auto rc = ComputeCCHalfLogisticCutoff(out.merged, cut_ss.GetResolutionCCTarget(), logger);
+                if (rc.d_cut) {
+                    effective_d_min = rc.d_cut;
+                    logger.Info("Auto resolution cutoff: {:.2f} A ({}; override with --scaling-high-resolution)",
+                                *rc.d_cut, rc.note);
+                } else {
+                    logger.Info("Auto resolution cutoff: none ({}); keeping the full resolution range", rc.note);
+                }
+            }
+            if (effective_d_min)
+                out.merged.erase(std::remove_if(out.merged.begin(), out.merged.end(),
+                    [&](const MergedReflection &m) { return std::isfinite(m.d) && m.d < *effective_d_min; }),
+                    out.merged.end());
+
             phase("Computing statistics");
-            out.statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data);
+            out.statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data, effective_d_min);
             result.error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0;
             logger.Info("Merge complete ({} unique reflections, {})", out.merged.size(), label);
             return out;

tools/rugnux_cli.cpp

@@ -4,6 +4,7 @@
 #include <algorithm>
 #include <atomic>
 #include <chrono>
+#include <cmath>
 #include <csignal>
 #include <getopt.h>
 #include <iostream>
@@ -24,6 +25,7 @@
 #include "../image_analysis/scale_merge/Merge.h"
 #include "../image_analysis/scale_merge/ScaleOnTheFly.h"
 #include "../image_analysis/scale_merge/RotationScaleMerge.h"
+#include "../image_analysis/scale_merge/ResolutionCutoff.h"
 #include "../image_analysis/scale_merge/TwinningAnalysis.h"
 #include "../rugnux/Rugnux.h"
 
@@ -79,7 +81,10 @@ void print_usage() {
     std::cout << "   --smooth-g[=deg]                rot3d: smooth per-frame scale G over a deg-degree rotation range (XDS DELPHI-like) before the combine (default: 5 for rot3d; 0 = off)" << std::endl;
     std::cout << "   -A, --anomalous                 Anomalous mode (don't merge Friedel pairs)" << std::endl;
     std::cout << "   -B, --refine-bfactor            Refine per image B-factor" << std::endl;
-    std::cout << "   --scaling-high-resolution <num> High resolution limit for scaling/merging (default: no limit)" << std::endl;
+    std::cout << "   --scaling-high-resolution <num> High resolution limit for scaling/merging (manual override; default: no limit)" << std::endl;
+    std::cout << "   --resolution-cutoff <txt>       Automatic high-resolution cutoff for the written reflections + reported shells: cc-logistic|off (default: cc-logistic; ignored when --scaling-high-resolution is set)" << std::endl;
+    std::cout << "   --resolution-cc-target <num>    CC1/2 target defining the cc-logistic fall-off (default: 0.30)" << std::endl;
+    std::cout << "   --resolution-shells <num>       Number of resolution shells in the reported statistics table (default: 10)" << std::endl;
     std::cout << "   --min-partiality <num>          Minimum partiality to accept reflection (default: 0.02)" << std::endl;
     std::cout << "   --capture-uncertainty <num>     rot3d: systematic sigma ~num*(1-captured_fraction)*I on under-captured fulls (default: 1.0 for rot3d, 0 otherwise)" << std::endl;
     std::cout << "   --mosaicity <num>               Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl;
@@ -123,6 +128,9 @@ enum {
     OPT_MIN_IMAGE_CC,
     OPT_SCALING_ITERATIONS,
     OPT_SCALING_HIGH_RESOLUTION,
+    OPT_RESOLUTION_CUTOFF,
+    OPT_RESOLUTION_CC_TARGET,
+    OPT_RESOLUTION_SHELLS,
     OPT_SCALING_OUTPUT,
     OPT_SINGLE_PASS_ROTATION,
     OPT_REDO_ROTATION_SPOTS,
@@ -213,6 +221,9 @@ static option long_options[] = {
     {"min-image-cc", required_argument, nullptr, OPT_MIN_IMAGE_CC},
     {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS},
     {"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION},
+    {"resolution-cutoff", required_argument, nullptr, OPT_RESOLUTION_CUTOFF},
+    {"resolution-cc-target", required_argument, nullptr, OPT_RESOLUTION_CC_TARGET},
+    {"resolution-shells", required_argument, nullptr, OPT_RESOLUTION_SHELLS},
     {"scaling-output", required_argument, nullptr, OPT_SCALING_OUTPUT},
     {"bandwidth", required_argument, nullptr, OPT_BANDWIDTH},
     {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS},
@@ -448,6 +459,9 @@ int main(int argc, char **argv) {
 
     float d_min_spot_finding = 1.5;
     std::optional<float> d_min_scale_merge;
+    std::optional<ResolutionCutoffMethod> resolution_cutoff_method;   // --resolution-cutoff cc-logistic|off
+    std::optional<double> resolution_cc_target;                      // --resolution-cc-target
+    std::optional<int> report_shell_count;                           // --resolution-shells
     std::optional<std::string> integration_radius_arg; // "r1" or "r1,r2,r3"
     std::optional<IntegratorMode> integrator_mode;     // --integrator boxsum|gaussian|empirical
     std::optional<double> outlier_reject_nsigma;       // merge per-observation outlier rejection
@@ -683,6 +697,27 @@ int main(int argc, char **argv) {
         case OPT_SCALING_HIGH_RESOLUTION:
             d_min_scale_merge = atof(optarg);
             break;
+        case OPT_RESOLUTION_CUTOFF:
+            if (strcmp(optarg, "cc-logistic") == 0)
+                resolution_cutoff_method = ResolutionCutoffMethod::CCHalfLogistic;
+            else if (strcmp(optarg, "off") == 0)
+                resolution_cutoff_method = ResolutionCutoffMethod::Off;
+            else {
+                logger.Error("Invalid --resolution-cutoff value: {} (expected cc-logistic|off)", optarg);
+                print_usage();
+                exit(EXIT_FAILURE);
+            }
+            break;
+        case OPT_RESOLUTION_CC_TARGET:
+            resolution_cc_target = parse_double_arg(optarg, "--resolution-cc-target", logger);
+            break;
+        case OPT_RESOLUTION_SHELLS:
+            report_shell_count = atoi(optarg);
+            if (report_shell_count.value() < 1) {
+                logger.Error("Invalid --resolution-shells value: {} (must be >= 1)", report_shell_count.value());
+                exit(EXIT_FAILURE);
+            }
+            break;
         case OPT_SCALING_OUTPUT:
             if (strcmp(optarg, "mtz") == 0) {
                 intensity_format = IntensityFormat::MTZ;
@@ -882,6 +917,9 @@ int main(int argc, char **argv) {
         ScalingSettings scaling_settings;
         if (d_min_scale_merge)
             scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value());
+        if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method);
+        if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target);
+        if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count);
         scaling_settings.MergeFriedel(!anomalous_mode);
         scaling_settings.RefineB(refine_bfactor);
         scaling_settings.MinPartiality(min_partiality);
@@ -933,7 +971,29 @@ int main(int argc, char **argv) {
             for (size_t i = 0; i < reflections.size(); ++i)
                 merge_engine.AddImage(reflections[i], static_cast<int64_t>(i));
             merged_reflections = merge_engine.ExportReflections();
-            merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data);
+
+            // Automatic high-resolution cutoff (post-merge), matching the full-analysis path: a manual
+            // --scaling-high-resolution wins, otherwise trim the written reflections + reported shells
+            // to the CC1/2 fall-off. (Rotation is cut inside RotationScaleMerge above.)
+            const auto &cut_ss = experiment.GetScalingSettings();
+            std::optional<double> effective_d_min = cut_ss.GetHighResolutionLimit_A();
+            if (!effective_d_min
+                && cut_ss.GetResolutionCutoff() == ResolutionCutoffMethod::CCHalfLogistic) {
+                const auto rc = ComputeCCHalfLogisticCutoff(merged_reflections, cut_ss.GetResolutionCCTarget(),
+                                                            logger);
+                if (rc.d_cut) {
+                    effective_d_min = rc.d_cut;
+                    logger.Info("Auto resolution cutoff: {:.2f} A ({}; override with --scaling-high-resolution)",
+                                *rc.d_cut, rc.note);
+                }
+            }
+            if (effective_d_min)
+                merged_reflections.erase(std::remove_if(merged_reflections.begin(), merged_reflections.end(),
+                    [&](const MergedReflection &m) { return std::isfinite(m.d) && m.d < *effective_d_min; }),
+                    merged_reflections.end());
+
+            merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data,
+                                                        effective_d_min);
             error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0;
         }
 
@@ -1123,6 +1183,9 @@ int main(int argc, char **argv) {
     scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(rotation_indexing ? SMOOTH_G_DEFAULT_DEG : 0.0));
     if (d_min_scale_merge)
         scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value());
+    if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method);
+    if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target);
+    if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count);
     scaling_settings.MergeFriedel(!anomalous_mode);
     scaling_settings.RefineB(refine_bfactor);
     scaling_settings.MinPartiality(min_partiality);