From fb399480e5c93808aa93cdb3ff9ee6fab888d735 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 12:03:21 +0200 Subject: [PATCH] rugnux: add opt-in stills partiality (--still-partiality) Serial stills currently treat every reflection as a full (partiality hardcoded to 1). Add an opt-in Gaussian excitation-error partiality set at prediction time (CPU + CUDA): p = exp(-dist_ewald^2 / (2*sigma^2)), sigma^2 = profile_radius^2 + (bandwidth_sigma*|recip_z|)^2, with sigma = the per-image profile radius (ewald_dist_cutoff/2), so an edge-of-acceptance reflection keeps p ~ exp(-2). Off by default; the merge weight (~p^2) then down-weights far-from-Ewald partials instead of trusting them as fulls. Validated: helps medium/strong stills (LOV R-free 0.336->0.329, lyso8 0.433->0.410, lowers the systematic error-model b in both) but HURTS weak OCP (dividing by a small, uncertain p amplifies orientation error -> high-res noise, resolution collapse), so it is left opt-in. A static forward p explains only ~6-10% of the partiality scatter; the full win needs per-image post-refinement (future work), for which this is the prediction-side groundwork. Co-Authored-By: Claude Opus 4.8 (1M context) --- common/BraggIntegrationSettings.cpp | 9 +++++++++ common/BraggIntegrationSettings.h | 3 +++ image_analysis/IndexAndRefine.cpp | 6 +++++- .../bragg_prediction/BraggPrediction.cpp | 9 ++++++++- .../bragg_prediction/BraggPrediction.h | 10 ++++++++-- .../bragg_prediction/BraggPredictionGPU.cu | 17 ++++++++++++++--- .../bragg_prediction/BraggPredictionGPU.h | 2 ++ rugnux/rugnux_cli.cpp | 14 ++++++++++++++ 8 files changed, 63 insertions(+), 7 deletions(-) diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index d69a0028..3436852f 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -96,3 +96,12 @@ float BraggIntegrationSettings::GetDMinLimit_A() const { float BraggIntegrationSettings::GetMinimumSigmaInRegardsToI() const { return minimum_sigma_in_regards_to_i; } + +BraggIntegrationSettings &BraggIntegrationSettings::StillPartiality(bool input) { + still_partiality = input; + return *this; +} + +bool BraggIntegrationSettings::GetStillPartiality() const { + return still_partiality; +} diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index bd53b7db..e0598756 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -20,6 +20,7 @@ class BraggIntegrationSettings { float d_min_limit_A = 1.0; std::optional fixed_profile_radius; float minimum_sigma_in_regards_to_i = 0.02; + bool still_partiality = false; // experimental stills excitation-error partiality (rugnux --still-partiality) public: BraggIntegrationSettings& R1(float input); @@ -28,6 +29,7 @@ public: BraggIntegrationSettings& DMinLimit_A(float input); BraggIntegrationSettings& FixedProfileRadius_recipA(std::optional input); BraggIntegrationSettings& Integrator(IntegratorMode input); + BraggIntegrationSettings& StillPartiality(bool input); [[nodiscard]] IntegratorMode GetIntegrator() const; @@ -38,4 +40,5 @@ public: [[nodiscard]] float GetDMinLimit_A() const; [[nodiscard]] float GetMinimumSigmaInRegardsToI() const; + [[nodiscard]] bool GetStillPartiality() const; }; diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index b0de29fa..8d542eac 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -324,7 +324,11 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, .wedge_deg = std::fabs(wedge_deg), .mosaicity_deg = std::fabs(mos_deg), // FWHM -> sigma; 0 when monochromatic, leaving the prediction unchanged. - .bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f + .bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f, + // Experimental stills partiality (off by default): sigma = ewald_dist_cutoff/2 = the per-image + // profile radius, so a reflection at the acceptance edge (dist_ewald ~ 2*sigma) keeps p ~ exp(-2). + .still_partiality = experiment.GetBraggIntegrationSettings().GetStillPartiality(), + .profile_radius_recipA = ewald_dist_cutoff * 0.5f }; // Predict, then integrate with the selected integrator (box-sum or profile-fit). diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp index 701da2b9..e19d0249 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.cpp +++ b/image_analysis/bragg_prediction/BraggPrediction.cpp @@ -156,6 +156,13 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal continue; float d = 1.0f / sqrtf(recip_sq); + float partiality = 1.0f; + if (settings.still_partiality && settings.profile_radius_recipA > 0.0f) { + const float sig_bw = settings.bandwidth_sigma * std::fabs(recip_z); + const float sigma2 = settings.profile_radius_recipA * settings.profile_radius_recipA + + sig_bw * sig_bw; + partiality = std::exp(-0.5f * dist_ewald_sphere * dist_ewald_sphere / sigma2); + } reflections[i] = Reflection{ .h = h, .k = k, @@ -168,7 +175,7 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal .d = d, .dist_ewald = dist_ewald_sphere, .rlp = 1.0, - .partiality = 1.0, + .partiality = partiality, .zeta = 1.0, .image_scale_corr = 1.0 }; diff --git a/image_analysis/bragg_prediction/BraggPrediction.h b/image_analysis/bragg_prediction/BraggPrediction.h index a2a5b52a..5406d23e 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.h +++ b/image_analysis/bragg_prediction/BraggPrediction.h @@ -21,9 +21,15 @@ struct BraggPredictionSettings { // Relative X-ray bandwidth Δλ/λ expressed as a Gaussian sigma (0 = monochromatic). // When > 0 the Ewald-shell acceptance is thickened radially per reflection by // σ_bw = |recip_z|·bandwidth_sigma (= bλ/2d²), so the 1/d² pink-beam smear no - // longer clips high-resolution reflections. Must stay the last member so existing - // designated initializers remain valid. + // longer clips high-resolution reflections. float bandwidth_sigma = 0.0f; + // Experimental stills partiality (rugnux --still-partiality). When still_partiality is set and + // profile_radius_recipA > 0, each reflection gets a Gaussian excitation-error partiality + // p = exp(-dist_ewald^2 / (2*sigma^2)), sigma^2 = profile_radius_recipA^2 + (bandwidth_sigma*|recip_z|)^2, + // instead of the fixed 1.0 (full). Off by default. Keep these two as the trailing members so the + // existing designated initializers (which stop at bandwidth_sigma) remain valid. + bool still_partiality = false; + float profile_radius_recipA = 0.0f; }; class BraggPrediction { diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.cu b/image_analysis/bragg_prediction/BraggPredictionGPU.cu index 615eb794..8e172867 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.cu @@ -127,7 +127,13 @@ namespace { out.d = 1.0f / sqrtf(recip_sq); out.dist_ewald = dist_ewald; out.rlp = 1.0f; - out.partiality = 1.0f; + float partiality = 1.0f; + if (C.still_partiality && C.profile_radius_recipA > 0.0f) { + const float sig_bw = C.bandwidth_sigma * fabsf(recip_z); + const float sigma2 = C.profile_radius_recipA * C.profile_radius_recipA + sig_bw * sig_bw; + partiality = expf(-0.5f * dist_ewald * dist_ewald / sigma2); + } + out.partiality = partiality; out.zeta = 1.0f; out.image_scale_corr = 1.0f; return true; @@ -158,7 +164,9 @@ namespace { float high_res_A, float ewald_dist_cutoff, char centering, - float bandwidth_sigma) { + float bandwidth_sigma, + bool still_partiality, + float profile_radius_recipA) { KernelConsts kc{}; auto geom = experiment.GetDiffractionGeometry(); kc.det_width_pxl = static_cast(experiment.GetXPixelsNum()); @@ -171,6 +179,8 @@ namespace { kc.one_over_wavelength = 1.0f / geom.GetWavelength_A(); kc.ewald_cutoff = ewald_dist_cutoff; kc.bandwidth_sigma = bandwidth_sigma; + kc.still_partiality = still_partiality; + kc.profile_radius_recipA = profile_radius_recipA; kc.Astar = lattice.Astar(); kc.Bstar = lattice.Bstar(); kc.Cstar = lattice.Cstar(); @@ -193,7 +203,8 @@ int BraggPredictionGPU::Calc(const DiffractionExperiment &experiment, const BraggPredictionSettings &settings) { // Build constants on host KernelConsts hK = BuildKernelConsts(experiment, lattice, settings.high_res_A, settings.ewald_dist_cutoff, - settings.centering, settings.bandwidth_sigma); + settings.centering, settings.bandwidth_sigma, + settings.still_partiality, settings.profile_radius_recipA); cudaMemcpyAsync(dK, &hK, sizeof(KernelConsts), cudaMemcpyHostToDevice, stream); cudaMemsetAsync(d_count, 0, sizeof(int), stream); diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.h b/image_analysis/bragg_prediction/BraggPredictionGPU.h index ba0fc1f5..a7d2032d 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.h +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.h @@ -19,6 +19,8 @@ struct KernelConsts { float one_over_dmax_sq; float ewald_cutoff; float bandwidth_sigma; // relative Δλ/λ (sigma); 0 = monochromatic + bool still_partiality; // experimental stills excitation-error partiality (off => p = 1, fulls) + float profile_radius_recipA; // Gaussian sigma [1/A] for the stills partiality Coord Astar, Bstar, Cstar, S0; float rot[9]; char centering; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 40621736..42446fec 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -109,6 +109,7 @@ void print_usage() { std::cout << " --bandwidth Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM); default from file or 0" << std::endl; 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 << " --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 << " -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; @@ -153,6 +154,7 @@ enum { OPT_MODEL, OPT_DUMP_OBSERVATIONS, OPT_INTEGRATOR, + OPT_STILL_PARTIALITY, OPT_SCALE_FULLS, OPT_CAPTURE_UNCERTAINTY, OPT_MIN_CAPTURED_FRACTION, @@ -246,6 +248,7 @@ static option long_options[] = { {"bandwidth", required_argument, nullptr, OPT_BANDWIDTH}, {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS}, {"integrator", required_argument, nullptr, OPT_INTEGRATOR}, + {"still-partiality", no_argument, nullptr, OPT_STILL_PARTIALITY}, {"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}, @@ -489,6 +492,7 @@ int main(int argc, char **argv) { std::optional report_shell_count; // --resolution-shells std::optional integration_radius_arg; // "r1" or "r1,r2,r3" std::optional integrator_mode; // --integrator boxsum|gaussian|empirical + bool still_partiality_flag = false; // --still-partiality (experimental stills partiality) std::optional outlier_reject_nsigma; // merge per-observation outlier rejection std::optional delta_cchalf_nsigma; // per-crystal CC1/2-delta rejection @@ -726,6 +730,9 @@ 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; + break; case OPT_REJECT_OUTLIERS: outlier_reject_nsigma = parse_double_arg(optarg, "--reject-outliers", logger); break; @@ -1331,6 +1338,13 @@ int main(int argc, char **argv) { : "profile (empirical)"); } + 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;