Merge branch 'd-cbf-gaps' into pool-C
This commit is contained in:
@@ -2,6 +2,8 @@ ADD_LIBRARY(JFJochReader STATIC
|
||||
JFJochReader.cpp JFJochReader.h
|
||||
JFJochHDF5Reader.cpp
|
||||
JFJochHDF5Reader.h
|
||||
SweepLayout.cpp
|
||||
SweepLayout.h
|
||||
MiniCBF.cpp
|
||||
MiniCBF.h
|
||||
JFJochCBFReader.cpp
|
||||
|
||||
+33
-15
@@ -17,6 +17,7 @@
|
||||
#include "../common/JFJochException.h"
|
||||
#include "../common/Logger.h"
|
||||
#include "../common/JFJochMath.h"
|
||||
#include "SweepLayout.h"
|
||||
|
||||
namespace {
|
||||
|
||||
@@ -286,32 +287,33 @@ void JFJochCBFReader::ReadFiles(const std::string &path) {
|
||||
std::chrono::duration_cast<std::chrono::nanoseconds>(
|
||||
std::chrono::duration<double>(header0_.exposure_s)));
|
||||
|
||||
// The rotation angle of every image, from its own header. Parsing one costs a few hundred
|
||||
// microseconds - it is two dozen regular expressions - so on a sweep of several thousand frames
|
||||
// this is seconds of startup before a single image is read, and the files are independent.
|
||||
std::vector<double> angles(files_.size());
|
||||
// Where every image sits on the spindle, and how the instrument stood, from its own header.
|
||||
// Parsing one costs a few hundred microseconds - it is two dozen regular expressions - so on a
|
||||
// sweep of several thousand frames this is seconds of startup before a single image is read, and
|
||||
// the files are independent.
|
||||
std::vector<sweep::Frame> frames(files_.size());
|
||||
{
|
||||
const size_t nthreads = std::min<size_t>(std::max(1u, std::thread::hardware_concurrency()), 8);
|
||||
std::vector<std::future<void>> futures;
|
||||
for (size_t t = 0; t < nthreads; t++)
|
||||
futures.push_back(std::async(std::launch::async, [&, t] {
|
||||
for (size_t i = t; i < files_.size(); i += nthreads)
|
||||
angles[i] = minicbf::ReadHeader(files_[i]).start_angle_deg;
|
||||
for (size_t i = t; i < files_.size(); i += nthreads) {
|
||||
const auto h = minicbf::ReadHeader(files_[i]);
|
||||
frames[i] = {files_[i], h.start_angle_deg, h.angle_increment_deg, h.distance_m,
|
||||
h.beam_x_px, h.beam_y_px, h.wavelength_A};
|
||||
}
|
||||
}));
|
||||
for (auto &f : futures)
|
||||
f.get();
|
||||
}
|
||||
|
||||
double increment = header0_.angle_increment_deg;
|
||||
if (files_.size() > 1) {
|
||||
// Prefer the measured step over the header's nominal one, and unwrap a sweep that passes 360.
|
||||
double d = angles[1] - angles[0];
|
||||
if (d < -180.0) d += 360.0;
|
||||
if (std::abs(d) > 1e-6) increment = d;
|
||||
}
|
||||
// The sweep the headers describe, which is not always the files laid out end to end: a deposited
|
||||
// series can be missing frames, and those are gaps in the rotation rather than images to close up.
|
||||
const auto layout = sweep::Place(frames, "CBFReader");
|
||||
files_ = layout.files;
|
||||
dataset_->experiment.Goniometer(GoniometerAxis(header0_.axis_name,
|
||||
static_cast<float>(angles.front()),
|
||||
static_cast<float>(increment),
|
||||
static_cast<float>(layout.start_deg),
|
||||
static_cast<float>(layout.increment_deg),
|
||||
RotationAxis(header0_), {}));
|
||||
|
||||
dataset_->error_value = -1;
|
||||
@@ -345,6 +347,9 @@ CompressedImage JFJochCBFReader::DecodeInto(int64_t image_number, Buffer &buffer
|
||||
if (image_number < 0 || static_cast<size_t>(image_number) >= files_.size())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"Image number out of range");
|
||||
if (files_[image_number].empty())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"No image at this point of the sweep");
|
||||
|
||||
const size_t npixel = static_cast<size_t>(header0_.nx) * static_cast<size_t>(header0_.ny);
|
||||
buffer.resize(npixel * sizeof(int32_t));
|
||||
@@ -371,6 +376,9 @@ bool JFJochCBFReader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &dataset,
|
||||
if (!dataset)
|
||||
return false;
|
||||
|
||||
if (!HasImage(image_number))
|
||||
return false;
|
||||
|
||||
// The image must outlive this call, so it is decoded straight into the caller's buffer - the same
|
||||
// thing the argument is for on the HDF5 path - and message.image only points at it.
|
||||
std::vector<uint8_t> scratch;
|
||||
@@ -379,7 +387,17 @@ bool JFJochCBFReader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &dataset,
|
||||
return true;
|
||||
}
|
||||
|
||||
// A slot the series has no file for is a missing image, not an error: every image loop in the
|
||||
// pipeline already treats "nothing to read" as a frame to pass over, which is exactly what a gap in
|
||||
// a deposited sweep is.
|
||||
bool JFJochCBFReader::HasImage(int64_t image_number) const {
|
||||
return image_number >= 0 && static_cast<size_t>(image_number) < files_.size()
|
||||
&& !files_[image_number].empty();
|
||||
}
|
||||
|
||||
bool JFJochCBFReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) {
|
||||
if (!HasImage(image_number))
|
||||
return false;
|
||||
image.image = DecodeInto(image_number, image.image_buffer, image.read_buffer);
|
||||
return true;
|
||||
}
|
||||
|
||||
@@ -50,6 +50,12 @@ public:
|
||||
[[nodiscard]] uint64_t GetNumberOfImages() const override;
|
||||
void Close() override;
|
||||
|
||||
|
||||
// Whether the sweep really has an image at this point. A deposited series can be missing frames,
|
||||
// and they are kept as gaps so that every image keeps the spindle angle its own header states;
|
||||
// ReadRawImage hands such a slot back as "nothing to read".
|
||||
[[nodiscard]] bool HasImage(int64_t image_number) const;
|
||||
|
||||
bool ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) override;
|
||||
[[nodiscard]] std::vector<SpotToSave> ReadSpots(int64_t image) const override;
|
||||
};
|
||||
|
||||
@@ -10,6 +10,7 @@
|
||||
#include "../common/JFJochException.h"
|
||||
#include "../common/JFJochMath.h"
|
||||
#include "../common/Logger.h"
|
||||
#include "SweepLayout.h"
|
||||
|
||||
namespace {
|
||||
|
||||
@@ -106,32 +107,34 @@ void JFJochMarCCDReader::ReadFiles(const std::string &path) {
|
||||
"default frame time is kept. Nothing in the geometry or the "
|
||||
"merge depends on it.", files_[0], header0_.exposure_s);
|
||||
|
||||
// The rotation angle of every image, from its own header. Reading one costs a 4 kB read, so on
|
||||
// a sweep of several thousand frames this is worth spreading over the cores, as the CBF path
|
||||
// does for the same reason.
|
||||
std::vector<double> angles(files_.size());
|
||||
// Where every image sits on the spindle, and how the instrument stood, from its own header.
|
||||
// Reading one costs a 4 kB read, so on a sweep of several thousand frames this is worth
|
||||
// spreading over the cores, as the CBF path does for the same reason.
|
||||
std::vector<sweep::Frame> frames(files_.size());
|
||||
{
|
||||
const size_t nthreads = std::min<size_t>(std::max(1u, std::thread::hardware_concurrency()), 8);
|
||||
std::vector<std::future<void>> futures;
|
||||
for (size_t t = 0; t < nthreads; t++)
|
||||
futures.push_back(std::async(std::launch::async, [&, t] {
|
||||
for (size_t i = t; i < files_.size(); i += nthreads)
|
||||
angles[i] = marccd::ReadHeader(files_[i]).start_angle_deg;
|
||||
for (size_t i = t; i < files_.size(); i += nthreads) {
|
||||
const auto h = marccd::ReadHeader(files_[i]);
|
||||
frames[i] = {files_[i], h.start_angle_deg, h.angle_increment_deg, h.distance_m,
|
||||
h.beam_x_px, h.beam_y_px, h.wavelength_A};
|
||||
}
|
||||
}));
|
||||
for (auto &f : futures)
|
||||
f.get();
|
||||
}
|
||||
|
||||
double increment = header0_.angle_increment_deg;
|
||||
if (files_.size() > 1) {
|
||||
// Prefer the measured step over the header's nominal one, and unwrap a sweep that passes 360.
|
||||
double d = angles[1] - angles[0];
|
||||
if (d < -180.0) d += 360.0;
|
||||
if (std::abs(d) > 1e-6) increment = d;
|
||||
}
|
||||
// The sweep the headers describe, which is not always the files laid out end to end: a deposited
|
||||
// series can be missing frames, and those are gaps in the rotation rather than images to close
|
||||
// up. A folder of screening shots taken at scattered angles is refused here by name instead of
|
||||
// failing later as a lattice nobody can explain.
|
||||
const auto layout = sweep::Place(frames, "MarCCDReader");
|
||||
files_ = layout.files;
|
||||
dataset_->experiment.Goniometer(GoniometerAxis(header0_.axis_name,
|
||||
static_cast<float>(angles.front()),
|
||||
static_cast<float>(increment),
|
||||
static_cast<float>(layout.start_deg),
|
||||
static_cast<float>(layout.increment_deg),
|
||||
ASSUMED_BASE_AXIS, {}));
|
||||
|
||||
dataset_->error_value = -1;
|
||||
@@ -159,6 +162,9 @@ CompressedImage JFJochMarCCDReader::DecodeInto(int64_t image_number, Buffer &buf
|
||||
if (image_number < 0 || static_cast<size_t>(image_number) >= files_.size())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"Image number out of range");
|
||||
if (files_[image_number].empty())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"No image at this point of the sweep");
|
||||
|
||||
const size_t npixel = static_cast<size_t>(header0_.nx) * static_cast<size_t>(header0_.ny);
|
||||
buffer.resize(npixel * sizeof(int32_t));
|
||||
@@ -183,13 +189,26 @@ bool JFJochMarCCDReader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &datas
|
||||
if (!dataset)
|
||||
return false;
|
||||
|
||||
if (!HasImage(image_number))
|
||||
return false;
|
||||
|
||||
std::vector<uint8_t> scratch;
|
||||
message.image = DecodeInto(image_number, buffer, scratch);
|
||||
message.number = image_number;
|
||||
return true;
|
||||
}
|
||||
|
||||
// A slot the series has no file for is a missing image, not an error: every image loop in the
|
||||
// pipeline already treats "nothing to read" as a frame to pass over, which is exactly what a gap in
|
||||
// a deposited sweep is.
|
||||
bool JFJochMarCCDReader::HasImage(int64_t image_number) const {
|
||||
return image_number >= 0 && static_cast<size_t>(image_number) < files_.size()
|
||||
&& !files_[image_number].empty();
|
||||
}
|
||||
|
||||
bool JFJochMarCCDReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) {
|
||||
if (!HasImage(image_number))
|
||||
return false;
|
||||
image.image = DecodeInto(image_number, image.image_buffer, image.read_buffer);
|
||||
return true;
|
||||
}
|
||||
|
||||
@@ -52,6 +52,12 @@ public:
|
||||
[[nodiscard]] uint64_t GetNumberOfImages() const override;
|
||||
void Close() override;
|
||||
|
||||
|
||||
// Whether the sweep really has an image at this point. A deposited series can be missing frames,
|
||||
// and they are kept as gaps so that every image keeps the spindle angle its own header states;
|
||||
// ReadRawImage hands such a slot back as "nothing to read".
|
||||
[[nodiscard]] bool HasImage(int64_t image_number) const;
|
||||
|
||||
bool ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) override;
|
||||
[[nodiscard]] std::vector<SpotToSave> ReadSpots(int64_t image) const override;
|
||||
};
|
||||
|
||||
+34
-15
@@ -10,6 +10,7 @@
|
||||
#include "../common/JFJochException.h"
|
||||
#include "../common/JFJochMath.h"
|
||||
#include "../common/Logger.h"
|
||||
#include "SweepLayout.h"
|
||||
|
||||
namespace {
|
||||
|
||||
@@ -107,32 +108,34 @@ void JFJochSMVReader::ReadFiles(const std::string &path) {
|
||||
"default frame time is kept. Nothing in the geometry or the "
|
||||
"merge depends on it.", files_[0], header0_.exposure_s);
|
||||
|
||||
// The rotation angle of every image, from its own header. Reading one costs a 4 kB read, so on
|
||||
// a sweep of several thousand frames this is worth spreading over the cores, as the CBF path
|
||||
// does for the same reason.
|
||||
std::vector<double> angles(files_.size());
|
||||
// Where every image sits on the spindle, and how the instrument stood, from its own header.
|
||||
// Reading one costs a 4 kB read, so on a sweep of several thousand frames this is worth
|
||||
// spreading over the cores, as the CBF path does for the same reason.
|
||||
std::vector<sweep::Frame> frames(files_.size());
|
||||
{
|
||||
const size_t nthreads = std::min<size_t>(std::max(1u, std::thread::hardware_concurrency()), 8);
|
||||
std::vector<std::future<void>> futures;
|
||||
for (size_t t = 0; t < nthreads; t++)
|
||||
futures.push_back(std::async(std::launch::async, [&, t] {
|
||||
for (size_t i = t; i < files_.size(); i += nthreads)
|
||||
angles[i] = smv::ReadHeader(files_[i]).start_angle_deg;
|
||||
for (size_t i = t; i < files_.size(); i += nthreads) {
|
||||
const auto h = smv::ReadHeader(files_[i]);
|
||||
frames[i] = {files_[i], h.start_angle_deg, h.angle_increment_deg, h.distance_m,
|
||||
h.beam_x_px, h.beam_y_px, h.wavelength_A};
|
||||
}
|
||||
}));
|
||||
for (auto &f : futures)
|
||||
f.get();
|
||||
}
|
||||
|
||||
double increment = header0_.angle_increment_deg;
|
||||
if (files_.size() > 1) {
|
||||
// Prefer the measured step over the header's nominal one, and unwrap a sweep that passes 360.
|
||||
double d = angles[1] - angles[0];
|
||||
if (d < -180.0) d += 360.0;
|
||||
if (std::abs(d) > 1e-6) increment = d;
|
||||
}
|
||||
// The sweep the headers describe, which is not always the files laid out end to end: a deposited
|
||||
// series can be missing frames, and those are gaps in the rotation rather than images to close
|
||||
// up. A folder of screening shots taken at scattered angles is refused here by name instead of
|
||||
// failing later as a lattice nobody can explain.
|
||||
const auto layout = sweep::Place(frames, "SMVReader");
|
||||
files_ = layout.files;
|
||||
dataset_->experiment.Goniometer(GoniometerAxis(header0_.axis_name,
|
||||
static_cast<float>(angles.front()),
|
||||
static_cast<float>(increment),
|
||||
static_cast<float>(layout.start_deg),
|
||||
static_cast<float>(layout.increment_deg),
|
||||
ASSUMED_BASE_AXIS, {}));
|
||||
|
||||
dataset_->error_value = -1;
|
||||
@@ -160,6 +163,9 @@ CompressedImage JFJochSMVReader::DecodeInto(int64_t image_number, Buffer &buffer
|
||||
if (image_number < 0 || static_cast<size_t>(image_number) >= files_.size())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"Image number out of range");
|
||||
if (files_[image_number].empty())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"No image at this point of the sweep");
|
||||
|
||||
const size_t npixel = static_cast<size_t>(header0_.nx) * static_cast<size_t>(header0_.ny);
|
||||
buffer.resize(npixel * sizeof(int32_t));
|
||||
@@ -184,13 +190,26 @@ bool JFJochSMVReader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &dataset,
|
||||
if (!dataset)
|
||||
return false;
|
||||
|
||||
if (!HasImage(image_number))
|
||||
return false;
|
||||
|
||||
std::vector<uint8_t> scratch;
|
||||
message.image = DecodeInto(image_number, buffer, scratch);
|
||||
message.number = image_number;
|
||||
return true;
|
||||
}
|
||||
|
||||
// A slot the series has no file for is a missing image, not an error: every image loop in the
|
||||
// pipeline already treats "nothing to read" as a frame to pass over, which is exactly what a gap in
|
||||
// a deposited sweep is.
|
||||
bool JFJochSMVReader::HasImage(int64_t image_number) const {
|
||||
return image_number >= 0 && static_cast<size_t>(image_number) < files_.size()
|
||||
&& !files_[image_number].empty();
|
||||
}
|
||||
|
||||
bool JFJochSMVReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) {
|
||||
if (!HasImage(image_number))
|
||||
return false;
|
||||
image.image = DecodeInto(image_number, image.image_buffer, image.read_buffer);
|
||||
return true;
|
||||
}
|
||||
|
||||
@@ -52,6 +52,12 @@ public:
|
||||
[[nodiscard]] uint64_t GetNumberOfImages() const override;
|
||||
void Close() override;
|
||||
|
||||
|
||||
// Whether the sweep really has an image at this point. A deposited series can be missing frames,
|
||||
// and they are kept as gaps so that every image keeps the spindle angle its own header states;
|
||||
// ReadRawImage hands such a slot back as "nothing to read".
|
||||
[[nodiscard]] bool HasImage(int64_t image_number) const;
|
||||
|
||||
bool ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) override;
|
||||
[[nodiscard]] std::vector<SpotToSave> ReadSpots(int64_t image) const override;
|
||||
};
|
||||
|
||||
@@ -0,0 +1,224 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include "SweepLayout.h"
|
||||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <cctype>
|
||||
#include <filesystem>
|
||||
#include <optional>
|
||||
|
||||
#include "../common/JFJochException.h"
|
||||
#include "../common/Logger.h"
|
||||
|
||||
namespace {
|
||||
|
||||
// The difference between two recorded angles, brought into (-180, 180]. A sweep that runs past 360
|
||||
// starts over at 0 in some writers' headers and keeps counting in others', and this reads both the
|
||||
// same way. Safe because no gap seen in a deposited series comes close to half a turn - the widest
|
||||
// in the corpus here is 0.7 degrees - so a folded step is the step and not an aliased one.
|
||||
double Fold(double d) {
|
||||
while (d <= -180.0) d += 360.0;
|
||||
while (d > 180.0) d -= 360.0;
|
||||
return d;
|
||||
}
|
||||
|
||||
std::string Name(const std::string &path) {
|
||||
return std::filesystem::path(path).filename().string();
|
||||
}
|
||||
|
||||
// The frame number a file name ends in. Every one-file-per-image format here numbers its frames that
|
||||
// way, and the numbering is a second, independent statement of how long the series should be: a
|
||||
// directory holding fewer files than its own numbering spans is a directory that was not unpacked
|
||||
// whole, which is worth saying out loud - one staged sweep here was short by 690 frames for months,
|
||||
// and the only sign of it was a resolution nobody could explain.
|
||||
std::optional<int64_t> TrailingNumber(const std::string &path) {
|
||||
const std::string name = Name(path);
|
||||
size_t end = name.find_last_of('.');
|
||||
if (end == std::string::npos)
|
||||
end = name.size();
|
||||
// ".cbf.gz" and the like: step back over as many trailing extensions as there are.
|
||||
while (end > 0 && !std::isdigit(static_cast<unsigned char>(name[end - 1]))) {
|
||||
const size_t dot = name.find_last_of('.', end - 1);
|
||||
if (dot == std::string::npos)
|
||||
return {};
|
||||
end = dot;
|
||||
}
|
||||
size_t begin = end;
|
||||
while (begin > 0 && std::isdigit(static_cast<unsigned char>(name[begin - 1])))
|
||||
begin--;
|
||||
if (begin == end)
|
||||
return {};
|
||||
return std::stoll(name.substr(begin, end - begin));
|
||||
}
|
||||
|
||||
// The first few frames of a list, by name, for a message a user has to act on. All of them would be
|
||||
// hundreds of lines on the series this exists for.
|
||||
std::string NameSome(const std::vector<std::string> &paths) {
|
||||
const size_t show = std::min<size_t>(paths.size(), 5);
|
||||
std::string out;
|
||||
for (size_t i = 0; i < show; i++)
|
||||
out += (i ? ", " : "") + Name(paths[i]);
|
||||
if (paths.size() > show)
|
||||
out += fmt::format(" and {} more", paths.size() - show);
|
||||
return out;
|
||||
}
|
||||
|
||||
// Whether two readings of the same instrument setting are the same reading. Relative, because what
|
||||
// counts as the same distance depends on the distance; the bounds are far wider than a read-back
|
||||
// jitters and far narrower than a real move.
|
||||
bool Same(double a, double b, double rel_tol, double abs_tol) {
|
||||
return std::abs(a - b) <= std::max(abs_tol, rel_tol * std::max(std::abs(a), std::abs(b)));
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
namespace sweep {
|
||||
|
||||
Layout Place(const std::vector<Frame> &frames, const std::string &logger_name) {
|
||||
if (frames.empty())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"No images to place on a sweep");
|
||||
Logger logger(logger_name);
|
||||
|
||||
// ---- one sweep, or several sets of images that happen to share a name?
|
||||
//
|
||||
// Everything below takes the geometry from the first frame and applies it to all of them, so a
|
||||
// series whose headers disagree about where the detector was is not a sweep at all and must not
|
||||
// be silently averaged into one. A folder of screening shots is the case that matters: it fails
|
||||
// deep inside indexing, as a lattice nobody can explain, when it should fail here by name.
|
||||
const Frame &f0 = frames[0];
|
||||
std::vector<std::string> bad_distance, bad_beam, bad_wavelength, bad_increment;
|
||||
for (const auto &f : frames) {
|
||||
if (!Same(f.distance_m, f0.distance_m, 0.005, 1e-6))
|
||||
bad_distance.push_back(f.path);
|
||||
if (std::abs(f.beam_x_px - f0.beam_x_px) > 2.0 || std::abs(f.beam_y_px - f0.beam_y_px) > 2.0)
|
||||
bad_beam.push_back(f.path);
|
||||
if (!Same(f.wavelength_A, f0.wavelength_A, 0.001, 1e-9))
|
||||
bad_wavelength.push_back(f.path);
|
||||
if (!Same(f.increment_deg, f0.increment_deg, 0.01, 1e-6))
|
||||
bad_increment.push_back(f.path);
|
||||
}
|
||||
const auto refuse = [&](const char *what, const std::vector<std::string> &who, double first) {
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
fmt::format("The images named here are not one sweep: {} differ(s) from "
|
||||
"{} ({:g}) in {}. Process the sweeps separately.",
|
||||
what, Name(f0.path), first, NameSome(who)));
|
||||
};
|
||||
if (!bad_distance.empty()) refuse("the detector distance", bad_distance, f0.distance_m);
|
||||
if (!bad_beam.empty()) refuse("the beam centre", bad_beam, f0.beam_x_px);
|
||||
if (!bad_wavelength.empty()) refuse("the wavelength", bad_wavelength, f0.wavelength_A);
|
||||
if (!bad_increment.empty()) refuse("the oscillation width", bad_increment, f0.increment_deg);
|
||||
|
||||
// ---- the angles, unwrapped so a sweep that passes 360 keeps counting
|
||||
std::vector<double> angle(frames.size());
|
||||
angle[0] = f0.angle_deg;
|
||||
for (size_t i = 1; i < frames.size(); i++)
|
||||
angle[i] = angle[i - 1] + Fold(frames[i].angle_deg - angle[i - 1]);
|
||||
|
||||
// ---- one rotation step
|
||||
//
|
||||
// The step is the smallest move between two frames that ARE adjacent in the series, signed with
|
||||
// the way it went. Taking the FIRST pair instead - which is what this code used to do - reads a
|
||||
// gap as the step and compresses the whole sweep by however much is missing. The header's own
|
||||
// Angle_increment is not used for this: it is the oscillation WIDTH, which a series with
|
||||
// overlapping or spaced wedges does not step by.
|
||||
//
|
||||
// A difference smaller than half the oscillation width is not a step but jitter in the recorded
|
||||
// angle: no instrument slices finer than it exposes, so wedges overlapping twofold would be a
|
||||
// read-back wobble, and taking one as the step would spread the sweep over millions of slots.
|
||||
const double too_fine = 0.5 * std::abs(f0.increment_deg);
|
||||
double step = 0;
|
||||
for (size_t i = 1; i < frames.size(); i++) {
|
||||
const double d = Fold(angle[i] - angle[i - 1]);
|
||||
if (std::abs(d) > std::max(too_fine, 1e-6) && (step == 0 || std::abs(d) < std::abs(step)))
|
||||
step = d;
|
||||
}
|
||||
|
||||
Layout out;
|
||||
// A series that never turns: a grid scan, a set of stills, or a single image. There is no sweep
|
||||
// to place anything on, so the files are the slots and the header's nominal increment stands.
|
||||
if (step == 0) {
|
||||
out.files.reserve(frames.size());
|
||||
for (const auto &f : frames)
|
||||
out.files.push_back(f.path);
|
||||
out.start_deg = f0.angle_deg;
|
||||
out.increment_deg = f0.increment_deg;
|
||||
out.present = frames.size();
|
||||
return out;
|
||||
}
|
||||
|
||||
// ---- every frame on that step, or this is not a rotation series
|
||||
std::vector<int64_t> slot(frames.size());
|
||||
std::vector<std::string> off_grid;
|
||||
for (size_t i = 0; i < frames.size(); i++) {
|
||||
const double k = (angle[i] - angle[0]) / step;
|
||||
slot[i] = std::llround(k);
|
||||
if (std::abs(k - static_cast<double>(slot[i])) > 0.25)
|
||||
off_grid.push_back(frames[i].path);
|
||||
}
|
||||
if (!off_grid.empty())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
fmt::format("The images named here do not lie on one rotation series: "
|
||||
"{} starts at {:.4f} deg, the series steps by {:.4f} deg, "
|
||||
"and {} sit(s) off that step. Screening images taken at "
|
||||
"scattered angles are not a sweep.",
|
||||
Name(f0.path), angle[0], step, NameSome(off_grid)));
|
||||
|
||||
// ---- the slots
|
||||
//
|
||||
// Numbered from the frame that comes FIRST on the spindle, which is not always the first file:
|
||||
// slot 0 is where the goniometer's start angle is, and the files were only ever sorted by name.
|
||||
const int64_t first = *std::min_element(slot.begin(), slot.end());
|
||||
const int64_t last = *std::max_element(slot.begin(), slot.end());
|
||||
out.files.assign(static_cast<size_t>(last - first) + 1, std::string());
|
||||
std::vector<std::string> duplicates;
|
||||
for (size_t i = 0; i < frames.size(); i++) {
|
||||
std::string &at = out.files[static_cast<size_t>(slot[i] - first)];
|
||||
if (!at.empty())
|
||||
duplicates.push_back(frames[i].path);
|
||||
at = frames[i].path;
|
||||
}
|
||||
if (!duplicates.empty())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
fmt::format("The images named here repeat an angle already taken by "
|
||||
"another image of the series: {}. Two sweeps of the same "
|
||||
"crystal have to be processed separately.",
|
||||
NameSome(duplicates)));
|
||||
|
||||
out.start_deg = angle[0] + step * static_cast<double>(first);
|
||||
out.increment_deg = step;
|
||||
out.present = frames.size();
|
||||
|
||||
if (out.present < out.files.size())
|
||||
logger.Warning("{} of the {} images the sweep spans are present; the {} missing ones are "
|
||||
"left as gaps, so every image keeps the spindle angle its own header states "
|
||||
"({:.2f} to {:.2f} deg). The merge will be that much less complete.",
|
||||
out.present, out.files.size(), out.files.size() - out.present,
|
||||
out.start_deg,
|
||||
out.start_deg + step * static_cast<double>(out.files.size() - 1));
|
||||
|
||||
// The numbering says the same thing a second way, and says it about the ends of the series too,
|
||||
// which the angles cannot: a sweep missing its first and last frames still spans only the angles
|
||||
// that are there. Where the two disagree with the file count, the directory is short.
|
||||
int64_t lo = 0, hi = 0;
|
||||
bool numbered = true;
|
||||
for (size_t i = 0; i < frames.size() && numbered; i++) {
|
||||
const auto n = TrailingNumber(frames[i].path);
|
||||
if (!n.has_value())
|
||||
numbered = false;
|
||||
else if (i == 0)
|
||||
lo = hi = *n;
|
||||
else
|
||||
lo = std::min(lo, *n), hi = std::max(hi, *n);
|
||||
}
|
||||
if (numbered && hi - lo + 1 > static_cast<int64_t>(out.present))
|
||||
logger.Warning("The file numbering runs {}..{}, which is {} frames, but the directory holds "
|
||||
"{}: {} are not there. If this series should be complete, it was not unpacked "
|
||||
"or copied whole - check the source.",
|
||||
lo, hi, hi - lo + 1, out.present, hi - lo + 1 - static_cast<int64_t>(out.present));
|
||||
return out;
|
||||
}
|
||||
|
||||
} // namespace sweep
|
||||
@@ -0,0 +1,60 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
// Where the frames of a one-file-per-image sweep sit on the spindle.
|
||||
//
|
||||
// A deposited series of CBF / marCCD / SMV files is not always complete: frames go missing between
|
||||
// the beamline and the archive, and the numbers that are left are scattered through the range. Every
|
||||
// one of these formats writes each image's own start angle in its own header, so the sweep is fully
|
||||
// recoverable - but only if the frames are placed at the angles they state rather than laid out
|
||||
// end to end. Reading a gapped series as contiguous compresses the sweep: one deposited 179.8 degree
|
||||
// series of 1108 files out of 1800 came out as 111 degrees, and every frame past the first gap was
|
||||
// analysed at the wrong spindle angle.
|
||||
//
|
||||
// So the sweep is a grid of SLOTS, one per rotation step over the range the headers span, and a frame
|
||||
// occupies the slot its own angle puts it in. A slot with no file is a missing image, not a compressed
|
||||
// out one: the reader hands it back as "nothing to read", which every image loop already handles, and
|
||||
// the goniometer's start + increment * image_number is then the true angle of every image.
|
||||
namespace sweep {
|
||||
|
||||
// What one frame's own header says about where it sits and how the instrument stood while it was
|
||||
// taken. The last four are not used to place the frame - they are what the whole series has to agree
|
||||
// on for it to be one sweep at all.
|
||||
struct Frame {
|
||||
std::string path;
|
||||
double angle_deg = 0;
|
||||
double increment_deg = 0;
|
||||
double distance_m = 0;
|
||||
double beam_x_px = 0;
|
||||
double beam_y_px = 0;
|
||||
double wavelength_A = 0;
|
||||
};
|
||||
|
||||
struct Layout {
|
||||
std::vector<std::string> files; // one entry per slot; empty where the series has no frame
|
||||
double start_deg = 0; // angle of slot 0
|
||||
double increment_deg = 0; // one slot, signed with the direction the sweep turns
|
||||
size_t present = 0; // files that are really there (files.size() - the gaps)
|
||||
};
|
||||
|
||||
// Place the frames on the sweep they describe. frames must be in collection order (the order the
|
||||
// file names sort in) and must not be empty.
|
||||
//
|
||||
// Throws JFJochException, naming the frames, when the series is not one sweep: when the headers
|
||||
// disagree about the instrument, when two frames claim the same angle, or when the angles do not sit
|
||||
// on a single rotation step - a folder of screening images taken at scattered angles is not a sweep,
|
||||
// and silently averaging it into one is how such a set comes out as an indexing failure instead of a
|
||||
// clear refusal.
|
||||
//
|
||||
// A series that does not turn at all - one angle repeated, a grid scan or a set of stills - is left
|
||||
// exactly as it came, one slot per file, with the header's nominal increment.
|
||||
//
|
||||
// logger_name is the reader's own logger, so the gap report names the format the user passed.
|
||||
Layout Place(const std::vector<Frame> &frames, const std::string &logger_name);
|
||||
|
||||
} // namespace sweep
|
||||
@@ -50,6 +50,7 @@ ADD_EXECUTABLE(jfjoch_test
|
||||
ModuleSummationTest.cpp
|
||||
ZMQMetadataSocketTest.cpp
|
||||
JFJochReaderTest.cpp
|
||||
SweepLayoutTest.cpp
|
||||
ShadowFinderTest.cpp
|
||||
ParallelForTest.cpp
|
||||
RugnuxTest.cpp
|
||||
|
||||
@@ -0,0 +1,96 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <catch2/catch_all.hpp>
|
||||
|
||||
#include "../common/JFJochException.h"
|
||||
#include "../reader/SweepLayout.h"
|
||||
|
||||
namespace {
|
||||
|
||||
// A series of frames all taken at the same instrument setting, at the angles given.
|
||||
std::vector<sweep::Frame> Series(const std::vector<double> &angles, double increment) {
|
||||
std::vector<sweep::Frame> out;
|
||||
out.reserve(angles.size());
|
||||
for (size_t i = 0; i < angles.size(); i++)
|
||||
out.push_back({"f" + std::to_string(i) + ".cbf", angles[i], increment, 0.2, 1000, 1000, 1.0});
|
||||
return out;
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
TEST_CASE("SweepLayout_Contiguous") {
|
||||
std::vector<double> angles;
|
||||
for (int i = 0; i < 100; i++)
|
||||
angles.push_back(20.0 + 0.1 * i);
|
||||
const auto l = sweep::Place(Series(angles, 0.1), "test");
|
||||
CHECK(l.files.size() == 100);
|
||||
CHECK(l.present == 100);
|
||||
CHECK(l.start_deg == Catch::Approx(20.0));
|
||||
CHECK(l.increment_deg == Catch::Approx(0.1));
|
||||
for (size_t i = 0; i < l.files.size(); i++)
|
||||
CHECK(l.files[i] == "f" + std::to_string(i) + ".cbf");
|
||||
}
|
||||
|
||||
TEST_CASE("SweepLayout_Gapped") {
|
||||
// Frames 3, 4 and 7 of a ten-frame sweep never made it into the archive. The sweep is still ten
|
||||
// steps wide and every frame keeps its own angle.
|
||||
const auto l = sweep::Place(Series({0.0, 0.1, 0.2, 0.5, 0.6, 0.8, 0.9}, 0.1), "test");
|
||||
REQUIRE(l.files.size() == 10);
|
||||
CHECK(l.present == 7);
|
||||
CHECK(l.increment_deg == Catch::Approx(0.1));
|
||||
CHECK(l.files[2] == "f2.cbf");
|
||||
CHECK(l.files[3].empty());
|
||||
CHECK(l.files[4].empty());
|
||||
CHECK(l.files[5] == "f3.cbf");
|
||||
CHECK(l.files[7].empty());
|
||||
CHECK(l.files[9] == "f6.cbf");
|
||||
}
|
||||
|
||||
TEST_CASE("SweepLayout_PastFullTurn") {
|
||||
// A writer that starts over at 0 rather than counting on past 360.
|
||||
const auto l = sweep::Place(Series({359.7, 359.8, 359.9, 0.0, 0.1}, 0.1), "test");
|
||||
CHECK(l.files.size() == 5);
|
||||
CHECK(l.present == 5);
|
||||
CHECK(l.increment_deg == Catch::Approx(0.1));
|
||||
}
|
||||
|
||||
TEST_CASE("SweepLayout_Reversed") {
|
||||
const auto l = sweep::Place(Series({10.0, 9.5, 9.0, 8.0}, 0.5), "test");
|
||||
REQUIRE(l.files.size() == 5);
|
||||
CHECK(l.increment_deg == Catch::Approx(-0.5));
|
||||
CHECK(l.files[3].empty());
|
||||
CHECK(l.files[4] == "f3.cbf");
|
||||
}
|
||||
|
||||
TEST_CASE("SweepLayout_Stills") {
|
||||
// Nothing turns: the files are the images, and the header's nominal increment stands.
|
||||
const auto l = sweep::Place(Series({45.0, 45.0, 45.0}, 0.1), "test");
|
||||
CHECK(l.files.size() == 3);
|
||||
CHECK(l.present == 3);
|
||||
CHECK(l.start_deg == Catch::Approx(45.0));
|
||||
CHECK(l.increment_deg == Catch::Approx(0.1));
|
||||
}
|
||||
|
||||
TEST_CASE("SweepLayout_ScreeningImagesRefused") {
|
||||
// Five shots at scattered angles are not a sweep, and must not be laid out as one.
|
||||
CHECK_THROWS_AS(sweep::Place(Series({0.0, 90.0, 45.0, 300.0, 270.0}, 0.5), "test"),
|
||||
JFJochException);
|
||||
}
|
||||
|
||||
TEST_CASE("SweepLayout_MovedDetectorRefused") {
|
||||
auto frames = Series({0.0, 0.1, 0.2, 0.3}, 0.1);
|
||||
frames[2].distance_m = 0.3;
|
||||
CHECK_THROWS_AS(sweep::Place(frames, "test"), JFJochException);
|
||||
}
|
||||
|
||||
TEST_CASE("SweepLayout_SecondWavelengthRefused") {
|
||||
auto frames = Series({0.0, 0.1, 0.2, 0.3}, 0.1);
|
||||
frames[3].wavelength_A = 1.9;
|
||||
CHECK_THROWS_AS(sweep::Place(frames, "test"), JFJochException);
|
||||
}
|
||||
|
||||
TEST_CASE("SweepLayout_RepeatedAngleRefused") {
|
||||
// Two sweeps of the same crystal concatenated: the second covers angles the first already has.
|
||||
CHECK_THROWS_AS(sweep::Place(Series({0.0, 0.1, 0.2, 0.1, 0.2}, 0.1), "test"), JFJochException);
|
||||
}
|
||||
Reference in New Issue
Block a user