stills: optional detector-plane modulation (flat-field) correction

The stills scaling/merging path (ScaleOnTheFly + MergeOnTheFly) had no
correction surfaces, whereas rotation fits decay/absorption/modulation inside
RotationScaleMerge. Serial stills hammer the same detector regions shot after
shot, so a detector-plane (flat-field) systematic is at least as well
determined there as in rotation - it just had no home.

Add MergeOnTheFly::RefineModulation: a 16x16 multiplicative surface over the
predicted detector position, fit against the merged reference and folded into
each reflection's image_scale_corr before the error model and the merge. It
mirrors RotationScaleMerge::ApplyCellSurface - alternating per-cell fit with a
den-weighted geometric-mean gauge (never drifts the overall scale) and
Tikhonov pull to 1 - and is cross-validated the same way: fit the surface on
even images, score the held-out odd equivalents by a sigma-independent
R-meas-like agreement (a fractional metric a sigma-reshaping surface cannot
game), and apply the full-data surface only if the held-out gain clears a
margin. A no-op when the systematic is absent or the data too sparse (< 8 obs
per cell), so it is safe to leave off by default and opt in.

Wired through ScalingSettings::StillsModulation and rugnux --stills-modulation
(default off), in both the full-analysis and --scale paths. On a 17 MP
JUNGFRAU serial set the surface cross-validates with a large held-out gain,
confirming a real detector-plane systematic; whether it improves the merged
data quality is under evaluation.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

common/ScalingSettings.cpp

@@ -151,6 +151,15 @@ bool ScalingSettings::GetCorrectionSurfaces() const {
     return correction_surfaces;
 }
 
