diff --git a/reader/CMakeLists.txt b/reader/CMakeLists.txt index fe0f0d7ec..1a07866cc 100644 --- a/reader/CMakeLists.txt +++ b/reader/CMakeLists.txt @@ -20,6 +20,10 @@ ADD_LIBRARY(JFJochReader STATIC SMV.h JFJochSMVReader.cpp JFJochSMVReader.h + XTC.cpp + XTC.h + JFJochXTCReader.cpp + JFJochXTCReader.h HDF5ImageLocator.cpp JFJochStartMessageReader.cpp JFJochStartMessageReader.h diff --git a/reader/JFJochXTCReader.cpp b/reader/JFJochXTCReader.cpp new file mode 100644 index 000000000..baea3d85d --- /dev/null +++ b/reader/JFJochXTCReader.cpp @@ -0,0 +1,186 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "JFJochXTCReader.h" + +#include +#include +#include + +#include + +#include "../common/JFJochException.h" +#include "../common/Logger.h" + +namespace { + +using Vec = std::array; + +double Dot(const Vec &a, const Vec &b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; } + +Vec Cross(const Vec &a, const Vec &b) { + return {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]}; +} + +Vec ReadVec(const nlohmann::json &j) { return {j.at(0).get(), j.at(1).get(), j.at(2).get()}; } + +// A DIALS laboratory vector in the internal frame. DIALS uses the imgCIF laboratory frame (z from +// the sample to the source, y up), which differs from the internal one by a half turn about x. +Coord DialsToInternal(const Vec &v) { + return {static_cast(v[0]), static_cast(-v[1]), static_cast(-v[2])}; +} + +// The DIALS experiment list beside the streams, or in the directory above them. +std::string FindExpt(const std::string &stream) { + const auto dir = std::filesystem::absolute(stream).parent_path(); + for (const auto &d: {dir, dir.parent_path()}) { + std::error_code ec; + for (const auto &e: std::filesystem::directory_iterator(d, ec)) + if (e.is_regular_file() && e.path().extension() == ".expt") + return e.path().string(); + } + return {}; +} + +} // namespace + +bool JFJochXTCReader::CanRead(const std::string &path) { + return xtc::CanRead(path); +} + +void JFJochXTCReader::ReadGeometry(const std::string &expt) { + std::ifstream in(expt); + const auto j = nlohmann::json::parse(in); + const auto &beam = j.at("beam").at(0); + const auto &panels = j.at("detector").at(0).at("panels"); + if (panels.size() != 1) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + expt + " describes a detector of " + std::to_string(panels.size()) + + " panels; only a single-panel detector is read from XTC"); + const auto &panel = panels.at(0); + + const double wavelength_A = beam.at("wavelength").get(); + // DIALS states the direction from the sample to the source; the beam travels the other way. + const Vec direction = ReadVec(beam.at("direction")); + const Vec b = {-direction[0], -direction[1], -direction[2]}; + const Vec fast = ReadVec(panel.at("fast_axis")), slow = ReadVec(panel.at("slow_axis")); + const Vec origin = ReadVec(panel.at("origin")); // mm, corner of the first pixel + const double pixel_mm = panel.at("pixel_size").at(0).get(); + + // Where the beam meets the detector plane, in pixels along fast and slow, and how far along the + // detector normal the plane is. + const Vec normal = Cross(fast, slow); + const double t = Dot(origin, normal) / Dot(b, normal); + const Vec hit = {t * b[0] - origin[0], t * b[1] - origin[1], t * b[2] - origin[2]}; + if (std::fabs(Dot(b, normal)) < 0.9999) + Logger("XTCReader").Warning("{}: the detector is tilted by {:.2f} deg against the beam; it is read " + "as square-on", expt, std::acos(std::fabs(Dot(b, normal))) * 180.0 / M_PI); + + auto &experiment = dataset_->experiment; + const auto &f0 = frames_.front(); + DetectorSetup detector = DetDECTRIS(f0.width, f0.height, "XTC camera", {}); + detector.PixelSize_um(static_cast(pixel_mm * 1000.0)); + // A CCD: the phosphor converts at the surface, so there is no depth to correct parallax for. + detector.SensorThickness_um(0); + detector.SaturationLimit(65535); + detector.BitDepthImage(32); + const auto orientation = DetectorOrientation::Match(DialsToInternal(fast), DialsToInternal(slow)); + if (!orientation.has_value()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + expt + ": the panel's fast and slow axes are not a square-on layout this reader supports"); + detector.ImageOrientation(orientation.value()); + detector.MinFrameTime(std::chrono::microseconds(0)); + detector.MinCountTime(std::chrono::microseconds(0)); + detector.ReadOutTime(std::chrono::nanoseconds(0)); + experiment.Detector(detector); + + experiment.BeamX_pxl(static_cast(Dot(hit, fast) / pixel_mm)); + experiment.BeamY_pxl(static_cast(Dot(hit, slow) / pixel_mm)); + experiment.DetectorDistance_mm(static_cast(std::fabs(Dot(origin, normal)))); + experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / static_cast(wavelength_A)); + + // The polarization degree is 2f - 1 of the fraction f DIALS states; positive where the plane + // contains the lab x axis (normal along y), as in the SMV reader. + if (beam.contains("polarization_fraction") && beam.contains("polarization_normal")) { + const Coord n = DialsToInternal(ReadVec(beam.at("polarization_normal"))); + const float degree = 2.0f * beam.at("polarization_fraction").get() - 1.0f; + experiment.PolarizationFactor(std::fabs(n.x) > std::fabs(n.y) ? -degree : degree); + } +} + +void JFJochXTCReader::ReadFiles(const std::string &path) { + frames_.clear(); + for (const auto &stream: xtc::CollectStreams(path)) { + auto f = xtc::IndexFrames(stream); + frames_.insert(frames_.end(), f.begin(), f.end()); + } + if (frames_.empty()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "No camera images found in the XTC of " + path); + for (const auto &f: frames_) + if (f.width != frames_.front().width || f.height != frames_.front().height) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + f.file + ": the images are not all the same size"); + + const std::string expt = FindExpt(frames_.front().file); + if (expt.empty()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "XTC carries no detector geometry: put the DIALS experiment list (.expt) of " + "the run beside the XTC files, or in the directory above them"); + Logger("XTCReader").Info("{} images in XTC, geometry from {}", frames_.size(), expt); + + dataset_ = std::make_shared(); + dataset_->experiment = default_experiment; + ReadGeometry(expt); + dataset_->error_value = -1; + dataset_->experiment.ImagesPerTrigger(static_cast(frames_.size())); + dataset_->pixel_mask = std::make_shared(frames_.front().width, frames_.front().height); + + SetStartMessage(dataset_); +} + +uint64_t JFJochXTCReader::GetNumberOfImages() const { + return frames_.size(); +} + +void JFJochXTCReader::Close() { + frames_.clear(); + dataset_.reset(); +} + +template +CompressedImage JFJochXTCReader::DecodeInto(int64_t image_number, Buffer &buffer, + std::vector &scratch) const { + if (image_number < 0 || static_cast(image_number) >= frames_.size()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Image number out of range"); + const auto &f = frames_[image_number]; + buffer.resize(static_cast(f.width) * f.height * sizeof(int32_t)); + xtc::ReadFrame(f, reinterpret_cast(buffer.data()), scratch); + return CompressedImage(buffer.data(), buffer.size(), f.width, f.height, + CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION); +} + +bool JFJochXTCReader::LoadImage_i(std::shared_ptr &dataset, + DataMessage &message, + std::vector &buffer, + int64_t image_number, + bool update_dataset) { + (void) update_dataset; + if (!dataset || image_number < 0 || static_cast(image_number) >= frames_.size()) + return false; + std::vector scratch; + message.image = DecodeInto(image_number, buffer, scratch); + message.number = image_number; + return true; +} + +bool JFJochXTCReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) { + if (image_number < 0 || static_cast(image_number) >= frames_.size()) + return false; + image.image = DecodeInto(image_number, image.image_buffer, image.read_buffer); + return true; +} + +std::vector JFJochXTCReader::ReadSpots(int64_t) const { + return {}; // raw XTC stores no analysis results +} diff --git a/reader/JFJochXTCReader.h b/reader/JFJochXTCReader.h new file mode 100644 index 000000000..6b22eb4ef --- /dev/null +++ b/reader/JFJochXTCReader.h @@ -0,0 +1,53 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "JFJochReader.h" +#include "XTC.h" + +// Reads serial-crystallography stills straight from LCLS-I XTC files: one camera image per event, +// from a single-panel detector such as the Rayonix CCDs at MFX and XPP. A multi-panel detector +// (CSPAD, ePix, JUNGFRAU at LCLS) is stored as raw per-module data with its own calibration, and is +// not read here. +// +// XTC records the events and not the experiment: there is no distance, beam centre or wavelength +// in it that this reader can use. Those come from a DIALS experiment list (.expt) beside the XTC +// files, or in their parent directory - the geometry a deposition ships for cctbx.xfel, and what +// cctbx.xfel itself reads them with. Without one the run is refused rather than guessed. +// +// The images are stills: no goniometer. Each is decoded on demand. The pixels are passed on as +// stored - the CCD's constant offset of a few counts is left in, as it is in an SMV image, and the +// pre-scan finds the defective pixels. +class JFJochXTCReader : public JFJochReader { + std::vector frames_; + std::shared_ptr dataset_; + + bool LoadImage_i(std::shared_ptr &dataset, + DataMessage &message, + std::vector &buffer, + int64_t image_number, + bool update_dataset) override; + + template + CompressedImage DecodeInto(int64_t image_number, Buffer &buffer, std::vector &scratch) const; + + void ReadGeometry(const std::string &expt); + +public: + ~JFJochXTCReader() override = default; + + static bool CanRead(const std::string &path); + + // path is one XTC stream, or a directory whose streams are all read. + void ReadFiles(const std::string &path); + + [[nodiscard]] uint64_t GetNumberOfImages() const override; + void Close() override; + + bool ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) override; + [[nodiscard]] std::vector ReadSpots(int64_t image) const override; +}; diff --git a/reader/XTC.cpp b/reader/XTC.cpp new file mode 100644 index 000000000..03dd0de4c --- /dev/null +++ b/reader/XTC.cpp @@ -0,0 +1,149 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "XTC.h" + +#include +#include +#include +#include +#include + +#include "../common/JFJochException.h" + +namespace { + +constexpr size_t DGRAM_HEADER = 40; // clock time 8, timestamp 8, env 4, Xtc header 20 +constexpr size_t XTC_HEADER = 20; +constexpr size_t FRAME_HEADER = 16; // Camera::FrameV1: width, height, depth, offset +constexpr uint32_t ID_XTC = 1; +constexpr uint32_t ID_FRAME = 2; +constexpr uint32_t TRANSITION_L1ACCEPT = 12; +constexpr uint32_t TRANSITION_LAST = 12; + +struct FileCloser { + void operator()(FILE *f) const { if (f) fclose(f); } +}; +using File = std::unique_ptr; + +File Open(const std::string &path) { + File f(fopen(path.c_str(), "rb")); + if (!f) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Cannot open " + path); + return f; +} + +uint32_t Word(const uint8_t *p, size_t i) { + uint32_t v; + memcpy(&v, p + 4 * i, 4); + return v; +} + +// The transition a datagram header names, from the top bits of its timestamp's first word. +uint32_t Transition(const uint8_t *dgram) { + return (Word(dgram, 2) >> 24) & 0xf; +} + +// A first datagram is a Configure with a container as its outermost Xtc, whose extent covers at +// least its own header. +bool PlausibleFirstDatagram(const uint8_t *h) { + return Transition(h) <= TRANSITION_LAST && (Word(h, 8) & 0xffff) == ID_XTC && Word(h, 9) >= XTC_HEADER; +} + +// Walks the Xtcs in [begin, end) of the file, looking for camera frames; the largest found is kept. +void FindFrame(FILE *f, uint64_t begin, uint64_t end, xtc::Frame &best) { + uint64_t pos = begin; + while (pos + XTC_HEADER <= end) { + uint8_t h[XTC_HEADER]; + if (fseeko(f, static_cast(pos), SEEK_SET) != 0 || fread(h, 1, XTC_HEADER, f) != XTC_HEADER) + return; + const uint32_t type = Word(h, 3) & 0xffff; + const uint32_t extent = Word(h, 4); + if (extent < XTC_HEADER || pos + extent > end) + return; + if (type == ID_XTC) + FindFrame(f, pos + XTC_HEADER, pos + extent, best); + else if (type == ID_FRAME && extent >= XTC_HEADER + FRAME_HEADER) { + uint8_t fh[FRAME_HEADER]; + if (fread(fh, 1, FRAME_HEADER, f) == FRAME_HEADER) { + const uint32_t width = Word(fh, 0), height = Word(fh, 1), depth = Word(fh, 2); + const uint64_t npixel = static_cast(width) * height; + if (depth > 8 && depth <= 16 && XTC_HEADER + FRAME_HEADER + 2 * npixel <= extent + && npixel > static_cast(best.width) * best.height) + best = {"", pos + XTC_HEADER + FRAME_HEADER, width, height}; + } + } + pos += extent; + } +} + +} // namespace + +namespace xtc { + +bool CanRead(const std::string &path) { + std::error_code ec; + if (std::filesystem::is_directory(path, ec)) + return !CollectStreams(path).empty(); + if (!std::filesystem::is_regular_file(path, ec) || !path.ends_with(".xtc") || path.ends_with(".smd.xtc")) + return false; + File f(fopen(path.c_str(), "rb")); + uint8_t h[DGRAM_HEADER]; + return f && fread(h, 1, DGRAM_HEADER, f.get()) == DGRAM_HEADER && PlausibleFirstDatagram(h); +} + +std::vector CollectStreams(const std::string &path) { + std::error_code ec; + if (!std::filesystem::is_directory(path, ec)) + return {path}; + std::vector out; + for (const auto &e: std::filesystem::directory_iterator(path, ec)) { + const std::string p = e.path().string(); + if (e.is_regular_file() && p.ends_with(".xtc") && !p.ends_with(".smd.xtc")) + out.push_back(p); + } + std::sort(out.begin(), out.end()); + return out; +} + +std::vector IndexFrames(const std::string &file) { + File f = Open(file); + std::vector out; + uint64_t pos = 0; + uint8_t h[DGRAM_HEADER]; + while (fseeko(f.get(), static_cast(pos), SEEK_SET) == 0 + && fread(h, 1, DGRAM_HEADER, f.get()) == DGRAM_HEADER) { + const uint32_t extent = Word(h, 9); + if (extent < XTC_HEADER) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + file + ": a datagram's extent is shorter than its own header"); + const uint64_t payload = pos + DGRAM_HEADER, end = pos + DGRAM_HEADER - XTC_HEADER + extent; + if (Transition(h) == TRANSITION_L1ACCEPT) { + Frame frame; + FindFrame(f.get(), payload, end, frame); + if (frame.width > 0) { + frame.file = file; + out.push_back(frame); + } + } + pos = end; + } + return out; +} + +void ReadFrame(const Frame &frame, int32_t *out, std::vector &scratch) { + const size_t npixel = static_cast(frame.width) * frame.height; + scratch.resize(2 * npixel); + File f = Open(frame.file); + if (fseeko(f.get(), static_cast(frame.offset), SEEK_SET) != 0 + || fread(scratch.data(), 1, scratch.size(), f.get()) != scratch.size()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + frame.file + ": an image ends past the end of the file"); + for (size_t i = 0; i < npixel; i++) { + uint16_t v; + memcpy(&v, scratch.data() + 2 * i, 2); + out[i] = v; + } +} + +} // namespace xtc diff --git a/reader/XTC.h b/reader/XTC.h new file mode 100644 index 000000000..8e2dcff6a --- /dev/null +++ b/reader/XTC.h @@ -0,0 +1,50 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +// The LCLS-I XTC container, as far as reading camera frames out of it needs: what the LCLS data +// system (pdsdata) wrote at LCLS until 2023, and what serial-crystallography depositions from that +// era (Rayonix at MFX and XPP) are archived as. +// +// A file is a sequence of datagrams. Each is a 40-byte header - the event's clock time and timestamp +// (16 bytes, whose top bits name the transition: Configure, BeginRun, ..., L1Accept for an event), +// an environment word, and the outermost Xtc header - followed by the Xtc's payload. An Xtc header is +// five 32-bit words: damage, the source (two words), the type and version, and the extent, which +// counts the header itself. A container (type Id_Xtc) holds further Xtcs one after another; anything +// else is a leaf holding one data object. A camera image is a leaf of type Id_Frame (Camera::FrameV1): +// width, height, bit depth and offset, 32 bits each, then the pixels row by row. Everything is +// little-endian and uncompressed. +// +// A run is written as several streams (-s00, -s01, ...) by different recorders, each holding a +// share of the events. Their order across streams is the event time, which this reader does not +// reconstruct: stills are independent, so the streams are read one after another. +namespace xtc { + +// One camera image: where its pixels start in which file. +struct Frame { + std::string file; + uint64_t offset = 0; // of the first pixel + uint32_t width = 0; + uint32_t height = 0; +}; + +// True if the path is an XTC file (or a directory holding one) - by its first datagram, not its name. +bool CanRead(const std::string &path); + +// The event streams of the path: the file itself, or every XTC file in a directory, by name. The +// "smalldata" copies (.smd.xtc) carry the events without their images and are left out. +std::vector CollectStreams(const std::string &path); + +// Every L1Accept of the file that holds a 16-bit camera image, in file order. Where an event holds +// images from more than one camera, the largest is taken - the diffraction detector. +std::vector IndexFrames(const std::string &file); + +// Reads one image as 32-bit pixels into out (width*height values), through scratch. +void ReadFrame(const Frame &frame, int32_t *out, std::vector &scratch); + +} // namespace xtc diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 956ff62a8..09d2d711d 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -29,6 +29,7 @@ #include "../reader/JFJochMarCCDReader.h" #include "../reader/JFJochSMVReader.h" #include "../reader/JFJochCrystFELReader.h" +#include "../reader/JFJochXTCReader.h" #include "../reader/JFJochStartMessageReader.h" #include "BrokerFeed.h" #include "../writer/HDF5NXmx.h" @@ -1730,6 +1731,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { SimulatedLive marccd_reader; SimulatedLive smv_reader; SimulatedLive crystfel_reader; + SimulatedLive xtc_reader; JFJochStartMessageReader live_reader; std::unique_ptr broker_feed; JFJochReader *reader_ptr = nullptr; @@ -1737,11 +1739,14 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { // first plausible member cannot be read throws here - and an input the program cannot read has to // come out as the message below, not as an unhandled exception with nothing printed at all. bool input_is_cbf = false, input_is_marccd = false, input_is_smv = false, input_is_crystfel = false; + bool input_is_xtc = false; try { input_is_crystfel = !crystfel_geometry.empty(); input_is_cbf = !input_is_crystfel && JFJochCBFReader::CanRead(input_file); input_is_marccd = !input_is_cbf && JFJochMarCCDReader::CanRead(input_file); input_is_smv = !input_is_cbf && !input_is_marccd && JFJochSMVReader::CanRead(input_file); + input_is_xtc = !input_is_crystfel && !input_is_cbf && !input_is_marccd && !input_is_smv + && JFJochXTCReader::CanRead(input_file); if (input_is_crystfel) { crystfel_reader.ReadFiles(input_file, crystfel_geometry); reader_ptr = &crystfel_reader; @@ -1762,11 +1767,15 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { reader_ptr = &smv_reader; logger.Info("Read {} SMV images of the sweep containing {}", smv_reader.GetNumberOfImages(), input_file); + } else if (input_is_xtc) { + xtc_reader.ReadFiles(input_file); + reader_ptr = &xtc_reader; + logger.Info("Read {} XTC still images from {}", xtc_reader.GetNumberOfImages(), input_file); } else if (std::filesystem::is_directory(input_file)) { // A directory is never an HDF5 master, so handing it to the HDF5 reader only buries the // real answer - that no one-file-per-image sweep was found in it - under a page of HDF5 // stack output. - logger.Error("{} holds no CBF, marCCD or SMV sweep this program can read (a directory is " + logger.Error("{} holds no CBF, marCCD, SMV or XTC data this program can read (a directory is " "read as the sweep of frames it holds)", input_file); exit(EXIT_FAILURE); } else if (!broker_url.empty() && !std::filesystem::exists(input_file)) { @@ -1936,7 +1945,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { "runs no space-group search. Use --mode mx, which makes the P1 merge the " "ledger is folded from."); // Re-scaling reads reflections a previous run integrated, which only a _process.h5 holds. - if (input_is_cbf || input_is_marccd || input_is_smv || input_is_crystfel) { + if (input_is_cbf || input_is_marccd || input_is_smv || input_is_crystfel || input_is_xtc) { logger.Error("--mode scale needs the integrated reflections in a _process.h5; " "a raw image sweep holds none. Run --mode mx on it first."); exit(EXIT_FAILURE); @@ -3167,7 +3176,8 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { config.detect_beam_stop = detect_beam_stop; config.beam_center_check = beam_center_check; config.adaptive_integration_radius = adaptive_integration_radius; - config.counts_in_adu = input_is_smv || input_is_marccd || (input_is_crystfel && !crystfel_reader.CountsInPhotons()); + config.counts_in_adu = input_is_smv || input_is_marccd || input_is_xtc + || (input_is_crystfel && !crystfel_reader.CountsInPhotons()); config.rotation_postrefine_geometry = rotation_postrefine_geometry; // Measure the spot budget from the data unless the user pinned it. config.measure_spot_budget = !max_spot_count_override.has_value();