diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index 288a447e..08ff80c6 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -160,6 +160,15 @@ bool ScalingSettings::GetStillsModulation() const { return stills_modulation; } +ScalingSettings &ScalingSettings::StillsPartialityRefine(bool input) { + stills_partiality_refine = input; + return *this; +} + +bool ScalingSettings::GetStillsPartialityRefine() const { + return stills_partiality_refine; +} + ScalingSettings &ScalingSettings::ExpectedVarianceMerge(bool input) { expected_variance_merge = input; return *this; diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index c541cc61..f215709f 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -51,6 +51,13 @@ class ScalingSettings { // rotation path fits its own modulation via RotationScaleMerge). Enabled by rugnux --stills-modulation. bool stills_modulation = false; + // Physical partiality post-refinement for the STILLS merge (StillsPartialityRefine): refine a per-crystal + // orientation tilt against the running merge, recompute each reflection's partiality from the refined + // geometry (angular Ewald-proximity model), and re-scale/merge - the "full model" for stills. ON by + // default (helps mono stills, neutral on pink beam, tames weak data via a soft prior). rugnux + // --simple-stills turns it OFF, reverting to treating every reflection as a full (p = 1, single pass). + bool stills_partiality_refine = true; + // Expected-variance merge weighting for the STILLS merge (MergeOnTheFly). When combining a reflection's // redundant observations by inverse variance, rebuild the Poisson signal part of each observation's // variance at the reflection's EXPECTED instead of the observation's own intensity. Weighting by an @@ -101,6 +108,7 @@ public: ScalingSettings& AbsorptionIter(int input); ScalingSettings& CorrectionSurfaces(bool input); ScalingSettings& StillsModulation(bool input); + ScalingSettings& StillsPartialityRefine(bool input); ScalingSettings& ExpectedVarianceMerge(bool input); ScalingSettings& SmoothGDegrees(double input); ScalingSettings& RelativeBDegrees(double input); @@ -141,6 +149,7 @@ public: [[nodiscard]] int GetAbsorptionIter() const; [[nodiscard]] bool GetCorrectionSurfaces() const; [[nodiscard]] bool GetStillsModulation() const; + [[nodiscard]] bool GetStillsPartialityRefine() const; [[nodiscard]] bool GetExpectedVarianceMerge() const; [[nodiscard]] double GetSmoothGDegrees() const; [[nodiscard]] double GetRelativeBDegrees() const; diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index 1c172a2f..338b1290 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -7,6 +7,8 @@ ADD_LIBRARY(JFJochScaleMerge Merge.h ScaleOnTheFly.cpp ScaleOnTheFly.h + StillsPartialityRefine.cpp + StillsPartialityRefine.h RotationScaleMerge.cpp RotationScaleMerge.h ResolutionCutoff.cpp diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp new file mode 100644 index 00000000..03f60189 --- /dev/null +++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp @@ -0,0 +1,329 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "StillsPartialityRefine.h" + +#include +#include +#include +#include +#include + +#include +#include + +#include "Merge.h" + +namespace { + constexpr size_t MIN_FIT_REFLECTIONS = 20; + constexpr double kRadToDeg = 180.0 / 3.14159265358979323846; + + double SafeInv(double x, double fallback) { + if (!std::isfinite(x) || x == 0.0) + return fallback; + return 1.0 / x; + } + + // One accepted reflection reduced to the physical partiality fit: the base reciprocal vector q (crystal + // frame at the stored per-image orientation), the reference full intensity, the measured intensity, the + // Lorentz factor and the weight. dist_ewald / partiality are recomputed from q under the refined tilt. + struct FitObs { + double qx, qy, qz; + double Iref; + double Iobs; + double lp; // 1 / rlp + double weight; // 1 / sigma + }; + + double VecLen(double x, double y, double z) { return std::sqrt(x * x + y * y + z * z); } + + // Analytic partiality for a reflection whose base reciprocal vector is q, tilted by (psi_x, psi_y). + // Mirrors BraggPrediction: dist_ewald = |S| - 1/lambda with S = q_rot + S0, and + // p = exp(-dist_ewald^2 / 2 sigma^2), sigma^2 = (gamma0 + gamma_e*d*)^2 + (bw*|q_z|)^2. + double ComputeP(double qx, double qy, double qz, + double psi_x, double psi_y, + double s0x, double s0y, double s0z, double inv_lambda, + double gamma0, double gamma_e, double bw) { + double aa[3] = {psi_x, psi_y, 0.0}; + double q[3] = {qx, qy, qz}; + double qr[3]; + ceres::AngleAxisRotatePoint(aa, q, qr); + const double Sx = qr[0] + s0x, Sy = qr[1] + s0y, Sz = qr[2] + s0z; + const double de = std::sqrt(Sx * Sx + Sy * Sy + Sz * Sz) - inv_lambda; + const double dstar = std::sqrt(qr[0] * qr[0] + qr[1] * qr[1] + qr[2] * qr[2]); + const double sig = gamma0 + gamma_e * dstar; + const double sbw = bw * std::fabs(qr[2]); + const double sig2 = sig * sig + sbw * sbw; + if (!(sig2 > 0.0)) + return 1.0; + return std::exp(-0.5 * de * de / sig2); + } + + // Robust per-crystal scale G (linear in G given the model coefficients), identical objective to + // ScaleOnTheFly::SolveScaleIRLS: minimise sum Cauchy_k( w (G*coeff - Iobs) ) over G >= 0. + double SolveScaleIRLS(const std::vector &coeff, const std::vector &Iobs, + const std::vector &weight, double robust_k) { + auto weighted_scale = [&](auto robust_weight) { + double num = 0.0, den = 0.0; + for (size_t i = 0; i < coeff.size(); ++i) { + const double rw = robust_weight(i); + const double w2 = weight[i] * weight[i]; + num += rw * w2 * coeff[i] * Iobs[i]; + den += rw * w2 * coeff[i] * coeff[i]; + } + return den > 0.0 ? num / den : NAN; + }; + + double G = weighted_scale([](size_t) { return 1.0; }); + if (!std::isfinite(G)) + return 1.0; + G = std::max(0.0, G); + + const double k2 = robust_k * robust_k; + for (int iter = 0; iter < 30; ++iter) { + const double G_prev = G; + const double G_next = weighted_scale([&](size_t i) { + const double res = weight[i] * (G * coeff[i] - Iobs[i]); + return 1.0 / (1.0 + res * res / k2); + }); + if (!std::isfinite(G_next)) + break; + G = std::max(0.0, G_next); + if (std::abs(G - G_prev) <= 1e-7 * std::max(G, 1.0)) + break; + } + return G; + } + + // Ceres residual: refine the orientation tilt (psi_x, psi_y) holding the scale G fixed. The tilt + // rotates the base reciprocal vector q; partiality follows analytically. Residual is the intensity + // mismatch weighted by 1/sigma, exactly matching ScaleOnTheFly's intensity-space objective. + struct PsiResidual { + double qx, qy, qz; + double s0x, s0y, s0z, inv_lambda; + double gamma0, gamma_e, bw; + double G, lp, Iref, Iobs, weight; + + template + bool operator()(const T *const psi, T *residual) const { + T q[3] = {T(qx), T(qy), T(qz)}; + T aa[3] = {psi[0], psi[1], T(0.0)}; + T qr[3]; + ceres::AngleAxisRotatePoint(aa, q, qr); + const T Sx = qr[0] + T(s0x), Sy = qr[1] + T(s0y), Sz = qr[2] + T(s0z); + const T de = ceres::sqrt(Sx * Sx + Sy * Sy + Sz * Sz) - T(inv_lambda); + const T dstar = ceres::sqrt(qr[0] * qr[0] + qr[1] * qr[1] + qr[2] * qr[2]); + const T sig = T(gamma0) + T(gamma_e) * dstar; + const T sbw = T(bw) * ceres::abs(qr[2]); + const T sig2 = sig * sig + sbw * sbw; + const T p = ceres::exp(T(-0.5) * de * de / sig2); + residual[0] = T(weight) * (T(G) * p * T(lp) * T(Iref) - T(Iobs)); + return true; + } + }; + + // Gaussian prior N(0, sigma_prior^2) on the tilt. The data residuals above are (model-obs)/sigma, a + // proper chi^2, so the MAP prior residual is simply dpsi/sigma_prior - no scale calibration needed. + struct PsiPrior { + double inv_sigma; + template + bool operator()(const T *const psi, T *residual) const { + residual[0] = T(inv_sigma) * psi[0]; + residual[1] = T(inv_sigma) * psi[1]; + return true; + } + }; +} + +StillsPartialityRefine::StillsPartialityRefine(const DiffractionExperiment &x) + : StillsPartialityRefine(x, Settings{}) {} + +StillsPartialityRefine::StillsPartialityRefine(const DiffractionExperiment &x, Settings settings) + : experiment_(x), + settings_(settings), + hkl_key_generator_(x.GetScalingSettings().GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)), + d_min_limit_(x.GetScalingSettings().GetHighResolutionLimit_A()), + bandwidth_sigma_(x.GetBandwidthFWHM().value_or(0.0f) / 2.3548f) {} + +double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, + const std::map &reference) const { + if (outcome.reflections.empty()) + return 0.0; + + const Coord Astar = outcome.latt.Astar(); + const Coord Bstar = outcome.latt.Bstar(); + const Coord Cstar = outcome.latt.Cstar(); + const Coord S0 = outcome.geom.GetScatteringVector(); + const double inv_lambda = 1.0 / outcome.geom.GetWavelength_A(); + const double bw = bandwidth_sigma_; + const double gamma0 = 0.0; // width is purely angular: sigma(d*) = gamma_e * d* (set per crystal below) + + auto base_q = [&](const Reflection &r) { + return Astar * static_cast(r.h) + Bstar * static_cast(r.k) + + Cstar * static_cast(r.l); + }; + + // Collect the reflections that constrain the fit (accepted, non-ice, finite, present in the reference). + std::vector obs; + obs.reserve(outcome.reflections.size()); + double sum_ang2 = 0.0; // RMS angular excitation error (dist_ewald / d*) -> per-crystal mosaic width + size_t n_de = 0; + for (const Reflection &r: outcome.reflections) { + if (r.on_ice_ring || !AcceptReflection(r, d_min_limit_)) + continue; + if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f) + continue; + const auto it = reference.find(hkl_key_generator_(r)); + if (it == reference.end() || !std::isfinite(it->second)) + continue; + + const Coord q = base_q(r); + obs.push_back(FitObs{ + .qx = q.x, .qy = q.y, .qz = q.z, + .Iref = it->second, + .Iobs = static_cast(r.I), + .lp = SafeInv(r.rlp, 1.0), + .weight = SafeInv(r.sigma, 1.0), + }); + + // Angular excitation error delta_psi = dist_ewald / d* at the stored orientation (psi = 0). Using + // the ANGULAR distance (not the linear reciprocal-space distance) makes the partiality width + // resolution-clean: a fixed mosaic angle smears high-resolution rlps more in reciprocal space, so a + // constant linear width computes p too small at high resolution and over-divides those shells. + const double dstar = VecLen(q.x, q.y, q.z); + const double de0 = VecLen(q.x + S0.x, q.y + S0.y, q.z + S0.z) - inv_lambda; + if (dstar > 1e-9) { + const double dpsi = de0 / dstar; + sum_ang2 += dpsi * dpsi; + ++n_de; + } + } + + if (obs.size() < MIN_FIT_REFLECTIONS || n_de == 0) + return 0.0; + + // Per-crystal angular mosaic width from the RMS angular excitation error. sigma(d*) = gamma_e * d* + // (gamma0 = 0), i.e. p = exp(-0.5 (delta_psi / gamma_e)^2) is a Gaussian in the angular distance from + // the Ewald sphere - the physical mosaic/divergence model, independent of resolution. A positive + // settings_.gamma_e overrides the per-crystal estimate with a shared (pooled) width. + const double gamma_e_ang = std::max(std::sqrt(sum_ang2 / static_cast(n_de)), 1e-9); + const double gamma_e = settings_.gamma_e > 0.0 ? settings_.gamma_e : gamma_e_ang; + + double psi[2] = {0.0, 0.0}; + double G = 1.0; + + const bool refine_tilt = obs.size() >= settings_.min_reflections; + const int inner = refine_tilt ? settings_.inner_iterations : 1; + + for (int it = 0; it < inner; ++it) { + // (1) Solve G given the current partialities. + std::vector coeff(obs.size()), Iobs(obs.size()), weight(obs.size()); + for (size_t j = 0; j < obs.size(); ++j) { + const double p = ComputeP(obs[j].qx, obs[j].qy, obs[j].qz, psi[0], psi[1], + S0.x, S0.y, S0.z, inv_lambda, gamma0, gamma_e, bw); + coeff[j] = p * obs[j].lp * obs[j].Iref; + Iobs[j] = obs[j].Iobs; + weight[j] = obs[j].weight; + } + G = SolveScaleIRLS(coeff, Iobs, weight, settings_.robust_k); + if (!(G > 0.0) || !std::isfinite(G)) + return 0.0; + + if (!refine_tilt) + break; + + // (2) Refine the tilt holding G fixed. + ceres::Problem problem; + for (const auto &o: obs) { + auto *cost = new ceres::AutoDiffCostFunction(new PsiResidual{ + .qx = o.qx, .qy = o.qy, .qz = o.qz, + .s0x = S0.x, .s0y = S0.y, .s0z = S0.z, .inv_lambda = inv_lambda, + .gamma0 = gamma0, .gamma_e = gamma_e, .bw = bw, + .G = G, .lp = o.lp, .Iref = o.Iref, .Iobs = o.Iobs, .weight = o.weight}); + problem.AddResidualBlock(cost, new ceres::CauchyLoss(settings_.robust_k), psi); + } + if (settings_.prior_sigma_deg > 0.0) { + const double inv_sigma = kRadToDeg / settings_.prior_sigma_deg; // 1 / sigma_prior (rad) + problem.AddResidualBlock(new ceres::AutoDiffCostFunction( + new PsiPrior{inv_sigma}), nullptr, psi); + } + problem.SetParameterLowerBound(psi, 0, -settings_.max_tilt_rad); + problem.SetParameterUpperBound(psi, 0, settings_.max_tilt_rad); + problem.SetParameterLowerBound(psi, 1, -settings_.max_tilt_rad); + problem.SetParameterUpperBound(psi, 1, settings_.max_tilt_rad); + + ceres::Solver::Options options; + options.linear_solver_type = ceres::DENSE_QR; + options.minimizer_progress_to_stdout = false; + options.num_threads = 1; + options.max_num_iterations = 25; + ceres::Solver::Summary summary; + ceres::Solve(options, &problem, &summary); + } + + // Write the refined partiality + scale correction onto every reflection of the crystal (not only the + // fit subset), so the merge sees a consistent model. image_scale_corr = rlp / (partiality * G). + for (auto &r: outcome.reflections) { + const Coord q = base_q(r); + const double p = ComputeP(q.x, q.y, q.z, psi[0], psi[1], S0.x, S0.y, S0.z, inv_lambda, + gamma0, gamma_e, bw); + r.partiality = static_cast(p); + const double denom = p * G; + r.image_scale_corr = (std::isfinite(r.rlp) && denom > 0.0) + ? static_cast(r.rlp / denom) + : NAN; + } + outcome.image_scale_g = static_cast(G); + + const double tilt_deg = std::sqrt(psi[0] * psi[0] + psi[1] * psi[1]) * kRadToDeg; + return tilt_deg; +} + +double StillsPartialityRefine::Run(std::vector &outcomes, size_t nthreads) const { + if (nthreads == 0) + nthreads = std::thread::hardware_concurrency(); + nthreads = std::max(1, nthreads); + + double last_mean_tilt = 0.0; + + for (int outer = 0; outer < settings_.outer_iterations; ++outer) { + // Reference full intensities from the current corrections. + const std::vector merged = MergeAll(experiment_, outcomes, false); + std::map reference; + for (const auto &m: merged) + reference[hkl_key_generator_(m)] = m.I; + + std::atomic tilt_sum{0.0}; + std::atomic tilt_n{0}; + std::atomic next{0}; + + auto worker = [&]() { + size_t i = next.fetch_add(1); + while (i < outcomes.size()) { + const double t = RefineOne(outcomes[i], reference); + if (t > 0.0) { + double prev = tilt_sum.load(); + while (!tilt_sum.compare_exchange_weak(prev, prev + t)) {} + tilt_n.fetch_add(1); + } + i = next.fetch_add(1); + } + }; + + const size_t nt = std::min(nthreads, std::max(1, outcomes.size())); + if (nt <= 1) { + worker(); + } else { + std::vector> futures; + futures.reserve(nt); + for (size_t t = 0; t < nt; ++t) + futures.emplace_back(std::async(std::launch::async, worker)); + for (auto &f: futures) + f.get(); + } + + last_mean_tilt = tilt_n > 0 ? tilt_sum.load() / static_cast(tilt_n.load()) : 0.0; + } + + return last_mean_tilt; +} diff --git a/image_analysis/scale_merge/StillsPartialityRefine.h b/image_analysis/scale_merge/StillsPartialityRefine.h new file mode 100644 index 00000000..108f7cac --- /dev/null +++ b/image_analysis/scale_merge/StillsPartialityRefine.h @@ -0,0 +1,59 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "../../common/DiffractionExperiment.h" +#include "../IntegrationOutcome.h" +#include "HKLKey.h" + +// Experimental physical partiality post-refinement for STILLS (env JFJOCH_STILL_POSTREFINE). +// +// The default stills partiality is a frozen scalar-sigma Gaussian (or p == 1): p is set once at +// prediction and never optimised, and any attempt to free a per-image sigma jointly with the per-image +// scale G collapses (within one still the excitation-error spread is narrow, so a scalar sigma is +// degenerate with G). This class instead parametrises partiality by the crystal ORIENTATION - a small +// tilt (dpsi_x, dpsi_y) about the two axes perpendicular to the beam - shared by all of a crystal's +// reflections. A tilt moves each reflection's excitation error by an amount that depends on where the +// reflection sits on the pattern (one side of the Ewald sphere approaches, the opposite recedes), so it +// reshapes the SPATIAL pattern of partialities in a way a single G cannot mimic. That breaks the +// degeneracy that killed the scalar-sigma fit. +// +// Because the integrated intensity I_obs is fixed, refining the tilt only recomputes the partiality +// analytically from the stored per-image lattice/geometry (q = A*.h + B*.k + C*.l, then +// dist_ewald = | |q + S0| - 1/lambda |) - NO pixel re-integration. The loop is: merge -> per-crystal +// refine tilt (G profiled out by the same robust IRLS ScaleOnTheFly uses) -> recompute p and +// image_scale_corr -> re-merge, iterated a few times. Mutates each reflection's `partiality` and +// `image_scale_corr` in place; the existing MergeOnTheFly then consumes the improved corrections. +class StillsPartialityRefine { +public: + struct Settings { + int outer_iterations = 2; // merge <-> refine cycles + int inner_iterations = 3; // (solve G) <-> (refine tilt) alternations per crystal + size_t min_reflections = 40; // skip tilt refinement below this (anti-overfit on sparse crystals) + double max_tilt_rad = 0.0175; // hard bound on |dpsi| (~1 deg); indexing already refined orientation + double prior_sigma_deg = 0.02; // soft prior pulling dpsi toward 0 (0 = off); tames weak-data overfit, + // inert on strong data (well-supported tilts overcome it) + double robust_k = 3.0; // Cauchy loss scale (sigma units) + double gamma_e = 0.0; // angular width sigma(d*) = gamma_e*d* (0 = estimate per crystal from data) + }; + + explicit StillsPartialityRefine(const DiffractionExperiment &x); + StillsPartialityRefine(const DiffractionExperiment &x, Settings settings); + + // Refine all crystals in place. Returns the mean |dpsi| applied (degrees), for diagnostics. + double Run(std::vector &outcomes, size_t nthreads = 0) const; + +private: + const DiffractionExperiment experiment_; + const Settings settings_; + const HKLKeyGenerator hkl_key_generator_; + const std::optional d_min_limit_; + const float bandwidth_sigma_; + + // Refine one crystal against the reference map; returns |dpsi| in degrees (0 if skipped). + double RefineOne(IntegrationOutcome &outcome, + const std::map &reference) const; +}; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 8502d569..172a6d76 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -40,6 +40,7 @@ #include "../image_analysis/lattice_search/LatticeSearch.h" #include "../image_analysis/scale_merge/TwinningAnalysis.h" #include "../image_analysis/scale_merge/HKLKey.h" +#include "../image_analysis/scale_merge/StillsPartialityRefine.h" #include "../image_analysis/WriteReflections.h" #include "../image_analysis/bragg_integration/CalcISigma.h" #include "../common/Definitions.h" @@ -970,6 +971,15 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), false); indexer->ScaleAllImages(merge_result); } + // Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal + // orientation tilt against the running merge and recompute each reflection's partiality + scale + // correction (no re-integration). Never for the P1 search pass. + if (!for_search && experiment_.GetScalingSettings().GetStillsPartialityRefine()) { + phase("Partiality post-refine (" + label + ")"); + StillsPartialityRefine refiner(experiment_); + const double mean_tilt = refiner.Run(indexer->GetIntegrationOutcome(), config_.nthreads); + logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt); + } const std::vector &merge_input = indexer->GetIntegrationOutcome(); phase("Merging"); diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 3a0071a5..4d8450a6 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -106,8 +106,6 @@ std::string RugnuxCommandLine(const ProcessConfig &config, std::ostringstream radii; radii << bragg.GetR1() << "," << bragg.GetR2() << "," << bragg.GetR3(); add("--integration-radius", radii.str()); - if (bragg.GetStillPartiality()) - args.emplace_back("--still-partiality"); // Background trim defaults to 0.10 in the CLI, so emit it whenever the GUI value differs (a // custom fraction, or 0 when the box is unchecked) to reproduce the GUI's choice faithfully. if (bragg.GetBackgroundTrimFraction() != 0.10f) @@ -153,6 +151,8 @@ std::string RugnuxCommandLine(const ProcessConfig &config, add("--partiality-uncertainty", num(sc.GetPartialityUncertaintyCoeff())); if (sc.GetStillsModulation()) args.emplace_back("--stills-modulation"); + if (!sc.GetStillsPartialityRefine()) + args.emplace_back("--simple-stills"); if (!sc.GetExpectedVarianceMerge()) args.emplace_back("--no-expected-variance-merge"); // When merging, the CLI skips the large _process.h5 unless asked; emit the flag when it is diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index e9fcfbc4..eafc7967 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -27,6 +27,7 @@ #include "../image_analysis/scale_merge/Merge.h" #include "../image_analysis/scale_merge/RfreeFlags.h" #include "../image_analysis/scale_merge/ScaleOnTheFly.h" +#include "../image_analysis/scale_merge/StillsPartialityRefine.h" #include "../image_analysis/scale_merge/RotationScaleMerge.h" #include "../image_analysis/scale_merge/ResolutionCutoff.h" #include "../image_analysis/scale_merge/TwinningAnalysis.h" @@ -103,7 +104,7 @@ void print_usage() { std::cout << " --resolution-shells Number of resolution shells in the reported statistics table (default: 10)" << std::endl; std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << std::endl; std::cout << " --capture-uncertainty rot3d: systematic sigma ~num*(1-captured_fraction)*I on under-captured fulls (default: 1.0 for rot3d, 0 otherwise)" << std::endl; - std::cout << " --partiality-uncertainty stills: extra merge sigma ~num*(1-partiality)* on partials (use with --still-partiality; auto-gated to error-model b>1 / ISa<1; default 0, ~2.5 recommended)" << std::endl; + std::cout << " --partiality-uncertainty stills: extra merge sigma ~num*(1-partiality)* on partials (auto-gated to error-model b>1 / ISa<1; default 0, ~2.5 recommended)" << std::endl; std::cout << " --min-captured-fraction rot3d: drop a combined full whose rocking curve was captured below this fraction (edge-of-sweep truncated fulls) (default: 0.7 for rotation, 0 otherwise; 0 = off)" << std::endl; std::cout << " --mosaicity Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl; std::cout << " --reject-outliers Per-observation merge outlier rejection, N sigma from the per-reflection median (default: 6 for rot3d, XDS/DIALS-style; 0 = off)" << std::endl; @@ -120,7 +121,7 @@ void print_usage() { std::cout << " --integration-radius Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl; std::cout << " --background-trim Monochromatic (rotation + still): symmetric trimmed-mean fraction for the background ring (0<=f<0.5, default 0.10; 0 = plain mean). Removes the high-side bias that over-subtracts weak high-angle spots (broadband data keep the sigma-clip instead)" << std::endl; std::cout << " --integrator Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl; - std::cout << " --still-partiality Experimental: weight stills reflections by a Gaussian excitation-error partiality exp(-dist_ewald^2/2sigma^2) instead of treating each as a full" << std::endl; + std::cout << " --simple-stills stills: treat every reflection as a full (p=1, single-pass scale/merge); disables the default physical partiality post-refinement" << std::endl; std::cout << " -q, --azim-q-spacing Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl; std::cout << " --azim-min-q Azimuthal-integration minimum Q (1/A)" << std::endl; std::cout << " --azim-max-q Azimuthal-integration maximum Q (1/A)" << std::endl; @@ -167,7 +168,7 @@ enum { OPT_MODEL, OPT_DUMP_OBSERVATIONS, OPT_INTEGRATOR, - OPT_STILL_PARTIALITY, + OPT_SIMPLE_STILLS, OPT_SCALE_FULLS, OPT_CAPTURE_UNCERTAINTY, OPT_PARTIALITY_UNCERTAINTY, @@ -271,7 +272,7 @@ static option long_options[] = { {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS}, {"background-trim", required_argument, nullptr, OPT_BACKGROUND_TRIM}, {"integrator", required_argument, nullptr, OPT_INTEGRATOR}, - {"still-partiality", no_argument, nullptr, OPT_STILL_PARTIALITY}, + {"simple-stills", no_argument, nullptr, OPT_SIMPLE_STILLS}, {"detect-ice-rings", optional_argument, nullptr, OPT_DETECT_ICE_RINGS}, {"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS}, {"reject-delta-cchalf", required_argument, nullptr, OPT_REJECT_DELTA_CCHALF}, @@ -547,7 +548,7 @@ int main(int argc, char **argv) { std::optional integration_radius_arg; // "r1" or "r1,r2,r3" std::optional background_trim_arg; // --background-trim: background-ring trimmed-mean fraction std::optional integrator_mode; // --integrator boxsum|gaussian|empirical - bool still_partiality_flag = false; // --still-partiality (experimental stills partiality) + bool simple_stills_flag = false; // --simple-stills: disable the default stills partiality post-refinement std::optional outlier_reject_nsigma; // merge per-observation outlier rejection std::optional delta_cchalf_nsigma; // per-crystal CC1/2-delta rejection @@ -829,8 +830,8 @@ int main(int argc, char **argv) { else if (strcmp(optarg, "empirical") == 0) integrator_mode = IntegratorMode::ProfileEmpirical; else { logger.Error("--integrator expects boxsum|gaussian|empirical"); return 1; } break; - case OPT_STILL_PARTIALITY: - still_partiality_flag = true; + case OPT_SIMPLE_STILLS: + simple_stills_flag = true; break; case OPT_REJECT_OUTLIERS: outlier_reject_nsigma = parse_double_arg(optarg, "--reject-outliers", logger); @@ -1081,6 +1082,7 @@ int main(int argc, char **argv) { (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.StillsPartialityRefine(!simple_stills_flag); scaling_settings.ExpectedVarianceMerge(!no_expected_variance_merge); scaling_settings.OutlierRejectNsigma( outlier_reject_nsigma.value_or( @@ -1137,6 +1139,14 @@ int main(int argc, char **argv) { // the data's own merge. The per-image scale G is the exact one-pass solution, so one pass // (iterating a self-rebuilt reference only re-fits the freshly-scaled noise on weak stills). ScaleOnTheFly(experiment, MergeAll(experiment, reflections)).Scale(reflections, nthreads); + // Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal + // orientation tilt against the merge and recompute each reflection's partiality + scale correction + // (no re-integration), then merge with the improved corrections. + if (experiment.GetScalingSettings().GetStillsPartialityRefine()) { + StillsPartialityRefine refiner(experiment); + const double mean_tilt = refiner.Run(reflections, nthreads); + logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt); + } MergeOnTheFly merge_engine(experiment); merge_engine.ReferenceCell(experiment.GetUnitCell()); // Optional detector-plane modulation (flat-field) correction, folded into each reflection's @@ -1394,6 +1404,7 @@ int main(int argc, char **argv) { if (no_scaling_corrections) scaling_settings.CorrectionSurfaces(false); scaling_settings.StillsModulation(stills_modulation_flag); + scaling_settings.StillsPartialityRefine(!simple_stills_flag); scaling_settings.ExpectedVarianceMerge(!no_expected_variance_merge); if (d_min_scale_merge) scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value()); @@ -1469,13 +1480,6 @@ int main(int argc, char **argv) { *background_trim_arg); } - if (still_partiality_flag) { - BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); - bis.StillPartiality(true); - experiment.ImportBraggIntegrationSettings(bis); - logger.Info("Stills partiality enabled (experimental Gaussian excitation-error weighting)"); - } - SpotFindingSettings spot_settings; spot_settings.enable = true; spot_settings.indexing = true; diff --git a/viewer/widgets/JFJochViewerSettingsDock.cpp b/viewer/widgets/JFJochViewerSettingsDock.cpp index ee2087aa..fea318b9 100644 --- a/viewer/widgets/JFJochViewerSettingsDock.cpp +++ b/viewer/widgets/JFJochViewerSettingsDock.cpp @@ -446,11 +446,6 @@ QWidget *JFJochViewerSettingsDock::BuildBraggSection() { auto *r3 = new NumberLineEdit(1.0f, 30.0f, bragg_.GetR3(), 1, "px", this); auto *radii = new QHBoxLayout(); radii->addWidget(r1); radii->addWidget(r2); radii->addWidget(r3); - auto *stillPartiality = new QCheckBox("Stills partiality (experimental)", this); - stillPartiality->setChecked(bragg_.GetStillPartiality()); - stillPartiality->setToolTip("Experimental, stills only: weight each reflection by a Gaussian " - "excitation-error partiality exp(-d_ewald²/2σ²) instead of treating it as a " - "full. Pairs with \"Partiality uncertainty\" in Scaling."); // Background trim: replace the r2..r3 ring mean with a symmetric trimmed mean (drop the lowest and // highest fraction of ring pixels), which removes the high-side bias that over-subtracts weak // high-angle reflections. Checkbox + fraction; 0 = plain mean. @@ -468,7 +463,6 @@ QWidget *JFJochViewerSettingsDock::BuildBraggSection() { trimRow->addWidget(bkgTrimFrac, 1); form->addRow("", gaussian); form->addRow("Radii r1/r2/r3", radii); - form->addRow("", stillPartiality); form->addRow("", trimRow); section->setContentLayout(form); section->setExpanded(false); // folded on start (only geometry + unit cell start open) @@ -477,12 +471,10 @@ QWidget *JFJochViewerSettingsDock::BuildBraggSection() { bragg_.Integrator(gaussian->isChecked() ? IntegratorMode::ProfileGaussian : IntegratorMode::BoxSum); bragg_.R1(static_cast(r1->value())).R2(static_cast(r2->value())) .R3(static_cast(r3->value())); - bragg_.StillPartiality(stillPartiality->isChecked()); bragg_.BackgroundTrimFraction(bkgTrim->isChecked() ? static_cast(bkgTrimFrac->value()) : 0.0f); emit braggChanged(bragg_); }; connect(gaussian, &QCheckBox::toggled, this, [emitBragg] { emitBragg(); }); - connect(stillPartiality, &QCheckBox::toggled, this, [emitBragg] { emitBragg(); }); connect(bkgTrim, &QCheckBox::toggled, this, [bkgTrim, bkgTrimFrac, emitBragg] { bkgTrimFrac->setEnabled(bkgTrim->isChecked()); emitBragg(); }); connect(bkgTrimFrac, &NumberLineEdit::newValue, this, [emitBragg] { emitBragg(); }); @@ -512,20 +504,27 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() { modulation->setToolTip("Stills: fit a detector-plane modulation (flat-field) surface over where each " "reflection lands, cross-validated so it no-ops when the systematic is absent. " "For rotation, modulation is part of \"Correction surfaces\" above."); + auto *partRefine = new QCheckBox("Partiality post-refinement (stills)", this); + partRefine->setChecked(scaling_.GetStillsPartialityRefine()); + partRefine->setToolTip("Stills: refine a per-crystal orientation tilt against the running merge and " + "recompute each reflection's partiality (physical Ewald-proximity model), then " + "re-scale/merge. On by default; uncheck for the simple model (each reflection a " + "full, single-pass). No effect on rotation data."); auto *limitRes = new QCheckBox("High-resolution limit", this); limitRes->setChecked(scaling_.GetHighResolutionLimit_A().has_value()); auto *highRes = new NumberLineEdit(0.3f, 5.0f, scaling_.GetHighResolutionLimit_A().value_or(2.0), 1, "Å", this); highRes->setEnabled(limitRes->isChecked()); // Stills partiality-uncertainty merge term: adds a systematic sigma ~c*(1-partiality)* on partials, - // so strong low-partiality partials are not over-trusted. Pairs with "Stills partiality" (Bragg); the - // library auto-gates it to strong/medium data. 0 = off; ~2.5 recommended. + // so strong low-partiality partials are not over-trusted. Relevant once partials exist (partiality + // post-refinement on, the default); the library auto-gates it to strong/medium data. 0 = off; ~2.5 rec. const double part_unc = scaling_.GetPartialityUncertaintyCoeff(); auto *partUncertain = new QCheckBox("Partiality uncertainty", this); partUncertain->setChecked(part_unc > 0.0); partUncertain->setToolTip("Stills: add a systematic merge σ ~c·(1−partiality)·⟨I⟩ to partials so strong " - "low-partiality partials are not over-trusted. Use with \"Stills partiality\"; " - "auto-gated to strong/medium data. ~2.5 recommended; unchecked = off."); + "low-partiality partials are not over-trusted. Relevant with partiality " + "post-refinement (the default); auto-gated to strong/medium data. ~2.5 " + "recommended; unchecked = off."); auto *partUncertainCoeff = new NumberLineEdit(0.1f, 10.0f, part_unc > 0.0 ? part_unc : 2.5, 1, "", this); partUncertainCoeff->setEnabled(partUncertain->isChecked()); @@ -533,6 +532,7 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() { form->addRow("", refineB); form->addRow("", corrections); form->addRow("", modulation); + form->addRow("", partRefine); // Compact, and aligned with the checkboxes above: the limit checkbox + value sit together in the // field column (not as a row label, which would indent it differently). auto *resRow = new QHBoxLayout(); @@ -551,6 +551,7 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() { scaling_.RefineB(refineB->isChecked()); scaling_.CorrectionSurfaces(corrections->isChecked()); scaling_.StillsModulation(modulation->isChecked()); + scaling_.StillsPartialityRefine(partRefine->isChecked()); scaling_.HighResolutionLimit_A(limitRes->isChecked() ? std::optional(highRes->value()) : std::nullopt); scaling_.PartialityUncertaintyCoeff(partUncertain->isChecked() ? partUncertainCoeff->value() : 0.0); @@ -560,6 +561,7 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() { connect(refineB, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); }); connect(corrections, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); }); connect(modulation, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); }); + connect(partRefine, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); }); connect(limitRes, &QCheckBox::toggled, this, [emitScaling, highRes](bool on) { highRes->setEnabled(on); emitScaling(); }); connect(highRes, &NumberLineEdit::newValue, this, [emitScaling] { emitScaling(); });