rugnux: dataset Wilson B-factor estimate + robust per-image Wilson B

Add a dataset-wide isotropic Wilson B-factor estimate, the analogue of XDS's
"WILSON LINE ... B=" which we did not export. CalcGlobalWilsonB fits ln<I> vs
1/d^2 over the merged reflections (B = -2*slope), skipping the low-resolution
non-linear region (d > 4 A) and shells past the signal limit (<I/sigma> < 1)
so the estimate is insensitive to how far the merged data were carried. It is
diagnostic only - not fed back into scaling - and is written to the mmCIF
(_reflns.B_iso_Wilson_estimate), the printed merge statistics, and the log.

Also harden the per-image Wilson B (CalcWilsonBFactor): accept the fit only
when it is well-correlated and physically plausible (0 < B < 200 A^2), else
leave b_factor unset. A bad frame (an indexing glitch, too few reflections)
otherwise produced a wildly steep Wilson line and a B of several hundred A^2
that polluted the per-image plot; NaN is preferable to garbage.

Diagnostic-only: the merged intensities and every merge statistic are
byte-identical (verified baseline vs modified on the rotation battery).

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

image_analysis/WriteReflections.cpp

@@ -144,6 +144,10 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
         out << "_reflns.pdbx_Rrim_I_all           " << Fmt(ov.r_meas, 4) << "\n";
         out << "_reflns.pdbx_CC_half              " << Fmt(ov.cc_half, 4) << "\n";
         out << "_reflns.jfjoch_diffrn_ISa             " << CifStr(isa) << "   # asymptotic I/sigma (Diederichs)\n";
+        // Dataset-wide isotropic Wilson B-factor estimate (standard PDBx item), analogous to XDS's
+        // "WILSON LINE ... B=". Emitted only when the log-linear fit succeeded.
+        if (std::isfinite(statistics.wilson_b) && statistics.wilson_b > 0.0)
+            out << "_reflns.B_iso_Wilson_estimate     " << Fmt(statistics.wilson_b, 2) << "\n";
         // Twinning indicators (no standard mmCIF item; same jfjoch local prefix as ISa above).
         if (twinning.l_test_pairs > 0) {
             out << "_reflns.jfjoch_L_test_mean_abs_L      " << Fmt(twinning.mean_abs_l, 3)

image_analysis/bragg_integration/CalcISigma.cpp

@@ -1,10 +1,21 @@
 // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
 // SPDX-License-Identifier: GPL-3.0-only
 
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
 #include "CalcISigma.h"
 #include "Regression.h"
 #include "../../common/ResolutionShells.h"
 
+// Upper bound on a physically plausible macromolecular isotropic B-factor (A^2). A fit that lands
+// outside (0, WILSON_B_MAX) comes from a bad frame / bad dataset (too few or mis-indexed reflections,
+// a degenerate resolution range) rather than real Debye-Waller falloff, so it is rejected as
+// indeterminate rather than reported. Radiation-damaged / low-resolution data rarely exceeds ~120 A^2;
+// 200 leaves generous head-room while still excluding the hundreds-of-A^2 garbage.
+static constexpr float WILSON_B_MAX = 200.0f;
+
 void CalcISigma(DataMessage &msg) {
     CalcISigma(msg, msg.reflections);
 }
@@ -81,11 +92,80 @@ void CalcWilsonBFactor(DataMessage &msg,
 
     if (replace_b && valid_shells > 2) {
         auto reg_result = regression(shells_mean_one_over_d_square, log_I_mean, valid_shells);
+        const float b_est = -2.0f * reg_result.slope;
 
-        if (reg_result.r_square > 0.3)
-             msg.b_factor = -2.0f * reg_result.slope;
+        // Accept only a well-correlated, physically plausible fit. Occasional bad frames (an indexing
+        // glitch, too few reflections) produce a wildly steep Wilson line and a B of several hundred
+        // A^2 that pollutes the per-image plot; leaving b_factor unset (rendered as NaN) is better than
+        // emitting garbage.
+        if (reg_result.r_square > 0.3 && std::isfinite(b_est) && b_est > 0.0f && b_est < WILSON_B_MAX)
+             msg.b_factor = b_est;
     }
 
     msg.integration_B_logI = log_I_mean;
     msg.integration_B_one_over_d_square = shells_mean_one_over_d_square;
 }
+
+GlobalWilsonB CalcGlobalWilsonB(const std::vector<MergedReflection> &merged) {
+    GlobalWilsonB out;
+    // A dataset-wide estimate needs enough reflections to average the shell means; below this the
+    // per-image estimate is the only thing on offer and a global number would be meaningless.
+    if (merged.size() < 100)
+        return out;
+
+    float d_min = std::numeric_limits<float>::infinity(), d_max = 0.0f;
+    for (const auto &r : merged) {
+        if (std::isfinite(r.I) && r.d > 0.0f) {
+            d_min = std::min(d_min, r.d);
+            d_max = std::max(d_max, r.d);
+        }
+    }
+    if (!(d_min < d_max))
+        return out;
+
+    // Below ~4 A the Wilson plot is non-linear (bonding/solvent structure), so when the data extend to
+    // lower resolution than that, restrict the fit to d <= 4 A - the standard Wilson-B convention. If
+    // the whole dataset is coarser than 4 A, fall back to using all of it.
+    constexpr double WILSON_LOW_RES_LIMIT_A = 4.0;
+    const float d_low = (d_max > WILSON_LOW_RES_LIMIT_A && d_min < WILSON_LOW_RES_LIMIT_A)
+                        ? static_cast<float>(WILSON_LOW_RES_LIMIT_A) : d_max;
+
+    const int nshells = 20;
+    ResolutionShells shells(d_min, d_low, nshells);
+    std::vector<double> I_sum(nshells, 0.0), sig_sum(nshells, 0.0), count(nshells, 0.0);
+    for (const auto &r : merged) {
+        if (!std::isfinite(r.I) || r.d <= 0.0f)
+            continue;
+        auto s = shells.GetShell(r.d);          // reflections coarser than d_low fall outside -> skipped
+        if (s) {
+            I_sum[*s] += r.I;
+            if (std::isfinite(r.sigma) && r.sigma > 0.0f)
+                sig_sum[*s] += r.sigma;
+            ++count[*s];
+        }
+    }
+
+    const auto s2 = shells.GetShellMeanOneOverResSq();
+    std::vector<float> x, y;
+    for (int i = 0; i < nshells; ++i) {
+        // Skip empty / net-negative shells, and shells past the signal limit (mean I/sigma < 1). The
+        // latter keeps the fit out of the noise floor: without a resolution cut the weakest high-angle
+        // shells are background-residual-dominated and flatten the Wilson line, deflating B. This mirrors
+        // XDS/ctruncate fitting only over the meaningful range and makes the estimate insensitive to how
+        // far the merged data were carried.
+        if (count[i] > 0 && I_sum[i] > 0.0 && I_sum[i] > sig_sum[i]) {
+            x.push_back(s2[i]);
+            y.push_back(std::log(static_cast<float>(I_sum[i] / count[i])));
+        }
+    }
+    if (x.size() < 3)
+        return out;
+
+    const auto reg = regression(x, y, x.size());
+    const double b = -2.0 * reg.slope;          // <I> ~ exp(-2 B s^2), s^2 = 1/(4 d^2), x = 1/d^2
+    out.n_shells = static_cast<int>(x.size());
+    out.correlation = std::sqrt(std::clamp(static_cast<double>(reg.r_square), 0.0, 1.0));
+    if (std::isfinite(b) && b > 0.0 && b < WILSON_B_MAX)
+        out.b = b;
+    return out;
+}

image_analysis/bragg_integration/CalcISigma.h

@@ -3,6 +3,9 @@
 
 #pragma once
 
+#include <limits>
+#include <vector>
+
 #include "../../common/JFJochMessages.h"
 #include "../../common/Reflection.h"
 
@@ -12,3 +15,14 @@ void CalcWilsonBFactor(DataMessage &msg, bool replace_b = true);
 void CalcISigma(DataMessage &msg, const std::vector<Reflection> &reflections);
 void CalcWilsonBFactor(DataMessage &msg, const std::vector<Reflection> &reflections, bool replace_b = true);
 
+// Dataset-wide isotropic Wilson B-factor estimate from a log-linear fit of the shell-averaged merged
+// intensity against 1/d^2 (B = -2*slope). Robust because it averages over the whole dataset - the
+// per-image CalcWilsonBFactor above is a per-frame estimate and much noisier. Diagnostic only (like
+// XDS's "WILSON LINE ... B="), not used in scaling. Returns b = NaN when it cannot be determined.
+struct GlobalWilsonB {
+    double b = std::numeric_limits<double>::quiet_NaN();   // Wilson B-factor estimate (A^2)
+    double correlation = std::numeric_limits<double>::quiet_NaN();  // |correlation| of the log-linear fit
+    int n_shells = 0;                                      // resolution shells actually used
+};
+GlobalWilsonB CalcGlobalWilsonB(const std::vector<MergedReflection> &merged);
+

image_analysis/scale_merge/Merge.cpp

@@ -920,6 +920,9 @@ std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) {
     output << fmt::format("  {:>8s} ", "Overall");
     output << in.overall;
     output << std::endl;
+    if (std::isfinite(in.wilson_b) && in.wilson_b > 0.0)
+        output << fmt::format("  Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f})",
+                              in.wilson_b, in.wilson_b_correlation) << std::endl;
     output << std::endl;
     return output;
 }

