mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 85 Wiki Activity
Files
rc178
Jungfraujoch/reader/JFJochXFELReader.cpp
T
leonarski_fandClaude Opus 5.5 3d8728e3a1 JFJochXFELReader: the positional-read handles are held for the most recent files only 
The direct chunk read opened one HDF5ImageSource::RawFile per file and never released it, so the
MAX_OPEN_FILES bound on the HDF5 handles did not hold for that path: a dataset stored one frame per
file would keep a descriptor per file open. They now follow the same least-recently-used bound; a
reading thread still holding a released handle keeps it through its shared_ptr.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
2026-10-21 01:52:55 +02:00

865 lines
41 KiB
C++
Raw Permalink Blame History

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "JFJochXFELReader.h"
#include "CrystFELGeometry.h"
#include <algorithm>
#include <climits>
#include <cmath>
#include <cstring>
#include <filesystem>
#include <numeric>
#include <zlib.h>
#include "../common/JFJochException.h"
#include "../common/Logger.h"
#include "../writer/HDF5Objects.h"
namespace {
constexpr double JOULE_PER_EV = 1.602176634e-19;
constexpr double PI_D = 3.14159265358979323846;
const std::string CXI_DATA = "/entry_1/data_1/data";
const std::string CXI_DETECTOR = "/entry_1/instrument_1/detector_1/";
// The bits of a Cheetah pixel mask that mark a pixel unusable in every frame: invalid, hot, dead,
// shadowed, bad, missing and noisy (Cheetah's PIXEL_IS_* flags). The peak mask, saturation and the
// resolution limits are left alone.
constexpr uint32_t CHEETAH_BAD_BITS = 1u | 4u | 8u | 16u | 128u | 512u | 1024u;
std::string ReplaceEvent(std::string path, const std::string &group) {
const auto pos = path.find('%');
if (pos != std::string::npos)
path.replace(pos, 1, group);
return path;
}
// A running mean, and the spread, of a value the geometry takes from the image files.
struct Mean {
double sum = 0, min = INFINITY, max = -INFINITY;
size_t n = 0;
void Add(double v) { sum += v; n++; min = std::min(min, v); max = std::max(max, v); }
double Get() const { return sum / static_cast<double>(n); }
};
// Adds what one file states at path: one value per shot group (a path with %), a per-frame vector, or
// one scalar for the file.
void Accumulate(HDF5ReadOnlyFile &file, const std::string &path, const std::vector<std::string> &groups,
Mean &mean) {
if (path.find('%') != std::string::npos) {
for (const auto &g: groups)
mean.Add(HDF5DataSet(file, ReplaceEvent(path, g)).ReadScalar<double>());
} else if (file.GetDimension(path).empty())
mean.Add(HDF5DataSet(file, path).ReadScalar<double>());
else
for (const double v: file.ReadVector<double>(path))
mean.Add(v);
}
} // namespace
JFJochXFELReader::JFJochXFELReader() = default;
JFJochXFELReader::~JFJochXFELReader() {
Close();
}
// Whether the file is HDF5 at all, asked without HDF5 printing its error stack for one that is not (an
// XTC stream offered to every reader in turn). Under hdf5_mutex.
static bool IsHDF5(const std::string &path) {
htri_t r = -1;
H5E_BEGIN_TRY {
r = H5Fis_accessible(path.c_str(), H5P_DEFAULT);
} H5E_END_TRY
return r > 0;
}
// The HDF5 filter IDs registered for bitshuffle and for plain LZ4.
constexpr H5Z_filter_t BITSHUFFLE_FILTER_ID = 32008;
constexpr H5Z_filter_t LZ4_FILTER_ID = 32004;
std::shared_ptr<const HDF5ImageSource::RawFile> JFJochXFELReader::RawFileAt(size_t file) const {
if (raw_files_.size() != files_.size()) {
raw_files_.assign(files_.size(), nullptr);
user_block_.assign(files_.size(), 0);
raw_order_.clear();
}
if (!raw_files_[file]) {
const hid_t fcpl = H5Fget_create_plist(FileAt(file).GetID());
hsize_t user_block = 0;
if (fcpl >= 0) {
H5Pget_userblock(fcpl, &user_block);
H5Pclose(fcpl);
}
auto raw = std::make_shared<const HDF5ImageSource::RawFile>(files_[file]);
if (!raw->IsOpen())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Cannot open " + files_[file]);
raw_files_[file] = raw;
user_block_[file] = user_block;
raw_order_.push_back(file);
if (raw_order_.size() > MAX_OPEN_FILES) {
raw_files_[raw_order_.front()].reset();
raw_order_.erase(raw_order_.begin());
}
}
return raw_files_[file];
}
HDF5ReadOnlyFile &JFJochXFELReader::FileAt(size_t file) const {
auto it = std::find_if(open_files_.begin(), open_files_.end(), [&](const auto &e) { return e.first == file; });
if (it == open_files_.end()) {
if (open_files_.size() >= MAX_OPEN_FILES)
open_files_.erase(open_files_.begin());
open_files_.emplace_back(file, std::make_unique<HDF5ReadOnlyFile>(files_[file]));
} else if (it != open_files_.end() - 1)
std::rotate(it, it + 1, open_files_.end());
return *open_files_.back().second;
}
// The files given, in the order given; a directory is not read (the caller expands a glob).
static void CheckFiles(const std::vector<std::string> &files) {
if (files.empty())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No file to read");
}
void JFJochXFELReader::ReadFiles(const std::vector<std::string> &files, const std::string &geometry_file,
std::optional<float> gain) {
Close();
CheckFiles(files);
std::vector<std::string> ignored;
geom_ = crystfel::Read(geometry_file, ignored);
source_ = geometry_file;
if (!ignored.empty()) {
std::string list;
for (const auto &key: ignored)
list += (list.empty() ? "" : ", ") + key;
Logger("XFELReader").Info("{}: only the module geometry (and camera length) is taken; ignored: {}",
geometry_file, list);
}
// What the geometry file no longer states comes from the data files themselves, where they are one of
// the layouts this reader knows: a Cheetah CXI file (photon energy, camera length, the shared mask
// and the per-shot one) or a SACLA Cheetah file (photon energy per shot, camera length, pixel mask).
// Otherwise the photon energy is the command line's to give, and the camera length too where the
// geometry has none.
{
std::unique_lock ul(hdf5_mutex);
HDF5ReadOnlyFile file(files[0]);
if (file.Exists(CXI_DATA)) {
if (file.Exists("/entry_1/instrument_1/source_1/energy"))
geom_.photon_energy = xfel::Value{"/entry_1/instrument_1/source_1/energy", 0, 1.0 / JOULE_PER_EV};
if (file.Exists(CXI_DETECTOR + "mask_shared"))
geom_.masks.push_back({CXI_DETECTOR + "mask_shared", "", 0, CHEETAH_BAD_BITS});
if (file.Exists("/entry_1/data_1/mask"))
geom_.masks.push_back({"/entry_1/data_1/mask", "", 0, CHEETAH_BAD_BITS});
if (!geom_.panels[0].clen && file.Exists(CXI_DETECTOR + "distance"))
for (auto &p: geom_.panels)
p.clen = xfel::Value{CXI_DETECTOR + "distance", 0, 1.0};
} else if (file.Exists("/metadata/distance_in_mm")) {
geom_.photon_energy = xfel::Value{"/%/photon_energy_ev", 0, 1.0};
geom_.masks.push_back({"/metadata/pixelmask", "", 0, 0xffffffffu});
if (!geom_.panels[0].clen)
for (auto &p: geom_.panels)
p.clen = xfel::Value{"/metadata/distance_in_mm", 0, 1e-3};
}
}
for (auto &p: geom_.panels) {
if (!p.clen)
p.clen = xfel::Value{"", 0.1, 1.0}; // a placeholder: the distance is given on the command line
if (gain)
p.adu_per_photon = *gain;
}
Open(files);
Place();
}
// A Cheetah CXI file carries what one flat panel needs - the frames at /entry_1/data_1/data, the pixel
// size, the photon energy and a camera length, all in SI units - but not where the beam is. The panel is
// put with its middle on the beam, which is as good as no centre at all: the run measures it.
//
// Only a square frame is taken this way. A Cheetah file from a tiled detector (CSPAD, AGIPD, ePix,
// JUNGFRAU, MPCCD) holds its panels stacked in a slab, which read as one flat panel would place every
// pixel wrongly and say nothing; those slabs are not square, the single-panel CCDs are.
bool JFJochXFELReader::CanReadCXI(const std::string &path) {
try {
if (!std::filesystem::is_regular_file(path))
return false;
std::unique_lock ul(hdf5_mutex);
if (!IsHDF5(path))
return false;
HDF5ReadOnlyFile file(path);
if (!file.Exists(CXI_DATA))
return false;
const auto dim = file.GetDimension(CXI_DATA);
return dim.size() == 3 && dim[1] == dim[2]
&& file.Exists(CXI_DETECTOR + "x_pixel_size")
&& file.Exists("/entry_1/instrument_1/source_1/energy");
} catch (...) {
return false;
}
}
std::optional<std::string> JFJochXFELReader::CXIRefusal(const std::string &path) {
try {
if (!std::filesystem::is_regular_file(path) || CanReadCXI(path))
return std::nullopt;
std::unique_lock ul(hdf5_mutex);
if (!IsHDF5(path))
return std::nullopt;
HDF5ReadOnlyFile file(path);
if (!file.Exists(CXI_DATA))
return std::nullopt;
const auto dim = file.GetDimension(CXI_DATA);
if (dim.size() == 3 && dim[1] != dim[2])
return fmt::format("{} is a Cheetah CXI file whose {} x {} frame is the panels of a tiled detector "
"stacked in a slab; their layout is not in the file - give it with "
"--crystfel-geometry", path, dim[2], dim[1]);
return fmt::format("{} is a CXI file without the pixel size and photon energy of a single-panel "
"detector; give its geometry with --crystfel-geometry", path);
} catch (...) {
return std::nullopt;
}
}
void JFJochXFELReader::ReadCXI(const std::vector<std::string> &files, std::optional<float> gain) {
Close();
CheckFiles(files);
std::vector<hsize_t> dim;
bool has_mask = false;
{
std::unique_lock ul(hdf5_mutex);
HDF5ReadOnlyFile file(files[0]);
dim = file.GetDimension(CXI_DATA);
has_mask = file.Exists(CXI_DETECTOR + "mask_shared");
const auto pixel = file.ReadVector<double>(CXI_DETECTOR + "x_pixel_size");
xfel::Panel p;
p.name = "cxi";
p.max_ss = static_cast<int64_t>(dim[1]) - 1;
p.max_fs = static_cast<int64_t>(dim[2]) - 1;
p.corner_x = -0.5 * static_cast<double>(dim[2]);
p.corner_y = -0.5 * static_cast<double>(dim[1]);
p.res = 1.0 / pixel.at(0);
if (file.Exists(CXI_DETECTOR + "distance"))
p.clen = xfel::Value{CXI_DETECTOR + "distance", 0, 1.0};
else
p.clen = xfel::Value{"", 0.1, 1.0}; // a placeholder: the distance is given or found
if (gain)
p.adu_per_photon = *gain;
geom_ = {};
geom_.panels.push_back(p);
geom_.data = CXI_DATA;
geom_.photon_energy = xfel::Value{"/entry_1/instrument_1/source_1/energy", 0, 1.0 / JOULE_PER_EV};
if (has_mask)
geom_.masks.push_back({CXI_DETECTOR + "mask_shared", "", 0, CHEETAH_BAD_BITS});
}
source_ = "the CXI metadata";
Open(files);
Place();
}
// A SACLA Cheetah file: one group per shot (tag-N) with the frame and the shot's photon energy, and
// /metadata stating where every sensor of the MPCCD lies (posx/posy in um, the in-plane angle - in
// degrees, whatever its name says - and the pixel size), the camera length and a pixel mask that
// already holds the sensors' edges and the rows behind each amplifier's shield. The sensors are stacked
// one after the other along the frame's slow axis.
bool JFJochXFELReader::CanReadPanelTable(const std::string &path) {
try {
if (!std::filesystem::is_regular_file(path))
return false;
std::unique_lock ul(hdf5_mutex);
if (!IsHDF5(path))
return false;
HDF5ReadOnlyFile file(path);
return file.Exists("/metadata/posx_in_um") && file.Exists("/metadata/posy_in_um")
&& file.Exists("/metadata/angle_in_rad") && file.Exists("/metadata/pixelsizex_in_um")
&& file.Exists("/metadata/distance_in_mm");
} catch (...) {
return false;
}
}
void JFJochXFELReader::ReadPanelTable(const std::vector<std::string> &files, std::optional<float> gain) {
Close();
CheckFiles(files);
struct Metrology {
std::vector<double> posx, posy, angle, pixel;
double distance_mm = 0;
bool operator==(const Metrology &) const = default;
};
auto read_metrology = [](HDF5ReadOnlyFile &file) {
Metrology m;
m.posx = file.ReadVector<double>("/metadata/posx_in_um");
m.posy = file.ReadVector<double>("/metadata/posy_in_um");
m.angle = file.ReadVector<double>("/metadata/angle_in_rad");
m.pixel = file.ReadVector<double>("/metadata/pixelsizex_in_um");
m.distance_mm = HDF5DataSet(file, "/metadata/distance_in_mm").ReadScalar<double>();
return m;
};
Metrology m;
std::string first_group;
{
std::unique_lock ul(hdf5_mutex);
HDF5ReadOnlyFile file(files[0]);
m = read_metrology(file);
// The sensors of one run are placed once: every file has to state the same metrology.
for (size_t f = 1; f < files.size(); f++) {
HDF5ReadOnlyFile other(files[f]);
if (!(read_metrology(other) == m))
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
files[f] + " places the sensors or the detector differently from "
+ files[0] + "; read the runs separately");
}
for (const auto &g: file.FindLeafs("/"))
if (file.Exists("/" + g + "/data")) {
first_group = g;
break;
}
if (first_group.empty())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No shot groups in " + files[0]);
}
const size_t nsensor = m.posx.size();
if (nsensor == 0 || m.posy.size() != nsensor || m.angle.size() != nsensor || m.pixel.size() != nsensor)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
files[0] + ": /metadata does not state the same number of positions, angles and "
"pixel sizes");
if (!(m.distance_mm > 0))
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, files[0] + " states no camera length");
std::vector<hsize_t> dim;
{
std::unique_lock ul(hdf5_mutex);
HDF5ReadOnlyFile file(files[0]);
dim = file.GetDimension("/" + first_group + "/data");
}
if (dim.size() != 2 || dim[0] % nsensor != 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
files[0] + ": the frame does not hold " + std::to_string(nsensor) + " sensors stacked");
const auto rows = static_cast<int64_t>(dim[0] / nsensor);
geom_ = {};
geom_.data = "/%/data";
geom_.photon_energy = xfel::Value{"/%/photon_energy_ev", 0, 1.0};
geom_.masks.push_back({"/metadata/pixelmask", "", 0, 0xffffffffu});
// The placement the SACLA pipeline itself writes into its CrystFEL geometry: the corner of the
// sensor's first pixel at (-posx, posy) and its fast axis at the stated angle.
for (size_t i = 0; i < nsensor; i++) {
const double a = m.angle[i] * PI_D / 180.0;
xfel::Panel p;
p.name = "q" + std::to_string(i + 1);
p.min_ss = static_cast<int64_t>(i) * rows;
p.max_ss = p.min_ss + rows - 1;
p.max_fs = static_cast<int64_t>(dim[1]) - 1;
p.corner_x = -m.posx[i] / m.pixel[i];
p.corner_y = m.posy[i] / m.pixel[i];
p.fs_x = std::cos(a);
p.fs_y = std::sin(a);
p.ss_x = -std::sin(a);
p.ss_y = std::cos(a);
p.res = 1e6 / m.pixel[i];
p.clen = xfel::Value{"", m.distance_mm * 1e-3, 1.0};
if (gain)
p.adu_per_photon = *gain;
geom_.panels.push_back(p);
}
source_ = "the /metadata panel table";
Open(files);
Place();
}
JFJochXFELReader::Layout JFJochXFELReader::ReadLayout(HDF5ReadOnlyFile &file, const Frame &frame) {
Layout l;
HDF5DataSet dataset(file, frame.dataset);
const auto dim = file.GetDimension(frame.dataset);
if (dim.size() < 2 || dim.size() > 3)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
frame.dataset + " is not a stack of frames or one frame");
l.height = dim[dim.size() - 2];
l.width = dim[dim.size() - 1];
const hid_t dcpl = H5Dget_create_plist(dataset.GetID());
const int rank = static_cast<int>(dim.size());
std::vector<hsize_t> chunk(rank, 0);
bool ok = H5Pget_layout(dcpl) == H5D_CHUNKED && H5Pget_chunk(dcpl, rank, chunk.data()) == rank
&& chunk[rank - 1] == l.width && chunk[rank - 2] == l.height && (rank == 2 || chunk[0] == 1);
bool deflate = false, bitshuffle = false;
for (int i = 0; ok && i < H5Pget_nfilters(dcpl); i++) {
unsigned int flags;
size_t nelements = 0;
const H5Z_filter_t filter = H5Pget_filter2(dcpl, i, &flags, &nelements, nullptr, 0, nullptr, nullptr);
if (filter == H5Z_FILTER_DEFLATE)
deflate = true;
else if (filter == H5Z_FILTER_SHUFFLE)
l.shuffled = true;
else if (filter == BITSHUFFLE_FILTER_ID || filter == LZ4_FILTER_ID)
bitshuffle = true;
else
ok = false;
}
H5Pclose(dcpl);
// Bitshuffle/LZ4 (or zstd, or plain LZ4) is the one filter of its chunk: CompressedImage decodes it.
if (ok && bitshuffle && !deflate && !l.shuffled)
l.algorithm = HDF5Dcpl(dataset).GetCompression();
const hid_t type = H5Dget_type(dataset.GetID());
l.stored_size = H5Tget_size(type);
l.stored_float = H5Tget_class(type) == H5T_FLOAT;
l.stored_signed = H5Tget_sign(type) == H5T_SGN_2;
ok &= ((H5Tget_class(type) == H5T_INTEGER && (l.stored_size == 2 || l.stored_size == 4))
|| (l.stored_float && l.stored_size == 4))
&& H5Tget_order(type) == H5T_ORDER_LE;
H5Tclose(type);
l.deflated = deflate;
l.direct_chunks = ok && (deflate || (bitshuffle && l.algorithm != CompressionAlgorithm::NO_COMPRESSION));
return l;
}
void JFJochXFELReader::Open(const std::vector<std::string> &files) {
files_ = files;
geom_.masks.insert(geom_.masks.end(), extra_masks_.begin(), extra_masks_.end());
// Every frame of every file, in file order, with the mean of every value taken from the files.
Mean energy, clen;
const auto &p0 = geom_.panels.at(0);
const auto pct = geom_.data.find('%');
for (const auto &p: geom_.panels)
if (p.res != p0.res || p.clen->path != p0.clen->path || p.clen->value != p0.clen->value)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
source_ + ": panels with different camera lengths or pixel sizes cannot be put "
"on one flat detector");
const xfel::Value *energy_value = geom_.photon_energy ? &*geom_.photon_energy
: (geom_.wavelength ? &*geom_.wavelength : nullptr);
std::unique_lock ul(hdf5_mutex);
for (size_t f = 0; f < files_.size(); f++) {
HDF5ReadOnlyFile file(files_[f]);
const size_t first = frames_.size();
std::vector<std::string> groups;
if (pct != std::string::npos) {
const std::string parent = geom_.data.substr(0, pct);
for (const auto &group: file.FindLeafs(parent.empty() ? "/" : parent)) {
const auto dataset = ReplaceEvent(geom_.data, group);
if (file.Exists(dataset)) {
frames_.push_back({f, dataset, -1, group});
groups.push_back(group);
}
}
} else if (file.Exists(geom_.data)) {
const auto dim = file.GetDimension(geom_.data);
if (dim.size() == 3) {
for (hsize_t i = 0; i < dim[0]; i++)
frames_.push_back({f, geom_.data, static_cast<int64_t>(i), ""});
} else if (dim.size() == 2)
frames_.push_back({f, geom_.data, -1, ""});
}
if (frames_.size() == first)
continue;
// Every file has to store its frames the way the first does: the decoding is set up once.
const Layout l = ReadLayout(file, frames_[first]);
if (first == 0)
layout_ = l;
else if (!(l == layout_))
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
files_[f] + " stores its frames differently (size, type or compression) from "
+ files_[0] + "; read them separately");
if (energy_value != nullptr && !energy_value->path.empty())
Accumulate(file, energy_value->path, groups, energy);
if (!p0.clen->path.empty())
Accumulate(file, p0.clen->path, groups, clen);
}
if (frames_.empty())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"No frames at " + geom_.data + " in the files given");
frame_height_ = layout_.height;
frame_width_ = layout_.width;
clen_m_ = p0.clen->path.empty() ? p0.clen->value : clen.Get() * p0.clen->scale;
double coffset_sum = 0;
for (const auto &p: geom_.panels)
coffset_sum += p.coffset;
distance_mm_ = (clen_m_ + coffset_sum / static_cast<double>(geom_.panels.size())) * 1000.0;
// Nothing states the photon energy: it is the command line's to give (--wavelength).
photon_energy_eV_ = 0;
if (energy_value != nullptr) {
const double v = energy_value->path.empty() ? energy_value->value : energy.Get() * energy_value->scale;
photon_energy_eV_ = geom_.photon_energy ? v : WVL_1A_IN_KEV * 1000.0 / v;
if (!(photon_energy_eV_ > 0) || !std::isfinite(photon_energy_eV_))
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
source_ + ": no usable photon energy in the files given");
}
if (energy.n > 1)
Logger("XFELReader").Info("Photon energy {:.1f} eV, the mean of {} values from {:.1f} to {:.1f} (the run "
"is read at one wavelength)", photon_energy_eV_, energy.n,
geom_.photon_energy ? energy.min * energy_value->scale
: WVL_1A_IN_KEV * 1000.0 / (energy.max * energy_value->scale),
geom_.photon_energy ? energy.max * energy_value->scale
: WVL_1A_IN_KEV * 1000.0 / (energy.min * energy_value->scale));
}
void JFJochXFELReader::Place() {
std::vector<uint32_t> file_mask(frame_width_ * frame_height_, 0);
std::vector<uint8_t> masked_by_file(frame_width_ * frame_height_, 0);
{
std::unique_lock ul(hdf5_mutex);
frame_masks_.clear();
for (const auto &m: geom_.masks) {
std::vector<uint32_t> v;
if (m.file.empty() && FileAt(0).GetDimension(m.data).size() == 3) {
frame_masks_.push_back(m); // one plane per frame: applied as each frame is read
continue;
}
if (!m.file.empty()) {
HDF5ReadOnlyFile mask_file(m.file);
v = mask_file.ReadVector<uint32_t>(m.data, {0, 0}, {frame_height_, frame_width_});
} else
v = FileAt(0).ReadVector<uint32_t>(m.data, {0, 0}, {frame_height_, frame_width_});
for (size_t i = 0; i < v.size(); i++)
if (m.Bad(v[i]))
masked_by_file[i] = 1;
}
}
// CrystFEL's laboratory frame and the flat image's (column, row, beam) are both right-handed with
// the beam along z, so the point (X, Y) lies at column beam_x + X, row beam_y + Y of the detector plane
// - and CrystFEL's horizontal x stays the image's horizontal, where the polarisation correction expects
// it. A pixel centre is at its index. A panel set nearer or further than the plane (its own coffset)
// is projected onto it along the ray from the sample, so each pixel keeps its scattering direction.
auto position = [this](const xfel::Panel &p, double fs, double ss) {
const double scale = distance_mm_ / ((clen_m_ + p.coffset) * 1000.0);
return std::pair{scale * (p.corner_x + fs * p.fs_x + ss * p.ss_x),
scale * (p.corner_y + fs * p.fs_y + ss * p.ss_y)};
};
// Each panel goes onto the image whole, as a block: its axes turned to the nearest image axis and its
// first pixel put at a whole pixel. Where it really lies - the fraction of a pixel, the small turn - is
// its module's placement (DetectorModules), so no pixel is moved on its own.
struct Block {
int64_t fx = 1, fy = 0, sx = 0, sy = 1; // image step per +1 fs / +1 ss
double x0 = 0, y0 = 0; // plane position of the first pixel's centre, beam at 0
};
auto snap = [](double x, double y, int64_t &ix, int64_t &iy) {
if (std::fabs(x) >= std::fabs(y)) { ix = x > 0 ? 1 : -1; iy = 0; }
else { ix = 0; iy = y > 0 ? 1 : -1; }
};
std::vector<Block> blocks;
double min_x = 1e9, min_y = 1e9, max_x = -1e9, max_y = -1e9;
for (const auto &p: geom_.panels) {
Block bl;
snap(p.fs_x, p.fs_y, bl.fx, bl.fy);
snap(p.ss_x, p.ss_y, bl.sx, bl.sy);
if (bl.fx * bl.sy - bl.fy * bl.sx == 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Panel " + p.name + " has its fast and slow axes along the same line");
std::tie(bl.x0, bl.y0) = position(p, 0.5, 0.5);
const double w = static_cast<double>(p.max_fs - p.min_fs), h = static_cast<double>(p.max_ss - p.min_ss);
for (double fs: {0.0, w})
for (double ss: {0.0, h}) {
const double x = bl.x0 + fs * static_cast<double>(bl.fx) + ss * static_cast<double>(bl.sx);
const double y = bl.y0 + fs * static_cast<double>(bl.fy) + ss * static_cast<double>(bl.sy);
min_x = std::min(min_x, x); max_x = std::max(max_x, x);
min_y = std::min(min_y, y); max_y = std::max(max_y, y);
}
blocks.push_back(bl);
}
// The first panel's pixel centres on whole columns and rows, two pixels of margin all round.
const double beam_x = std::ceil(blocks[0].x0 - min_x) + 2 - blocks[0].x0;
const double beam_y = std::ceil(blocks[0].y0 - min_y) + 2 - blocks[0].y0;
width_ = static_cast<size_t>(std::ceil(beam_x + max_x)) + 3;
height_ = static_cast<size_t>(std::ceil(beam_y + max_y)) + 3;
// Everything starts as a gap; the panels' pixels are put on top, less the bad regions.
std::vector<uint32_t> mask(width_ * height_, 1u << PixelMask::ModuleGapPixelBit);
std::vector<uint8_t> taken(width_ * height_, 0);
destination_.assign(frame_width_ * frame_height_, -1);
adu_per_photon_.assign(frame_width_ * frame_height_, 1.0f);
counts_in_photons_ = true;
std::optional<double> saturation;
std::vector<DetectorModule> modules;
for (size_t n = 0; n < geom_.panels.size(); n++) {
const auto &p = geom_.panels[n];
const Block &bl = blocks[n];
if (p.max_fs >= static_cast<int64_t>(frame_width_) || p.max_ss >= static_cast<int64_t>(frame_height_))
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Panel " + p.name + " lies outside the stored frame");
double gain = 1.0;
if (p.adu_per_photon.has_value())
gain = p.adu_per_photon.value();
else if (p.adu_per_eV.has_value())
gain = p.adu_per_eV.value() * photon_energy_eV_;
else
counts_in_photons_ = false;
if (p.max_adu.has_value())
saturation = std::min(saturation.value_or(INT64_MAX), p.max_adu.value() / gain);
if (p.no_index)
continue;
// The image pixel the panel's first pixel lands on, and the block it fills.
const int64_t c0 = std::lround(beam_x + bl.x0), r0 = std::lround(beam_y + bl.y0);
const int64_t w = p.max_fs - p.min_fs, h = p.max_ss - p.min_ss;
const int64_t bx0 = c0 + std::min<int64_t>(0, bl.fx * w) + std::min<int64_t>(0, bl.sx * h);
const int64_t by0 = r0 + std::min<int64_t>(0, bl.fy * w) + std::min<int64_t>(0, bl.sy * h);
// Its placement: the plane position of image pixel (bx0, by0), and the plane steps of +1 column
// and +1 row - the panel's own axes, scaled by its projection, taken through the block's turn.
const double scale = distance_mm_ / ((clen_m_ + p.coffset) * 1000.0);
DetectorModule m;
m.x0 = bx0;
m.y0 = by0;
m.width = std::abs(bl.fx * w + bl.sx * h) + 1;
m.height = std::abs(bl.fy * w + bl.sy * h) + 1;
// image (dc, dr) = fs (fx, fy) + ss (sx, sy); the inverse of a signed permutation is its transpose.
auto plane_of = [&](double dc, double dr) {
const double fs = dc * static_cast<double>(bl.fx) + dr * static_cast<double>(bl.fy);
const double ss = dc * static_cast<double>(bl.sx) + dr * static_cast<double>(bl.sy);
const auto [x, y] = position(p, fs + 0.5, ss + 0.5);
return std::pair{beam_x + x, beam_y + y};
};
std::tie(m.ox, m.oy) = plane_of(static_cast<double>(bx0 - c0), static_cast<double>(by0 - r0));
m.fx = scale * (static_cast<double>(bl.fx) * p.fs_x + static_cast<double>(bl.sx) * p.ss_x);
m.fy = scale * (static_cast<double>(bl.fx) * p.fs_y + static_cast<double>(bl.sx) * p.ss_y);
m.sx = scale * (static_cast<double>(bl.fy) * p.fs_x + static_cast<double>(bl.sy) * p.ss_x);
m.sy = scale * (static_cast<double>(bl.fy) * p.fs_y + static_cast<double>(bl.sy) * p.ss_y);
modules.push_back(m);
const int64_t e = p.mask_edge_pixels;
for (int64_t ss = p.min_ss; ss <= p.max_ss; ss++) {
for (int64_t fs = p.min_fs; fs <= p.max_fs; fs++) {
const int64_t col = c0 + (fs - p.min_fs) * bl.fx + (ss - p.min_ss) * bl.sx;
const int64_t row = r0 + (fs - p.min_fs) * bl.fy + (ss - p.min_ss) * bl.sy;
const size_t dst = static_cast<size_t>(row) * width_ + static_cast<size_t>(col);
if (taken[dst])
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Panel " + p.name + " overlaps another on the image");
taken[dst] = 1;
if (fs < p.min_fs + e || fs > p.max_fs - e || ss < p.min_ss + e || ss > p.max_ss - e)
continue;
const auto [x, y] = position(p, fs - p.min_fs + 0.5, ss - p.min_ss + 0.5);
bool bad = false;
for (const auto &b: geom_.bad) {
if (b.in_lab)
bad |= x >= b.min_x && x <= b.max_x && y >= b.min_y && y <= b.max_y;
else if (b.panel.empty() || b.panel == p.name)
bad |= fs >= b.min_fs && fs <= b.max_fs && ss >= b.min_ss && ss <= b.max_ss;
}
if (bad)
continue;
const size_t src = ss * frame_width_ + fs;
destination_[src] = static_cast<int64_t>(dst);
adu_per_photon_[src] = static_cast<float>(gain);
mask[dst] = masked_by_file[src] ? (1u << PixelMask::ErrorPixelBit) : 0;
}
}
}
dataset_ = std::make_shared<JFJochReaderDataset>();
dataset_->experiment = default_experiment;
DetectorSetup detector = DetDECTRIS(static_cast<int64_t>(width_), static_cast<int64_t>(height_),
"placed by " + source_, {});
detector.PixelSize_um(static_cast<float>(1e6 / geom_.panels[0].res));
// Nothing a parallax or sensor-efficiency correction needs is stated in these files.
detector.SensorThickness_um(0);
// Where nothing states a limit, the stored integer type still has one - a Cheetah file keeps its
// frames as int16, and a value clipped at 32767 is no measurement either.
if (!saturation.has_value() && !layout_.stored_float && layout_.stored_size > 0) {
const double type_max = std::exp2(8.0 * static_cast<double>(layout_.stored_size)
- (layout_.stored_signed ? 1 : 0)) - 1;
saturation = type_max / *std::max_element(adu_per_photon_.begin(), adu_per_photon_.end());
}
if (saturation.has_value())
detector.SaturationLimit(static_cast<int64_t>(std::floor(saturation.value())));
detector.BitDepthImage(32);
detector.MinFrameTime(std::chrono::microseconds(0));
detector.MinCountTime(std::chrono::microseconds(0));
detector.ReadOutTime(std::chrono::nanoseconds(0));
dataset_->experiment.Detector(detector);
// Dark-subtracted counts go below zero; the images are handed out as signed int32.
dataset_->experiment.BitDepthImage(32);
dataset_->experiment.PixelSigned(true);
dataset_->experiment.BeamX_pxl(static_cast<float>(beam_x));
dataset_->experiment.BeamY_pxl(static_cast<float>(beam_y));
dataset_->experiment.DetectorDistance_mm(static_cast<float>(distance_mm_));
if (photon_energy_eV_ > 0)
dataset_->experiment.IncidentEnergy_keV(static_cast<float>(photon_energy_eV_ / 1000.0));
dataset_->experiment.ImagesPerTrigger(static_cast<int64_t>(frames_.size()));
// Each panel is a module, placed where the layout puts it; one panel square with the image needs none.
const DetectorModules placed(modules);
if (modules.size() >= 2 || placed.Moved())
dataset_->experiment.Modules(placed);
dataset_->pixel_mask = std::make_shared<const PixelMask>(mask);
Logger("XFELReader").Info("{} frames in {} files, {} panels placed by {}; {} x {} image, beam at ({:.2f}, "
"{:.2f}) px, {:.3f} mm, {:.1f} eV, {}{}", frames_.size(), files_.size(),
geom_.panels.size(), source_, width_, height_, beam_x, beam_y, distance_mm_,
photon_energy_eV_,
counts_in_photons_ ? "converted to photons" : "in detector units (no gain known)",
layout_.direct_chunks ? ", chunks decoded in parallel" : "");
SetStartMessage(dataset_);
}
void JFJochXFELReader::ApplyGain(float adu_per_photon) {
for (auto &p: geom_.panels)
p.adu_per_photon = adu_per_photon;
Place();
}
bool JFJochXFELReader::HasPhotonEnergy() const {
return photon_energy_eV_ > 0;
}
bool JFJochXFELReader::CountsInPhotons() const {
return counts_in_photons_;
}
uint64_t JFJochXFELReader::GetNumberOfImages() const {
return frames_.size();
}
void JFJochXFELReader::Close() {
std::unique_lock ul(hdf5_mutex);
open_files_.clear();
raw_files_.clear();
user_block_.clear();
raw_order_.clear();
files_.clear();
frames_.clear();
dataset_.reset();
}
template <class Buffer>
void JFJochXFELReader::DecodeInto(int64_t image_number, Buffer &buffer) const {
const auto &frame = frames_.at(image_number);
std::vector<float> stored(frame_width_ * frame_height_);
if (layout_.direct_chunks) {
// Under the lock only where the chunk is (a lookup in its index); the read and the decoding not.
std::shared_ptr<const HDF5ImageSource::RawFile> raw_file;
uint64_t address = 0, chunk_size = 0;
{
std::unique_lock ul(hdf5_mutex);
HDF5DataSet dataset(FileAt(frame.file), frame.dataset);
const hsize_t coord_3d[3] = {static_cast<hsize_t>(frame.index), 0, 0};
const hsize_t coord_2d[2] = {0, 0};
unsigned filter_mask = 0;
haddr_t addr = HADDR_UNDEF;
hsize_t size = 0;
if (H5Dget_chunk_info_by_coord(dataset.GetID(), frame.index >= 0 ? coord_3d : coord_2d, &filter_mask,
&addr, &size) < 0 || addr == HADDR_UNDEF || size == 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"No stored chunk for " + frame.dataset + " of " + files_[frame.file]);
raw_file = RawFileAt(frame.file);
address = user_block_[frame.file] + addr;
chunk_size = size;
}
std::vector<uint8_t> compressed(chunk_size);
raw_file->ReadAt(compressed.data(), compressed.size(), address);
const size_t n = stored.size();
const size_t size = layout_.stored_size;
std::vector<uint8_t> raw(n * size);
if (layout_.deflated) {
uLongf raw_size = raw.size();
if (uncompress(raw.data(), &raw_size, compressed.data(), compressed.size()) != Z_OK
|| raw_size != raw.size())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Cannot inflate " + frame.dataset + " of " + files_[frame.file]);
} else {
const CompressedImage image(compressed.data(), compressed.size(), frame_width_, frame_height_,
CalcImageMode(size, layout_.stored_float, layout_.stored_signed),
layout_.algorithm);
std::vector<uint8_t> decoded;
std::memcpy(raw.data(), image.GetUncompressedPtr(decoded), raw.size());
}
if (layout_.shuffled) {
std::vector<uint8_t> tmp(raw.size());
for (size_t i = 0; i < n; i++)
for (size_t b = 0; b < size; b++)
tmp[i * size + b] = raw[b * n + i];
raw.swap(tmp);
}
for (size_t i = 0; i < n; i++) {
const uint8_t *v = raw.data() + i * size;
if (layout_.stored_float) {
float x;
std::memcpy(&x, v, 4);
stored[i] = x;
} else if (size == 2) {
uint16_t x;
std::memcpy(&x, v, 2);
stored[i] = layout_.stored_signed ? static_cast<float>(static_cast<int16_t>(x)) : x;
} else {
uint32_t x;
std::memcpy(&x, v, 4);
stored[i] = layout_.stored_signed ? static_cast<float>(static_cast<int32_t>(x))
: static_cast<float>(x);
}
}
} else {
std::unique_lock ul(hdf5_mutex);
HDF5DataSet dataset(FileAt(frame.file), frame.dataset);
if (frame.index >= 0)
dataset.ReadVector(stored, {static_cast<hsize_t>(frame.index), 0, 0}, {1, frame_height_, frame_width_});
else
dataset.ReadVector(stored, {0, 0}, {frame_height_, frame_width_});
}
// A pixel a per-frame mask marks in this frame is handed out as no measurement.
std::vector<uint8_t> masked;
if (!frame_masks_.empty() && frame.index >= 0) {
masked.assign(stored.size(), 0);
std::unique_lock ul(hdf5_mutex);
for (const auto &m: frame_masks_) {
const auto v = FileAt(frame.file).ReadVector<uint32_t>(m.data, {static_cast<hsize_t>(frame.index), 0, 0},
{1, frame_height_, frame_width_});
for (size_t i = 0; i < v.size(); i++)
if (m.Bad(v[i]))
masked[i] = 1;
}
}
// And so is one whose value the layout flags as no measurement.
const auto flagged = [this](float v) {
return (geom_.flag_lessthan && v < *geom_.flag_lessthan) || (geom_.flag_morethan && v > *geom_.flag_morethan)
|| (geom_.flag_equal && v == *geom_.flag_equal);
};
buffer.assign(width_ * height_ * sizeof(int32_t), 0);
auto out = reinterpret_cast<int32_t *>(buffer.data());
for (size_t i = 0; i < stored.size(); i++)
if (destination_[i] >= 0)
out[destination_[i]] = (!masked.empty() && masked[i]) || flagged(stored[i])
? INT32_MIN
: static_cast<int32_t>(std::lround(stored[i] / adu_per_photon_[i]));
}
bool JFJochXFELReader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &dataset,
DataMessage &message,
std::vector<uint8_t> &buffer,
int64_t image_number,
bool update_dataset) {
(void) update_dataset;
if (!dataset || image_number < 0 || static_cast<size_t>(image_number) >= frames_.size())
return false;
DecodeInto(image_number, buffer);
message.image = CompressedImage(buffer.data(), buffer.size(), width_, height_,
CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION);
message.number = image_number;
return true;
}
bool JFJochXFELReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) {
if (image_number < 0 || static_cast<size_t>(image_number) >= frames_.size())
return false;
DecodeInto(image_number, image.image_buffer);
image.image = CompressedImage(image.image_buffer.data(), image.image_buffer.size(), width_, height_,
CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION);
return true;
}
std::vector<SpotToSave> JFJochXFELReader::ReadSpots(int64_t) const {
return {};
}
Reference in new issue View Git Blame Copy Permalink