+ScalingSettings &ScalingSettings::StillsModulation(bool input) {
+    stills_modulation = input;
+    return *this;
+}
+
+bool ScalingSettings::GetStillsModulation() const {
+    return stills_modulation;
+}
+
 ScalingSettings &ScalingSettings::SmoothGDegrees(double input) {
     if (input < 0)
         throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Smooth-G range must be non-negative");

common/ScalingSettings.h

@@ -46,6 +46,10 @@ class ScalingSettings {
     bool correction_surfaces = true;
     // Absorption-surface refinement iteration count (used when correction_surfaces is on).
     int absorption_iter = 3;
+    // Detector-plane modulation (flat-field) correction for the STILLS merge path (ScaleOnTheFly +
+    // MergeOnTheFly), which otherwise has no correction surfaces. Off by default (experimental; the
+    // rotation path fits its own modulation via RotationScaleMerge). Enabled by rugnux --stills-modulation.
+    bool stills_modulation = false;
 
     // Smooth the per-frame scale G across frames (centered moving average of log G) before the rot3d
     // combine, so a rocking event's partials share a consistent scale. Given as a ROTATION RANGE in
@@ -80,6 +84,7 @@ public:
     ScalingSettings& ScaleFulls(bool input);
     ScalingSettings& AbsorptionIter(int input);
     ScalingSettings& CorrectionSurfaces(bool input);
+    ScalingSettings& StillsModulation(bool input);
     ScalingSettings& SmoothGDegrees(double input);
 
     ScalingSettings& RfreeFraction(double input);
@@ -117,6 +122,7 @@ public:
     [[nodiscard]] bool GetScaleFulls() const;
     [[nodiscard]] int GetAbsorptionIter() const;
     [[nodiscard]] bool GetCorrectionSurfaces() const;
+    [[nodiscard]] bool GetStillsModulation() const;
     [[nodiscard]] double GetSmoothGDegrees() const;
 
     [[nodiscard]] double GetRfreeFraction() const;

docs/CHANGELOG.md

@@ -23,6 +23,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as
 * rugnux: Write anomalous data as a standard CCP4 anomalous MTZ - one row per reflection with `IMEAN`, `I(+)`/`I(-)` (and the matching `F`/`F(+)`/`F(-)` amplitudes) - instead of two `IMEAN` rows per Bijvoet pair, so `aimless`/`mtz2sca`/ANODE read it directly. The non-anomalous MTZ output is unchanged.
 * rugnux: Always write the merged reflections as both `<prefix>.mtz` and `<prefix>.cif` (each has its uses downstream - MTZ for the CCP4/phenix tools, mmCIF for deposition). The `--scaling-output` format selector and the plain-text `.hkl` output are removed.
 * rugnux: Add a detector-plane **modulation** (flat-field) correction surface to rotation scaling - a smooth multiplicative factor over the detector position where a reflection lands, fitted against the merged reference (symmetry-equivalents land at different positions as the crystal rotates) and cross-validated. It corrects detector-response/geometric systematics that inflate R-meas; on the rotation battery it lowers R-meas by several to tens of percent (e.g. lysoC 23->16%, lyso_2 47->28%, EcwtAL500 53->28%) with CC1/2 held or improved and the anomalous signal preserved. On by default with decay/absorption (`--no-scaling-corrections` disables all). The correction-surface cross-validation now scores a sigma-independent R-meas-like agreement instead of a studentized chi^2, so a surface can no longer pass by reshaping the sigmas (which also hardens the absorption surface against a mis-indexed-data regression).
+* rugnux: Add an optional **stills detector-plane modulation** (flat-field) correction (`--stills-modulation`, default off). The stills scaling/merging path (ScaleOnTheFly + MergeOnTheFly) previously had no correction surfaces - unlike rotation, which fits decay/absorption/modulation inside RotationScaleMerge. This fits the same 16x16 detector-frame multiplicative surface over where each reflection lands (predicted x,y), against the merged reference, cross-validated (fit even images, score the held-out odd equivalents by a sigma-independent R-meas-like metric) so it is a no-op when the systematic is absent, and folds it into each reflection's scale before the error model and merge. Applies in both full analysis and `--scale`. Serial data hammers the same detector regions, so a flat-field systematic is well-determined; on a 17 MP JUNGFRAU serial set the surface cross-validates with a large held-out gain.
 * Bragg integration: The local background under each Bragg spot is now a **symmetric trimmed mean** of the background ring rather than the plain mean, on by default for monochromatic/rotation data (`--background-trim <f>`, default 0.10; 0 restores the plain mean). The plain ring mean reads high because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided contaminants, so it over-subtracts - and since a weak intensity is a small difference of large numbers, that bias is fractionally largest at the resolution edge. Dropping the lowest and highest fraction of ring pixels removes it while leaving a clean Poisson ring essentially unchanged. On the rotation battery `<I/sigma>` improved on every crystal (median +50%) and resolution-edge R-meas fell several-fold (e.g. lyso_ref 1.0 A 108->43%); the high-resolution CC1/2 is rescued where the plain mean had collapsed it (e.g. a thaumatin case regains ~0.5 A) at a small CC1/2 cost in already-clean shells. Implemented identically in the CPU and GPU integrators; stills keep their existing high-side sigma-clip (§9.2 of CPU_DATA_ANALYSIS.md). The dead, never-read `Reflection.completeness` field was removed.
 
 ### 1.0.0-rc.159

image_analysis/scale_merge/Merge.cpp

@@ -131,6 +131,137 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id
     }
 }
 
