leonarski_f
91c714f1a4
Build Packages / build:rpm (rocky8_nocuda) (push) Failing after 9m9s
Build Packages / build:rpm (rocky9_nocuda) (push) Failing after 10m59s
Build Packages / build:rpm (ubuntu2404_nocuda) (push) Failing after 8m49s
Build Packages / build:rpm (ubuntu2204_nocuda) (push) Failing after 9m36s
Build Packages / build:rpm (rocky8_sls9) (push) Failing after 10m11s
Build Packages / build:rpm (rocky9_sls9) (push) Failing after 11m6s
Build Packages / build:rpm (rocky8) (push) Failing after 8m55s
Build Packages / build:rpm (rocky9) (push) Failing after 9m6s
Build Packages / Generate python client (push) Successful in 41s
Build Packages / Build documentation (push) Successful in 1m17s
Build Packages / Create release (push) Has been skipped
Build Packages / XDS test (durin plugin) (push) Successful in 9m9s
Build Packages / build:rpm (ubuntu2204) (push) Failing after 11m51s
Build Packages / build:rpm (ubuntu2404) (push) Failing after 12m46s
Build Packages / XDS test (neggia plugin) (push) Successful in 8m46s
Build Packages / XDS test (JFJoch plugin) (push) Successful in 9m32s
Build Packages / DIALS test (push) Successful in 16m2s
Build Packages / Unit tests (push) Failing after 36m36s
1133 lines
49 KiB
C++
1133 lines
49 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <cmath>
|
|
#include "JFJochHDF5Reader.h"
|
|
#include "spdlog/fmt/fmt.h"
|
|
#include "../image_analysis/bragg_integration/CalcISigma.h"
|
|
#include "../image_analysis/spot_finding/SpotUtils.h"
|
|
#include "../common/GridScanSettings.h"
|
|
|
|
inline std::pair<gemmi::CrystalSystem, char> parse_niggli_class(int64_t val) {
|
|
switch (val) {
|
|
case 1: return {gemmi::CrystalSystem::Cubic, 'F'};
|
|
case 2: return {gemmi::CrystalSystem::Trigonal, 'R'};
|
|
case 3: return {gemmi::CrystalSystem::Cubic, 'P'};
|
|
case 5: return {gemmi::CrystalSystem::Cubic, 'I'};
|
|
case 4: return {gemmi::CrystalSystem::Trigonal, 'R'};
|
|
case 6: return {gemmi::CrystalSystem::Tetragonal, 'I'};
|
|
case 7: return {gemmi::CrystalSystem::Tetragonal, 'I'};
|
|
case 9: return {gemmi::CrystalSystem::Trigonal, 'R'};
|
|
case 10: return {gemmi::CrystalSystem::Monoclinic, 'C'};
|
|
case 11: return {gemmi::CrystalSystem::Tetragonal, 'P'};
|
|
case 12: return {gemmi::CrystalSystem::Hexagonal, 'P'};
|
|
case 13: return {gemmi::CrystalSystem::Orthorhombic, 'C'};
|
|
case 15: return {gemmi::CrystalSystem::Tetragonal, 'I'};
|
|
case 16: return {gemmi::CrystalSystem::Orthorhombic, 'F'};
|
|
case 14: return {gemmi::CrystalSystem::Monoclinic, 'C'};
|
|
case 17: return {gemmi::CrystalSystem::Monoclinic, 'C'};
|
|
case 18: return {gemmi::CrystalSystem::Tetragonal, 'I'};
|
|
case 19: return {gemmi::CrystalSystem::Orthorhombic, 'I'};
|
|
case 20: return {gemmi::CrystalSystem::Monoclinic, 'C'};
|
|
case 21: return {gemmi::CrystalSystem::Tetragonal, 'P'};
|
|
case 22: return {gemmi::CrystalSystem::Hexagonal, 'P'};
|
|
case 23: return {gemmi::CrystalSystem::Orthorhombic, 'C'};
|
|
case 24: return {gemmi::CrystalSystem::Trigonal, 'R'};
|
|
case 25: return {gemmi::CrystalSystem::Monoclinic, 'C'};
|
|
case 26: return {gemmi::CrystalSystem::Orthorhombic, 'F'};
|
|
case 27: return {gemmi::CrystalSystem::Monoclinic, 'C'};
|
|
case 28: return {gemmi::CrystalSystem::Monoclinic, 'C'};
|
|
case 29: return {gemmi::CrystalSystem::Monoclinic, 'C'};
|
|
case 30: return {gemmi::CrystalSystem::Monoclinic, 'C'};
|
|
case 31: return {gemmi::CrystalSystem::Triclinic, 'P'};
|
|
case 32: return {gemmi::CrystalSystem::Orthorhombic, 'P'};
|
|
case 40: return {gemmi::CrystalSystem::Orthorhombic, 'C'};
|
|
case 35: return {gemmi::CrystalSystem::Monoclinic, 'P'};
|
|
case 36: return {gemmi::CrystalSystem::Orthorhombic, 'C'};
|
|
case 33: return {gemmi::CrystalSystem::Monoclinic, 'P'};
|
|
case 38: return {gemmi::CrystalSystem::Orthorhombic, 'C'};
|
|
case 34: return {gemmi::CrystalSystem::Monoclinic, 'P'};
|
|
case 42: return {gemmi::CrystalSystem::Orthorhombic, 'I'};
|
|
case 41: return {gemmi::CrystalSystem::Monoclinic, 'C'};
|
|
case 37: return {gemmi::CrystalSystem::Monoclinic, 'C'};
|
|
case 39: return {gemmi::CrystalSystem::Monoclinic, 'C'};
|
|
case 44: return {gemmi::CrystalSystem::Triclinic, 'P'};
|
|
default: return {gemmi::CrystalSystem::Triclinic, 'P'};
|
|
}
|
|
}
|
|
|
|
std::vector<hsize_t> GetDimension(HDF5Object &object, const std::string &path) {
|
|
const auto dim = object.GetDimension(path);
|
|
if (dim.size() != 3)
|
|
throw JFJochException(JFJochExceptionCategory::HDF5, "Wrong dimension of /entry/data/data");
|
|
return dim;
|
|
}
|
|
|
|
std::vector<HDF5VirtualDatasetMapping> ReadVDSImageMappings(HDF5Object &file,
|
|
const std::string &dataset_name) {
|
|
HDF5DataSet dataset(file, dataset_name);
|
|
HDF5Dcpl dcpl(dataset);
|
|
auto mappings = dcpl.GetVirtualMappings();
|
|
|
|
if (mappings.empty())
|
|
throw JFJochException(JFJochExceptionCategory::HDF5,
|
|
dataset_name + " is not a virtual dataset");
|
|
|
|
for (const auto &mapping: mappings) {
|
|
if (mapping.dataset.empty())
|
|
throw JFJochException(JFJochExceptionCategory::HDF5,
|
|
"VDS mapping has empty source dataset name");
|
|
if (mapping.virtual_start.size() != 3)
|
|
throw JFJochException(JFJochExceptionCategory::HDF5,
|
|
"Only 3D image VDS mappings are supported");
|
|
}
|
|
|
|
return mappings;
|
|
}
|
|
|
|
std::string ResolveRelativeToMaster(const std::string &directory,
|
|
const std::string &filename) {
|
|
std::filesystem::path path(filename);
|
|
if (path.is_absolute() || directory.empty())
|
|
return filename;
|
|
|
|
return (std::filesystem::path(directory) / path).string();
|
|
}
|
|
|
|
template<class T>
|
|
void JFJochHDF5Reader::ReadVector(std::vector<T> &v,
|
|
HDF5Object &file,
|
|
const std::string &dataset_name,
|
|
size_t image0,
|
|
size_t nimages) {
|
|
try {
|
|
auto tmp = file.ReadOptVector<T>(dataset_name);
|
|
if (tmp.size() <= nimages) {
|
|
v.resize(image0 + nimages);
|
|
for (int i = 0; i < tmp.size(); i++)
|
|
v[image0 + i] = tmp[i];
|
|
}
|
|
} catch (JFJochException &e) {
|
|
}
|
|
}
|
|
|
|
std::string removeSuffix(const std::string &s, const std::string &suffix) {
|
|
if (s.ends_with(suffix))
|
|
return s.substr(0, s.size() - suffix.size());
|
|
|
|
return s;
|
|
}
|
|
|
|
std::string dataset_name(const std::string &path) {
|
|
std::string file = path;
|
|
int pos = file.rfind('/');
|
|
if (pos != std::string::npos)
|
|
file = file.substr(pos + 1);
|
|
file = removeSuffix(file, "_master.h5");
|
|
// If previous suffix was not found, try removing this one
|
|
file = removeSuffix(file, ".h5");
|
|
return file;
|
|
}
|
|
|
|
bool ReadReflectionsFromGroup(HDF5Object &file,
|
|
const std::string &image_group_name,
|
|
DataMessage &message) {
|
|
if (!file.Exists("/entry/reflections") || !file.Exists(image_group_name))
|
|
return false;
|
|
|
|
auto h = file.ReadOptVector<int32_t>(image_group_name + "/h");
|
|
auto k = file.ReadOptVector<int32_t>(image_group_name + "/k");
|
|
auto l = file.ReadOptVector<int32_t>(image_group_name + "/l");
|
|
auto predicted_x = file.ReadOptVector<float>(image_group_name + "/predicted_x");
|
|
auto predicted_y = file.ReadOptVector<float>(image_group_name + "/predicted_y");
|
|
|
|
auto d = file.ReadOptVector<float>(image_group_name + "/d");
|
|
auto int_sum = file.ReadOptVector<float>(image_group_name + "/int_sum");
|
|
auto int_err = file.ReadOptVector<float>(image_group_name + "/int_err");
|
|
auto bkg = file.ReadOptVector<float>(image_group_name + "/background_mean");
|
|
auto lp = file.ReadOptVector<float>(image_group_name + "/lp");
|
|
auto partiality = file.ReadOptVector<float>(image_group_name + "/partiality");
|
|
auto phi = file.ReadOptVector<float>(image_group_name + "/delta_phi");
|
|
auto zeta = file.ReadOptVector<float>(image_group_name + "/zeta");
|
|
|
|
if (h.size() != l.size() || h.size() != k.size() || h.size() != d.size()
|
|
|| h.size() != predicted_x.size() || h.size() != predicted_y.size()
|
|
|| h.size() != int_sum.size() || h.size() != int_err.size() || h.size() != bkg.size())
|
|
throw JFJochException(JFJochExceptionCategory::HDF5, "Wrong size of reflections dataset");
|
|
|
|
for (size_t i = 0; i < h.size(); i++) {
|
|
float lp_val = 0.0;
|
|
if (lp.size() > i && lp[i] != 0.0f)
|
|
lp_val = 1.0f / lp[i];
|
|
|
|
float partiality_val = -1.0f;
|
|
if (partiality.size() > i && partiality[i] >= 0.0f)
|
|
partiality_val = partiality[i];
|
|
float delta_phi_val = NAN;
|
|
if (phi.size() > i)
|
|
delta_phi_val = phi[i];
|
|
float zeta_val = NAN;
|
|
if (zeta.size() > i)
|
|
zeta_val = zeta[i];
|
|
|
|
Reflection r{
|
|
.h = h.at(i),
|
|
.k = k.at(i),
|
|
.l = l.at(i),
|
|
.delta_phi_deg = delta_phi_val,
|
|
.predicted_x = predicted_x.at(i),
|
|
.predicted_y = predicted_y.at(i),
|
|
.d = d.at(i),
|
|
.I = int_sum.at(i),
|
|
.bkg = bkg.at(i),
|
|
.sigma = int_err.at(i),
|
|
.rlp = lp_val,
|
|
.partiality = partiality_val,
|
|
.zeta = zeta_val
|
|
};
|
|
message.reflections.emplace_back(r);
|
|
}
|
|
|
|
CalcISigma(message);
|
|
CalcWilsonBFactor(message, !message.b_factor.has_value());
|
|
return true;
|
|
}
|
|
|
|
template<class T>
|
|
std::optional<T> ReadElementMasterFirst(HDF5Object &master_file,
|
|
HDF5Object &source_file,
|
|
const std::string &path,
|
|
hsize_t master_image,
|
|
hsize_t source_image) {
|
|
if (master_file.Exists(path))
|
|
return master_file.ReadElement<T>(path, master_image);
|
|
if (source_file.Exists(path))
|
|
return source_file.ReadElement<T>(path, source_image);
|
|
return {};
|
|
}
|
|
|
|
template<class T>
|
|
std::vector<T> ReadVectorMasterFirst(HDF5Object &master_file,
|
|
HDF5Object &source_file,
|
|
const std::string &path,
|
|
const std::vector<hsize_t> &master_start,
|
|
const std::vector<hsize_t> &source_start,
|
|
const std::vector<hsize_t> &size) {
|
|
if (master_file.Exists(path))
|
|
return master_file.ReadOptVector<T>(path, master_start, size);
|
|
if (source_file.Exists(path))
|
|
return source_file.ReadOptVector<T>(path, source_start, size);
|
|
return {};
|
|
}
|
|
|
|
void JFJochHDF5Reader::ReadFile(const std::string &filename) {
|
|
std::unique_lock ul(hdf5_mutex);
|
|
try {
|
|
auto dataset = std::make_shared<JFJochReaderDataset>();
|
|
master_file = std::make_shared<HDF5ReadOnlyFile>(filename);
|
|
master_filename = filename;
|
|
dataset->experiment = default_experiment;
|
|
|
|
std::filesystem::path master_path(filename);
|
|
master_file_directory = master_path.parent_path().string();
|
|
vds_data_mappings.clear();
|
|
|
|
dataset->arm_date = master_file->GetString("/entry/start_time");
|
|
|
|
dataset->experiment.FilePrefix(dataset_name(filename));
|
|
|
|
// JFJochReader is always using int32_t
|
|
dataset->experiment.BitDepthImage(32);
|
|
dataset->experiment.PixelSigned(true);
|
|
|
|
size_t image_size_x = 0;
|
|
size_t image_size_y = 0;
|
|
|
|
if (master_file->Exists("/entry/data/data")) {
|
|
HDF5DataSet data_dataset(*master_file, "/entry/data/data");
|
|
HDF5Dcpl dcpl(data_dataset);
|
|
data_layout = dcpl.GetLayout();
|
|
|
|
auto dim = GetDimension(*master_file, "/entry/data/data");
|
|
number_of_images = dim[0];
|
|
image_size_y = dim[1];
|
|
image_size_x = dim[2];
|
|
|
|
images_per_file = number_of_images;
|
|
if (data_layout == HDF5DataSetLayout::VIRTUAL)
|
|
vds_data_mappings = ReadVDSImageMappings(*master_file, "/entry/data/data");
|
|
|
|
if (master_file->Exists("/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image"))
|
|
dataset->efficiency = master_file->ReadVector<float>(
|
|
"/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image");
|
|
else
|
|
dataset->efficiency = std::vector<float>(number_of_images, 1.0);
|
|
|
|
if (master_file->Exists("/entry/roi"))
|
|
dataset->roi = master_file->FindLeafs("/entry/roi");
|
|
|
|
for (const auto &s: dataset->roi) {
|
|
dataset->roi_max.emplace_back(master_file->ReadVector<int64_t>("/entry/roi/" + s + "/max"));
|
|
dataset->roi_sum.emplace_back(master_file->ReadVector<int64_t>("/entry/roi/" + s + "/sum"));
|
|
dataset->roi_sum_sq.emplace_back(master_file->ReadVector<int64_t>("/entry/roi/" + s + "/sum_sq"));
|
|
dataset->roi_npixel.emplace_back(master_file->ReadVector<int64_t>("/entry/roi/" + s + "/npixel"));
|
|
dataset->roi_x.emplace_back(master_file->ReadVector<float>("/entry/roi/" + s + "/x"));
|
|
dataset->roi_y.emplace_back(master_file->ReadVector<float>("/entry/roi/" + s + "/y"));
|
|
}
|
|
|
|
if (master_file->Exists("/entry/MX")) {
|
|
if (master_file->Exists("/entry/MX/peakCountUnfiltered"))
|
|
dataset->spot_count = master_file->ReadOptVector<float>("/entry/MX/peakCountUnfiltered");
|
|
else
|
|
dataset->spot_count = master_file->ReadOptVector<float>("/entry/MX/nPeaks");
|
|
|
|
dataset->spot_count_low_res = master_file->ReadOptVector<float>("/entry/MX/peakCountLowRes");
|
|
dataset->spot_count_indexed = master_file->ReadOptVector<float>("/entry/MX/peakCountIndexed");
|
|
dataset->spot_count_ice_rings = master_file->ReadOptVector<float>("/entry/MX/peakCountIceRingRes");
|
|
|
|
dataset->indexing_result = master_file->ReadOptVector<float>("/entry/MX/imageIndexed");
|
|
dataset->bkg_estimate = master_file->ReadOptVector<float>("/entry/MX/bkgEstimate");
|
|
dataset->resolution_estimate = master_file->ReadOptVector<float>("/entry/MX/resolutionEstimate");
|
|
dataset->profile_radius = master_file->ReadOptVector<float>("/entry/MX/profileRadius");
|
|
dataset->mosaicity_deg = master_file->ReadOptVector<float>("/entry/MX/mosaicity");
|
|
dataset->b_factor = master_file->ReadOptVector<float>("/entry/MX/bFactor");
|
|
dataset->scale_factor = master_file->ReadOptVector<float>("/entry/MX/imageScaleFactor");
|
|
dataset->integrated_reflections = master_file->ReadOptVector<float>("/entry/MX/integratedReflections");
|
|
}
|
|
if (master_file->Exists("/entry/image"))
|
|
dataset->max_value = master_file->ReadOptVector<int64_t>("/entry/image/max_value");
|
|
|
|
format = FileWriterFormat::NXmxVDS;
|
|
} else if (master_file->Exists("/entry/data/data_000001")) {
|
|
format = FileWriterFormat::NXmxLegacy;
|
|
data_layout = HDF5DataSetLayout::CONTIGUOUS;
|
|
|
|
legacy_format_files.clear();
|
|
|
|
image_size_x = master_file->GetInt("/entry/instrument/detector/detectorSpecific/x_pixels_in_detector");
|
|
image_size_y = master_file->GetInt("/entry/instrument/detector/detectorSpecific/y_pixels_in_detector");
|
|
|
|
//size_t expected_images = master_file->GetInt("/entry/instrument/detector/detectorSpecific/nimages");
|
|
|
|
images_per_file = 0;
|
|
number_of_images = 0;
|
|
uint32_t nfiles = 0;
|
|
|
|
std::filesystem::path file_path(filename);
|
|
std::filesystem::path directory = file_path.parent_path();
|
|
|
|
while (true) {
|
|
std::string dname = fmt::format("/entry/data/data_{:06d}", nfiles + 1);
|
|
if (!master_file->Exists(dname))
|
|
break;
|
|
|
|
size_t fimages = 0;
|
|
|
|
try {
|
|
auto fname = master_file->GetLinkedFileName(dname);
|
|
|
|
if (!directory.empty())
|
|
fname = fmt::format("{}/{}", directory.string(), fname);
|
|
|
|
HDF5ReadOnlyFile data_file(fname);
|
|
|
|
fimages = GetDimension(data_file, "/entry/data/data")[0];
|
|
|
|
legacy_format_files.push_back(fname);
|
|
|
|
if (nfiles == 0 && data_file.Exists("/entry/roi"))
|
|
dataset->roi = data_file.FindLeafs("/entry/roi");
|
|
|
|
dataset->roi_max.resize(dataset->roi.size());
|
|
dataset->roi_npixel.resize(dataset->roi.size());
|
|
dataset->roi_sum.resize(dataset->roi.size());
|
|
dataset->roi_sum_sq.resize(dataset->roi.size());
|
|
dataset->roi_x.resize(dataset->roi.size());
|
|
dataset->roi_y.resize(dataset->roi.size());
|
|
|
|
for (int i = 0; i < dataset->roi.size(); i++) {
|
|
auto roi_name = dataset->roi[i];
|
|
ReadVector(dataset->roi_max.at(i),
|
|
data_file, "/entry/roi/" + roi_name + "/max",
|
|
number_of_images, fimages);
|
|
ReadVector(dataset->roi_npixel.at(i),
|
|
data_file, "/entry/roi/" + roi_name + "/npixel",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->roi_sum.at(i),
|
|
data_file, "/entry/roi/" + roi_name + "/sum",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->roi_sum_sq.at(i),
|
|
data_file, "/entry/roi/" + roi_name + "/sum_sq",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->roi_x.at(i),
|
|
data_file, "/entry/roi/" + roi_name + "/x",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->roi_y.at(i),
|
|
data_file, "/entry/roi/" + roi_name + "/y",
|
|
number_of_images, fimages);
|
|
}
|
|
|
|
if (data_file.Exists("/entry/detector")) {
|
|
ReadVector(dataset->efficiency,
|
|
data_file, "/entry/detector/data_collection_efficiency_image",
|
|
number_of_images, fimages);
|
|
}
|
|
|
|
if (data_file.Exists("/entry/MX")) {
|
|
if (data_file.Exists("/entry/MX/peakCountUnfiltered"))
|
|
ReadVector(dataset->spot_count,
|
|
data_file, "/entry/MX/peakCountUnfiltered",
|
|
number_of_images, fimages);
|
|
else
|
|
ReadVector(dataset->spot_count,
|
|
data_file, "/entry/MX/nPeaks",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->spot_count_ice_rings,
|
|
data_file, "/entry/MX/peakCountIceRingRes",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->spot_count_low_res,
|
|
data_file, "/entry/MX/peakCountLowRes",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->spot_count_indexed,
|
|
data_file, "/entry/MX/peakCountIndexed",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->indexing_result,
|
|
data_file, "/entry/MX/imageIndexed",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->bkg_estimate,
|
|
data_file, "/entry/MX/bkgEstimate",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->profile_radius,
|
|
data_file, "/entry/MX/profileRadius",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->mosaicity_deg,
|
|
data_file, "/entry/MX/mosaicity",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->b_factor,
|
|
data_file, "/entry/MX/bFactor",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->resolution_estimate,
|
|
data_file, "/entry/MX/resolutionEstimate",
|
|
number_of_images, fimages);
|
|
}
|
|
|
|
if (data_file.Exists("/entry/image")) {
|
|
ReadVector(dataset->max_value,
|
|
data_file, "/entry/image/max_value",
|
|
number_of_images, fimages);
|
|
}
|
|
} catch (JFJochException &e) {
|
|
}
|
|
|
|
if (nfiles == 0)
|
|
images_per_file = fimages;
|
|
number_of_images += fimages;
|
|
nfiles++;
|
|
}
|
|
} else {
|
|
image_size_x = master_file->GetInt("/entry/instrument/detector/detectorSpecific/x_pixels_in_detector");
|
|
image_size_y = master_file->GetInt("/entry/instrument/detector/detectorSpecific/y_pixels_in_detector");
|
|
number_of_images = 0;
|
|
}
|
|
|
|
if (master_file->Exists("/entry/MX")) {
|
|
auto indexing = master_file->GetString("/entry/MX/indexing_algorithm", "none");
|
|
if (indexing == "fft")
|
|
dataset->experiment.IndexingAlgorithm(IndexingAlgorithmEnum::FFT);
|
|
else if (indexing == "ffbidx")
|
|
dataset->experiment.IndexingAlgorithm(IndexingAlgorithmEnum::FFBIDX);
|
|
}
|
|
|
|
auto ring_current_A = master_file->GetOptFloat("/entry/source/current");
|
|
if (ring_current_A) dataset->experiment.RingCurrent_mA(ring_current_A.value() * 1000.0);
|
|
|
|
dataset->experiment.DetectIceRings(
|
|
master_file->GetOptBool("/entry/instrument/detector/detectorSpecific/detect_ice_rings").value_or(false));
|
|
dataset->experiment.PoniRot1_rad(
|
|
master_file->GetOptFloat("/entry/instrument/detector/transformations/rot1").value_or(0.0));
|
|
dataset->experiment.PoniRot2_rad(
|
|
master_file->GetOptFloat("/entry/instrument/detector/transformations/rot2").value_or(0.0));
|
|
dataset->experiment.PoniRot3_rad(
|
|
master_file->GetOptFloat("/entry/instrument/detector/transformations/rot3").value_or(0.0));
|
|
dataset->experiment.SampleTemperature_K(master_file->GetOptFloat("/entry/sample/temperature"));
|
|
|
|
dataset->experiment.BeamX_pxl(master_file->GetFloat("/entry/instrument/detector/beam_center_x"));
|
|
dataset->experiment.BeamY_pxl(master_file->GetFloat("/entry/instrument/detector/beam_center_y"));
|
|
|
|
float det_distance = master_file->GetFloat("/entry/instrument/detector/distance");
|
|
if (det_distance < 0.001)
|
|
det_distance = 0.1; // Set to 100 mm, if det distance is less than 1 mm
|
|
dataset->experiment.DetectorDistance_mm(det_distance * 1000.0);
|
|
|
|
dataset->experiment.IncidentEnergy_keV(
|
|
WVL_1A_IN_KEV / master_file->GetFloat("/entry/instrument/beam/incident_wavelength"));
|
|
|
|
dataset->error_value = master_file->GetOptInt("/entry/instrument/detector/error_value");
|
|
|
|
dataset->jfjoch_release = master_file->GetString("/entry/instrument/detector/detectorSpecific/jfjoch_release");
|
|
|
|
InstrumentMetadata metadata;
|
|
metadata.InstrumentName(master_file->GetString("/entry/instrument/name"));
|
|
metadata.SourceName(master_file->GetString("/entry/source/name"));
|
|
dataset->experiment.ImportInstrumentMetadata(metadata);
|
|
|
|
if (master_file->Exists("/entry/sample/transformations")) {
|
|
if (master_file->Exists("/entry/sample/transformations/omega")) {
|
|
auto omega = ReadAxis(master_file.get(), "omega");
|
|
dataset->experiment.Goniometer(omega);
|
|
} else if (master_file->Exists("/entry/sample/grid_scan")) {
|
|
GridScanSettings grid(
|
|
master_file->GetInt("/entry/sample/grid_scan/n_fast"),
|
|
master_file->GetFloat("/entry/sample/grid_scan/step_x") * 1e6f,
|
|
master_file->GetFloat("/entry/sample/grid_scan/step_y") * 1e6f,
|
|
master_file->GetOptBool("/entry/sample/grid_scan/snake_scan").value_or(false),
|
|
master_file->GetOptBool("/entry/sample/grid_scan/vertical_scan").value_or(false)
|
|
);
|
|
grid.ImageNum(number_of_images);
|
|
dataset->experiment.GridScan(grid);
|
|
}
|
|
}
|
|
|
|
auto tmp = master_file->ReadOptVector<float>("/entry/sample/unit_cell");
|
|
if (tmp.size() == 6)
|
|
dataset->experiment.SetUnitCell(UnitCell{
|
|
.a = tmp[0],
|
|
.b = tmp[1],
|
|
.c = tmp[2],
|
|
.alpha = tmp[3],
|
|
.beta = tmp[4],
|
|
.gamma = tmp[5]
|
|
});
|
|
dataset->experiment.SpaceGroupNumber(master_file->GetOptInt("/entry/sample/space_group_number"));
|
|
dataset->experiment.SampleName(master_file->GetString("/entry/sample/name"));
|
|
|
|
|
|
if (master_file->Exists("/entry/instrument/attenuator"))
|
|
dataset->experiment.AttenuatorTransmission(
|
|
master_file->GetOptFloat("/entry/instrument/attenuator/attenuator_transmission"));
|
|
dataset->experiment.TotalFlux(master_file->GetOptFloat("/entry/instrument/beam/total_flux"));
|
|
|
|
if (master_file->Exists("/entry/azint") && master_file->Exists("/entry/azint/bin_to_q")) {
|
|
HDF5DataSet bin_to_q_dataset(*master_file, "/entry/azint/bin_to_q");
|
|
HDF5DataSpace bin_to_q_dataspace(bin_to_q_dataset);
|
|
auto dim = bin_to_q_dataspace.GetDimensions();
|
|
|
|
if (dim.size() == 1) {
|
|
dataset->azimuthal_bins = 0;
|
|
dataset->q_bins = dim[0];
|
|
bin_to_q_dataset.ReadVector(dataset->az_int_bin_to_q);
|
|
} else if (dim.size() == 2) {
|
|
dataset->azimuthal_bins = dim[0];
|
|
dataset->q_bins = dim[1];
|
|
dataset->az_int_bin_to_q.resize(dim[0] * dim[1]);
|
|
bin_to_q_dataset.ReadVector(dataset->az_int_bin_to_q, {0, 0}, dim);
|
|
} else
|
|
throw JFJochException(JFJochExceptionCategory::HDF5, "Wrong dimension of /entry/azint/image dataset");
|
|
if (master_file->Exists("/entry/azint/bin_to_phi")) {
|
|
HDF5DataSet bin_to_phi_dataset(*master_file, "/entry/azint/bin_to_phi");
|
|
if (dataset->q_bins > 0) {
|
|
dataset->az_int_bin_to_phi.resize(dim[0] * dim[1]);
|
|
bin_to_phi_dataset.ReadVector(dataset->az_int_bin_to_phi, {0, 0}, dim);
|
|
} else {
|
|
bin_to_phi_dataset.ReadVector(dataset->az_int_bin_to_phi);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Read fluorescence spectrum if present
|
|
if (master_file->Exists("/entry/instrument/fluorescence")) {
|
|
auto energy = master_file->ReadOptVector<float>("/entry/instrument/fluorescence/energy");
|
|
auto data = master_file->ReadOptVector<float>("/entry/instrument/fluorescence/data");
|
|
if (!energy.empty() && energy.size() == data.size())
|
|
dataset->experiment.FluorescenceSpectrum(XrayFluorescenceSpectrum(energy, data));
|
|
}
|
|
|
|
auto detector_name = master_file->GetString("/entry/instrument/detector/description");
|
|
|
|
DetectorSetup detector = DetDECTRIS(image_size_x, image_size_y, detector_name, {});
|
|
detector.PixelSize_um(master_file->GetFloat("/entry/instrument/detector/x_pixel_size") * 1e6);
|
|
detector.SaturationLimit(master_file->GetInt("/entry/instrument/detector/saturation_value"));
|
|
detector.MinFrameTime(std::chrono::microseconds(0));
|
|
detector.MinCountTime(std::chrono::microseconds(0));
|
|
detector.ReadOutTime(std::chrono::nanoseconds(0));
|
|
dataset->experiment.Detector(detector);
|
|
|
|
dataset->experiment.FrameTime(
|
|
std::chrono::duration_cast<std::chrono::nanoseconds>(
|
|
std::chrono::duration<float>(
|
|
master_file->GetFloat("/entry/instrument/detector/frame_time"))),
|
|
std::chrono::duration_cast<std::chrono::nanoseconds>(
|
|
std::chrono::duration<float>(
|
|
master_file->GetFloat("/entry/instrument/detector/count_time")))
|
|
);
|
|
|
|
if (master_file->Exists("/entry/instrument/detector/calibration")) {
|
|
dataset->calibration_data = master_file->FindLeafs("/entry/instrument/detector/calibration");
|
|
std::sort(dataset->calibration_data.begin(), dataset->calibration_data.end());
|
|
}
|
|
|
|
if (image_size_x * image_size_y > 0) {
|
|
auto mask_tmp = master_file->ReadOptVector<uint32_t>(
|
|
"/entry/instrument/detector/pixel_mask",
|
|
{0, 0},
|
|
{image_size_y, image_size_x}
|
|
);
|
|
if (mask_tmp.empty())
|
|
mask_tmp = std::vector<uint32_t>(image_size_x * image_size_y);
|
|
dataset->pixel_mask = PixelMask(mask_tmp);
|
|
}
|
|
dataset->experiment.ImagesPerTrigger(number_of_images);
|
|
SetStartMessage(dataset);
|
|
} catch (const std::exception &e) {
|
|
master_file = {};
|
|
number_of_images = 0;
|
|
SetStartMessage({});
|
|
throw;
|
|
}
|
|
}
|
|
|
|
uint64_t JFJochHDF5Reader::GetNumberOfImages() const {
|
|
std::unique_lock ul(hdf5_mutex);
|
|
return number_of_images;
|
|
}
|
|
|
|
CompressedImage JFJochHDF5Reader::LoadImageDataset(std::vector<uint8_t> &tmp, HDF5Object &file, hsize_t number) {
|
|
std::vector<hsize_t> start = {static_cast<hsize_t>(number), 0, 0};
|
|
|
|
const std::string dataset_name = "/entry/data/data";
|
|
|
|
HDF5DataSet dataset(file, dataset_name);
|
|
HDF5DataSpace dataspace(dataset);
|
|
HDF5DataType datatype(dataset);
|
|
HDF5Dcpl dcpl(dataset);
|
|
|
|
if (dataspace.GetNumOfDimensions() != 3)
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"/entry/data/data dataset must be 3D");
|
|
|
|
auto dim = dataspace.GetDimensions();
|
|
|
|
CompressionAlgorithm algorithm = CompressionAlgorithm::NO_COMPRESSION;
|
|
auto chunk_size = dcpl.GetChunking();
|
|
|
|
if ((chunk_size.size() == 3) && (chunk_size[0] == 1) && (chunk_size[1] == dim[1]) && (chunk_size[2] == dim[2])) {
|
|
dataset.ReadDirectChunk(tmp, start);
|
|
algorithm = dcpl.GetCompression();
|
|
} else {
|
|
dataset.ReadVectorToU8(tmp, start, {1, dim[1], dim[2]});
|
|
algorithm = CompressionAlgorithm::NO_COMPRESSION;
|
|
}
|
|
|
|
if (datatype.IsFloat())
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Float datasets not supported at this time");
|
|
|
|
return {
|
|
tmp, dim[2], dim[1],
|
|
CalcImageMode(datatype.GetElemSize(), datatype.IsFloat(), datatype.IsSigned()),
|
|
algorithm
|
|
};
|
|
}
|
|
|
|
std::pair<std::shared_ptr<HDF5ReadOnlyFile>, uint32_t> JFJochHDF5Reader::GetImageLocation(int64_t image_number) {
|
|
if (image_number >= number_of_images || image_number < 0)
|
|
throw JFJochException(JFJochExceptionCategory::HDF5, "Image out of bounds");
|
|
|
|
uint32_t image_id;
|
|
std::shared_ptr<HDF5ReadOnlyFile> data_file;
|
|
|
|
if (format == FileWriterFormat::NXmxLegacy) {
|
|
uint32_t file_id = image_number / images_per_file;
|
|
image_id = image_number % images_per_file;
|
|
data_file = std::make_shared<HDF5ReadOnlyFile>(legacy_format_files.at(file_id));
|
|
} else if (format == FileWriterFormat::NXmxVDS && data_layout == HDF5DataSetLayout::VIRTUAL) {
|
|
const auto image = static_cast<hsize_t>(image_number);
|
|
|
|
for (const auto &mapping: vds_data_mappings) {
|
|
if (!mapping.ContainsVirtualImage(image))
|
|
continue;
|
|
|
|
image_id = static_cast<uint32_t>(mapping.SourceImage(image));
|
|
data_file = std::make_shared<HDF5ReadOnlyFile>(
|
|
ResolveRelativeToMaster(master_file_directory, mapping.filename)
|
|
);
|
|
return {std::move(data_file), image_id};
|
|
}
|
|
|
|
throw JFJochException(JFJochExceptionCategory::HDF5,
|
|
"Image not covered by /entry/data/data VDS mappings");
|
|
} else {
|
|
image_id = image_number;
|
|
data_file = master_file;
|
|
}
|
|
|
|
return {std::move(data_file), image_id};
|
|
}
|
|
|
|
std::shared_ptr<JFJochReaderRawImage> JFJochHDF5Reader::GetRawImage(int64_t image_number) {
|
|
std::unique_lock ul(hdf5_mutex);
|
|
|
|
if (!master_file)
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Cannot load image if file not loaded");
|
|
|
|
auto [source_file, image_id] = GetImageLocation(image_number);
|
|
auto ret = std::make_shared<JFJochReaderRawImage>();
|
|
ret->image = LoadImageDataset(ret->image_buffer, *source_file, image_id);
|
|
return ret;
|
|
}
|
|
|
|
bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &dataset,
|
|
DataMessage &message,
|
|
std::vector<uint8_t> &buffer,
|
|
int64_t image_number,
|
|
bool update_dataset) {
|
|
std::unique_lock ul(hdf5_mutex);
|
|
|
|
if (!dataset)
|
|
return {};
|
|
|
|
if (!master_file)
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Cannot load image if file not loaded");
|
|
|
|
auto [source_file, image_id] = GetImageLocation(image_number);
|
|
|
|
message.image = LoadImageDataset(buffer, *source_file, image_id);
|
|
message.number = image_number;
|
|
|
|
const auto master_image = static_cast<hsize_t>(image_number);
|
|
const auto source_image = static_cast<hsize_t>(image_id);
|
|
|
|
auto spot_count_opt = ReadElementMasterFirst<uint32_t>(*master_file,
|
|
*source_file,
|
|
"/entry/MX/nPeaks",
|
|
master_image,
|
|
source_image);
|
|
|
|
if (spot_count_opt.has_value() && spot_count_opt.value() > 0) {
|
|
size_t spot_count = spot_count_opt.value();
|
|
|
|
auto spot_x = ReadVectorMasterFirst<float>(
|
|
*master_file,
|
|
*source_file,
|
|
"/entry/MX/peakXPosRaw",
|
|
{master_image, 0},
|
|
{source_image, 0},
|
|
{1, spot_count}
|
|
);
|
|
auto spot_y = ReadVectorMasterFirst<float>(
|
|
*master_file,
|
|
*source_file,
|
|
"/entry/MX/peakYPosRaw",
|
|
{master_image, 0},
|
|
{source_image, 0},
|
|
{1, spot_count}
|
|
);
|
|
auto spot_intensity = ReadVectorMasterFirst<float>(
|
|
*master_file,
|
|
*source_file,
|
|
"/entry/MX/peakTotalIntensity",
|
|
{master_image, 0},
|
|
{source_image, 0},
|
|
{1, spot_count}
|
|
);
|
|
|
|
if (spot_x.size() < spot_count || spot_y.size() < spot_count || spot_intensity.size() < spot_count)
|
|
throw JFJochException(JFJochExceptionCategory::HDF5, "Wrong size of spot dataset");
|
|
|
|
auto spot_indexed = ReadVectorMasterFirst<uint8_t>(
|
|
*master_file,
|
|
*source_file,
|
|
"/entry/MX/peakIndexed",
|
|
{master_image, 0},
|
|
{source_image, 0},
|
|
{1, spot_count}
|
|
);
|
|
|
|
auto spot_ice = ReadVectorMasterFirst<uint8_t>(
|
|
*master_file,
|
|
*source_file,
|
|
"/entry/MX/peakIceRingRes",
|
|
{master_image, 0},
|
|
{source_image, 0},
|
|
{1, spot_count}
|
|
);
|
|
|
|
auto spot_h = ReadVectorMasterFirst<int32_t>(
|
|
*master_file,
|
|
*source_file,
|
|
"/entry/MX/peakH",
|
|
{master_image, 0},
|
|
{source_image, 0},
|
|
{1, spot_count}
|
|
);
|
|
|
|
auto spot_k = ReadVectorMasterFirst<int32_t>(
|
|
*master_file,
|
|
*source_file,
|
|
"/entry/MX/peakK",
|
|
{master_image, 0},
|
|
{source_image, 0},
|
|
{1, spot_count}
|
|
);
|
|
|
|
auto spot_l = ReadVectorMasterFirst<int32_t>(
|
|
*master_file,
|
|
*source_file,
|
|
"/entry/MX/peakL",
|
|
{master_image, 0},
|
|
{source_image, 0},
|
|
{1, spot_count}
|
|
);
|
|
|
|
auto spot_dist_ewald_sphere = ReadVectorMasterFirst<float>(
|
|
*master_file,
|
|
*source_file,
|
|
"/entry/MX/peakDistEwaldSphere",
|
|
{master_image, 0},
|
|
{source_image, 0},
|
|
{1, spot_count}
|
|
);
|
|
|
|
auto geom = dataset->experiment.GetDiffractionGeometry();
|
|
for (int i = 0; i < spot_count; i++) {
|
|
auto x = spot_x.at(i);
|
|
auto y = spot_y.at(i);
|
|
auto d = geom.PxlToRes(x, y);
|
|
|
|
SpotToSave s{
|
|
.x = x,
|
|
.y = y,
|
|
.intensity = spot_intensity.at(i),
|
|
.d_A = d,
|
|
.image = image_number
|
|
};
|
|
if (spot_indexed.size() > i)
|
|
s.indexed = (spot_indexed.at(i) != 0);
|
|
if (spot_h.size() > i)
|
|
s.h = spot_h.at(i);
|
|
if (spot_k.size() > i)
|
|
s.k = spot_k.at(i);
|
|
if (spot_l.size() > i)
|
|
s.l = spot_l.at(i);
|
|
if (spot_dist_ewald_sphere.size() > i)
|
|
s.dist_ewald_sphere = spot_dist_ewald_sphere.at(i);
|
|
if (spot_ice.size() > i)
|
|
s.ice_ring = (spot_ice.at(i) != 0);
|
|
message.spots.emplace_back(s);
|
|
}
|
|
|
|
if (auto v = ReadElementMasterFirst<uint32_t>(*master_file, *source_file,
|
|
"/entry/MX/peakCountUnfiltered",
|
|
master_image, source_image); v)
|
|
message.spot_count = v;
|
|
else
|
|
message.spot_count = spot_count_opt;
|
|
|
|
message.spot_count_ice_rings = ReadElementMasterFirst<uint32_t>(
|
|
*master_file, *source_file, "/entry/MX/peakCountIceRingRes", master_image, source_image);
|
|
message.spot_count_low_res = ReadElementMasterFirst<uint32_t>(
|
|
*master_file, *source_file, "/entry/MX/peakCountLowRes", master_image, source_image);
|
|
message.spot_count_indexed = ReadElementMasterFirst<uint32_t>(
|
|
*master_file, *source_file, "/entry/MX/peakCountIndexed", master_image, source_image);
|
|
|
|
GenerateSpotPlot(message, 1.5);
|
|
}
|
|
|
|
if (!dataset->az_int_bin_to_q.empty()) {
|
|
if (dataset->azimuthal_bins == 0) {
|
|
message.az_int_profile = ReadVectorMasterFirst<float>(
|
|
*master_file,
|
|
*source_file,
|
|
"/entry/azint/image",
|
|
{master_image, 0},
|
|
{source_image, 0},
|
|
{1, dataset->az_int_bin_to_q.size()}
|
|
);
|
|
} else {
|
|
message.az_int_profile = ReadVectorMasterFirst<float>(
|
|
*master_file,
|
|
*source_file,
|
|
"/entry/azint/image",
|
|
{master_image, 0, 0},
|
|
{source_image, 0, 0},
|
|
{1, dataset->azimuthal_bins, dataset->q_bins}
|
|
);
|
|
}
|
|
}
|
|
if (dataset->integrated_reflections.size() > image_number)
|
|
message.integrated_reflections = static_cast<int64_t>(std::lround(dataset->integrated_reflections.at(image_number)));
|
|
if (dataset->resolution_estimate.size() > image_number)
|
|
message.resolution_estimate = dataset->resolution_estimate[image_number];
|
|
if (dataset->indexing_result.size() > image_number)
|
|
message.indexing_result = dataset->indexing_result[image_number];
|
|
if (dataset->bkg_estimate.size() > image_number)
|
|
message.bkg_estimate = dataset->bkg_estimate[image_number];
|
|
if (dataset->efficiency.size() > image_number)
|
|
message.image_collection_efficiency = dataset->efficiency[image_number];
|
|
if (dataset->profile_radius.size() > image_number)
|
|
message.profile_radius = dataset->profile_radius[image_number];
|
|
if (dataset->mosaicity_deg.size() > image_number)
|
|
message.mosaicity_deg = dataset->mosaicity_deg[image_number];
|
|
if (dataset->b_factor.size() > image_number)
|
|
message.b_factor = dataset->b_factor[image_number];
|
|
|
|
if (dataset->indexing_result.size() > image_number
|
|
&& dataset->indexing_result[image_number] != 0
|
|
&& (master_file->Exists("/entry/MX/latticeIndexed") ||
|
|
source_file->Exists("/entry/MX/latticeIndexed"))) {
|
|
std::vector<float> tmp = ReadVectorMasterFirst<float>(
|
|
*master_file,
|
|
*source_file,
|
|
"/entry/MX/latticeIndexed",
|
|
{master_image, 0},
|
|
{source_image, 0},
|
|
{1, 9}
|
|
);
|
|
|
|
if (tmp.size() == 9)
|
|
message.indexing_lattice = CrystalLattice(tmp);
|
|
|
|
std::optional<uint32_t> niggli_opt;
|
|
if (master_file->Exists("/entry/MX/niggli_class"))
|
|
niggli_opt = master_file->ReadElement<uint32_t>("/entry/MX/niggli_class", image_number);
|
|
else if (master_file->Exists("/entry/MX/niggliClass"))
|
|
niggli_opt = master_file->ReadElement<uint32_t>("/entry/MX/niggliClass", image_number);
|
|
else if (source_file->Exists("/entry/MX/niggli_class"))
|
|
niggli_opt = source_file->ReadElement<uint32_t>("/entry/MX/niggli_class", image_id);
|
|
else if (source_file->Exists("/entry/MX/niggliClass"))
|
|
niggli_opt = source_file->ReadElement<uint32_t>("/entry/MX/niggliClass", image_id);
|
|
|
|
if (niggli_opt) {
|
|
auto symm_info = parse_niggli_class(niggli_opt.value());
|
|
|
|
message.lattice_type = LatticeMessage{
|
|
.centering = symm_info.second,
|
|
.niggli_class = static_cast<int64_t>(niggli_opt.value()),
|
|
.crystal_system = symm_info.first,
|
|
};
|
|
}
|
|
}
|
|
|
|
const std::string master_reflection_group_name = fmt::format("/entry/reflections/image_{:06d}", image_number);
|
|
const std::string source_reflection_group_name = fmt::format("/entry/reflections/image_{:06d}", image_id);
|
|
|
|
if (!ReadReflectionsFromGroup(*master_file, master_reflection_group_name, message))
|
|
ReadReflectionsFromGroup(*source_file, source_reflection_group_name, message);
|
|
|
|
return true;
|
|
}
|
|
|
|
void JFJochHDF5Reader::Close() {
|
|
std::unique_lock ul(hdf5_mutex);
|
|
master_file = {};
|
|
number_of_images = 0;
|
|
legacy_format_files.clear();
|
|
vds_data_mappings.clear();
|
|
master_file_directory.clear();
|
|
master_filename.clear();
|
|
SetStartMessage({});
|
|
}
|
|
|
|
std::optional<GoniometerAxis> JFJochHDF5Reader::ReadAxis(HDF5Object *file, const std::string &name) {
|
|
std::string dname = "/entry/sample/transformations/" + name;
|
|
|
|
if (!file->Exists(dname))
|
|
return {};
|
|
|
|
|
|
HDF5DataSet dataset(*file, dname);
|
|
std::vector<double> angle;
|
|
dataset.ReadVector(angle);
|
|
|
|
if (angle.size() < 2)
|
|
return {};
|
|
|
|
std::vector<double> end = file->ReadOptVector<double>(dname + "_end");
|
|
|
|
double start = angle[0];
|
|
double incr = angle[1] - angle[0];
|
|
|
|
if (dataset.ReadAttrStr("transformation_type") != "rotation")
|
|
return {};
|
|
|
|
std::vector<double> axis_vec = dataset.ReadAttrVec("vector");
|
|
if (axis_vec.size() != 3)
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
dname + " Vector must have 3 elements");
|
|
|
|
Coord axis(axis_vec[0], axis_vec[1], axis_vec[2]);
|
|
GoniometerAxis g_axis(name, start, incr, axis, {});
|
|
if (!end.empty())
|
|
g_axis.ScreeningWedge(end[0] - angle[0]);
|
|
|
|
return g_axis;
|
|
}
|
|
|
|
CompressedImage JFJochHDF5Reader::ReadCalibration(std::vector<uint8_t> &tmp, const std::string &name) const {
|
|
std::vector<hsize_t> start = {0, 0};
|
|
if (!master_file)
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Master file not loaded");
|
|
if (!master_file->Exists("/entry/instrument/detector/calibration/" + name))
|
|
throw JFJochException(JFJochExceptionCategory::HDF5, "Calibration dataset not found");
|
|
|
|
HDF5DataSet dataset(*master_file, "/entry/instrument/detector/calibration/" + name);
|
|
HDF5DataSpace dataspace(dataset);
|
|
HDF5DataType datatype(dataset);
|
|
HDF5Dcpl dcpl(dataset);
|
|
|
|
if (dataspace.GetNumOfDimensions() != 2)
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Calibration dataset must be 2D");
|
|
|
|
auto dim = dataspace.GetDimensions();
|
|
|
|
CompressionAlgorithm algorithm = CompressionAlgorithm::NO_COMPRESSION;
|
|
dataset.ReadVectorToU8(tmp, start, {dim[0], dim[1]});
|
|
algorithm = CompressionAlgorithm::NO_COMPRESSION;
|
|
|
|
return {
|
|
tmp, dim[1], dim[0],
|
|
CalcImageMode(datatype.GetElemSize(), datatype.IsFloat(), datatype.IsSigned()),
|
|
algorithm
|
|
};
|
|
}
|
|
|
|
void AppendOrExtendSourceMapping(std::vector<HDF5DataSourceMessage> &ret,
|
|
const std::string &filename,
|
|
const std::string &dataset,
|
|
uint64_t source_first_image,
|
|
uint64_t virtual_first_image,
|
|
uint64_t image_count) {
|
|
if (image_count == 0)
|
|
return;
|
|
|
|
if (!ret.empty()) {
|
|
auto &last = ret.back();
|
|
if (last.filename == filename
|
|
&& last.dataset == dataset
|
|
&& last.source_first_image + last.image_count == source_first_image
|
|
&& last.virtual_first_image + last.image_count == virtual_first_image) {
|
|
last.image_count += image_count;
|
|
return;
|
|
}
|
|
}
|
|
|
|
ret.push_back(HDF5DataSourceMessage{
|
|
.filename = filename,
|
|
.dataset = dataset,
|
|
.source_first_image = source_first_image,
|
|
.virtual_first_image = virtual_first_image,
|
|
.image_count = image_count
|
|
});
|
|
};
|
|
|
|
std::vector<HDF5DataSourceMessage> JFJochHDF5Reader::GetHDF5DataSource(
|
|
uint64_t first_image,
|
|
std::optional<uint64_t> image_count
|
|
) const {
|
|
std::unique_lock ul(hdf5_mutex);
|
|
|
|
if (!master_file)
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Cannot generate HDF5 source mapping if file not loaded");
|
|
|
|
if (first_image > number_of_images)
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"First image outside dataset range");
|
|
|
|
const uint64_t requested_count = image_count.value_or(number_of_images - first_image);
|
|
if (first_image + requested_count > number_of_images)
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Requested image range outside dataset range");
|
|
|
|
std::vector<HDF5DataSourceMessage> ret;
|
|
if (requested_count == 0)
|
|
return ret;
|
|
|
|
// Integrated / contiguous source: link directly to original master file.
|
|
if (format == FileWriterFormat::NXmxVDS && data_layout != HDF5DataSetLayout::VIRTUAL) {
|
|
AppendOrExtendSourceMapping(ret,
|
|
master_filename,
|
|
"/entry/data/data",
|
|
first_image,
|
|
0,
|
|
requested_count);
|
|
return ret;
|
|
}
|
|
|
|
// VDS source: expand VDS mappings to original source files, not to the VDS master.
|
|
if (format == FileWriterFormat::NXmxVDS && data_layout == HDF5DataSetLayout::VIRTUAL) {
|
|
for (uint64_t local_image = 0; local_image < requested_count; ++local_image) {
|
|
const hsize_t virtual_image = first_image + local_image;
|
|
|
|
bool found = false;
|
|
for (const auto &mapping: vds_data_mappings) {
|
|
if (!mapping.ContainsVirtualImage(virtual_image))
|
|
continue;
|
|
|
|
const uint64_t source_image = mapping.SourceImage(virtual_image);
|
|
const std::string filename = ResolveRelativeToMaster(master_file_directory, mapping.filename);
|
|
const std::string dataset = mapping.dataset.empty() ? "/entry/data/data" : mapping.dataset;
|
|
|
|
AppendOrExtendSourceMapping(ret,
|
|
filename,
|
|
dataset,
|
|
source_image,
|
|
local_image,
|
|
1);
|
|
found = true;
|
|
break;
|
|
}
|
|
|
|
if (!found)
|
|
throw JFJochException(JFJochExceptionCategory::HDF5,
|
|
"Image not covered by /entry/data/data VDS mappings");
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
// Legacy source: link directly to linked data files.
|
|
if (format == FileWriterFormat::NXmxLegacy) {
|
|
if (images_per_file == 0)
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Cannot generate HDF5 source mapping: images_per_file is zero");
|
|
|
|
for (uint64_t local_image = 0; local_image < requested_count; ++local_image) {
|
|
const uint64_t source_global_image = first_image + local_image;
|
|
const uint64_t file_id = source_global_image / images_per_file;
|
|
const uint64_t source_image = source_global_image % images_per_file;
|
|
|
|
if (file_id >= legacy_format_files.size())
|
|
throw JFJochException(JFJochExceptionCategory::HDF5,
|
|
"Legacy image source file missing");
|
|
|
|
AppendOrExtendSourceMapping(ret,
|
|
legacy_format_files.at(file_id),
|
|
"/entry/data/data",
|
|
source_image,
|
|
local_image,
|
|
1);
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Unsupported HDF5 file layout for source mapping");
|
|
} |