diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index 6af20bf7..171cbfb4 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -121,6 +121,7 @@ ADD_LIBRARY(JFJochCommon STATIC DarkMaskSettings.h TopPixels.cpp TopPixels.h + hkl_key.h ) TARGET_LINK_LIBRARIES(JFJochCommon JFJochLogger Compression JFCalibration gemmi Threads::Threads -lrt ) diff --git a/common/hkl_key.h b/common/hkl_key.h new file mode 100644 index 00000000..fdfd4e36 --- /dev/null +++ b/common/hkl_key.h @@ -0,0 +1,39 @@ +// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +# pragma once + +#include +#include "JFJochException.h" + +constexpr int64_t bias = 512; + +inline uint64_t encode(int64_t v) { + const int64_t tmp = v + bias; + if (tmp < 0 || tmp >= 2 * bias) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Value out of range"); + return static_cast(tmp); +} + +inline uint64_t hkl_key(const int64_t h, const int64_t k, const int64_t l) { + const uint64_t uh = encode(h); + const uint64_t uk = encode(k); + const uint64_t ul = encode(l); + return uh | (uk << 11) | (ul << 22); +} + +inline int64_t decode(uint64_t v) { + if (v >= 2 * bias) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Value out of range"); + return static_cast(v) - bias; +} + +inline void hkl_from_key(uint64_t key, int64_t& h, int64_t& k, int64_t& l) { + constexpr uint64_t mask = (1ULL << 11) - 1; + const uint64_t uh = (key >> 0) & mask; + const uint64_t uk = (key >> 11) & mask; + const uint64_t ul = (key >> 22) & mask; + h = decode(uh); + k = decode(uk); + l = decode(ul); +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 41156158..248d7dac 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -63,6 +63,7 @@ ADD_EXECUTABLE(jfjoch_test TimeTest.cpp RotationIndexerTest.cpp TopPixelsTest.cpp + HKLKeyTest.cpp ) target_link_libraries(jfjoch_test Catch2WithMain JFJochBroker JFJochReceiver JFJochReader JFJochWriter JFJochImageAnalysis JFJochCommon JFJochHLSSimulation JFJochPreview) diff --git a/tests/HKLKeyTest.cpp b/tests/HKLKeyTest.cpp new file mode 100644 index 00000000..c147bd5e --- /dev/null +++ b/tests/HKLKeyTest.cpp @@ -0,0 +1,40 @@ +// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include "../common/hkl_key.h" + +TEST_CASE("HKL key round-trip", "[hkl_key]") { + struct Case { int64_t h, k, l; }; + const Case cases[] = { + {0, 0, 0}, + {1, -2, 3}, + {-511, 0, 511}, + {512 - 1, -(512 - 1), 7}, + {-128, 255, -7} + }; + + for (const auto& c : cases) { + const uint64_t key = hkl_key(c.h, c.k, c.l); + int64_t h = 0, k = 0, l = 0; + hkl_from_key(key, h, k, l); + + CHECK(h == c.h); + CHECK(k == c.k); + CHECK(l == c.l); + } +} + +TEST_CASE("HKL key boundaries", "[hkl_key]") { + const int64_t min = -512; + const int64_t max = 511; + + const uint64_t key = hkl_key(min, 0, max); + int64_t h = 0, k = 0, l = 0; + hkl_from_key(key, h, k, l); + + CHECK(h == min); + CHECK(k == 0); + CHECK(l == max); +} \ No newline at end of file diff --git a/tools/XdsIntegrateParser.cpp b/tools/XdsIntegrateParser.cpp index 7db549be..f5f7f363 100644 --- a/tools/XdsIntegrateParser.cpp +++ b/tools/XdsIntegrateParser.cpp @@ -8,6 +8,7 @@ #include #include "../common/ResolutionShells.h" +#include "../common/hkl_key.h" IntegrateMap ParseXdsIntegrateHkl(const std::string& filename) { std::ifstream in(filename); @@ -47,7 +48,7 @@ IntegrateMap ParseXdsIntegrateHkl(const std::string& filename) { entry.tail.push_back(v); } - result[hkl_key_16(-h, k, l)].push_back(std::move(entry)); + result[hkl_key(-h, k, l)].push_back(std::move(entry)); } return result; diff --git a/tools/XdsIntegrateParser.h b/tools/XdsIntegrateParser.h index 776e220b..ea5c19f0 100644 --- a/tools/XdsIntegrateParser.h +++ b/tools/XdsIntegrateParser.h @@ -12,21 +12,12 @@ #include "../common/CrystalLattice.h" - struct CcHalfByResolutionResult { std::vector cc; std::vector pairs; std::vector shell_mean_one_over_d2; }; -inline uint64_t hkl_key_16(int64_t h, int64_t k, int64_t l) { - const uint16_t bias = 512; // maps -512..512 -> 0..1024 - uint64_t uh = static_cast(std::clamp(h + bias, int64_t(0), int64_t(1024))); - uint64_t uk = static_cast(std::clamp(k + bias, int64_t(0), int64_t(1024))); - uint64_t ul = static_cast(std::clamp(l + bias, int64_t(0), int64_t(1024))); - return uh | (uk << 16) | (ul << 32); -} - struct HKLData { int64_t h ,k,l; double I = 0.0;