366 lines
15 KiB
C++
366 lines
15 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "JFJochHDF5Reader.h"
|
|
#include "../include/spdlog/fmt/fmt.h"
|
|
|
|
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;
|
|
}
|
|
|
|
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 < nimages; i++)
|
|
v[image0 + i] = tmp[i];
|
|
}
|
|
} catch (JFJochException &e) {}
|
|
}
|
|
|
|
std::string removeSuffix(const std::string& s, const std::string& suffix)
|
|
{
|
|
if (s.rfind(suffix) == s.size() - suffix.size()) {
|
|
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;
|
|
}
|
|
|
|
void JFJochHDF5Reader::ReadFile(const std::string& filename) {
|
|
std::unique_lock ul(master_file_mutex);
|
|
try {
|
|
auto dataset = std::make_shared<JFJochReaderDataset>();
|
|
master_file = std::make_unique<HDF5ReadOnlyFile>(filename);
|
|
|
|
std::filesystem::path fsPath(filename);
|
|
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")) {
|
|
legacy_format = false;
|
|
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;
|
|
|
|
dataset->efficiency = master_file->ReadVector<float>(
|
|
"/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image");
|
|
|
|
dataset->spot_count = master_file->ReadOptVector<uint32_t>("/entry/MX/nPeaks");
|
|
dataset->indexing_result = master_file->ReadOptVector<uint8_t>("/entry/MX/imageIndexed");
|
|
dataset->bkg_estimate = master_file->ReadOptVector<float>("/entry/MX/bkgEstimate");
|
|
dataset->resolution_estimate = master_file->ReadOptVector<float>("/entry/MX/resolutionEstimate");
|
|
} else if (master_file->Exists("/entry/data/data_000001")) {
|
|
legacy_format = true;
|
|
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);
|
|
|
|
ReadVector(dataset->efficiency,
|
|
data_file, "/entry/detector/data_collection_efficiency_image",
|
|
number_of_images, fimages);
|
|
|
|
ReadVector(dataset->spot_count,
|
|
data_file, "/entry/MX/nPeaks",
|
|
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);
|
|
} 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;
|
|
}
|
|
|
|
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/jfjoch_release");
|
|
|
|
InstrumentMetadata metadata;
|
|
metadata.InstrumentName(master_file->GetString("/entry/instrument/name"));
|
|
metadata.SourceName(master_file->GetString("/entry/source/name"));
|
|
dataset->experiment.ImportInstrumentMetadata(metadata);
|
|
|
|
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").value_or(0));
|
|
dataset->experiment.SampleName(master_file->GetString("/entry/sample/name"));
|
|
|
|
|
|
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"));
|
|
}
|
|
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/image/max_value"))
|
|
dataset->max_value = master_file->ReadVector<int64_t>("/entry/image/max_value");
|
|
|
|
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");
|
|
bin_to_q_dataset.ReadVector(dataset->az_int_bin_to_q);
|
|
}
|
|
|
|
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::microseconds (0));
|
|
dataset->experiment.Detector(detector);
|
|
|
|
dataset->experiment.FrameTime(
|
|
std::chrono::microseconds(std::lround(master_file->GetFloat("/entry/instrument/detector/frame_time") * 1e6)),
|
|
std::chrono::microseconds(std::lround(master_file->GetFloat("/entry/instrument/detector/count_time") * 1e6)));
|
|
|
|
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(master_file_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};
|
|
|
|
HDF5DataSet dataset(file, "/entry/data/data");
|
|
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[1], dim[2],
|
|
CalcImageMode(datatype.GetElemSize(), datatype.IsFloat(), datatype.IsSigned()),
|
|
algorithm};
|
|
}
|
|
|
|
std::shared_ptr<JFJochReaderImage> JFJochHDF5Reader::LoadImageInternal(int64_t image_number) {
|
|
std::unique_lock ul(master_file_mutex);
|
|
|
|
if (!master_file)
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Cannot load image if file not loaded");
|
|
|
|
if (image_number >= number_of_images)
|
|
throw JFJochException(JFJochExceptionCategory::HDF5, "Image out of bounds");
|
|
|
|
std::unique_ptr<HDF5ReadOnlyFile> tmp_data_file;
|
|
uint32_t image_id;
|
|
HDF5Object *source_file;
|
|
|
|
if (legacy_format) {
|
|
uint32_t file_id = image_number / images_per_file;
|
|
image_id = image_number % images_per_file;
|
|
tmp_data_file = std::make_unique<HDF5ReadOnlyFile>(legacy_format_files.at(file_id));
|
|
source_file = tmp_data_file.get();
|
|
} else {
|
|
image_id = image_number;
|
|
source_file = master_file.get();
|
|
}
|
|
|
|
DataMessage message;
|
|
message.image = LoadImageDataset(buffer, *source_file, image_id);
|
|
message.number = image_number;
|
|
auto dataset = GetStartMessage();
|
|
|
|
if (dataset->indexing_result.size() > image_number && dataset->indexing_result[image_number] != 0
|
|
&& source_file->Exists("/entry/MX") && source_file->Exists("/entry/MX/latticeIndexed")) {
|
|
|
|
std::vector<float> tmp = source_file->ReadVector<float>(
|
|
"/entry/MX/latticeIndexed",
|
|
{(hsize_t) image_id, 0},
|
|
{1, 9}
|
|
);
|
|
message.indexing_lattice = CrystalLattice(tmp);
|
|
}
|
|
|
|
if ((dataset->spot_count.size() > image_number)
|
|
&& (dataset->spot_count[image_number] > 0)) {
|
|
auto spot_count = dataset->spot_count[image_number];
|
|
auto spot_x = source_file->ReadVector<float>(
|
|
"/entry/MX/peakXPosRaw",
|
|
{(hsize_t) image_id, 0},
|
|
{1, spot_count}
|
|
);
|
|
auto spot_y = source_file->ReadVector<float>(
|
|
"/entry/MX/peakYPosRaw",
|
|
{(hsize_t) image_id, 0},
|
|
{1, spot_count}
|
|
);
|
|
auto spot_intensity = source_file->ReadVector<float>(
|
|
"/entry/MX/peakTotalIntensity",
|
|
{(hsize_t) image_id, 0},
|
|
{1, spot_count}
|
|
);
|
|
auto spot_indexed = source_file->ReadVector<uint8_t>(
|
|
"/entry/MX/peakIndexed",
|
|
{(hsize_t) image_id, 0},
|
|
{1, spot_count}
|
|
);
|
|
for (int i = 0; i < spot_count; i++) {
|
|
message.spots.emplace_back(SpotToSave{
|
|
.x = spot_x.at(i),
|
|
.y = spot_y.at(i),
|
|
.intensity = spot_intensity.at(i),
|
|
.indexed = (spot_indexed.at(i) != 0)
|
|
});
|
|
}
|
|
}
|
|
|
|
if (!dataset->az_int_bin_to_q.empty()) {
|
|
message.az_int_profile = source_file->ReadOptVector<float>(
|
|
"/entry/azint/image",
|
|
{(hsize_t) image_id, 0},
|
|
{1, dataset->az_int_bin_to_q.size()}
|
|
);
|
|
}
|
|
|
|
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];
|
|
|
|
return std::make_shared<JFJochReaderImage>(message, dataset);
|
|
}
|
|
|
|
void JFJochHDF5Reader::Close() {
|
|
std::unique_lock ul(master_file_mutex);
|
|
master_file = {};
|
|
number_of_images = 0;
|
|
legacy_format_files.clear();
|
|
SetStartMessage({});
|
|
} |