diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 424ec962..f43016e3 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -39,7 +39,9 @@ ADD_EXECUTABLE(jfjoch_indexing_test jfjoch_indexing_test.cpp) TARGET_LINK_LIBRARIES(jfjoch_indexing_test JFJochReader JFJochImageAnalysis JFJochWriter) INSTALL(TARGETS jfjoch_indexing_test RUNTIME COMPONENT jfjoch) -ADD_EXECUTABLE(jfjoch_extract_hkl jfjoch_extract_hkl.cpp) +ADD_EXECUTABLE(jfjoch_extract_hkl jfjoch_extract_hkl.cpp + XdsIntegrateParser.cpp + XdsIntegrateParser.h) TARGET_LINK_LIBRARIES(jfjoch_extract_hkl JFJochReader) INSTALL(TARGETS jfjoch_extract_hkl RUNTIME COMPONENT viewer) diff --git a/tools/XdsIntegrateParser.cpp b/tools/XdsIntegrateParser.cpp new file mode 100644 index 00000000..8d94dbdb --- /dev/null +++ b/tools/XdsIntegrateParser.cpp @@ -0,0 +1,48 @@ +// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "XdsIntegrateParser.h" + +#include +#include +#include + +IntegrateMap ParseXdsIntegrateHkl(const std::string& filename) { + std::ifstream in(filename); + if (!in.is_open()) { + throw std::runtime_error("Failed to open INTEGRATE.HKL file: " + filename); + } + + IntegrateMap result; + std::string line; + + while (std::getline(in, line)) { + if (line.empty()) continue; + char c0 = line.front(); + if (c0 == '!' || c0 == '#') continue; + + std::istringstream iss(line); + int32_t h, k, l; + double I, sigma; + + if (!(iss >> h >> k >> l >> I >> sigma)) { + continue; + } + + HKLData entry; + entry.h = h; + entry.k = k; + entry.l = l; + entry.I = I; + entry.sigma = sigma; + + double v = 0.0; + while (iss >> v) { + entry.tail.push_back(v); + } + + result[hkl_key_16(h, k, l)].push_back(std::move(entry)); + } + + return result; +} \ No newline at end of file diff --git a/tools/XdsIntegrateParser.h b/tools/XdsIntegrateParser.h new file mode 100644 index 00000000..1d5e4c7e --- /dev/null +++ b/tools/XdsIntegrateParser.h @@ -0,0 +1,31 @@ +// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + + +#pragma once + +#include +#include +#include +#include +#include + +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; + double sigma = 0.0; + float last_image = 0; + std::vector tail; // any remaining numeric columns +}; + +using IntegrateMap = std::unordered_map>; + +IntegrateMap ParseXdsIntegrateHkl(const std::string& filename); \ No newline at end of file diff --git a/tools/jfjoch_extract_hkl.cpp b/tools/jfjoch_extract_hkl.cpp index c6dd9865..126db9fc 100644 --- a/tools/jfjoch_extract_hkl.cpp +++ b/tools/jfjoch_extract_hkl.cpp @@ -3,38 +3,27 @@ #include +#include "XdsIntegrateParser.h" #include "../reader/JFJochHDF5Reader.h" #include "../common/print_license.h" #include "../common/Logger.h" - +#include "XdsIntegrateParser.h" void print_usage(Logger &logger) { logger.Info("Usage ./jfjoch_extract_hkl {} "); logger.Info("Options:"); - logger.Info(" -I Number of images"); - logger.Info(" -o Output filename"); - logger.Info(" -x Max dist from Ewald sphere"); - logger.Info(" -R Sum same HKL across neighboring images (image +/- 1)"); + logger.Info(" -I Number of images"); + logger.Info(" -o Output filename"); + logger.Info(" -x Max dist from Ewald sphere"); + logger.Info(" -R{} Sum same HKL across neighboring images (image +/- 1) - optional reference INTEGRATE.HKL file name"); } -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 { - double I_sum = 0.0; - float first_image = 0; - int64_t h ,k,l; -}; int main(int argc, char **argv) { int64_t image_number = 0; float max_dist_ewald_sphere = 1.0; bool sum_neighboring = false; + std::string ref_file; std::string output_filename = "out.hkl"; @@ -43,7 +32,7 @@ int main(int argc, char **argv) { Logger logger("jfjoch_extract_hkl"); logger.Verbose(true); int opt; - while ((opt = getopt(argc, argv, "I:x:o:R")) != -1) { + while ((opt = getopt(argc, argv, "I:x:o:R::")) != -1) { switch (opt) { case 'I': image_number = atol(optarg); @@ -56,6 +45,8 @@ int main(int argc, char **argv) { break; case 'R': sum_neighboring = true; + if (optarg) + ref_file = std::string(optarg); break; default: /* '?' */ print_usage(logger); @@ -75,37 +66,51 @@ int main(int argc, char **argv) { if (dataset) { if (sum_neighboring) { - std::unordered_map agg; + IntegrateMap agg; int64_t total_images = reader.GetNumberOfImages(); - auto geom = dataset->Dataset().experiment.GetDiffractionGeometry(); + + IntegrateMap xds_result; + if (!ref_file.empty()) + xds_result = ParseXdsIntegrateHkl(ref_file); for (int i = 0; i < total_images; ++i) { auto dataset = reader.LoadImage(i); if (!dataset) continue; for (const auto &r: dataset->ImageData().reflections) { int64_t key = hkl_key_16(r.h, r.k, r.l); - - if (agg.find(key) == agg.end()) { - agg[key] = HKLData{ - .I_sum = r.I, - .first_image = r.image_number, - .h = r.h, .k = r.k, .l = r.l - }; - } else { - // Handle situation of too many images apart - if (std::fabs(agg[key].first_image - r.image_number) > 30.0) - continue; - - agg[key].I_sum += r.I; + auto it = agg.find(key); + HKLData data{ + .h = r.h, + .k = r.k, + .l = r.l, + .I = r.I, + .last_image = r.image_number + }; + bool found = false; + if (it != agg.end()) { + for (auto &val: it->second) { + if (val.last_image == r.image_number - 1) { + val.I += r.I; + val.last_image = r.image_number; + found = true; + break; + } + } } + if (!found) + agg[key].push_back(data); } } std::fstream f(output_filename, std::ios::out); - for (const auto &val: agg | std::views::values) { - double I = val.I_sum; - - f << val.h << " " << val.k << " " << val.l << " " << I << std::endl; + for (const auto &[key, val]: agg) { + if (!val.empty()) { + f << val[0].h << " " << val[0].k << " " << val[0].l << " " << val[0].I; + auto xds_it = xds_result.find(key); + if (xds_it != xds_result.end() && !xds_it->second.empty()) + f << " " << xds_it->second[0].I; + f << std::endl; + } } } else { auto geom = dataset->Dataset().experiment.GetDiffractionGeometry();