+double MergeOnTheFly::RefineModulation(std::vector<IntegrationOutcome> &outcomes) {
+    // A minimum held-out generalizing gain (fraction of the held-out scatter) before the surface is
+    // applied - a margin, so a noise-level "improvement" never engages the correction. Matches the
+    // rotation ApplyCellSurface gate.
+    constexpr double CV_MIN_RELATIVE_GAIN = 0.02;
+    constexpr int NB = 16;
+    const int ncell = NB * NB;
+
+    // One scaled observation reduced to what the surface fit needs. `parity` (image index & 1) drives the
+    // even/odd cross-validation split - the stills analogue of the rotation frame parity.
+    struct MObs { double I, sigma, corr; float px, py; int32_t group; int parity; int cell; };
+    std::vector<MObs> obs;
+
+    // Accept exactly what AddImage merges (systematic absence, scale/resolution/ice/partiality filters),
+    // and additionally require a finite detector position. Returns the dense ASU-group id or -1.
+    std::unordered_map<uint64_t, int> group_of;
+    auto accept = [&](const Reflection &r, MObs &out) -> bool {
+        if (generator.IsSystematicallyAbsent(r)) return false;
+        if (r.image_scale_corr <= 0.0f || !std::isfinite(r.image_scale_corr)) return false;
+        if (!AcceptReflection(r, high_resolution_limit)) return false;
+        if (exclude_ice_rings && r.on_ice_ring) return false;
+        if (IsMaskedRing(r)) return false;
+        if (r.partiality < min_partiality) return false;
+        if (!std::isfinite(r.predicted_x) || !std::isfinite(r.predicted_y)) return false;
+        const float I_corr = r.I * r.image_scale_corr, sigma_corr = r.sigma * r.image_scale_corr;
+        if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) return false;
+        const uint64_t key = generator(r).pack();
+        auto [it, inserted] = group_of.try_emplace(key, static_cast<int>(group_of.size()));
+        out.I = r.I; out.sigma = r.sigma; out.corr = r.image_scale_corr;
+        out.px = r.predicted_x; out.py = r.predicted_y; out.group = it->second;
+        return true;
+    };
+
+    float pxmin = std::numeric_limits<float>::infinity(), pxmax = -pxmin, pymin = pxmin, pymax = -pxmin;
+    for (size_t i = 0; i < outcomes.size(); ++i) {
+        const int parity = static_cast<int>(i & 1);
+        for (const auto &r : outcomes[i].reflections) {
+            MObs m{};
+            if (!accept(r, m)) continue;
+            m.parity = parity;
+            pxmin = std::min(pxmin, m.px); pxmax = std::max(pxmax, m.px);
+            pymin = std::min(pymin, m.py); pymax = std::max(pymax, m.py);
+            obs.push_back(m);
+        }
+    }
+    const int n_groups = static_cast<int>(group_of.size());
+    if (obs.size() < static_cast<size_t>(8 * ncell) || !(pxmax > pxmin) || !(pymax > pymin))
+        return 0.0;   // too sparse to over-determine a 16x16 surface, or degenerate detector footprint
+
+    const float sx = NB / (pxmax - pxmin), sy = NB / (pymax - pymin);
+    for (auto &m : obs) {
+        const int ix = std::clamp(static_cast<int>((m.px - pxmin) * sx), 0, NB - 1);
+        const int iy = std::clamp(static_cast<int>((m.py - pymin) * sy), 0, NB - 1);
+        m.cell = ix * NB + iy;
+    }
+
+    // Fit the per-cell factor over {parity subset} (parity < 0 = all obs), n_iter alternating rounds against
+    // that subset's own inverse-variance reference: Tikhonov pull to 1, gauge-fixed to a den-weighted
+    // geometric mean of 1 so it never drifts the overall scale. Mirrors RotationScaleMerge::ApplyCellSurface.
+    constexpr int N_ITER = 3;
+    auto fit_surface = [&](int parity) -> std::vector<double> {
+        std::vector<double> A(ncell, 1.0);
+        for (int it = 0; it < N_ITER; ++it) {
+            std::vector<double> sw(n_groups, 0.0), swI(n_groups, 0.0);
+            for (const auto &o : obs) {
+                if (parity >= 0 && o.parity != parity) continue;
+                const double a = A[o.cell], sc = o.sigma * o.corr * a, w = 1.0 / (sc * sc);
+                sw[o.group] += w; swI[o.group] += w * o.I * o.corr * a;
+            }
+            std::vector<double> num(ncell, 0.0), den(ncell, 0.0);
+            for (const auto &o : obs) {
+                if ((parity >= 0 && o.parity != parity) || sw[o.group] <= 0.0) continue;
+                const double Iref = swI[o.group] / sw[o.group], a = A[o.cell];
+                const double Is = o.I * o.corr * a, sc = o.sigma * o.corr * a;
+                if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0) || !(sc > 0.0)) continue;
+                const double w = 1.0 / (sc * sc);
+                num[o.cell] += w * Is * Iref; den[o.cell] += w * Is * Is;
+            }
+            std::vector<double> dsorted = den;
+            std::nth_element(dsorted.begin(), dsorted.begin() + dsorted.size() / 2, dsorted.end());
+            const double lambda = 0.1 * std::max(1e-30, dsorted[dsorted.size() / 2]);
+            double logsum = 0.0, wsum = 0.0;
+            std::vector<double> upd(ncell, 1.0);
+            for (int c = 0; c < ncell; ++c) upd[c] = (num[c] + lambda) / (den[c] + lambda);
+            for (int c = 0; c < ncell; ++c) if (den[c] > 0.0) { logsum += den[c] * std::log(upd[c]); wsum += den[c]; }
+            const double gm = wsum > 0.0 ? std::exp(logsum / wsum) : 1.0;
+            for (int c = 0; c < ncell; ++c) A[c] = std::clamp(A[c] * upd[c] / gm, 0.25, 4.0);
+        }
+        return A;
+    };
+    // Sigma-independent (R-meas-like) agreement of the held-out equivalents: sum|Is - Iref| / sum|Iref|.
+    // A fractional metric cannot be gamed by a surface that merely reshapes sigma via corr.
+    auto score = [&](int parity, const std::vector<double> &A) -> double {
+        std::vector<double> sw(n_groups, 0.0), swI(n_groups, 0.0);
+        for (const auto &o : obs) {
+            if (o.parity != parity) continue;
+            const double a = A[o.cell], Is = o.I * o.corr * a, sc = o.sigma * o.corr * a, w = 1.0 / (sc * sc);
+            sw[o.group] += w; swI[o.group] += w * Is;
+        }
+        double num = 0.0, den = 0.0;
+        for (const auto &o : obs) {
+            if (o.parity != parity || sw[o.group] <= 0.0) continue;
+            const double a = A[o.cell], Is = o.I * o.corr * a, Iref = swI[o.group] / sw[o.group];
+            if (!std::isfinite(Iref) || Iref <= 0.0) continue;
+            num += std::abs(Is - Iref); den += Iref;
+        }
+        return den > 0.0 ? num / den : 0.0;
+    };
+
+    // Cross-validate: fit on even images, score the held-out odd equivalents (and vice versa). Apply the
+    // full-data surface only if the held-out agreement improves by a clear margin.
+    const std::vector<double> ident(ncell, 1.0);
+    const std::vector<double> A_even = fit_surface(0), A_odd = fit_surface(1);
+    const double base = score(1, ident) + score(0, ident);
+    const double gain = base - (score(1, A_even) + score(0, A_odd));
+    if (!(gain > CV_MIN_RELATIVE_GAIN * base))
+        return 0.0;   // not cross-validated: the correction stays a no-op (caller logs)
+    const std::vector<double> A = fit_surface(-1);
+
+    // Fold the surface into each accepted reflection's image_scale_corr (recompute its cell deterministically).
+    for (auto &outcome : outcomes)
+        for (auto &r : outcome.reflections) {
+            MObs m{};
+            if (!accept(r, m)) continue;
+            const int ix = std::clamp(static_cast<int>((m.px - pxmin) * sx), 0, NB - 1);
+            const int iy = std::clamp(static_cast<int>((m.py - pymin) * sy), 0, NB - 1);
+            r.image_scale_corr = static_cast<float>(r.image_scale_corr * A[ix * NB + iy]);
+        }
+    return gain / std::max(base, 1e-30);   // held-out gain fraction (caller logs)
+}
+
 float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl_key, float partiality) const {
     if (!error_model_active)
         return sigma_corr;

image_analysis/scale_merge/Merge.h

@@ -147,6 +147,16 @@ public:
     // order (or threading) of AddImage calls - not a draw from a shared RNG in call order.
     void AddImage(const IntegrationOutcome& outcome, int64_t image_id, bool cc_mask = false);
 
+    // Detector-plane modulation (flat-field) correction for the STILLS path: fit a smooth multiplicative
+    // factor over where a reflection lands on the detector (predicted x,y) against the merged reference and
+    // fold it into each reflection's image_scale_corr, so the following AddImage merge (and the error model)
+    // see the corrected scale. A 16x16 grid, cross-validated (fit even images, score the held-out odd
+    // equivalents by a sigma-independent R-meas-like metric) so it is a no-op when the systematic is absent
+    // or the data too sparse. The stills analogue of RotationScaleMerge::RefineModulation - the stills
+    // ScaleOnTheFly/MergeOnTheFly path otherwise has no correction surfaces. Mutates `outcomes`; call after
+    // scaling and before RefineErrorModel/AddImage. Returns the applied held-out gain fraction (0 = no-op).
+    double RefineModulation(std::vector<IntegrationOutcome> &outcomes);
+
     // Per-crystal CC1/2-delta rejection (CrystFEL deltaCChalf): returns a per-image flag
     // marking images whose removal would raise CC1/2 by a low-side outlier amount
     // (deltaCChalf < mean - nsigma*stddev). Skip the flagged images when merging.

rugnux/Rugnux.cpp

@@ -789,6 +789,17 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) {
             merge_engine.MaskIceRings(masked_ice_rings, config_.spot_finding.ice_ring_width_Q_recipA);
             if (result.consensus_cell.has_value())
                 merge_engine.ReferenceCell(*result.consensus_cell);
+            // Optional detector-plane modulation (flat-field) correction for stills, folded into each
+            // reflection's scale before the error model and merge (never for the P1 search pass). The
+            // rotation path fits its own modulation inside RotationScaleMerge.
+            if (experiment_.GetScalingSettings().GetStillsModulation() && !for_search) {
+                const double mod_gain = merge_engine.RefineModulation(indexer->GetIntegrationOutcome());
+                if (mod_gain > 0.0)
+                    logger.Info("Stills modulation: detector-frame 16x16 surface applied "
+                                "(cross-validated, held-out gain {:.1f}%)", 100.0 * mod_gain);
+                else
+                    logger.Info("Stills modulation: no cross-validated gain (skipped)");
+            }
             merge_engine.RefineErrorModel(merge_input);
             if (merge_engine.ErrorModelActive())
                 logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", merge_engine.ErrorModelA(),

rugnux/rugnux_cli.cpp

@@ -88,6 +88,7 @@ void print_usage() {
     std::cout << "   --write-process-h5              Also write the (large) _process.h5 when merging (default: only .mtz/.cif when merging)" << std::endl;
     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 << "   --no-scaling-corrections        rot3d: disable the (default-on) decay + absorption correction surfaces fitted on the fulls after scale-fulls" << std::endl;
+    std::cout << "   --stills-modulation             stills: fit a detector-plane modulation (flat-field) surface over the merged reflections (cross-validated; experimental, default 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 (stills only)" << std::endl;
     std::cout << "   --scaling-high-resolution <num> High resolution limit for scaling/merging (manual override; default: no limit)" << std::endl;
@@ -167,6 +168,7 @@ enum {
     OPT_MOSAICITY,
     OPT_SMOOTH_G,
     OPT_NO_SCALING_CORRECTIONS,
+    OPT_STILLS_MODULATION,
     OPT_DETECT_ICE_RINGS,
     OPT_NO_SCALE_FULLS,
     OPT_WRITE_PROCESS_H5,
@@ -212,6 +214,7 @@ static option long_options[] = {
     {"write-process-h5", no_argument, nullptr, OPT_WRITE_PROCESS_H5},
     {"smooth-g", optional_argument, nullptr, OPT_SMOOTH_G},
     {"no-scaling-corrections", no_argument, nullptr, OPT_NO_SCALING_CORRECTIONS},
+    {"stills-modulation", no_argument, nullptr, OPT_STILLS_MODULATION},
     {"refine", required_argument, nullptr, 'r'},
 
     {"two-pass-rotation", optional_argument, nullptr, 'R'},
@@ -491,6 +494,7 @@ int main(int argc, char **argv) {
     std::optional<float> beam_x, beam_y, detector_distance_mm, wavelength_A, rot1_rad, rot2_rad, polarization_factor;
     std::optional<double> smooth_g_deg_arg;   // --smooth-g[=deg]; default 5 deg for rot3d, 0 (off) otherwise
     bool no_scaling_corrections = false;   // --no-scaling-corrections: disable rot3d decay+absorption surfaces
+    bool stills_modulation_flag = false;   // --stills-modulation: detector-plane flat-field surface for stills
     bool anomalous_mode = false;
     std::optional<int64_t> space_group_number;
     std::optional<UnitCell> fixed_reference_unit_cell;
@@ -769,6 +773,9 @@ int main(int argc, char **argv) {
         case OPT_SMOOTH_G:
             smooth_g_deg_arg = optarg ? parse_double_arg(optarg, "--smooth-g", logger) : SMOOTH_G_DEFAULT_DEG;
             break;
+        case OPT_STILLS_MODULATION:
+            stills_modulation_flag = true;
+            break;
         case OPT_NO_SCALING_CORRECTIONS:
             no_scaling_corrections = true;
             break;
@@ -1050,6 +1057,7 @@ int main(int argc, char **argv) {
         scaling_settings.MinCapturedFraction(min_captured_fraction_arg.value_or(
             (experiment.GetGoniometer().has_value() && !force_still) ? 0.7 : 0.0));
         scaling_settings.MinCCForImage(min_image_cc / 100.0);   // --min-image-cc is percent; the setting is a fraction
+        scaling_settings.StillsModulation(stills_modulation_flag);
         scaling_settings.OutlierRejectNsigma(
             outlier_reject_nsigma.value_or(
                 (experiment.GetGoniometer().has_value() && !force_still) ? REJECT_OUTLIERS_DEFAULT_NSIGMA : 0.0));
@@ -1096,6 +1104,16 @@ int main(int argc, char **argv) {
             ScaleOnTheFly(experiment, MergeAll(experiment, reflections)).Scale(reflections, nthreads);
             MergeOnTheFly merge_engine(experiment);
             merge_engine.ReferenceCell(experiment.GetUnitCell());
+            // Optional detector-plane modulation (flat-field) correction, folded into each reflection's
+            // scale before the error model + merge (mirrors the full-analysis stills path in Rugnux.cpp).
+            if (experiment.GetScalingSettings().GetStillsModulation()) {
+                const double mod_gain = merge_engine.RefineModulation(reflections);
+                if (mod_gain > 0.0)
+                    logger.Info("Stills modulation: detector-frame 16x16 surface applied "
+                                "(cross-validated, held-out gain {:.1f}%)", 100.0 * mod_gain);
+                else
+                    logger.Info("Stills modulation: no cross-validated gain (skipped)");
+            }
             // Fit the (a, b) error model from symmetry-mate scatter before merging, exactly as the full
             // pipeline does (Rugnux.cpp). Without this the offline --scale merge would use the identity
             // model and produce much worse stills intensities (no (b*I)^2 systematic term, no sigma floor).
@@ -1339,6 +1357,7 @@ int main(int argc, char **argv) {
     scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(rotation_indexing ? SMOOTH_G_DEFAULT_DEG : 0.0));
     if (no_scaling_corrections)
         scaling_settings.CorrectionSurfaces(false);
+    scaling_settings.StillsModulation(stills_modulation_flag);
     if (d_min_scale_merge)
         scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value());
     if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method);