diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 4b94d76f8..b1e31d853 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,10 @@ # Changelog ## 1.0.0 +### 1.0.0-rc.171 + +* rugnux reads marCCD images natively: naming any frame of a Rayonix or mar Mosaic sweep processes the whole sweep, with no conversion step. + ### 1.0.0-rc.170 * Fixed a `jfjoch_broker` crash during indexing: sorting no longer misbehaves on non-finite values, and GPU FFT indexer kernel launches are now error-checked. diff --git a/docs/RUGNUX_FORMATS.md b/docs/RUGNUX_FORMATS.md index 11f3309b0..6102d3e8d 100644 --- a/docs/RUGNUX_FORMATS.md +++ b/docs/RUGNUX_FORMATS.md @@ -1,10 +1,10 @@ # What rugnux reads Which data `rugnux` opens, before anything is typed. The short answer: an HDF5 master (NXmx or -DECTRIS, from any facility) or a PILATUS miniCBF sweep — nothing else is read, so any other format -has to be converted to one of these two first. +DECTRIS, from any facility), a PILATUS miniCBF sweep or a marCCD sweep — nothing else is read, so +any other format has to be converted to one of these three first. -**Input** is either an HDF5 master file or a directory of PILATUS miniCBF frames. One input is +**Input** is an HDF5 master file, or a directory of PILATUS miniCBF or marCCD frames. One input is **one sweep of one crystal** — rugnux does not combine sweeps or crystals in a run; process each sweep to its own `_unmerged.mtz` and merge them downstream (see [Taking the data onward](RUGNUX_INTEGRATION.md#taking-the-data-onward)). @@ -26,6 +26,16 @@ sweep to its own `_unmerged.mtz` and merge them downstream header, including the imgCIF axis table where the header carries one (see [Detector geometry](DETECTOR_GEOMETRY.md)). A raw CBF carries no analysis results, so `--mode scale` — which re-scales the reflections stored in a `_process.h5` — does not accept one. +* **marCCD sweep** — what Rayonix MX-series and mar Mosaic detectors write, and what a decade of + deposited CCD data is archived as: an uncompressed TIFF with the instrument header in the gap + before the pixels. One frame per file, read natively. Naming a frame or its directory selects the + sweep exactly as for miniCBF, under either naming scheme these detectors use — a numbered stem + (`xtal_1_00042.mccd`) or the frame number as the file extension (`D1.042`). The distance, beam + centre, pixel size, wavelength and the circle that turned come from the header; the rotation + axis' direction does not, because the format has nowhere to state it, so the run settles its sign + from the data as it does for a miniCBF that carries no axis table. A CCD frame marks no untrusted + pixels, so the sweep starts with nothing masked. Like a raw CBF, it carries no analysis results + and `--mode scale` does not accept one. Spots are always found by `rugnux` itself, including for the two-pass rotation first pass — the spot lists a dataset may already carry were found online, at the acquisition's threshold and with diff --git a/reader/CMakeLists.txt b/reader/CMakeLists.txt index 5fb20d52a..3f336aae5 100644 --- a/reader/CMakeLists.txt +++ b/reader/CMakeLists.txt @@ -6,6 +6,10 @@ ADD_LIBRARY(JFJochReader STATIC MiniCBF.h JFJochCBFReader.cpp JFJochCBFReader.h + MarCCD.cpp + MarCCD.h + JFJochMarCCDReader.cpp + JFJochMarCCDReader.h HDF5ImageLocator.cpp HDF5ImageLocator.h HDF5ImageSource.cpp @@ -22,6 +26,8 @@ ADD_LIBRARY(JFJochReader STATIC # NB: JFJochHttpReader (the viewer's live-broker HTTP client) lives in viewer/ and links libcurl; # it is deliberately NOT part of this always-built library, so the broker/writer never pull in a # TLS/Kerberos stack. +# tiff: the marCCD reader's files are ordinary uncompressed TIFFs, so libtiff - which the build +# already fetches for JFJochPreview in every mode - reads the pixels and there is no new dependency. TARGET_LINK_LIBRARIES(JFJochReader JFJochImageAnalysis JFJochAPI JFJochCommon JFJochZMQ JFJochLogger - JFJochHDF5Wrappers CBORStream2FrameSerialize + JFJochHDF5Wrappers CBORStream2FrameSerialize tiff ${CMAKE_DL_LIBS}) diff --git a/reader/JFJochMarCCDReader.cpp b/reader/JFJochMarCCDReader.cpp new file mode 100644 index 000000000..5e3cdc902 --- /dev/null +++ b/reader/JFJochMarCCDReader.cpp @@ -0,0 +1,190 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "JFJochMarCCDReader.h" + +#include +#include +#include + +#include "../common/JFJochException.h" +#include "../common/JFJochMath.h" +#include "../common/Logger.h" + +namespace { + +// The base rotation axis in the internal frame (x along increasing detector column, y along +// increasing row, z along the beam). A marCCD header names the circle that turned but never states +// a direction for it, so this is the convention an NXmx master writes for the same instruments, and +// a file that needs the other sign is settled from the data by the run's axis-sign rescue - the +// same arrangement JFJochCBFReader makes for a miniCBF that states no axis table. +const Coord ASSUMED_BASE_AXIS(-1.0f, 0.0f, 0.0f); + +} // namespace + +bool JFJochMarCCDReader::CanRead(const std::string &path) { + return marccd::CanRead(path); +} + +void JFJochMarCCDReader::ReadFiles(const std::string &path) { + files_ = marccd::CollectSweep(path); + if (files_.empty()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "No marCCD images found for " + path); + + header0_ = marccd::ReadHeader(files_[0]); + // A pixel size is what makes the file an IMAGE: every resolution, every scattering vector and + // the beam centre in millimetres scale by it, and a default of 0 collapses all of them without + // a word. + if (!(header0_.pixel_x_m > 0.0) || !(header0_.pixel_y_m > 0.0)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + files_[0] + " states no pixel size in its marCCD header"); + if (!(header0_.wavelength_A > 0.0)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + files_[0] + " states no wavelength in its marCCD header"); + if (!(header0_.distance_m > 0.0)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + files_[0] + " states no detector distance in its marCCD header"); + + dataset_ = std::make_shared(); + dataset_->experiment = default_experiment; + + DetectorSetup detector = DetDECTRIS(header0_.nx, header0_.ny, + header0_.detector.empty() ? "marCCD" : header0_.detector, {}); + // Not rounded to whole micrometres, as the miniCBF path can afford to be: a PILATUS pixel is + // exactly 172 um, but these are 73.242 um, and rounding that to 73 is a 0.33% scale error on + // every cell edge the run reports. + detector.PixelSize_um(static_cast(header0_.pixel_x_m * 1e6)); + // A CCD has no sensor thickness worth correcting for: the phosphor converts at the surface and + // the fibre optic carries light, not X-rays, so the parallax correction a silicon sensor needs + // does not apply. Left at zero, which is what the geometry means by "no depth". + detector.SensorThickness_um(0); + if (header0_.saturated_value > 0) + detector.SaturationLimit(SaturationLimitFromValue(header0_.saturated_value)); + else + Logger("MarCCDReader").Warning("{} states no saturated value, so no pixel will be called " + "saturated; if this detector overloads, its strongest " + "reflections will be integrated as if they were valid.", + files_[0]); + // Images are handed out as signed 32-bit whatever the file stored, so that is the depth the + // rest of the code must see; the real overflow is the header's saturated value, set above. + detector.BitDepthImage(32); + detector.MinFrameTime(std::chrono::microseconds(0)); + detector.MinCountTime(std::chrono::microseconds(0)); + detector.ReadOutTime(std::chrono::nanoseconds(0)); + dataset_->experiment.Detector(detector); + + dataset_->experiment.BeamX_pxl(static_cast(header0_.beam_x_px)); + dataset_->experiment.BeamY_pxl(static_cast(header0_.beam_y_px)); + dataset_->experiment.DetectorDistance_mm(static_cast(header0_.distance_m * 1000.0)); + + // A detector swung out on a 2theta arm. The arm turns the detector about the sample and so + // carries the square-on geometry with it: the header's distance stays the distance along the + // detector normal and the beam centre stays the point of normal incidence, which is exactly + // what the PONI convention wants, so the swing is a PONI rotation and nothing else changes. + // The arm turns about the same axis as the spindle on the geometries these headers describe. + if (header0_.two_theta_deg != 0.0) { + float rot1 = 0, rot2 = 0, rot3 = 0; + PoniAnglesFromMatrix(RotMatrix(static_cast(header0_.two_theta_deg * PI / 180.0), + ASSUMED_BASE_AXIS), rot1, rot2, rot3); + dataset_->experiment.PoniRot1_rad(rot1).PoniRot2_rad(rot2).PoniRot3_rad(rot3); + } + + dataset_->experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / static_cast(header0_.wavelength_A)); + dataset_->experiment.FrameTime( + std::chrono::duration_cast( + std::chrono::duration(header0_.exposure_s)), + std::chrono::duration_cast( + std::chrono::duration(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()); + { + 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 (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; + } + dataset_->experiment.Goniometer(GoniometerAxis(header0_.axis_name, + static_cast(angles.front()), + static_cast(increment), + ASSUMED_BASE_AXIS, {})); + + dataset_->error_value = -1; + dataset_->experiment.ImagesPerTrigger(static_cast(files_.size())); + // A CCD frame stores no untrusted-pixel marker - every value is a real reading, and the + // detector has no module gaps - so the sweep starts with nothing masked. + dataset_->pixel_mask = std::make_shared(static_cast(header0_.nx), + static_cast(header0_.ny)); + + SetStartMessage(dataset_); +} + +uint64_t JFJochMarCCDReader::GetNumberOfImages() const { + return files_.size(); +} + +void JFJochMarCCDReader::Close() { + files_.clear(); + dataset_.reset(); +} + +template +CompressedImage JFJochMarCCDReader::DecodeInto(int64_t image_number, Buffer &buffer, + std::vector &scratch) const { + if (image_number < 0 || static_cast(image_number) >= files_.size()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Image number out of range"); + + const size_t npixel = static_cast(header0_.nx) * static_cast(header0_.ny); + buffer.resize(npixel * sizeof(int32_t)); + + const auto h = marccd::ReadInto(files_[image_number], + reinterpret_cast(buffer.data()), npixel, scratch); + if (h.nx != header0_.nx || h.ny != header0_.ny) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "marCCD image size differs from the first image of the sweep"); + + return CompressedImage(buffer.data(), buffer.size(), + static_cast(header0_.nx), static_cast(header0_.ny), + CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION); +} + +bool JFJochMarCCDReader::LoadImage_i(std::shared_ptr &dataset, + DataMessage &message, + std::vector &buffer, + int64_t image_number, + bool update_dataset) { + (void) update_dataset; + if (!dataset) + return false; + + std::vector scratch; + message.image = DecodeInto(image_number, buffer, scratch); + message.number = image_number; + return true; +} + +bool JFJochMarCCDReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) { + image.image = DecodeInto(image_number, image.image_buffer, image.read_buffer); + return true; +} + +std::vector JFJochMarCCDReader::ReadSpots(int64_t) const { + return {}; // a raw marCCD file stores no analysis results +} diff --git a/reader/JFJochMarCCDReader.h b/reader/JFJochMarCCDReader.h new file mode 100644 index 000000000..536993f81 --- /dev/null +++ b/reader/JFJochMarCCDReader.h @@ -0,0 +1,57 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "JFJochReader.h" +#include "MarCCD.h" + +// Reads a rotation sweep straight from a directory of marCCD files - what Rayonix MX-series and mar +// Mosaic detectors write, and what a decade of deposited CCD data is archived as - with no +// conversion step. +// +// The same shape as JFJochCBFReader: the sweep's geometry comes from the first file's header, the +// rotation angle of each image from its own, and images are decoded on demand so any number of +// workers can read at once. A raw marCCD file carries no analysis results, so there are no spots +// and no reflections here either. +// +// The one thing a CCD needs that a pixel detector does not: its pixels are unsigned 16-bit with no +// negative marker for untrusted ones, so the mask is what the header calls saturated and nothing +// else. A CCD's point spread and read-out noise are left to the analysis, which measures the +// background from the image rather than assuming a detector. +class JFJochMarCCDReader : public JFJochReader { + std::vector files_; + std::shared_ptr dataset_; + marccd::Header header0_; + + bool LoadImage_i(std::shared_ptr &dataset, + DataMessage &message, + std::vector &buffer, + int64_t image_number, + bool update_dataset) override; + + // Decodes one image into the caller's byte buffer and returns the image that points at it. The + // stored 16-bit strip is read through scratch, which the caller keeps between frames. + template + CompressedImage DecodeInto(int64_t image_number, Buffer &buffer, std::vector &scratch) const; + +public: + ~JFJochMarCCDReader() override = default; + + // True if the path names something this reader can open: a marCCD file, or a directory holding + // at least one. Cheap - it reads a few kB at most. + static bool CanRead(const std::string &path); + + // path is a directory of marCCD frames, or one frame inside the sweep to take the whole sweep + // from. + 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/MarCCD.cpp b/reader/MarCCD.cpp new file mode 100644 index 000000000..825a50f54 --- /dev/null +++ b/reader/MarCCD.cpp @@ -0,0 +1,322 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "MarCCD.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../common/JFJochException.h" + +namespace marccd { + +namespace { + +// The instrument header, as marCCD writes it: 3072 bytes of little-endian int32 at fixed offsets. +// Only the fields below are read; the rest are detector housekeeping. Offsets are from the start of +// the block, which is where the "MMX" name sits. +constexpr size_t HEADER_BYTES = 3072; +constexpr size_t OFF_NAME = 4; // char[16], "MMX" +constexpr size_t OFF_NFAST = 80; +constexpr size_t OFF_NSLOW = 84; +constexpr size_t OFF_DEPTH = 88; // bytes per pixel +constexpr size_t OFF_SATURATED = 104; +// Goniostat block: distance and beam centre, then the driven circles as start/end pairs. +constexpr size_t OFF_XTAL_TO_DETECTOR = 640; // micrometres +constexpr size_t OFF_BEAM_X = 644; // 1/1000 pixel, along fast +constexpr size_t OFF_BEAM_Y = 648; // 1/1000 pixel, along slow +constexpr size_t OFF_EXPOSURE_TIME = 656; // milliseconds +constexpr size_t OFF_START_TWOTHETA = 668; // millidegrees, and the seven that follow it +constexpr size_t OFF_END_TWOTHETA = 700; // the same seven, at the end of the exposure +constexpr size_t OFF_ROTATION_RANGE = 736; // millidegrees +constexpr size_t OFF_PIXELSIZE_X = 772; // nanometres +constexpr size_t OFF_PIXELSIZE_Y = 776; +constexpr size_t OFF_SOURCE_WAVELENGTH = 908; // 1e-5 angstrom + +// The circles of the goniostat block, in the order it stores them. The scanned one is whichever +// pair of start/end angles differs, so this table is what turns that into a name. +constexpr const char *CIRCLE_NAMES[] = {"two_theta", "omega", "chi", "kappa", + "phi", "delta", "gamma"}; +constexpr size_t N_CIRCLES = sizeof(CIRCLE_NAMES) / sizeof(CIRCLE_NAMES[0]); + +int32_t I32(const std::vector &h, size_t off) { + int32_t v; + std::memcpy(&v, h.data() + off, sizeof(v)); + return v; // the header declares little-endian byte order and no writer has ever used the other +} + +// The first bytes of a file: the TIFF header, the instrument header, and nothing else. +std::vector ReadPrefix(const std::string &path, size_t bytes) { + std::ifstream f(path, std::ios::binary); + if (!f) + throw JFJochException(JFJochExceptionCategory::MockFileOpenError, "Cannot open " + path); + std::vector out(bytes); + f.read(reinterpret_cast(out.data()), static_cast(bytes)); + out.resize(static_cast(f.gcount())); + return out; +} + +// Where the instrument header begins. The file states it in a private TIFF tag; every writer seen +// puts it at 1024 and says so, but the tag is what the format defines, so it wins where present. +uint32_t HeaderOffset(const std::vector &prefix) { + if (prefix.size() < 8) + return DEFAULT_HEADER_OFFSET; + const bool little = prefix[0] == 'I' && prefix[1] == 'I'; + const auto u16 = [&](size_t o) -> uint32_t { + if (o + 2 > prefix.size()) return 0; + return little ? (uint32_t(prefix[o]) | uint32_t(prefix[o + 1]) << 8) + : (uint32_t(prefix[o]) << 8 | uint32_t(prefix[o + 1])); + }; + const auto u32 = [&](size_t o) -> uint32_t { + if (o + 4 > prefix.size()) return 0; + return little ? (uint32_t(prefix[o]) | uint32_t(prefix[o + 1]) << 8 | + uint32_t(prefix[o + 2]) << 16 | uint32_t(prefix[o + 3]) << 24) + : (uint32_t(prefix[o]) << 24 | uint32_t(prefix[o + 1]) << 16 | + uint32_t(prefix[o + 2]) << 8 | uint32_t(prefix[o + 3])); + }; + const uint32_t ifd = u32(4); + const uint32_t n = u16(ifd); + for (uint32_t i = 0; i < n; i++) { + const size_t e = ifd + 2 + static_cast(i) * 12; + if (e + 12 > prefix.size()) + break; + if (u16(e) == TIFFTAG_MARCCD_HEADER_OFFSET) + return u32(e + 8); + } + return DEFAULT_HEADER_OFFSET; +} + +// The instrument header block, or an empty vector where the file does not carry one. +std::vector InstrumentHeader(const std::string &path) { + const uint32_t off = HeaderOffset(ReadPrefix(path, DEFAULT_HEADER_OFFSET)); + auto prefix = ReadPrefix(path, off + HEADER_BYTES); + if (prefix.size() < off + HEADER_BYTES) + return {}; + std::vector h(prefix.begin() + off, prefix.begin() + off + HEADER_BYTES); + if (std::memcmp(h.data() + OFF_NAME, "MMX", 4) != 0) + return {}; + return h; +} + +Header Parse(const std::vector &h, const std::string &path) { + Header out; + out.nx = I32(h, OFF_NFAST); + out.ny = I32(h, OFF_NSLOW); + if (out.nx <= 0 || out.ny <= 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + path + ": marCCD header states no image size"); + if (I32(h, OFF_DEPTH) != 2) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + path + ": only 16-bit marCCD images are supported"); + + out.pixel_x_m = I32(h, OFF_PIXELSIZE_X) * 1e-9; + out.pixel_y_m = I32(h, OFF_PIXELSIZE_Y) * 1e-9; + out.distance_m = I32(h, OFF_XTAL_TO_DETECTOR) * 1e-6; + // The beam centre is stored in thousandths of a pixel, x along the fast direction and y along + // the slow one, which is the same order the internal frame counts columns and rows in. + out.beam_x_px = I32(h, OFF_BEAM_X) * 1e-3; + out.beam_y_px = I32(h, OFF_BEAM_Y) * 1e-3; + out.wavelength_A = I32(h, OFF_SOURCE_WAVELENGTH) * 1e-5; + out.two_theta_deg = I32(h, OFF_START_TWOTHETA) * 1e-3; + out.exposure_s = I32(h, OFF_EXPOSURE_TIME) * 1e-3; + out.saturated_value = I32(h, OFF_SATURATED); + + // Which circle moved. The header also carries a rotation_axis index, but the start/end pair + // that differs is the same answer read off the angles themselves, and it stays right on a file + // whose index field was never filled in. + out.angle_increment_deg = I32(h, OFF_ROTATION_RANGE) * 1e-3; + for (size_t i = 0; i < N_CIRCLES; i++) { + const double start = I32(h, OFF_START_TWOTHETA + 4 * i) * 1e-3; + const double end = I32(h, OFF_END_TWOTHETA + 4 * i) * 1e-3; + if (start != end) { + out.axis_name = CIRCLE_NAMES[i]; + out.start_angle_deg = start; + if (out.angle_increment_deg == 0.0) + out.angle_increment_deg = end - start; + break; + } + } + return out; +} + +// The extensions a marCCD frame is written under. A frame whose extension is the frame number +// (D1.001) is named by the mar software itself and is just as common as a named one. +bool PlausibleExtension(const std::filesystem::path &p) { + std::string ext = p.extension().string(); + if (ext.empty()) + return false; + ext.erase(ext.begin()); // drop the dot + if (std::all_of(ext.begin(), ext.end(), [](unsigned char c) { return std::isdigit(c); })) + return true; + std::transform(ext.begin(), ext.end(), ext.begin(), + [](unsigned char c) { return std::tolower(c); }); + return ext == "mccd" || ext == "img" || ext == "tif" || ext == "tiff" || ext == "marccd"; +} + +// The sweep a file belongs to, as a template: everything before the LAST run of digits in the whole +// file name, the length of that run, and whatever follows it. That one rule covers both naming +// schemes in use - "xtal_1_00042.mccd" -> {"xtal_1_", 5, ".mccd"} and "D1.042" -> {"D1.", 3, ""} - +// because in both the frame number is the last number in the name. +struct Template { + std::string prefix; + size_t digits = 0; + std::string suffix; + + bool operator<(const Template &o) const { + return std::tie(prefix, digits, suffix) < std::tie(o.prefix, o.digits, o.suffix); + } + bool Matches(const std::string &name) const { + if (name.size() != prefix.size() + digits + suffix.size()) + return false; + if (name.compare(0, prefix.size(), prefix) != 0) + return false; + if (name.compare(prefix.size() + digits, suffix.size(), suffix) != 0) + return false; + for (size_t i = 0; i < digits; i++) + if (!std::isdigit(static_cast(name[prefix.size() + i]))) + return false; + return true; + } +}; + +std::optional