diff --git a/common/Definitions.h b/common/Definitions.h index f583919f..76fb17c2 100644 --- a/common/Definitions.h +++ b/common/Definitions.h @@ -70,3 +70,16 @@ inline bool IsOnIceRing(float d_A, float half_width_q_recipA) { return true; return false; } + +// Index into ICE_RING_RES_A of the hexagonal-ice ring resolution d sits on (within half_width in +// q = 2*pi/d), or -1 if none. Used to look up that ring's per-image strength for the merge gate. +inline int IceRingIndex(float d_A, float half_width_q_recipA) { + if (!(d_A > 0.0f)) + return -1; + constexpr float two_pi = 6.283185307f; + const float q = two_pi / d_A; + for (size_t i = 0; i < ICE_RING_RES_A.size(); ++i) + if (std::fabs(q - two_pi / ICE_RING_RES_A[i]) < half_width_q_recipA) + return static_cast(i); + return -1; +} diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index ff663c18..3fb697fa 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -15,6 +15,7 @@ #include "../../common/CorrelationCoefficient.h" #include "../../common/ResolutionShells.h" +#include "../../common/Definitions.h" #include "HKLKey.h" namespace { @@ -54,6 +55,13 @@ MergeOnTheFly &MergeOnTheFly::ReferenceCell(const std::optional &cell) return *this; } +bool MergeOnTheFly::IsMaskedRing(const Reflection &r) const { + if (masked_ice_rings.empty()) + return false; + const int ring = IceRingIndex(r.d, mask_ice_half_width_q); + return ring >= 0 && ring < static_cast(masked_ice_rings.size()) && masked_ice_rings[ring]; +} + void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id, bool cc_mask) { std::unique_lock ul(merged_mutex); @@ -72,6 +80,8 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id continue; if (exclude_ice_rings && r.on_ice_ring) continue; + if (IsMaskedRing(r)) + continue; if (r.partiality < min_partiality) continue; @@ -168,6 +178,8 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc continue; if (exclude_ice_rings && r.on_ice_ring) continue; + if (IsMaskedRing(r)) + continue; if (r.partiality < min_partiality) continue; const float I_corr = r.I * r.image_scale_corr; diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index ce2238ec..84d8af2e 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -73,6 +73,11 @@ class MergeOnTheFly { // see the ice-contaminated intensities. The final in-symmetry merge keeps them (for completeness). bool exclude_ice_rings = false; + // Ice rings (indices into ICE_RING_RES_A) to exclude from this merge because their merged CC1/2 + // collapsed relative to their resolution neighbours - ice decorrelated them. Empty = none masked. + std::vector masked_ice_rings; + float mask_ice_half_width_q = 0.03f; + HKLKeyGenerator generator; std::map accumulator; @@ -103,10 +108,14 @@ class MergeOnTheFly { size_t reject_count = 0; bool Mask(const IntegrationOutcome &outcome, bool cc_mask); + [[nodiscard]] bool IsMaskedRing(const Reflection &r) const; public: MergeOnTheFly(const DiffractionExperiment &x); MergeOnTheFly& ReferenceCell(const std::optional &cell); MergeOnTheFly& ExcludeIceRings(bool input) { exclude_ice_rings = input; return *this; } + MergeOnTheFly& MaskIceRings(std::vector masked, float half_width_q) { + masked_ice_rings = std::move(masked); mask_ice_half_width_q = half_width_q; return *this; + } // Fit the global error model from the spread of symmetry-equivalent observations. // Call once before merging; AddImage then applies it. diff --git a/process/JFJochProcess.cpp b/process/JFJochProcess.cpp index af1253ac..f2981917 100644 --- a/process/JFJochProcess.cpp +++ b/process/JFJochProcess.cpp @@ -35,6 +35,7 @@ #include "../image_analysis/scale_merge/HKLKey.h" #include "../image_analysis/WriteReflections.h" #include "../common/Definitions.h" +#include "../common/CorrelationCoefficient.h" #include #include #include @@ -566,6 +567,9 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { // Smoothing it more than once would compound the correction, so do it only on the first // pass. The no-reference path recomputes G from scratch each pass and re-smooths correctly. bool reference_g_smoothed = false; + // Ice rings masked from the merge by the CC1/2 test after the final merge; empty until then, + // at which point a re-merge applies them. + std::vector masked_ice_rings; auto scale_and_merge = [&](const std::string &label, bool for_search) -> ScaleMergeResult { // ScaleOnTheFly self-scaling is only for the no-reference path; with a reference each // image is already scaled against it during the per-image pass, so we merge directly. @@ -620,6 +624,7 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { // intensities and the error model, so the space-group search sees clean data; the final // in-symmetry merge keeps them so completeness is not lost. merge_engine.ExcludeIceRings(for_search); + 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); merge_engine.RefineErrorModel(merge_input); @@ -670,6 +675,39 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { } } + // Ice-ring CC1/2 mask: a hexagonal-ice ring whose merged half-set CC1/2 collapses well below its + // resolution shoulders has been decorrelated by ice (its Bragg is unrecoverable) - drop it and + // re-merge. Weak/absent rings track their neighbours and stay, so completeness on clean crystals + // is untouched. Uses the ring band (+/- ice width in q=2pi/d) vs the shoulders either side. + if (experiment_.IsDetectIceRings() && !sm.merged.empty()) { + constexpr float two_pi = 6.283185307f; + const float w = config_.spot_finding.ice_ring_width_Q_recipA; + std::vector mask(ICE_RING_RES_A.size(), 0); + for (size_t i = 0; i < ICE_RING_RES_A.size(); ++i) { + const float q_ring = two_pi / ICE_RING_RES_A[i]; + CorrelationCoefficient ring, shoulder; + size_t n_ring = 0, n_shoulder = 0; + for (const auto &m : sm.merged) { + if (!(m.d > 0.0f) || !std::isfinite(m.I_half[0]) || !std::isfinite(m.I_half[1])) + continue; + const float dq = std::fabs(two_pi / m.d - q_ring); + if (dq < w) { ring.Add(m.I_half[0], m.I_half[1]); ++n_ring; } + else if (dq < 3.0f * w) { shoulder.Add(m.I_half[0], m.I_half[1]); ++n_shoulder; } + } + if (n_ring >= 20 && n_shoulder >= 20 && shoulder.GetCC() > 0.5 + && ring.GetCC() < shoulder.GetCC() - 0.05) { + mask[i] = 1; + logger.Info("Ice-ring mask: {:.2f} A ring CC1/2 {:.3f} << shoulders {:.3f}; masked from merge", + ICE_RING_RES_A[i], ring.GetCC(), shoulder.GetCC()); + } + } + if (std::any_of(mask.begin(), mask.end(), [](char c) { return c != 0; })) { + masked_ice_rings = std::move(mask); + const auto final_sg = experiment_.GetGemmiSpaceGroup(); + sm = scale_and_merge(final_sg ? final_sg->short_name() : "P1", false); + } + } + const auto twin_sg_number = experiment_.GetSpaceGroupNumber(); const gemmi::SpaceGroup *twin_sg = twin_sg_number ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr;