From 5f78fc156f555ffc4ce61ec4038fb3137962f7d7 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 20 Sep 2026 10:10:29 +0200 Subject: [PATCH] Place CBF/marCCD/SMV frames on the sweep their own headers state A series of one file per image was laid out end to end: the rotation start came from the first file and the step from the difference between the first two, so a series with frames missing came out compressed - one 179.8 degree deposited sweep of 1108 files out of 1800 was read as 111 degrees, and every frame past the first gap was analysed at the wrong spindle angle. Indexing then found a lattice that took 4% of the validation spots, and two other gapped series aborted outright with "it is not this crystal's lattice". Every one of these formats writes each image's own start angle in its own header, so the sweep is fully recoverable. The new reader/SweepLayout places each frame at the slot its own angle puts it in and leaves a missing frame as a gap - a slot with no file, which ReadRawImage reports as nothing to read, which every image loop in the pipeline already passes over. The goniometer's start + increment * image_number is then the true angle of every image, and the sweep range, the per-10-degree delta-CC1/2 batches and the sweep-quality ledger all read the rotation the headers describe. The rotation step is the smallest move between two frames that really are adjacent, not the first pair. The three readers shared this code by duplication; it is now written once. The same place refuses what is not a sweep rather than averaging it into one: headers that disagree about the detector distance, the beam centre, the wavelength or the oscillation width, angles that do not sit on a single step (a folder of screening shots), or two frames claiming the same angle - each naming the frames. A series that does not turn at all is left exactly as it was. A directory holding fewer files than its own numbering spans is also reported, with both counts: that is the signal that a sweep was not unpacked or copied whole, which otherwise shows up only as a resolution nobody can explain. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT --- reader/CMakeLists.txt | 2 + reader/JFJochCBFReader.cpp | 48 +++++--- reader/JFJochCBFReader.h | 6 + reader/JFJochMarCCDReader.cpp | 49 +++++--- reader/JFJochMarCCDReader.h | 6 + reader/JFJochSMVReader.cpp | 49 +++++--- reader/JFJochSMVReader.h | 6 + reader/SweepLayout.cpp | 224 ++++++++++++++++++++++++++++++++++ reader/SweepLayout.h | 60 +++++++++ tests/CMakeLists.txt | 1 + tests/SweepLayoutTest.cpp | 96 +++++++++++++++ 11 files changed, 502 insertions(+), 45 deletions(-) create mode 100644 reader/SweepLayout.cpp create mode 100644 reader/SweepLayout.h create mode 100644 tests/SweepLayoutTest.cpp diff --git a/reader/CMakeLists.txt b/reader/CMakeLists.txt index 23029e7a7..06f304776 100644 --- a/reader/CMakeLists.txt +++ b/reader/CMakeLists.txt @@ -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 diff --git a/reader/JFJochCBFReader.cpp b/reader/JFJochCBFReader.cpp index f8dbddf88..c4390541c 100644 --- a/reader/JFJochCBFReader.cpp +++ b/reader/JFJochCBFReader.cpp @@ -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::duration(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 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 frames(files_.size()); { const size_t nthreads = std::min(std::max(1u, std::thread::hardware_concurrency()), 8); std::vector> 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(angles.front()), - static_cast(increment), + static_cast(layout.start_deg), + static_cast(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(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(header0_.nx) * static_cast(header0_.ny); buffer.resize(npixel * sizeof(int32_t)); @@ -371,6 +376,9 @@ bool JFJochCBFReader::LoadImage_i(std::shared_ptr &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 scratch; @@ -379,7 +387,17 @@ bool JFJochCBFReader::LoadImage_i(std::shared_ptr &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(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; } diff --git a/reader/JFJochCBFReader.h b/reader/JFJochCBFReader.h index e13555d6a..1ee73ab89 100644 --- a/reader/JFJochCBFReader.h +++ b/reader/JFJochCBFReader.h @@ -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 ReadSpots(int64_t image) const override; }; diff --git a/reader/JFJochMarCCDReader.cpp b/reader/JFJochMarCCDReader.cpp index 5e517cf25..b2821993c 100644 --- a/reader/JFJochMarCCDReader.cpp +++ b/reader/JFJochMarCCDReader.cpp @@ -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 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 frames(files_.size()); { const size_t nthreads = std::min(std::max(1u, std::thread::hardware_concurrency()), 8); std::vector> 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(angles.front()), - static_cast(increment), + static_cast(layout.start_deg), + static_cast(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(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(header0_.nx) * static_cast(header0_.ny); buffer.resize(npixel * sizeof(int32_t)); @@ -183,13 +189,26 @@ bool JFJochMarCCDReader::LoadImage_i(std::shared_ptr &datas if (!dataset) return false; + if (!HasImage(image_number)) + return false; + std::vector 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(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; } diff --git a/reader/JFJochMarCCDReader.h b/reader/JFJochMarCCDReader.h index 536993f81..886009162 100644 --- a/reader/JFJochMarCCDReader.h +++ b/reader/JFJochMarCCDReader.h @@ -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 ReadSpots(int64_t image) const override; }; diff --git a/reader/JFJochSMVReader.cpp b/reader/JFJochSMVReader.cpp index 51e0a5d6c..761ff287d 100644 --- a/reader/JFJochSMVReader.cpp +++ b/reader/JFJochSMVReader.cpp @@ -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 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 frames(files_.size()); { const size_t nthreads = std::min(std::max(1u, std::thread::hardware_concurrency()), 8); std::vector> 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(angles.front()), - static_cast(increment), + static_cast(layout.start_deg), + static_cast(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(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(header0_.nx) * static_cast(header0_.ny); buffer.resize(npixel * sizeof(int32_t)); @@ -184,13 +190,26 @@ bool JFJochSMVReader::LoadImage_i(std::shared_ptr &dataset, if (!dataset) return false; + if (!HasImage(image_number)) + return false; + std::vector 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(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; } diff --git a/reader/JFJochSMVReader.h b/reader/JFJochSMVReader.h index 02f12fb55..ec67a8187 100644 --- a/reader/JFJochSMVReader.h +++ b/reader/JFJochSMVReader.h @@ -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 ReadSpots(int64_t image) const override; }; diff --git a/reader/SweepLayout.cpp b/reader/SweepLayout.cpp new file mode 100644 index 000000000..8127f7741 --- /dev/null +++ b/reader/SweepLayout.cpp @@ -0,0 +1,224 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "SweepLayout.h" + +#include +#include +#include +#include +#include + +#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 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(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(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 &paths) { + const size_t show = std::min(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 &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 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 &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 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 slot(frames.size()); + std::vector 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(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(last - first) + 1, std::string()); + std::vector duplicates; + for (size_t i = 0; i < frames.size(); i++) { + std::string &at = out.files[static_cast(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(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(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(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(out.present)); + return out; +} + +} // namespace sweep diff --git a/reader/SweepLayout.h b/reader/SweepLayout.h new file mode 100644 index 000000000..31767f1c6 --- /dev/null +++ b/reader/SweepLayout.h @@ -0,0 +1,60 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +// 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 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 &frames, const std::string &logger_name); + +} // namespace sweep diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1b9ef6040..fa5452b7e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -50,6 +50,7 @@ ADD_EXECUTABLE(jfjoch_test ModuleSummationTest.cpp ZMQMetadataSocketTest.cpp JFJochReaderTest.cpp + SweepLayoutTest.cpp ShadowFinderTest.cpp ParallelForTest.cpp RugnuxTest.cpp diff --git a/tests/SweepLayoutTest.cpp b/tests/SweepLayoutTest.cpp new file mode 100644 index 000000000..9655afe01 --- /dev/null +++ b/tests/SweepLayoutTest.cpp @@ -0,0 +1,96 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#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 Series(const std::vector &angles, double increment) { + std::vector 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 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); +}