image_analysis/scale_merge/Merge.h

@@ -51,6 +51,12 @@ struct MergeStatisticsShell {
 struct MergeStatistics {
     std::vector<MergeStatisticsShell> shells;
     MergeStatisticsShell overall;
+
+    // Dataset-wide isotropic Wilson B-factor estimate (A^2) from the log-linear fit of the shell-mean
+    // merged intensity against 1/d^2 (CalcGlobalWilsonB) - the analogue of XDS's "WILSON LINE ... B=".
+    // Diagnostic only; not used in scaling. NaN when not determined.
+    double wilson_b = NAN;
+    double wilson_b_correlation = NAN;
 };
 
 

rugnux/Rugnux.cpp

@@ -40,6 +40,7 @@
 #include "../image_analysis/scale_merge/TwinningAnalysis.h"
 #include "../image_analysis/scale_merge/HKLKey.h"
 #include "../image_analysis/WriteReflections.h"
+#include "../image_analysis/bragg_integration/CalcISigma.h"
 #include "../common/Definitions.h"
 #include "../common/CorrelationCoefficient.h"
 #include <array>
@@ -1110,6 +1111,17 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) {
             }
         }
 
+        // Dataset-wide Wilson B-factor estimate (like XDS's WILSON LINE B). Diagnostic only - it is not
+        // fed back into scaling; it just lands in the printed statistics, the mmCIF, and the log.
+        {
+            const GlobalWilsonB wilson = CalcGlobalWilsonB(sm.merged);
+            sm.statistics.wilson_b = wilson.b;
+            sm.statistics.wilson_b_correlation = wilson.correlation;
+            if (std::isfinite(wilson.b) && wilson.b > 0.0)
+                logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells)",
+                            wilson.b, wilson.correlation, wilson.n_shells);
+        }
+
         stats_text << sm.statistics;
         result.merge_statistics_text = stats_text.str();
         result.has_merge_statistics = true;