mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 87 Wiki Activity
Files
d345deaecdcfd6cff1590098a4c80a7a998056c5
Jungfraujoch/reader/JFJochXFELReader.cpp
T
leonarski_fandClaude Opus 5.5 d345deaecd XFELReader: per-frame mask and photon conversion without a library call or a branch per pixel 
On a deflated Cheetah CXI with a per-frame mask (AGIPD 1M), half of the reader's time after the inflate
went to turning the stored frame into the image: the mask plane widened into a zeroed uint32 copy before
being tested, MarkBad testing it with a branch per pixel through a reference the store could alias, and
the photon conversion calling lroundf and branching on the mask, the finiteness and the clamp. Now the
mask is tested as stored (uint16 or uint32), and both loops are written so that the compiler vectorises
them: std::lround spelled out exactly (truncate, the exact fraction, a half or more away from zero) and
the selects done as arithmetic (Pick). Plain C++, no intrinsics.

Bit-identical: FNV hash of every image handed out (2000 frames) plus the pixel mask, before and after,
on a CXI with a CrystFEL geometry (with and without a gain), a panel-table HDF5 (with and without a
gain) and a bitshuffle/LZ4 int16 set: all identical. Reader alone, one thread: 14.1 -> 7.0 ms per
frame. rugnux on the 14,445-frame CXI run (16 workers, A/B back to back): geometry refinement sample
11.4 -> 8.9 s, image pass 27.6 -> 21.1 s, whole run 67.4 -> 58.2 s; geometry and refinement log lines
identical.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
2026-10-23 10:47:56 +02:00

