From e70b21f6de781a7530f52ee3000a72730111a23a Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Mon, 11 May 2026 08:29:41 +0200 Subject: [PATCH] HKLKey: Test properly, fix bug with merging Friedel pairs when no SG given --- image_analysis/scale_merge/HKLKey.cpp | 43 ++++++------- image_analysis/scale_merge/HKLKey.h | 13 ++++ image_analysis/scale_merge/Merge.cpp | 9 +-- image_analysis/scale_merge/ScaleOnTheFly.cpp | 7 ++- tests/CMakeLists.txt | 1 + tests/MergeScaleTest.cpp | 66 ++++++++++++++++++++ 6 files changed, 109 insertions(+), 30 deletions(-) create mode 100644 tests/MergeScaleTest.cpp diff --git a/image_analysis/scale_merge/HKLKey.cpp b/image_analysis/scale_merge/HKLKey.cpp index b8cc31cb..ffa28162 100644 --- a/image_analysis/scale_merge/HKLKey.cpp +++ b/image_analysis/scale_merge/HKLKey.cpp @@ -6,23 +6,32 @@ #include "HKLKey.h" #include "gemmi/symmetry.hpp" -HKLKey CanonicalHKL(int32_t h, int32_t k, int32_t l, bool merge_friedel, const std::optional &in_sg) { +HKLKeyGenerator::HKLKeyGenerator(bool merge_friedel, const std::optional &sg) +: merge_friedel(merge_friedel), sg(sg) {} + +HKLKey HKLKeyGenerator::operator()(const MergedReflection &r) const { + return operator()(r.h, r.k, r.l); +} + +HKLKey HKLKeyGenerator::operator()(const Reflection &r) const { + return operator()(r.h, r.k, r.l); +} + +HKLKey HKLKeyGenerator::operator()(int32_t h, int32_t k, int32_t l) const { HKLKey key{h, k, l, true}; - if (!in_sg.has_value()) { - if (!merge_friedel) { - const HKLKey neg{-h, -k, -l, true}; - if (std::tie(key.h, key.k, key.l) < std::tie(neg.h, neg.k, neg.l)) { - key.h = -key.h; - key.k = -key.k; - key.l = -key.l; - key.plus = false; - } + if (!sg.has_value()) { + const HKLKey neg{-h, -k, -l, true}; + if (std::tie(key.h, key.k, key.l) < std::tie(neg.h, neg.k, neg.l)) { + key.h = -key.h; + key.k = -key.k; + key.l = -key.l; + key.plus = merge_friedel; } } else { - gemmi::SpaceGroup sg = in_sg.value(); - const auto ops = sg.operations(); - const gemmi::ReciprocalAsu asu(&sg); + const auto sg_local = sg.value(); + const auto ops = sg_local.operations(); + const gemmi::ReciprocalAsu asu(&sg_local); const gemmi::Op::Miller in{h, k, l}; const auto [hkl, sign_plus] = asu.to_asu_sign(in, ops); @@ -35,14 +44,6 @@ HKLKey CanonicalHKL(int32_t h, int32_t k, int32_t l, bool merge_friedel, const s return key; } -HKLKey CanonicalHKL(const Reflection &r, bool merge_friedel, const std::optional &sg) { - return CanonicalHKL(r.h, r.k, r.l, merge_friedel, sg); -} - -HKLKey CanonicalHKL(const MergedReflection &r, bool merge_friedel, const std::optional &sg) { - return CanonicalHKL(r.h, r.k, r.l, merge_friedel, sg); -} - bool AcceptReflection(const Reflection &r, std::optional d_min_limit) { if (!std::isfinite(r.I)) return false; diff --git a/image_analysis/scale_merge/HKLKey.h b/image_analysis/scale_merge/HKLKey.h index c03bf092..93c5d505 100644 --- a/image_analysis/scale_merge/HKLKey.h +++ b/image_analysis/scale_merge/HKLKey.h @@ -16,6 +16,19 @@ struct HKLKey { bool operator<(const HKLKey &o) const { return std::tie(h, k, l, plus) < std::tie(o.h, o.k, o.l, o.plus); } + bool operator==(const HKLKey &o) const { + return h == o.h && k == o.k && l == o.l && plus == o.plus; + } +}; + +class HKLKeyGenerator { + bool merge_friedel; + std::optional sg; +public: + HKLKeyGenerator(bool merge_friedel, const std::optional &sg); + HKLKey operator()(const Reflection &r) const; + HKLKey operator()(const MergedReflection &r) const; + HKLKey operator()(int32_t h, int32_t k, int32_t l) const; }; HKLKey CanonicalHKL(const Reflection &r, bool merge_friedel, const std::optional &sg); diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index ddc56f16..cb68acd1 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -41,6 +41,8 @@ namespace { auto scaling_settings = x.GetScalingSettings(); + HKLKeyGenerator key_generator(scaling_settings.GetMergeFriedel(), x.GetGemmiSpaceGroup() ); + for (const auto &image: observations) { for (const auto &r: image) { @@ -49,12 +51,7 @@ namespace { if (!AcceptReflection(r, scaling_settings.GetHighResolutionLimit_A())) continue; - HKLKey key; - try { - key = CanonicalHKL(r, scaling_settings.GetMergeFriedel(), x.GetGemmiSpaceGroup()); - } catch (...) { - continue; - } + HKLKey key = key_generator(r); auto it = hkl_to_slot.find(key); if (it == hkl_to_slot.end()) { diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index bf40c0d4..91e47d6e 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -90,8 +90,9 @@ ScaleOnTheFly::ScaleOnTheFly(const std::vector &ref, const Dif model(x.GetPartialityModel()), s(x.GetScalingSettings()), rot_wedge_deg(GetWedge(x)) { + HKLKeyGenerator key_generator(s.GetMergeFriedel(), sg); for (const auto &r: ref) { - const auto key = CanonicalHKL(r, s.GetMergeFriedel(), sg); + const auto key = key_generator(r); reference_data[key] = r.I; } } @@ -126,9 +127,9 @@ ScaleOnTheFlyResult ScaleOnTheFly::Scale(std::vector &reflections, s }; size_t n_reflections = 0; - + HKLKeyGenerator key_generator(s.GetMergeFriedel(), sg); for (const auto &r: reflections) { - const HKLKey key = CanonicalHKL(r, s.GetMergeFriedel(), sg); + const HKLKey key = key_generator(r); if (!Accept(r)) continue; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index a0e1756b..1761d5ac 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -67,6 +67,7 @@ ADD_EXECUTABLE(jfjoch_test BraggIntegrate2DTest.cpp SearchSpaceGroupTest.cpp XDSPluginTest.cpp + MergeScaleTest.cpp ) target_link_libraries(jfjoch_test Catch2WithMain JFJochBroker JFJochReceiver JFJochReader JFJochWriter diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp new file mode 100644 index 00000000..b821b8d8 --- /dev/null +++ b/tests/MergeScaleTest.cpp @@ -0,0 +1,66 @@ +// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include "../image_analysis/scale_merge/HKLKey.h" + +TEST_CASE("HKLKey_NoSG_noMergeFriedel") { + HKLKeyGenerator hkl_key_gen(false, std::nullopt); + CHECK(hkl_key_gen(-1, -2, -3) != hkl_key_gen(1,2,3)); + CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3)); + CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3)); +} + +TEST_CASE("HKLKey_NoSG_MergeFriedel") { + HKLKeyGenerator hkl_key_gen(true, std::nullopt); + CHECK(hkl_key_gen(-1, -2, -3) == hkl_key_gen(1,2,3)); + CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3)); + CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3)); +} + +TEST_CASE("HKLKey_SG1_MergeFriedel") { + HKLKeyGenerator hkl_key_gen(true, *gemmi::find_spacegroup_by_number(1)); + CHECK(hkl_key_gen(-1, -2, -3) == hkl_key_gen(1,2,3)); + CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3)); + CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3)); +} + +TEST_CASE("HKLKey_SG1_NoMergeFriedel") { + HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(1)); + CHECK(hkl_key_gen(-1, -2, -3) != hkl_key_gen(1,2,3)); + CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3)); + CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3)); +} + +TEST_CASE("HKLKey_SG96_MergeFriedel") { + HKLKeyGenerator hkl_key_gen(true, *gemmi::find_spacegroup_by_number(96)); + CHECK(hkl_key_gen(-1, -2, -3) == hkl_key_gen(1,2,3)); + CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2,1,3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-1,-2,3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2,-1,3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(1,-2,-3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-1,2,-3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2,1,-3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2, -1, -3)); + + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2,-1,3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2, 1, 3)); +} + + +TEST_CASE("HKLKey_SG96_NoMergeFriedel") { + HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(96)); + CHECK(hkl_key_gen(-1, -2, -3) != hkl_key_gen(1,2,3)); + CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2,1,3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-1,-2,3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2,-1,3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(1,-2,-3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-1,2,-3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2,1,-3)); + CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2, -1, -3)); + + CHECK(hkl_key_gen(1,2,3) != hkl_key_gen(-2,-1,3)); + CHECK(hkl_key_gen(1,2,3) != hkl_key_gen(2, 1, 3)); +}