From a74b8767de7084250773c270702a12b79a31f5ef Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 30 Jun 2026 20:40:51 +0200 Subject: [PATCH] ice rings: add --detect-ice-rings and exclude ice-ring reflections from scaling jfjoch_process gains --detect-ice-rings, which (a) activates the existing spot-finder ice flagging (ice spots de-prioritised in indexing) and (b) drops reflections sitting on a hexagonal-ice powder ring from scaling/combine/merge/stats, via a new shared IsOnIceRing() helper over ICE_RING_RES_A using the spot-finder's q half-width. Their integrated intensity is contaminated by the strong, variable ice background, so leaving them in mis-scales the whole frame and inflates the error model. On EP0117 (ice-ring crystal): de-novo space-group determination recovers from P1 to P2 and CC1/2 improves (31->37%). Off by default; a no-op without the flag. This is the first, non-controversial step - the residual gap needs ice-aware background/integration (follow-up). Co-Authored-By: Claude Opus 4.8 (1M context) --- common/Definitions.h | 15 +++++++++++++++ process/JFJochProcess.cpp | 20 ++++++++++++++++++++ tools/jfjoch_process.cpp | 10 ++++++++++ 3 files changed, 45 insertions(+) diff --git a/common/Definitions.h b/common/Definitions.h index 8fd9eb83..f583919f 100644 --- a/common/Definitions.h +++ b/common/Definitions.h @@ -9,6 +9,7 @@ #include #include #include +#include constexpr float WVL_1A_IN_KEV = 12.39854f; @@ -55,3 +56,17 @@ constexpr float LAB6_CELL_A = 4.156468f; // https://journals.iucr.org/d/issues/2021/04/00/tz5104/index.html constexpr std::array ICE_RING_RES_A = {3.895, 3.661, 3.438, 2.667, 2.249, 2.068, 1.947, 1.916, 1.882, 1.719, 1.522}; + +// True when resolution d (Angstrom) sits within half_width of a hexagonal-ice powder ring, in the +// q = 2*pi/d units the spot-finder uses (ice_ring_width_Q_recipA). Used to drop ice-contaminated +// reflections from scaling/merging when ice-ring handling is enabled. +inline bool IsOnIceRing(float d_A, float half_width_q_recipA) { + if (!(d_A > 0.0f)) + return false; + constexpr float two_pi = 6.283185307f; + const float q = two_pi / d_A; + for (const float ice_d : ICE_RING_RES_A) + if (std::fabs(q - two_pi / ice_d) < half_width_q_recipA) + return true; + return false; +} diff --git a/process/JFJochProcess.cpp b/process/JFJochProcess.cpp index bbb1fc40..e0976638 100644 --- a/process/JFJochProcess.cpp +++ b/process/JFJochProcess.cpp @@ -34,6 +34,7 @@ #include "../image_analysis/scale_merge/Combine3D.h" #include "../image_analysis/scale_merge/HKLKey.h" #include "../image_analysis/WriteReflections.h" +#include "../common/Definitions.h" #include #include #include @@ -532,6 +533,25 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { }; phase("Scaling and merging"); + // Ice-ring handling (--detect-ice-rings): drop reflections sitting on a hexagonal-ice powder + // ring before scaling. Their integrated intensity is contaminated by the strong, variable ice + // background, so leaving them in mis-scales the whole frame (the per-image fit is dragged) and + // inflates the error model. Removed once here, so every scaling pass, the combine, the merge + // and the statistics all skip them. + if (experiment_.IsDetectIceRings()) { + const float ice_width = config_.spot_finding.ice_ring_width_Q_recipA; + size_t total = 0, dropped = 0; + for (auto &outcome : indexer->GetIntegrationOutcome()) { + total += outcome.reflections.size(); + const size_t before = outcome.reflections.size(); + std::erase_if(outcome.reflections, + [&](const Reflection &r) { return IsOnIceRing(r.d, ice_width); }); + dropped += before - outcome.reflections.size(); + } + logger.Info("Ice-ring exclusion: dropped {} of {} reflections on ice rings (half-width {:.3f} A^-1)", + dropped, total, ice_width); + } + // Scale the images and merge. Factored so it can run twice: first in P1 to give the // space-group search a merged dataset, then again in the determined space group so the // scaling sees symmetry equivalents and the final statistics are in the right symmetry. diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 0d439394..6357de37 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -37,6 +37,7 @@ void print_usage() { std::cout << " --spot-threshold Photon count threshold for spot finding (default: 10)" << std::endl; std::cout << " --spot-high-resolution High resolution limit for spot finding (default: 1.5)" << std::endl; std::cout << " --max-spots Max spot count (default: 250)" << std::endl; + std::cout << " --detect-ice-rings Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling/merging" << std::endl; std::cout << std::endl; std::cout << " Indexing" << std::endl; @@ -106,6 +107,7 @@ enum { OPT_MOSAICITY, OPT_SMOOTH_G, OPT_RECIPROCAL_PROFILE, + OPT_DETECT_ICE_RINGS, OPT_NO_SCALE_FULLS, OPT_WRITE_PROCESS_H5, OPT_ABSORPTION @@ -157,6 +159,7 @@ static option long_options[] = { {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS}, {"integrator", required_argument, nullptr, OPT_INTEGRATOR}, {"reciprocal-profile", no_argument, nullptr, OPT_RECIPROCAL_PROFILE}, + {"detect-ice-rings", no_argument, nullptr, OPT_DETECT_ICE_RINGS}, {"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS}, {"reject-delta-cchalf", required_argument, nullptr, OPT_REJECT_DELTA_CCHALF}, {nullptr, 0, nullptr, 0} @@ -297,6 +300,7 @@ int main(int argc, char **argv) { bool run_scaling = false; std::optional scale_fulls_arg; // --scale-fulls / --no-scale-fulls; default on for rot3d bool write_process_h5_flag = false; // --write-process-h5; also write _process.h5 when merging + bool detect_ice_rings_flag = false; // --detect-ice-rings; flag ice-ring spots and exclude them from scaling int absorption_iter = 0; // --absorption[=n]; rot3d absorption-surface iterations (0 = off) std::optional smooth_g_deg_arg; // --smooth-g[=deg]; default 5 deg for rot3d, 0 (off) otherwise bool anomalous_mode = false; @@ -535,6 +539,9 @@ int main(int argc, char **argv) { case OPT_NO_SCALE_FULLS: scale_fulls_arg = false; break; + case OPT_DETECT_ICE_RINGS: + detect_ice_rings_flag = true; + break; case OPT_WRITE_PROCESS_H5: write_process_h5_flag = true; break; @@ -751,6 +758,9 @@ int main(int argc, char **argv) { indexing_settings.GeomRefinementAlgorithm(refinement_algorithm); experiment.ImportIndexingSettings(indexing_settings); + if (detect_ice_rings_flag) + experiment.DetectIceRings(true); + // Scale-fulls refits the per-frame scale on the rot3d combined fulls; it is on by default for // rot3d (where it lifts ISa substantially) and off otherwise. --no-scale-fulls overrides. const bool scale_fulls = scale_fulls_arg.value_or(combine_3d);