1053 lines
51 KiB
C++
Raw 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 <type_traits>
#include <limits>
#include <cstring>
#include <filesystem>
#include <numeric>
#include <zlib.h>
#include "../common/JFJochException.h"
#include "../common/Logger.h"
#include "../writer/HDF5Objects.h"
#include "../compression/JFJochDecompress.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);
}
// b where use_b, otherwise a - as arithmetic rather than a branch, which the compiler would move the work
// computing a into, and then not vectorise the loop over the pixels.
int32_t Pick(bool use_b, int32_t a, int32_t b) {
const int32_t all = -static_cast<int32_t>(use_b);
return (a & ~all) | (b & all);
}
} // 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("{}: the module geometry, camera length and masks are taken; ignored: {}",
geometry_file, list);
}
Logger("XFELReader").Info("{}: {} mask array(s), {} bad region(s)", geometry_file, geom_.masks.size(),
geom_.bad.size());
// What the geometry file does not state 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).
// The file's own masks are taken only where the geometry names none: a geometry that does names the
// mask arrays with the bits that matter for this detector. 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]);
const bool file_masks = geom_.masks.empty();
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_masks && file.Exists(CXI_DETECTOR + "mask_shared"))
geom_.masks.push_back({CXI_DETECTOR + "mask_shared", "", 0, CHEETAH_BAD_BITS});
if (file_masks && 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};
if (file_masks)
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;
l.plain_type = ((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;
ok &= l.plain_type;
H5Tclose(type);
l.deflated = deflate;
l.direct_chunks = ok && (deflate || (bitshuffle && l.algorithm != CompressionAlgorithm::NO_COMPRESSION));
return l;
}
// A per-frame mask is read as the frames are where it is stored in plain unsigned integers - read as
// stored, they are what HDF5 would convert them to - of the frame's size, in every file the way the
// first file stores it. Any other is read by HDF5.
JFJochXFELReader::Layout JFJochXFELReader::FrameMaskLayout(const xfel::Mask &m) const {
Layout layout;
for (size_t i = 0; i < frames_.size(); i++) {
if (i > 0 && frames_[i].file == frames_[i - 1].file)
continue; // a file's first frame stands for all of its frames, as for the frames themselves
const Frame mask{frames_[i].file, ReplaceEvent(m.data, frames_[i].group), frames_[i].index, frames_[i].group};
if (!FileAt(mask.file).Exists(mask.dataset))
return {};
const Layout l = ReadLayout(FileAt(mask.file), mask);
if (i == 0)
layout = l;
else if (!(l == layout))
return {};
}
if (layout.stored_float || layout.stored_signed || layout.height != frame_height_ || layout.width != frame_width_)
return {};
return layout;
}
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();
frame_mask_layouts_.clear();
for (const auto &m: geom_.masks) {
std::vector<uint32_t> v;
if (m.file.empty() && (m.data.find('%') != std::string::npos
|| FileAt(0).GetDimension(m.data).size() == 3)) {
frame_masks_.push_back(m); // one per shot or one plane per frame: applied as each is read
frame_mask_layouts_.push_back(FrameMaskLayout(m));
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_});
MarkBad(m, v.data(), masked_by_file);
}
}
// CrystFEL's laboratory frame and the detector plane'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)};
};
// A detector of several panels is analysed as stored: the image is the stored frame, and each panel a
// module that maps its box of the frame onto the detector plane exactly where the layout puts it.
// A single panel goes onto an image of its own 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 then its module's placement, so no pixel is moved on its own.
const bool stored = geom_.panels.size() >= 2;
struct Block {
int64_t fx = 1, fy = 0, sx = 0, sy = 1; // image step per +1 fs / +1 ss (single panel)
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;
if (!stored) {
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 auto [x, y] = stored
? position(p, fs + 0.5, ss + 0.5)
: std::pair{bl.x0 + fs * static_cast<double>(bl.fx) + ss * static_cast<double>(bl.sx),
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 first pixel centre on a whole column and row, 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_ = stored ? frame_width_ : static_cast<size_t>(std::ceil(beam_x + max_x)) + 3;
height_ = stored ? frame_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);
blocks_.clear();
// One gain for every panel: what the run is given (ReadFiles, ApplyGain) is set on every panel alike.
const auto &p0 = geom_.panels[0];
counts_in_photons_ = p0.adu_per_photon.has_value() || p0.adu_per_eV.has_value();
const double gain = p0.adu_per_photon.has_value() ? p0.adu_per_photon.value()
: p0.adu_per_eV.has_value() ? p0.adu_per_eV.value() * photon_energy_eV_ : 1.0;
adu_per_photon_ = static_cast<float>(gain);
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");
if (p.max_adu.has_value())
saturation = std::min(saturation.value_or(INT64_MAX), p.max_adu.value() / gain);
if (p.no_index)
continue;
const int64_t w = p.max_fs - p.min_fs, h = p.max_ss - p.min_ss;
const double scale = distance_mm_ / ((clen_m_ + p.coffset) * 1000.0);
// The image pixel the panel's first pixel lands on (single panel).
const int64_t c0 = std::lround(beam_x + bl.x0), r0 = std::lround(beam_y + bl.y0);
DetectorModule m;
if (stored) {
// Its box in the stored frame, and the plane position and steps of its pixels: the panel's own axes,
// scaled by its projection.
m.x0 = p.min_fs;
m.y0 = p.min_ss;
m.width = w + 1;
m.height = h + 1;
m.ox = beam_x + bl.x0;
m.oy = beam_y + bl.y0;
m.fx = scale * p.fs_x;
m.fy = scale * p.fs_y;
m.sx = scale * p.ss_x;
m.sy = scale * p.ss_y;
} else {
// The block it fills, and 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 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);
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);
blocks_.push_back({static_cast<size_t>(p.min_ss) * frame_width_ + static_cast<size_t>(p.min_fs), w + 1,
h + 1, r0 * static_cast<int64_t>(width_) + c0,
bl.fy * static_cast<int64_t>(width_) + bl.fx,
bl.sy * static_cast<int64_t>(width_) + bl.sx});
}
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 = stored ? fs : c0 + (fs - p.min_fs) * bl.fx + (ss - p.min_ss) * bl.sx;
const int64_t row = stored ? ss : 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;
mask[dst] = masked_by_file[ss * frame_width_ + fs] ? (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 / adu_per_photon_;
}
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 (stored || 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_);
}
template <class T>
void JFJochXFELReader::MarkBad(const xfel::Mask &m, const T *v, std::vector<uint8_t> &bad) const {
int64_t min_fs = 0, max_fs = static_cast<int64_t>(frame_width_) - 1;
int64_t min_ss = 0, max_ss = static_cast<int64_t>(frame_height_) - 1;
for (const auto &p: geom_.panels)
if (!m.panel.empty() && p.name == m.panel) {
min_fs = p.min_fs; max_fs = p.max_fs;
min_ss = p.min_ss; max_ss = p.max_ss;
}
// The mask's bits in locals and plain pointers, a row at a time: a byte stored through bad could otherwise
// be m's or bad's own, every pixel would read them again, and the loop would not be vectorised.
const uint32_t good = m.good, bad_bits = m.bad;
const bool nonzero_good = m.nonzero_good;
const size_t width = max_fs - min_fs + 1;
for (int64_t ss = min_ss; ss <= max_ss; ss++) {
const size_t first = static_cast<size_t>(ss) * frame_width_ + min_fs;
const T *in = v + first;
uint8_t *out = bad.data() + first;
for (size_t i = 0; i < width; i++)
out[i] |= xfel::Mask::Bad(in[i], good, bad_bits, nonzero_good);
}
}
void JFJochXFELReader::ApplyGain(float adu_per_photon) {
for (auto &p: geom_.panels)
p.adu_per_photon = adu_per_photon;
Place();
}
void JFJochXFELReader::ClearGain() {
for (auto &p: geom_.panels)
p.adu_per_photon.reset();
Place();
}
bool JFJochXFELReader::StoredFloat() const {
return layout_.stored_float || !layout_.plain_type;
}
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>
bool JFJochXFELReader::ReadStored(const Frame &frame, const Layout &layout, Buffer &stored,
bool keep_compressed) const {
const std::vector<hsize_t> start = frame.index >= 0
? std::vector<hsize_t>{static_cast<hsize_t>(frame.index), 0, 0}
: std::vector<hsize_t>{0, 0};
const std::vector<hsize_t> count = frame.index >= 0 ? std::vector<hsize_t>{1, frame_height_, frame_width_}
: std::vector<hsize_t>{frame_height_, frame_width_};
// 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);
unsigned filter_mask = 0;
if (layout.direct_chunks) {
haddr_t addr = HADDR_UNDEF;
hsize_t size = 0;
if (H5Dget_chunk_info_by_coord(dataset.GetID(), start.data(), &filter_mask, &addr, &size) < 0
|| addr == HADDR_UNDEF || size == 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"No stored chunk for " + frame.dataset + " of " + files_[frame.file]);
address = addr;
chunk_size = size;
}
// A chunk stored with one of its filters skipped (its filter mask) is not the encoding decoded
// below: HDF5 reads it, as it reads a frame that is not stored as one chunk.
if (!layout.direct_chunks || filter_mask != 0) {
if (layout.plain_type)
dataset.ReadVectorToU8(stored, start, count);
else {
std::vector<float> v(frame_width_ * frame_height_);
dataset.ReadVector(v, start, count);
stored.resize(v.size() * sizeof(float));
std::memcpy(stored.data(), v.data(), stored.size());
}
return false;
}
raw_file = RawFileAt(frame.file);
address += user_block_[frame.file]; // HDF5's chunk addresses count from the end of the user block
}
if (keep_compressed && !layout.deflated) {
stored.resize(chunk_size);
raw_file->ReadAt(stored.data(), stored.size(), address);
return true;
}
thread_local std::vector<uint8_t> compressed, shuffled;
compressed.resize(chunk_size);
raw_file->ReadAt(compressed.data(), compressed.size(), address);
const size_t n = frame_width_ * frame_height_;
const size_t size = layout.stored_size;
if (layout.deflated) {
if (!layout.shuffled)
stored.resize(n * size);
else
shuffled.resize(n * size);
uint8_t *inflated = layout.shuffled ? shuffled.data() : stored.data();
uLongf inflated_size = n * size;
if (uncompress(inflated, &inflated_size, compressed.data(), compressed.size()) != Z_OK
|| inflated_size != n * size)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Cannot inflate " + frame.dataset + " of " + files_[frame.file]);
if (layout.shuffled) {
stored.resize(n * size);
for (size_t i = 0; i < n; i++)
for (size_t b = 0; b < size; b++)
stored[i * size + b] = shuffled[b * n + i];
}
} else {
stored.resize(n * size);
JFJochDecompressPtr(stored.data(), layout.algorithm, compressed.data(), compressed.size(), n, size);
}
return false;
}
template <class Out, class In, class Buffer, class Convert>
void JFJochXFELReader::PlaceInto(const In *in, const std::vector<uint8_t> &masked, Out err, Buffer &buffer,
Convert convert) const {
buffer.resize(width_ * height_ * sizeof(Out));
auto out = reinterpret_cast<Out *>(buffer.data());
if (blocks_.empty()) { // the stored frame, as stored
const size_t n = width_ * height_;
if (masked.empty())
for (size_t i = 0; i < n; i++)
out[i] = static_cast<Out>(convert(in[i]));
else {
const uint8_t *m = masked.data();
for (size_t i = 0; i < n; i++)
out[i] = static_cast<Out>(Pick(m[i] != 0, convert(in[i]), err));
}
return;
}
std::fill(out, out + width_ * height_, Out{0}); // between the panels
for (const auto &b: blocks_)
for (int64_t ss = 0; ss < b.height; ss++) {
const size_t src = b.src + static_cast<size_t>(ss) * frame_width_;
Out *dst = out + b.dst + ss * b.step_ss;
if (masked.empty() && b.step_fs == 1)
std::transform(in + src, in + src + b.width, dst, [&](In v) { return static_cast<Out>(convert(v)); });
else
for (int64_t fs = 0; fs < b.width; fs++)
dst[fs * b.step_fs] = !masked.empty() && masked[src + fs] ? err : static_cast<Out>(convert(in[src + fs]));
}
}
template <class Buffer>
CompressedImage JFJochXFELReader::DecodeInto(int64_t image_number, Buffer &buffer) const {
const auto &frame = frames_.at(image_number);
const CompressedImageMode stored_mode = CalcImageMode(layout_.stored_size, layout_.stored_float,
layout_.stored_signed);
// A frame of several panels in signed integers, read in its own units, goes up as stored, with no copy:
// a bitshuffle/LZ4 chunk still compressed, to be decoded where it is analysed (on the GPU, as
// Jungfraujoch's own images are), anything else decoded once. Its smallest value is no measurement
// (below), which the analysis takes it as anyway.
const bool as_stored = blocks_.empty() && layout_.plain_type && !layout_.stored_float && layout_.stored_signed
&& adu_per_photon_ == 1.0f;
if (as_stored && frame_masks_.empty()) {
const bool compressed = ReadStored(frame, layout_, buffer, true);
return CompressedImage(buffer.data(), buffer.size(), width_, height_, stored_mode,
compressed ? layout_.algorithm : CompressionAlgorithm::NO_COMPRESSION);
}
thread_local std::vector<uint8_t> stored;
ReadStored(frame, layout_, stored, false);
// 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()) {
masked.assign(frame_width_ * frame_height_, 0);
for (size_t i = 0; i < frame_masks_.size(); i++) {
const auto &m = frame_masks_[i];
const auto &mask_layout = frame_mask_layouts_[i];
const Frame mask{frame.file, ReplaceEvent(m.data, frame.group), frame.index, frame.group};
if (mask_layout.plain_type) {
thread_local std::vector<uint8_t> mask_stored;
ReadStored(mask, mask_layout, mask_stored, false);
if (mask_layout.stored_size == 2)
MarkBad(m, reinterpret_cast<const uint16_t *>(mask_stored.data()), masked);
else
MarkBad(m, reinterpret_cast<const uint32_t *>(mask_stored.data()), masked);
} else {
std::vector<uint32_t> v;
{
std::unique_lock ul(hdf5_mutex);
v = frame.index >= 0
? FileAt(frame.file).ReadVector<uint32_t>(mask.dataset, {static_cast<hsize_t>(frame.index), 0, 0},
{1, frame_height_, frame_width_})
: FileAt(frame.file).ReadVector<uint32_t>(mask.dataset, {0, 0}, {frame_height_, frame_width_});
}
MarkBad(m, v.data(), masked);
}
}
}
const float gain = adu_per_photon_;
// A float that is no number (NaN, inf) is no measurement; anything past int32 is clamped to it.
// Written without a branch or a library call, so that the loop over the pixels is vectorised. The clamp
// keeps the sign of an infinity; std::lround is spelled out: t is c truncated, c - t is exact, and a
// fraction of a half or more rounds away from zero.
const auto photons = [gain](auto v) {
const float x = static_cast<float>(v) / gain;
constexpr float limit = 2147483520.0f;
const float c = std::fabs(x) < limit ? x : std::copysign(limit, x);
const auto t = static_cast<int32_t>(c);
const float fraction = c - static_cast<float>(t);
const int32_t rounded = t + (fraction >= 0.5f) - (fraction <= -0.5f);
return Pick(!std::isfinite(x), rounded, INT32_MIN);
};
const auto as_is = [](auto v) { return v; };
const auto image = [&](CompressedImageMode mode) {
return CompressedImage(buffer.data(), buffer.size(), width_, height_, mode, CompressionAlgorithm::NO_COMPRESSION);
};
// A signed integer frame's smallest value is no measurement: what a converted file
// (tools/battery/still/agipd_to_int16.py) writes for a masked pixel or a gap.
const auto integers = [&](auto to_int32) {
const auto convert = [to_int32](auto v) -> int32_t {
if constexpr (std::is_signed_v<decltype(v)>)
return Pick(v == std::numeric_limits<decltype(v)>::min(), to_int32(v), INT32_MIN);
return to_int32(v);
};
const uint8_t *p = stored.data();
if (layout_.stored_size == 2 && layout_.stored_signed)
PlaceInto<int32_t>(reinterpret_cast<const int16_t *>(p), masked, INT32_MIN, buffer, convert);
else if (layout_.stored_size == 2)
PlaceInto<int32_t>(reinterpret_cast<const uint16_t *>(p), masked, INT32_MIN, buffer, convert);
else if (layout_.stored_signed)
PlaceInto<int32_t>(reinterpret_cast<const int32_t *>(p), masked, INT32_MIN, buffer, convert);
else
PlaceInto<int32_t>(reinterpret_cast<const uint32_t *>(p), masked, INT32_MIN, buffer, convert);
};
// Floats become photons, rounded. Integers stay integers: divided only by a gain other than one, and a
// 16-bit signed frame - as Cheetah writes - is handed out in its own type.
if (layout_.stored_float || !layout_.plain_type)
PlaceInto<int32_t>(reinterpret_cast<const float *>(stored.data()), masked, INT32_MIN, buffer, photons);
else if (gain != 1.0f)
integers(photons);
else if (layout_.stored_size == 2 && layout_.stored_signed) {
PlaceInto<int16_t>(reinterpret_cast<const int16_t *>(stored.data()), masked, INT16_MIN, buffer, as_is);
return image(CompressedImageMode::Int16);
} else
integers([](auto v) { return static_cast<int32_t>(v); });
return image(CompressedImageMode::Int32);
}
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;
message.image = DecodeInto(image_number, buffer);
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;
image.image = DecodeInto(image_number, image.image_buffer);
return true;
}
std::vector<SpotToSave> JFJochXFELReader::ReadSpots(int64_t) const {
return {};
}
Reference in new issue View Git Blame Copy Permalink