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); +}