mx/Jungfraujoch
0
0
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases 37 Wiki Activity

jfjoch_extract_hkl: Work in progress (early)

Browse Source
This commit is contained in:
leonarski_f
2026-01-21 10:17:37 +01:00
parent f14eca977d
commit 2000affffa
4 changed files with 125 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
tools/CMakeLists.txt
View File

@@ -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)

48
tools/XdsIntegrateParser.cpp Normal file
View File

@@ -0,0 +1,48 @@
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "XdsIntegrateParser.h"
#include <fstream>
#include <sstream>
#include <stdexcept>
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;
}

31
tools/XdsIntegrateParser.h Normal file
View File

@@ -0,0 +1,31 @@
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
#include <algorithm>
#include <string>
#include <vector>
#include <unordered_map>
#include <cstdint>
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<uint16_t>(std::clamp(h + bias, int64_t(0), int64_t(1024)));
uint64_t uk = static_cast<uint16_t>(std::clamp(k + bias, int64_t(0), int64_t(1024)));
uint64_t ul = static_cast<uint16_t>(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<double> tail; // any remaining numeric columns
};
using IntegrateMap = std::unordered_map<uint64_t, std::vector<HKLData>>;
IntegrateMap ParseXdsIntegrateHkl(const std::string& filename);

81
tools/jfjoch_extract_hkl.cpp
View File

@@ -3,38 +3,27 @@
#include <fstream>
#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 {<options>} <path to master file>");
logger.Info("Options:");
logger.Info(" -I<num> Number of images");
logger.Info(" -o<file> Output filename");
logger.Info(" -x<num> Max dist from Ewald sphere");
logger.Info(" -R Sum same HKL across neighboring images (image +/- 1)");
logger.Info(" -I<num> Number of images");
logger.Info(" -o<file> Output filename");
logger.Info(" -x<num> Max dist from Ewald sphere");
logger.Info(" -R{<file>} 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<uint16_t>(std::clamp(h + bias, int64_t(0), int64_t(1024)));
uint64_t uk = static_cast<uint16_t>(std::clamp(k + bias, int64_t(0), int64_t(1024)));
uint64_t ul = static_cast<uint16_t>(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<int64_t, HKLData> 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();