diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index 52ca4d79..a86492ae 100644 --- a/common/ScalingSettings.cpp +++ b/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"); diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index 958f6e31..d4c9ec6e 100644 --- a/common/ScalingSettings.h +++ b/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; diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index fa15475e..bed541e5 100644 --- a/docs/CHANGELOG.md +++ b/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 `.mtz` and `.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 `, 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 `` 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 diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 916e72e9..153c1f0f 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/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 &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 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 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(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::infinity(), pxmax = -pxmin, pymin = pxmin, pymax = -pxmin; + for (size_t i = 0; i < outcomes.size(); ++i) { + const int parity = static_cast(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(group_of.size()); + if (obs.size() < static_cast(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((m.px - pxmin) * sx), 0, NB - 1); + const int iy = std::clamp(static_cast((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 { + std::vector A(ncell, 1.0); + for (int it = 0; it < N_ITER; ++it) { + std::vector 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 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 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 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 &A) -> double { + std::vector 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 ident(ncell, 1.0); + const std::vector 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 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((m.px - pxmin) * sx), 0, NB - 1); + const int iy = std::clamp(static_cast((m.py - pymin) * sy), 0, NB - 1); + r.image_scale_corr = static_cast(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; diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 32f5c58a..bc468980 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/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 &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. diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 0e89b3db..d53331f0 100644 --- a/rugnux/Rugnux.cpp +++ b/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(), diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 1941a102..99430b3e 100644 --- a/rugnux/rugnux_cli.cpp +++ b/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 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 beam_x, beam_y, detector_distance_mm, wavelength_A, rot1_rad, rot2_rad, polarization_factor; std::optional 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 space_group_number; std::optional 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);