diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index e91172f1d..dcd7f7ddd 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -130,6 +130,9 @@ ADD_LIBRARY(JFJochCommon STATIC print_license.h DetectorGeometryModular.cpp DetectorGeometryModular.h + DetectorModules.cpp + DetectorModules.h + ModuleShift.h DetectorGeometryFixed.cpp DetectorGeometryFixed.h GridScanSettings.cpp diff --git a/common/DetectorModules.cpp b/common/DetectorModules.cpp new file mode 100644 index 000000000..9e06edbb1 --- /dev/null +++ b/common/DetectorModules.cpp @@ -0,0 +1,154 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "DetectorModules.h" +#include "DetectorGeometry.h" +#include "Definitions.h" + +#include +#include +#include + +DetectorModule DetectorModule::AsImaged(int64_t x0, int64_t y0, int64_t width, int64_t height) { + DetectorModule m; + m.x0 = x0; + m.y0 = y0; + m.width = width; + m.height = height; + m.ox = static_cast(x0); + m.oy = static_cast(y0); + return m; +} + +bool DetectorModule::ContainsImagePixel(double x, double y) const { + // A pixel's centre is at its index, so the module covers [x0 - 0.5, x0 + width - 0.5). + return x >= static_cast(x0) - 0.5 && x < static_cast(x0 + width) - 0.5 + && y >= static_cast(y0) - 0.5 && y < static_cast(y0 + height) - 0.5; +} + +std::pair DetectorModule::ImageToPlane(double x, double y) const { + const double u = x - static_cast(x0), v = y - static_cast(y0); + return {ox + u * fx + v * sx, oy + u * fy + v * sy}; +} + +std::pair DetectorModule::PlaneToImage(double px, double py) const { + // Inverse of the 2x2 [f s] + const double det = fx * sy - sx * fy; + const double dx = px - ox, dy = py - oy; + const double u = (sy * dx - sx * dy) / det; + const double v = (-fy * dx + fx * dy) / det; + return {static_cast(x0) + u, static_cast(y0) + v}; +} + +std::pair DetectorModule::PlaneCentre() const { + return ImageToPlane(static_cast(x0) + 0.5 * static_cast(width - 1), + static_cast(y0) + 0.5 * static_cast(height - 1)); +} + +void DetectorModule::Apply(const ModuleShift &shift) { + const double c = std::cos(shift.rot_rad), s = std::sin(shift.rot_rad); + const double ux = ox - shift.x0, uy = oy - shift.y0; + ox = shift.x0 + c * ux - s * uy + shift.dx; + oy = shift.y0 + s * ux + c * uy + shift.dy; + const double nfx = c * fx - s * fy, nfy = s * fx + c * fy; + const double nsx = c * sx - s * sy, nsy = s * sx + c * sy; + fx = nfx; fy = nfy; + sx = nsx; sy = nsy; +} + +DetectorModules::DetectorModules(std::vector in) : modules(std::move(in)) {} + +bool DetectorModules::empty() const { + return modules.empty(); +} + +size_t DetectorModules::size() const { + return modules.size(); +} + +const DetectorModule &DetectorModules::operator[](size_t i) const { + return modules.at(i); +} + +const std::vector &DetectorModules::Get() const { + return modules; +} + +bool DetectorModules::Moved() const { + for (const auto &m: modules) + if (m.ox != static_cast(m.x0) || m.oy != static_cast(m.y0) + || m.fx != 1.0 || m.fy != 0.0 || m.sx != 0.0 || m.sy != 1.0) + return true; + return false; +} + +int64_t DetectorModules::ModuleOfImagePixel(double x, double y) const { + for (size_t i = 0; i < modules.size(); i++) + if (modules[i].ContainsImagePixel(x, y)) + return static_cast(i); + return -1; +} + +int64_t DetectorModules::ModuleOfPlanePoint(double px, double py) const { + for (size_t i = 0; i < modules.size(); i++) { + const auto [x, y] = modules[i].PlaneToImage(px, py); + if (modules[i].ContainsImagePixel(x, y)) + return static_cast(i); + } + return -1; +} + +std::pair DetectorModules::ImageToPlane(double x, double y) const { + const int64_t m = ModuleOfImagePixel(x, y); + if (m < 0) + return {x, y}; + return modules[m].ImageToPlane(x, y); +} + +std::pair DetectorModules::PlaneToImage(double px, double py) const { + const int64_t m = ModuleOfPlanePoint(px, py); + if (m < 0) + return {px, py}; + return modules[m].PlaneToImage(px, py); +} + +void DetectorModules::Apply(const std::vector &shifts) { + for (auto &m: modules) { + const auto [cx, cy] = m.PlaneCentre(); + const double reach = 0.5 * static_cast(std::min(m.width, m.height)); + const ModuleShift *best = nullptr; + for (const auto &s: shifts) + if (std::hypot(s.x0 - cx, s.y0 - cy) < reach + && (!best || std::hypot(s.x0 - cx, s.y0 - cy) < std::hypot(best->x0 - cx, best->y0 - cy))) + best = &s; + if (best) + m.Apply(*best); + } +} + +std::vector DetectorModules::Centres() const { + std::vector out; + for (const auto &m: modules) { + ModuleShift s; + std::tie(s.x0, s.y0) = m.PlaneCentre(); + out.push_back(s); + } + return out; +} + +DetectorModules ModulesAsImaged(const DetectorGeometry &geom) { + if (geom.GetModulesNum() < 2) + return {}; + std::vector out; + for (int64_t m = 0; m < geom.GetModulesNum(); m++) { + const Coord f = geom.GetFastDirection(m), s = geom.GetSlowDirection(m); + const double cols = CONVERTED_MODULE_COLS - 1, lines = CONVERTED_MODULE_LINES - 1; + const double xa = static_cast(geom.GetX0(m)), ya = static_cast(geom.GetY0(m)); + const double xb = xa + cols * f.x + lines * s.x, yb = ya + cols * f.y + lines * s.y; + const auto x0 = static_cast(std::lround(std::min(xa, xb))); + const auto y0 = static_cast(std::lround(std::min(ya, yb))); + out.push_back(DetectorModule::AsImaged(x0, y0, std::lround(std::abs(xb - xa)) + 1, + std::lround(std::abs(yb - ya)) + 1)); + } + return DetectorModules(out); +} diff --git a/common/DetectorModules.h b/common/DetectorModules.h new file mode 100644 index 000000000..9b3e23d78 --- /dev/null +++ b/common/DetectorModules.h @@ -0,0 +1,75 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +#include "ModuleShift.h" + +class DetectorGeometry; + +// One detector module: where its pixels are in the image rugnux reads, and where they really are on the +// detector plane. The plane is the one DiffractionGeometry describes (beam centre, distance, tilt), in its +// pixel units, so a module placed exactly as the image has it maps every pixel onto itself. +// +// Image pixel (x, y) of the module, x0 <= x < x0 + width, y0 <= y < y0 + height, lies on the plane at +// (ox, oy) + (x - x0) (fx, fy) + (y - y0) (sx, sy) +// - the NXdetector_module description (data_origin, data_size, module_offset, fast_pixel_direction, +// slow_pixel_direction) taken in the plane. A pixel's centre is at its index, as everywhere in the image. +struct DetectorModule { + int64_t x0 = 0, y0 = 0, width = 0, height = 0; + double ox = 0, oy = 0; + double fx = 1, fy = 0; + double sx = 0, sy = 1; + + // The module where the image has it: every pixel maps onto itself. + static DetectorModule AsImaged(int64_t x0, int64_t y0, int64_t width, int64_t height); + + [[nodiscard]] bool ContainsImagePixel(double x, double y) const; + [[nodiscard]] std::pair ImageToPlane(double x, double y) const; + [[nodiscard]] std::pair PlaneToImage(double px, double py) const; + // Plane position of the module's centre. + [[nodiscard]] std::pair PlaneCentre() const; + // Moved by a refinement's rigid step, which is expressed on the plane (ModuleShift: turn about + // (x0, y0) of the shift, then shift). + void Apply(const ModuleShift &shift); +}; + +// The modules of one detector. Empty means the image is one flat detector: every mapping is the identity. +class DetectorModules { + std::vector modules; +public: + DetectorModules() = default; + explicit DetectorModules(std::vector modules); + + [[nodiscard]] bool empty() const; + [[nodiscard]] size_t size() const; + [[nodiscard]] const DetectorModule &operator[](size_t i) const; + [[nodiscard]] const std::vector &Get() const; + + // True where some module is placed other than as imaged. + [[nodiscard]] bool Moved() const; + + // Module of an image pixel, or -1 (a gap, or outside every module). + [[nodiscard]] int64_t ModuleOfImagePixel(double x, double y) const; + // Module whose placed footprint holds a plane point, or -1. + [[nodiscard]] int64_t ModuleOfPlanePoint(double px, double py) const; + + // Image -> plane; a point on no module is left as it is. + [[nodiscard]] std::pair ImageToPlane(double x, double y) const; + // Plane -> image; a point on no placed module is left as it is. + [[nodiscard]] std::pair PlaneToImage(double px, double py) const; + + // Each module moved by the shift whose centre is nearest its own plane centre (and within it). + void Apply(const std::vector &shifts); + // The modules as refinement starting points: their plane centres, no move. + [[nodiscard]] std::vector Centres() const; +}; + +// The modules of a modular (JUNGFRAU) detector where its assembled image puts them: one box per module, +// CONVERTED_MODULE_COLS x CONVERTED_MODULE_LINES turned the way the geometry turns it, placed as imaged. +// Empty for a detector of one module. +DetectorModules ModulesAsImaged(const DetectorGeometry &geom); diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index d01da5c51..db1ce112b 100644 --- a/common/DiffractionExperiment.cpp +++ b/common/DiffractionExperiment.cpp @@ -1448,6 +1448,15 @@ std::optional DiffractionExperiment::GetBandwidthFWHM() const { return dataset.GetBandwidthFWHM(); } +DiffractionExperiment &DiffractionExperiment::Modules(const DetectorModules &input) { + detector_modules = input; + return *this; +} + +const DetectorModules &DiffractionExperiment::GetModules() const { + return detector_modules; +} + DiffractionExperiment &DiffractionExperiment::SaveCalibration(const std::optional &input) { dataset.SaveCalibration(input); return *this; diff --git a/common/DiffractionExperiment.h b/common/DiffractionExperiment.h index 03515c17c..8d670e260 100644 --- a/common/DiffractionExperiment.h +++ b/common/DiffractionExperiment.h @@ -34,6 +34,7 @@ constexpr int64_t SaturationLimitFromValue(int64_t saturation_value) { return sa #include "AzimuthalIntegrationSettings.h" #include "FileWriterSettings.h" #include "DiffractionGeometry.h" +#include "DetectorModules.h" #include "CompressedImage.h" #include "IndexingSettings.h" #include "BraggIntegrationSettings.h" @@ -93,6 +94,10 @@ class DiffractionExperiment { bool cpu_summation; std::string detector_update_zmq_addr; + + // Where the detector's modules are on the plane the diffraction geometry describes (DetectorModules.h). + // Empty: the image is one flat detector. + DetectorModules detector_modules; public: // Public methods are atomic DiffractionExperiment(); @@ -421,6 +426,8 @@ public: [[nodiscard]] FileWriterFormat GetFileWriterFormat() const; DiffractionGeometry GetDiffractionGeometry() const; + DiffractionExperiment& Modules(const DetectorModules &input); + [[nodiscard]] const DetectorModules &GetModules() const; void CalcAzIntCorrRawCoord(float *output, size_t module_number) const; void CalcSpotFinderResolutionMap(float *data, size_t module_number) const; CompressedImageMode GetImageMode() const; diff --git a/common/ModuleShift.h b/common/ModuleShift.h new file mode 100644 index 000000000..cf4dcc907 --- /dev/null +++ b/common/ModuleShift.h @@ -0,0 +1,13 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +// A rigid move of one detector module on the detector plane (DetectorModules.h): a turn by rot_rad about +// the module's centre (x0, y0) followed by a shift (dx, dy), in the plane's pixels. A point the module +// put at p is moved to (x0, y0) + R(rot_rad) (p - (x0, y0)) + (dx, dy). +struct ModuleShift { + double x0 = 0, y0 = 0; + double dx = 0, dy = 0, rot_rad = 0; + int spots = 0; // indexed spots the refinement had on the module (output only) +}; diff --git a/docs/HDF5.md b/docs/HDF5.md index 2cdc4c73b..0bddb0b26 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -258,6 +258,24 @@ pixel-direction vectors carry the discrete image orientation (mirror in Y, multi beam); for a detector that this system assembled itself, they are the McStas form of the internal +x and +y, i.e. `(-1, 0, 0)` and `(0, -1, 0)`. +### `/entry/instrument/detector/module_` (NXdetector_module, one per module) + +The master files written here keep the single `module` group above. Rugnux writes the placement of a +detector's modules, where a stills run placed them, into `_modules.h5`: one group per module in +the NeXus `NXdetector_module` layout that NXmx uses, and reads the same groups back (`--module-geometry`, +or a master that states its modules one by one): + +| Field | Units / notes | +|-------|-------| +| `data_origin` | `[row, column]` of the module's first pixel in the image (slowest index first) | +| `data_size` | `[rows, columns]` of the module | +| `module_offset` | translation (m) from the image's pixel (0, 0) to the module's first pixel; `depends_on` the detector's `transformations/translation` | +| `fast_pixel_direction`, `slow_pixel_direction` | translation of one pixel (m) along the module's columns and rows; `depends_on` the module's `module_offset` | + +Vectors are in the McStas frame. A pixel (row i, column j) of a module is at +`translation + module_offset + j * fast_pixel_direction + i * slow_pixel_direction`, which is how +DIALS/dxtbx reads it (one panel per module). + ### `/entry/sample` (NXsample) | Field | Std | Units / notes | diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index ef58470c6..807579b96 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -194,7 +194,8 @@ One page per job, so the answer needed is near the top of a short page: - [Installing Rugnux](RUGNUX_INSTALL.md) — packages, the release archive, GPU drivers, building from source, hardware. - [What Rugnux reads](RUGNUX_FORMATS.md) — will it open your data: NXmx / EIGER masters, PILATUS - miniCBF, marCCD and SMV (ADSC, Rigaku d\*TREK) sweeps, one sweep per input. + miniCBF, marCCD and SMV (ADSC, Rigaku d\*TREK) sweeps, one sweep per input, and serial XFEL stills + from SACLA (Cheetah HDF5), Cheetah CXI single-panel files and LCLS-I XTC. - [Running Rugnux](RUGNUX_TUTORIAL.md) — a first run in detail, rotation, serial and small-molecule data, and every file a run writes. - [Rugnux with other programs](RUGNUX_INTEGRATION.md) — the reflection-file conventions, the diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 6c8d56d0b..a1a983dfc 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -492,6 +492,13 @@ Geometry overrides (defaults are taken from the input file; override them to rep | `--beam-y ` | Beam centre Y (pixel). `--beam-center-check` still runs and can adopt its measured centre ([§1.4](CPU_DATA_ANALYSIS_IMAGE.md)) | | `--detector-distance ` | Detector distance (mm) | | `--wavelength ` | Wavelength (Å) | +| `--pixel-size ` | Pixel size (µm), for a file that does not state it (an XTC stream of an unknown camera) | +| `--gain ` | Detector units per photon, for SACLA, CXI, XTC and `--crystfel-geometry` input: the frames are converted to photons with it. Default: for SACLA and CXI input, measured on the background of a few frames; XTC frames stay in detector units ([Serial XFEL inputs](RUGNUX_FORMATS.md#serial-xfel-inputs)) | +| `--saturation ` | Pixel value from which a pixel is overloaded, in the units the frames are read in. Replaces the file's limit and the one stills read off the frames | +| `--crystfel-geometry ` | Place the panels of an HDF5/CXI input by a CrystFEL geometry file instead of the file's own metadata: module geometry and camera length only (optional; [Serial XFEL inputs](RUGNUX_FORMATS.md#serial-xfel-inputs)) | +| `--module-geometry ` | Where the detector's modules are: an NXmx file with one NXdetector_module per module (what a run writes as `_modules.h5`) or a CrystFEL geometry file whose panels are laid on the image. Overrides the input's own | +| `--refine-modules[=off]` | Stills: place each module as a rigid block in the geometry bundle adjustment, where that predicts held-out spots better (default: on). Used for a foreign multi-panel detector, measured and reported only for a JUNGFRAU; the placement is written to `_modules.h5` | +| `--mask-good :` / `--mask-bad :` | HDF5/CXI stills: a mask array in another file, in the frames' stored layout; usable only where non-zero / unusable where non-zero. Repeatable ([Serial XFEL inputs](RUGNUX_FORMATS.md#serial-xfel-inputs)) | | `--rot1 ` | PONI detector rotation 1 (rad) | | `--rot2 ` | PONI detector rotation 2 (rad) | | `--rot3 ` | PONI detector rotation 3, about the beam (rad) — the one `rugnux --mode calibration` writes into its `.poni`, so a refined third rotation can be given straight back | diff --git a/docs/RUGNUX_FORMATS.md b/docs/RUGNUX_FORMATS.md index bc6dada36..9e5838fbc 100644 --- a/docs/RUGNUX_FORMATS.md +++ b/docs/RUGNUX_FORMATS.md @@ -10,7 +10,7 @@ sweep to its own `_unmerged.mtz` and merge them downstream (see [Taking the data onward](RUGNUX_INTEGRATION.md#taking-the-data-onward)). Stills are given as all the files of the dataset, which the shell expands from a glob and which are -read in the order given: `rugnux run0237_master.h5 run0238_master.h5` (or `rugnux *_master.h5`) for +read in the order given: `rugnux run1_master.h5 run2_master.h5` (or `rugnux *_master.h5`) for several Jungfraujoch/NXmx collections of one setup, `rugnux run*.h5` for SACLA, `rugnux *.cxi` for a Cheetah CXI set and `rugnux *.xtc` for LCLS streams. Several files must all be of one format; a directory is refused. @@ -65,12 +65,81 @@ are read so that a CCD sweep does not need a conversion step to be processed. On fields named above are read, and the headers vary by site and vintage, so no promise is made that every writer's dialect opens. The HDF5 and miniCBF paths are the maintained ones. -**XFEL data saved as individual panels is not read.** The detectors of serial XFEL endstations — -CSPAD, AGIPD, ePix and their kin — write each panel as its own array, with the panel positions kept -in a separate geometry file rather than in the image file. Rugnux has no reader for those per-panel -containers and no way to take an external geometry, so such data cannot be opened; it belongs in -the pipelines built around panel geometries, such as CrystFEL or cctbx.xfel. An XFEL dataset that -has already been assembled into single images in an NXmx master opens like any other HDF5 input. +### Serial XFEL inputs + +Stills from the XFEL serial-crystallography pipelines are read straight from the files those +pipelines write, **without a geometry file**. Each reader takes from the file what it states; what +the file does not state is either measured by rugnux or given on the command line. Every run is read +at **one wavelength** (the mean of the per-shot photon energies where the file states them), the +**beam centre is measured** (the file's centre, the background-isotropy measurement and the +point-symmetry capture are indexed on a sample of frames and the best is kept), and a SACLA or CXI +file that states **no gain** is converted to photons with the gain its background shows (`--gain` +overrides; XTC frames stay in detector units unless `--gain` is given). +The saturation level is read off the frames where the file states none. + +| Format | Read from the file | Measured by rugnux | Given on the command line | +|---|---|---|---| +| **Cheetah HDF5 with per-shot groups and a panel table** (`/tag-N/data`, sensors stacked in one frame, their placement in `/metadata` - as SACLA writes for the MPCCD) | the sensors' positions and in-plane angles, pixel size, nominal camera length and pixel mask (`/metadata`); each shot's photon energy | beam centre, gain, saturation; with a cell, the camera length (the stated one is nominal) and each sensor's placement | nothing; `--detector-distance` to give the camera length | +| **Cheetah CXI, single-panel detector** (`/entry_1/data_1/data`, square frames, e.g. a Rayonix at PAL-XFEL) | pixel size, photon energy per shot, camera length, `mask_shared` | beam centre, gain, saturation | `--detector-distance` where the stated camera length is wrong | +| **LCLS-I XTC, Rayonix** (Camera::FrameV1 events; one or more streams) | photon energy per event (electron-beam record), pixel size (Rayonix model and binning), the camera's pixel offset | beam centre, saturation | `--detector-distance` (not recorded in XTC); `--pixel-size`/`--wavelength` only where the XTC lacks them | + +The SACLA sensor positions in `/metadata` are the detector's nominal metrology, a few pixels from a +per-beamtime refined geometry; the stills geometry refinement places each sensor as a module (see +[Module geometry](#module-geometry)). + +What is **refused**, with a message rather than a wrong answer: + +* a Cheetah CXI file of a **tiled detector** (CSPAD, AGIPD, ePix, JUNGFRAU, MPCCD stored as a slab of + panels): the file does not say where the panels lie; +* **multi-panel XTC** (CSPAD, ePix, JUNGFRAU at LCLS), stored raw with per-module calibration, and an XTC + stream holding images of the same size from two cameras; +* **unconverted detector data** in general (raw ADC values with gain stages and pedestals); +* **per-shot wavelengths**: they are averaged, never followed shot by shot; +* files of one run that store their frames differently (size, type, compression) or place the sensors + differently. + +Damaged or compressed XTC images are left out, and so is a last event that a partly downloaded stream +ends inside. An XFEL dataset that has already been assembled into single images in an NXmx master +opens like any other HDF5 input. + +`--crystfel-geometry ` places the panels of an HDF5/CXI input by a CrystFEL geometry file. Only +the module geometry is taken from it (panel pixel ranges, fast/slow directions, corners, pixel size, +per-panel `coffset`, the frames' data path) and, optionally, the camera length; everything else it may +hold (gain, masks, bad regions, photon energy, value flags) is ignored, with a note in the log. DIALS +`.expt` files are not read. + +A mask array the data files do not hold themselves - e.g. a per-run pixel mask beside the frames - is +given with `--mask-good :` (a pixel is usable only where the array is non-zero) or +`--mask-bad :` (unusable where non-zero), in the frames' own layout, as stored. Both may be +given more than once, and add to the masks the data files state (the CXI `mask_shared` and per-shot +mask, the `/metadata/pixelmask` of a panel-table file). + +### Module geometry + +A detector of several modules is described module by module: where each module's pixels are in the image +rugnux reads, and where they really are on the detector plane. The modules come from an explicit +geometry, never from the image: + +* an NXmx file with one `NXdetector_module` group per module (`/entry/instrument/detector/module_0`, + `module_1`, ...: what a run writes as `_modules.h5`; the masters Jungfraujoch writes still state + the image as one module); +* the panels the XFEL readers know (the SACLA sensors from `/metadata`, the panels of a + `--crystfel-geometry` file); each panel is put on the image whole, at a whole pixel, and its fraction of + a pixel and small turn are its placement; +* `--module-geometry `: an NXmx file with such groups (what a run writes as `_modules.h5`) + or a CrystFEL geometry file whose panels are laid on the image (`fs` = column, `ss` = row, corners + relative to the beam). It replaces the input's own modules - the way to give the module layout of a + master written before the modules were recorded. + +Spots are mapped from the image onto the plane, and predictions from the plane back onto the image for +integration, so the image itself is never re-assembled. With a cell given, the stills geometry refinement +places each module with enough spots on its frames' lattices as a rigid block (shift and in-plane turn, +`--refine-modules`, on by default) where that predicts held-out spots better; modules move only relative +to one another, the beam centre, distance and crystal orientations carrying what all of them share. The +placement is used from there on for the panels of a foreign multi-panel detector (the sensors of a +panel-table file, the panels of a `--crystfel-geometry` file); for a JUNGFRAU (an NXmx master) it is +measured and reported only. Either way it is written to `_modules.h5`, which a later run reads +back with `--module-geometry`. 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/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 61b9884f2..0e50fbc00 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -4,6 +4,7 @@ #include "MXAnalysisWithoutFPGA.h" #include +#include #include #include "spot_finding/StrongPixelSet.h" @@ -205,14 +206,14 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, int64_t image_number) { if (count_crowding) crowding_engine->CountNeighbourCrowding(predicted, npredicted); - return bragg_engine->Run(*preprocessor_buffer, predicted, npredicted, image_number); + return IntegrateOnImage(predicted, npredicted, image_number); }; // The supercell probe integrates on the same engine, but what the engine counts is the run's // measurement of its own stencil, and the probe's reflections are not the run's: keep them out. const auto probe_integrate_fn = [this](const std::vector &predicted, size_t npredicted, int64_t image_number) { const BraggIntegrationCounts before = bragg_engine->Counts(); - auto ret = bragg_engine->Run(*preprocessor_buffer, predicted, npredicted, image_number); + auto ret = IntegrateOnImage(predicted, npredicted, image_number); const BraggIntegrationCounts after = bragg_engine->Counts(); probe_counts.predicted += after.predicted - before.predicted; probe_counts.bkg_starved += after.bkg_starved - before.bkg_starved; @@ -504,3 +505,31 @@ void MXAnalysisWithoutFPGA::UpdateMaskResolution(const SpotFindingSettings &sett spotFinder->SetResolutionMaskBits(*mask_resolution); adaptiveSpotFinder->SetResolutionMaskBits(*mask_resolution); } + +std::vector MXAnalysisWithoutFPGA::IntegrateOnImage(const std::vector &predicted, + size_t npredicted, int64_t image_number) { + const auto &modules = experiment.GetModules(); + if (modules.empty()) + return bragg_engine->Run(*preprocessor_buffer, predicted, npredicted, image_number); + + // Predictions are where the reflections fall on the detector plane; the engine reads the image, so + // each goes to where its module put those pixels, and comes back to the plane once integrated. + std::vector on_image(predicted.begin(), predicted.begin() + static_cast(npredicted)); + for (auto &r: on_image) { + const auto [x, y] = modules.PlaneToImage(r.predicted_x, r.predicted_y); + r.predicted_x = static_cast(x); + r.predicted_y = static_cast(y); + } + auto ret = bragg_engine->Run(*preprocessor_buffer, on_image, npredicted, image_number); + for (auto &r: ret) { + const auto [px, py] = modules.ImageToPlane(r.predicted_x, r.predicted_y); + r.predicted_x = static_cast(px); + r.predicted_y = static_cast(py); + if (std::isfinite(r.observed_x) && std::isfinite(r.observed_y)) { + const auto [ox, oy] = modules.ImageToPlane(r.observed_x, r.observed_y); + r.observed_x = static_cast(ox); + r.observed_y = static_cast(oy); + } + } + return ret; +} diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index b12d1c38d..3d3b6e802 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -63,6 +63,10 @@ class MXAnalysisWithoutFPGA { IndexAndRefine &indexer; std::unique_ptr prediction; std::unique_ptr bragg_engine; + // The engine run on reflections predicted on the detector plane: each read where its module put it on + // the image (DetectorModules), handed back on the plane. + std::vector IntegrateOnImage(const std::vector &predicted, size_t npredicted, + int64_t image_number); // What the supercell probe's integrations added to bragg_engine's counts (see BraggCounts). BraggIntegrationCounts probe_counts; // Counts the neighbour crowding of the integrated reflections at another stencil (CountCrowdingAt), diff --git a/image_analysis/geom_refinement/GeometryRefiner.cpp b/image_analysis/geom_refinement/GeometryRefiner.cpp index ded2b3dd5..642db86d3 100644 --- a/image_analysis/geom_refinement/GeometryRefiner.cpp +++ b/image_analysis/geom_refinement/GeometryRefiner.cpp @@ -130,6 +130,76 @@ struct CellLengthRegularizer { } }; +// Where a spot found at (x, y) on module m really is (ModuleShift). +void PlaceOnModule(const ModuleShift &m, const double move[3], double x, double y, double &px, double &py) { + const double ux = x - m.x0, uy = y - m.y0; + const double c = std::cos(move[2]), s = std::sin(move[2]); + px = m.x0 + c * ux - s * uy + move[0]; + py = m.y0 + s * ux + c * uy + move[1]; +} + +// XtalResidual for a spot on a module that is being placed: block `move` = (dx, dy, rot) of the module, +// (ux, uy) the spot's offset from the module centre. Moving the spot by R(rot) u - u + (dx, dy) shifts its +// offset from the beam exactly as moving the beam the other way does, so the beam handed to XtalResidual +// carries the move and everything else is the residual of an unmoved spot. +struct ModuleXtalResidual { + XtalResidual residual; + double ux, uy; + + template + bool operator()(const T *const beam, const T *const distance_mm, const T *const detector_rot, + const T *const rotation_axis, const T *const p0, const T *const p1, const T *const p2, + const T *const move, T *out) const { + using std::cos; + using std::sin; + const T c = cos(move[2]), s = sin(move[2]); + const T moved_beam[2] = {beam[0] - ((c - T(1)) * T(ux) - s * T(uy) + move[0]), + beam[1] - (s * T(ux) + (c - T(1)) * T(uy) + move[1])}; + return residual(moved_beam, distance_mm, detector_rot, rotation_axis, p0, p1, p2, out); + } +}; + +// What module moves cannot be told apart from: all modules shifted alike is the beam centre, all turned +// alike (about their own centres) is the crystal orientations, and all pushed outwards alike is the +// distance. The sums of those three over the moved modules are held at zero - the beam, the orientations +// and the distance then carry them - so the modules are left with what only they can explain. The sums are +// weighted by the spots each module holds: a module the frames barely reach is placed poorly, and an +// unweighted sum would hand its error to all the others. +// +// With the tilt refined, the two tilt-like patterns are held at zero as well: a detector tilted by a small +// angle about y moves a point (x, y) (from the beam) by (x^2, x y) / D times the angle, one about x by +// (x y, y^2) / D. Each is taken as the modules' moves projected on that field, in units of the lever. +struct ModuleGauge { + std::vector> outward; // unit vector from the beam to each module centre + std::vector> centre; // each module centre from the beam, over lever_px + std::vector share; // each module's fraction of the spots + double lever_px; // turns are weighed at this distance from the centre + double weight; + bool tilt = false; + + template + bool operator()(T const *const *move, T *out) const { + T sx(0), sy(0), sr(0), so(0), t1(0), t2(0); + for (size_t i = 0; i < outward.size(); ++i) { + const T w(share[i]); + sx += w * move[i][0]; + sy += w * move[i][1]; + sr += w * move[i][2] * T(lever_px); + so += w * (move[i][0] * T(outward[i][0]) + move[i][1] * T(outward[i][1])); + const double x = centre[i][0], y = centre[i][1]; + t1 += w * (move[i][0] * T(x * x) + move[i][1] * T(x * y)); + t2 += w * (move[i][0] * T(x * y) + move[i][1] * T(y * y)); + } + out[0] = T(weight) * sx; + out[1] = T(weight) * sy; + out[2] = T(weight) * sr; + out[3] = T(weight) * so; + out[4] = tilt ? T(weight) * t1 : T(0); + out[5] = tilt ? T(weight) * t2 : T(0); + return true; + } +}; + } // namespace GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_geom, @@ -170,6 +240,55 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge const double recip_per_px = pixel / (distance_mm * lambda); const double loss_scale = 3.0 * recip_per_px; + // The modules' moves (dx, dy, rot), starting where settings put them. Which of them are refined is + // decided once, on the indexed spots each holds, so the problem keeps its shape from round to round. + const auto &modules = settings.modules; + std::vector> move(modules.size()); + for (size_t m = 0; m < modules.size(); ++m) + move[m] = {modules[m].dx, modules[m].dy, modules[m].rot_rad}; + // A module counts the spots on it that its frame's lattice indexes where the module starts - not + // every spot found on it: a frame's other spots (noise, a second lattice, ice) would otherwise make a + // module the frames barely reach look well determined. + std::vector module_spots(modules.size(), 0); + for (const auto &f : frames) { + const Coord v0 = f.lattice.Vec0(), v1 = f.lattice.Vec1(), v2 = f.lattice.Vec2(); + for (const auto &s : f.spots) + if (s.module >= 0 && static_cast(s.module) < modules.size()) { + double x, y; + PlaceOnModule(modules[s.module], move[s.module].data(), s.x, s.y, x, y); + const Coord recip = nominal_geom.DetectorToRecip(x, y); + const double hf = recip * v0, kf = recip * v1, lf = recip * v2; + if (std::fabs(hf - std::round(hf)) < 0.15 && std::fabs(kf - std::round(kf)) < 0.15 + && std::fabs(lf - std::round(lf)) < 0.15) + module_spots[s.module]++; + } + } + std::vector refined(modules.size(), false); + int refined_count = 0; + for (size_t m = 0; m < modules.size(); ++m) + if (settings.refine_modules && module_spots[m] >= settings.min_module_spots) { + refined[m] = true; + refined_count++; + } + // One module alone has nothing to be placed against: its move is the beam and the orientation. + if (refined_count < 2) + std::fill(refined.begin(), refined.end(), false); + // A module too thin to be placed stays where it was given, whatever its start: a search on a few dozen + // spots finds noise. + if (settings.refine_modules) + for (size_t m = 0; m < modules.size(); ++m) + if (!refined[m]) + move[m] = {0, 0, 0}; + + auto place = [&](const GeomRefineSpot &s, double &x, double &y) { + if (s.module >= 0 && static_cast(s.module) < modules.size()) + PlaceOnModule(modules[s.module], move[s.module].data(), s.x, s.y, x, y); + else { + x = s.x; + y = s.y; + } + }; + // Tolerance schedule (fractional-hkl norm), tightening each round like the per-image optimizer. auto round_tol = [&](int r) { if (settings.rounds <= 1) return 0.2; @@ -185,7 +304,8 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge EffectiveCell(system, cell_len, cell_ang, eff_len, eff_alpha, eff_beta, eff_gamma); DiffractionGeometry cur = nominal_geom; - cur.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]).DetectorDistance_mm(distance_mm); + cur.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]).DetectorDistance_mm(distance_mm) + .PoniRot1_rad(static_cast(detector_rot[0])).PoniRot2_rad(static_cast(detector_rot[1])); const double tol = round_tol(round); const double tol_sq = tol * tol; @@ -201,19 +321,32 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge int used_in_frame = 0; for (const auto &s : frames[i].spots) { - const Coord recip = cur.DetectorToRecip(s.x, s.y); + double x, y; + place(s, x, y); + const Coord recip = cur.DetectorToRecip(x, y); const double hf = recip * v0, kf = recip * v1, lf = recip * v2; const double h = std::round(hf), k = std::round(kf), l = std::round(lf); const double dn = (hf - h) * (hf - h) + (kf - k) * (kf - k) + (lf - l) * (lf - l); if (dn > tol_sq) continue; - problem.AddResidualBlock( - new ceres::AutoDiffCostFunction( - new XtalResidual(s.x, s.y, lambda, pixel, cos_rot3, sin_rot3, 0.0, - h, k, l, system, orientation)), - new ceres::CauchyLoss(loss_scale), - beam, &distance_mm, detector_rot, rot_vec, orient[i].data(), cell_len, cell_ang); + if (s.module >= 0 && static_cast(s.module) < modules.size() && refined[s.module]) { + const auto &m = modules[s.module]; + problem.AddResidualBlock( + new ceres::AutoDiffCostFunction( + new ModuleXtalResidual{XtalResidual(s.x, s.y, lambda, pixel, cos_rot3, sin_rot3, + 0.0, h, k, l, system, orientation), + s.x - m.x0, s.y - m.y0}), + new ceres::CauchyLoss(loss_scale), + beam, &distance_mm, detector_rot, rot_vec, orient[i].data(), cell_len, cell_ang, + move[s.module].data()); + } else + problem.AddResidualBlock( + new ceres::AutoDiffCostFunction( + new XtalResidual(x, y, lambda, pixel, cos_rot3, sin_rot3, 0.0, + h, k, l, system, orientation)), + new ceres::CauchyLoss(loss_scale), + beam, &distance_mm, detector_rot, rot_vec, orient[i].data(), cell_len, cell_ang); ++used_in_frame; ++spots_used; } @@ -231,18 +364,65 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge // Cell-length regularizer (weight ~ coeff * sqrt(N_obs) / L0^2 -> effectively pins the cell). const double sqrtN = std::sqrt(static_cast(spots_used)); - CellLengthRegularizer *reg = new CellLengthRegularizer{ - settings.cell_reg_coeff * sqrtN / (cell_len0[0] * cell_len0[0]), - settings.cell_reg_coeff * sqrtN / (cell_len0[1] * cell_len0[1]), - settings.cell_reg_coeff * sqrtN / (cell_len0[2] * cell_len0[2]), - cell_len0[0], cell_len0[1], cell_len0[2]}; - problem.AddResidualBlock( - new ceres::AutoDiffCostFunction(reg), nullptr, cell_len); + if (settings.refine_beam_distance) { + CellLengthRegularizer *reg = new CellLengthRegularizer{ + settings.cell_reg_coeff * sqrtN / (cell_len0[0] * cell_len0[0]), + settings.cell_reg_coeff * sqrtN / (cell_len0[1] * cell_len0[1]), + settings.cell_reg_coeff * sqrtN / (cell_len0[2] * cell_len0[2]), + cell_len0[0], cell_len0[1], cell_len0[2]}; + problem.AddResidualBlock( + new ceres::AutoDiffCostFunction(reg), nullptr, cell_len); + } // Shared blocks go in the second elimination group; constant ones are stripped by Ceres. for (double *b : {beam, &distance_mm, detector_rot, rot_vec, cell_len, cell_ang}) ordering->AddElementToGroup(b, 1); + // The refined modules' moves: bounded to what a misplaced module plausibly is, and held to + // zero in sum by the gauge (ModuleGauge). A 1 px mean shift costs as much there as moving + // every spot by 1 px would, ten times over. + std::vector moved; + for (size_t m = 0; m < modules.size(); ++m) + if (refined[m] && problem.HasParameterBlock(move[m].data())) + moved.push_back(m); + if (moved.size() >= 2) { + auto *gauge = new ModuleGauge; + double lever = 0, total = 0; + for (size_t m : moved) + total += module_spots[m]; + for (size_t m : moved) { + const double ox = modules[m].x0 - nominal_geom.GetBeamX_pxl(); + const double oy = modules[m].y0 - nominal_geom.GetBeamY_pxl(); + const double r = std::max(1.0, std::hypot(ox, oy)); + gauge->outward.push_back({ox / r, oy / r}); + gauge->centre.push_back({ox, oy}); + gauge->share.push_back(module_spots[m] / total); + lever += r; + } + gauge->lever_px = lever / static_cast(moved.size()); + for (auto &c : gauge->centre) + c = {c[0] / gauge->lever_px, c[1] / gauge->lever_px}; + gauge->weight = 10.0 * recip_per_px * sqrtN; + gauge->tilt = settings.refine_tilt; + auto *cost = new ceres::DynamicAutoDiffCostFunction(gauge); + std::vector blocks; + for (size_t m : moved) { + cost->AddParameterBlock(3); + blocks.push_back(move[m].data()); + } + cost->SetNumResiduals(6); + problem.AddResidualBlock(cost, nullptr, blocks); + } + for (size_t m : moved) { + ordering->AddElementToGroup(move[m].data(), 1); + const double start[3] = {modules[m].dx, modules[m].dy, modules[m].rot_rad}; + const double range[3] = {12.0, 12.0, 1.5 * PI / 180.0}; + for (int j = 0; j < 3; ++j) { + problem.SetParameterLowerBound(move[m].data(), j, start[j] - range[j]); + problem.SetParameterUpperBound(move[m].data(), j, start[j] + range[j]); + } + } + // Beam + distance bounds around the nominal value; cell lengths loosely bounded (the // regularizer does the anchoring). problem.SetParameterLowerBound(beam, 0, nominal_geom.GetBeamX_pxl() - settings.beam_range_px); @@ -260,9 +440,23 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge if (!settings.refine_distance) problem.SetParameterBlockConstant(&distance_mm); - problem.SetParameterBlockConstant(detector_rot); // tilt not refined + if (settings.refine_tilt && settings.refine_beam_distance) { + const double tilt_range = 2.0 * PI / 180.0; + for (int j = 0; j < 2; ++j) { + problem.SetParameterLowerBound(detector_rot, j, (j == 0 ? nominal_geom.GetPoniRot1_rad() + : nominal_geom.GetPoniRot2_rad()) - tilt_range); + problem.SetParameterUpperBound(detector_rot, j, (j == 0 ? nominal_geom.GetPoniRot1_rad() + : nominal_geom.GetPoniRot2_rad()) + tilt_range); + } + } else + problem.SetParameterBlockConstant(detector_rot); problem.SetParameterBlockConstant(rot_vec); problem.SetParameterBlockConstant(cell_ang); // angles not refined in v1 + if (!settings.refine_beam_distance) { + problem.SetParameterBlockConstant(beam); + problem.SetParameterBlockConstant(&distance_mm); + problem.SetParameterBlockConstant(cell_len); + } ceres::Solver::Options options; options.linear_solver_type = ceres::DENSE_SCHUR; @@ -284,26 +478,52 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge // Final fit quality: median detector residual of the retained spots at the refined geometry. EffectiveCell(system, cell_len, cell_ang, eff_len, eff_alpha, eff_beta, eff_gamma); DiffractionGeometry cur = nominal_geom; - cur.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]).DetectorDistance_mm(distance_mm); + cur.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]).DetectorDistance_mm(distance_mm) + .PoniRot1_rad(static_cast(detector_rot[0])).PoniRot2_rad(static_cast(detector_rot[1])); std::vector residuals; + constexpr double CAP_PX = 2.0; + double capped_sum = 0; + size_t capped_n = 0; + std::vector module_sq(modules.size(), 0), module_d(modules.size(), 0); + std::vector module_n(modules.size(), 0); + result.predicted_px.resize(frames.size()); for (size_t i = 0; i < frames.size(); ++i) { const CrystalLattice latt = AngleAxisAndCellToLattice( orient[i].data(), eff_len, eff_alpha, eff_beta, eff_gamma); const Coord a = latt.Astar(), b = latt.Bstar(), c = latt.Cstar(); const Coord v0 = latt.Vec0(), v1 = latt.Vec1(), v2 = latt.Vec2(); for (const auto &s : frames[i].spots) { - const Coord recip = cur.DetectorToRecip(s.x, s.y); + double x, y; + place(s, x, y); + const Coord recip = cur.DetectorToRecip(x, y); const double hf = recip * v0, kf = recip * v1, lf = recip * v2; const double h = std::round(hf), k = std::round(kf), l = std::round(lf); const double dn = (hf - h) * (hf - h) + (kf - k) * (kf - k) + (lf - l) * (lf - l); - if (dn > 0.1 * 0.1) - continue; const Coord recip_pred = a * (float)h + b * (float)k + c * (float)l; const auto [px, py] = cur.RecipToDetector(recip_pred); + const double res = (std::isfinite(px) && std::isfinite(py)) ? std::hypot(px - x, py - y) : CAP_PX; + capped_sum += std::min(res, CAP_PX) * std::min(res, CAP_PX); + capped_n++; + // The prediction where the spot would be found: the spot's own module's move taken back. + float fx = NAN, fy = NAN; + if (std::isfinite(px) && std::isfinite(py)) { + fx = static_cast(s.x + (px - x)); + fy = static_cast(s.y + (py - y)); + } + result.predicted_px[i].push_back({fx, fy}); + if (s.module >= 0 && static_cast(s.module) < modules.size()) { + module_sq[s.module] += std::min(res, CAP_PX) * std::min(res, CAP_PX); + module_d[s.module] += 1.0 / std::max(1e-6f, recip.Length()); + module_n[s.module]++; + } + if (dn > 0.1 * 0.1) + continue; if (std::isfinite(px) && std::isfinite(py)) - residuals.push_back(std::hypot(px - s.x, py - s.y)); + residuals.push_back(res); } } + if (capped_n > 0) + result.capped_rms_px = std::sqrt(capped_sum / static_cast(capped_n)); if (!residuals.empty()) { std::nth_element(residuals.begin(), residuals.begin() + residuals.size() / 2, residuals.end()); result.median_residual_px = residuals[residuals.size() / 2]; @@ -314,9 +534,85 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge result.beam_y_px = beam[1]; result.distance_mm = distance_mm; result.cell = BuildUnitCell(system, cell_len, cell_ang); + result.rot1_rad = detector_rot[0]; + result.rot2_rad = detector_rot[1]; + result.module_capped_rms_px.resize(modules.size(), 0); + result.module_mean_d_A.resize(modules.size(), 0); + for (size_t m = 0; m < modules.size(); ++m) + if (module_n[m] > 0) { + result.module_capped_rms_px[m] = std::sqrt(module_sq[m] / static_cast(module_n[m])); + result.module_mean_d_A[m] = module_d[m] / static_cast(module_n[m]); + } + result.modules = modules; + for (size_t m = 0; m < modules.size(); ++m) { + result.modules[m].dx = move[m][0]; + result.modules[m].dy = move[m][1]; + result.modules[m].rot_rad = move[m][2]; + result.modules[m].spots = module_spots[m]; + } return result; } catch (...) { result.ok = false; return result; } } + +std::vector ScanModuleShifts(const DiffractionGeometry &geom, + const std::vector &frames, + const std::vector &modules, + double range_px) { + // A spot is on its lattice when every fractional index is within this of a whole number. + constexpr double ON_LATTICE = 0.15; + + struct Spot { float x, y; int frame; }; + std::vector> on_module(modules.size()); + for (size_t i = 0; i < frames.size(); i++) + for (const auto *list : {&frames[i].spots, &frames[i].unindexed}) + for (const auto &s : *list) + if (s.module >= 0 && static_cast(s.module) < modules.size()) + on_module[s.module].push_back({s.x, s.y, static_cast(i)}); + + std::vector> axes(frames.size()); + for (size_t i = 0; i < frames.size(); i++) + axes[i] = {frames[i].lattice.Vec0(), frames[i].lattice.Vec1(), frames[i].lattice.Vec2()}; + + auto score = [&](const std::vector &spots, const ModuleShift &m, double dx, double dy) { + int n = 0; + for (const auto &s : spots) { + double x, y; + const double move[3] = {m.dx + dx, m.dy + dy, m.rot_rad}; + PlaceOnModule(m, move, s.x, s.y, x, y); + const Coord recip = geom.DetectorToRecip(x, y); + bool on = true; + for (const auto &v : axes[s.frame]) { + const double f = recip * v; + on &= std::fabs(f - std::round(f)) < ON_LATTICE; + } + n += on; + } + return n; + }; + + std::vector out(modules.size()); + for (size_t m = 0; m < modules.size(); m++) { + auto &r = out[m]; + r.on_lattice_at_zero = r.on_lattice_at_best = score(on_module[m], modules[m], 0, 0); + // Coarse: whole-pixel steps over the range; then half-pixel steps around the best of them. + for (int pass = 0; pass < 2; pass++) { + const double step = pass == 0 ? 1.0 : 0.5; + const double half = pass == 0 ? range_px : 1.0; + const double cx = r.dx, cy = r.dy; + for (double dy = -half; dy <= half; dy += step) + for (double dx = -half; dx <= half; dx += step) { + const int n = score(on_module[m], modules[m], cx + dx, cy + dy); + if (n > r.on_lattice_at_best) { + r.on_lattice_at_best = n; + r.dx = cx + dx; + r.dy = cy + dy; + } + } + } + } + return out; +} + diff --git a/image_analysis/geom_refinement/GeometryRefiner.h b/image_analysis/geom_refinement/GeometryRefiner.h index ccd18cd4e..cb1d79aa1 100644 --- a/image_analysis/geom_refinement/GeometryRefiner.h +++ b/image_analysis/geom_refinement/GeometryRefiner.h @@ -3,11 +3,13 @@ #pragma once +#include #include #include #include "../../common/CrystalLattice.h" #include "../../common/DiffractionGeometry.h" +#include "../../common/ModuleShift.h" #include "../../common/UnitCell.h" #include "gemmi/symmetry.hpp" @@ -22,11 +24,14 @@ struct GeomRefineSpot { float x = 0, y = 0; // observed (raw) spot position, pixels int32_t h = 0, k = 0, l = 0; // integer Miller index assigned in pass 1 + int16_t module = -1; // detector module the spot lies on (GeometryRefinerSettings::modules); -1: none }; + struct GeomRefineFrame { CrystalLattice lattice; // per-frame indexed lattice (real space), for the orientation seed std::vector spots; // indexed spots on this frame + std::vector unindexed; // the frame's other spots (no index); see ScanModuleShifts int32_t ordinal = 0; // image this came from; ties are broken on it, so the // selected bundle does not depend on collection order }; @@ -43,6 +48,21 @@ struct GeometryRefinerSettings { double distance_range_mm = 8.0; // distance bound around the nominal value int min_spots_total = 100; // refuse to refine below this many indexed spots int num_threads = 1; + // false: beam, distance and cell are held where nominal_geom and input_cell put them and only the + // crystal orientations follow - the fit quality of a geometry decided elsewhere, on frames it was + // not fitted to. + bool refine_beam_distance = true; + // The detector's modules, with where each is placed to start from; empty: the image is one rigid + // detector. With refine_modules, every module that has at least min_module_spots spots on its frames' + // lattices where it starts is moved as a rigid block (shift and turn in the detector plane); the others + // go back to where they were given (no move), whatever their start. + std::vector modules; + bool refine_modules = false; + int min_module_spots = 500; + // Refine the detector tilt (PONI rot1, rot2) with the beam and distance. A tilt moves the modules far + // from the beam more than those near it, in a pattern module moves could otherwise take up: with it + // refined, the modules are also held to no tilt-like pattern of their own (ModuleGauge). + bool refine_tilt = false; }; struct GeometryRefinerResult { @@ -53,6 +73,19 @@ struct GeometryRefinerResult { int frames_used = 0; int spots_used = 0; double median_residual_px = 0; // final per-spot detector residual (fit quality) + // The same over EVERY spot of the frames, each spot's residual to its nearest predicted reflection + // capped at 2 px: unlike the median of the spots that fit, it does not improve by fitting fewer of them, + // so two geometries can be compared on the same spots. + double capped_rms_px = 0; + std::vector modules; // the modules' placement at the end (settings.modules if not refined) + double rot1_rad = 0, rot2_rad = 0; // detector tilt at the end (the nominal one where not refined) + // Per module (as settings.modules): rms of every spot's residual capped at 2 px, and the mean + // resolution of its spots - what a module's fit rests on. + std::vector module_capped_rms_px; + std::vector module_mean_d_A; + // Per frame, per spot (as frames[i].spots): where the refined geometry predicts the spot's reflection on + // the detector, as the spot would be found (before its module's move); NaN where it predicts none. + std::vector>> predicted_px; }; // nominal_geom supplies the starting beam/distance/tilt/wavelength/pixel size; input_cell is anchored by @@ -61,3 +94,21 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge const UnitCell &input_cell, const std::vector &frames, const GeometryRefinerSettings &settings); + +// A coarse search for where each module is, ahead of RefineGlobalGeometry: a module misplaced by more than a +// fraction of the spot spacing has its spots left unindexed (or indexed wrongly), so the refinement has +// nothing to pull it with. Each module is shifted over a grid of +-range_px, and at each shift every spot on +// it (indexed or not) is scored against its frame's lattice as indexed: the shift that puts the most spots on +// the lattices wins. The frames' lattices come from the other modules, so this works while most of the +// detector is roughly right. Turns are not searched. +struct ModuleScan { + double dx = 0, dy = 0; // the best shift found + int on_lattice_at_zero = 0; + int on_lattice_at_best = 0; +}; + +std::vector ScanModuleShifts(const DiffractionGeometry &geom, + const std::vector &frames, + const std::vector &modules, + double range_px); + diff --git a/image_analysis/spot_finding/SpotUtils.cpp b/image_analysis/spot_finding/SpotUtils.cpp index a844d7033..1ada9ee7f 100644 --- a/image_analysis/spot_finding/SpotUtils.cpp +++ b/image_analysis/spot_finding/SpotUtils.cpp @@ -372,9 +372,18 @@ void SpotAnalyze(const DiffractionExperiment &experiment, std::vector spots_out; + // A spot is found on the image; from here on it is where its module puts it on the detector plane. + const auto &modules = experiment.GetModules(); for (const auto &spot: spots) { - if (auto s = spot.Export(geom, output.number); s.has_value()) + if (auto s = spot.Export(geom, output.number); s.has_value()) { + if (!modules.empty()) { + const auto [px, py] = modules.ImageToPlane(s->x, s->y); + s->x = static_cast(px); + s->y = static_cast(py); + s->d_A = geom.PxlToRes(s->x, s->y); + } spots_out.push_back(s.value()); + } } if (spot_finding_settings.high_res_gap_Q_recipA.has_value()) diff --git a/reader/CMakeLists.txt b/reader/CMakeLists.txt index b9b3c75a3..20111de43 100644 --- a/reader/CMakeLists.txt +++ b/reader/CMakeLists.txt @@ -16,8 +16,9 @@ ADD_LIBRARY(JFJochReader STATIC JFJochMarCCDReader.h CrystFELGeometry.cpp CrystFELGeometry.h - JFJochCrystFELReader.cpp - JFJochCrystFELReader.h + JFJochXFELReader.cpp + JFJochXFELReader.h + DetectorLayout.h SMV.cpp SMV.h JFJochSMVReader.cpp @@ -26,6 +27,8 @@ ADD_LIBRARY(JFJochReader STATIC XTC.h JFJochXTCReader.cpp JFJochXTCReader.h + NXmxModules.cpp + NXmxModules.h HDF5ImageLocator.cpp JFJochStartMessageReader.cpp JFJochStartMessageReader.h diff --git a/reader/CrystFELGeometry.cpp b/reader/CrystFELGeometry.cpp index f377fbd17..f66e39e22 100644 --- a/reader/CrystFELGeometry.cpp +++ b/reader/CrystFELGeometry.cpp @@ -6,11 +6,16 @@ #include #include #include +#include +#include #include "../common/JFJochException.h" namespace crystfel { +using xfel::Value; +using xfel::Panel; + namespace { std::string Trim(const std::string &s) { @@ -20,8 +25,54 @@ std::string Trim(const std::string &s) { return s.substr(b, s.find_last_not_of(" \t\r") - b + 1); } -// "+0.000390x -1.000000y", "-x", "0.5x+0.5y": the in-plane direction of a panel axis. -void ParseDirection(const std::string &txt, double &x, double &y) { +[[noreturn]] void Refuse(const std::string &filename, const std::string &key, const std::string &why) { + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + filename + ": '" + key + "' " + why + "; this geometry is not read"); +} + +// "0.1 m", "70.5 mm", "/LCLS/detector_1/EncoderValue", "9.3 keV": a number or a path, with an optional +// unit after it. units maps each accepted unit to its factor into the unit held, and "" to the default. +Value ParseValue(const std::string &filename, const std::string &key, const std::string &text, + const std::map &units, double path_default) { + std::istringstream ss(text); + std::string token, unit, rest; + ss >> token >> unit >> rest; + if (token.empty() || !rest.empty()) + Refuse(filename, key, "has a value that is not a number or a path with one unit (" + text + ")"); + const bool is_path = token[0] == '/'; + double scale = is_path ? path_default : units.at(""); + if (!unit.empty()) { + const auto it = units.find(unit); + if (it == units.end() || unit.empty()) + Refuse(filename, key, "is stated in a unit not understood (" + unit + ")"); + scale = it->second; + } + Value v; + v.scale = scale; + if (is_path) + v.path = token; + else { + size_t used = 0; + double number = 0; + try { + number = std::stod(token, &used); + } catch (const std::exception &) { + used = 0; + } + if (used != token.size()) + Refuse(filename, key, "has a value that is not a number or a path (" + text + ")"); + v.value = number * scale; + v.scale = 1; + } + return v; +} + +const std::map LENGTH_UNITS = {{"", 1.0}, {"m", 1.0}, {"mm", 1e-3}}; + +// "+0.000390x -1.000000y", "-x", "0.5x+0.5y": the in-plane direction of a panel axis. A panel tilted +// out of the detector plane (a z component) is not read: it would be placed flat. +void ParseDirection(const std::string &filename, const std::string &key, const std::string &txt, + double &x, double &y) { x = 0; y = 0; static const std::regex term(R"(([+-]?)\s*([0-9.]*(?:[eE][+-]?[0-9]+)?)\s*([xyz]))"); @@ -33,26 +84,34 @@ void ParseDirection(const std::string &txt, double &x, double &y) { x = v; else if ((*it)[3] == "y") y = v; + else if (v != 0) + Refuse(filename, key, "tilts the panel out of the detector plane (" + txt + ")"); } } -// Sets one panel property; false if the key is not one. -bool SetPanel(Panel &p, const std::string &key, const std::string &value) { +// Keywords that only matter to refinement or peak search, read past without a word. +bool Bookkeeping(const std::string &key) { + return key.starts_with("rigid_group") || key.starts_with("group_") || key.starts_with("peak_list") + || key == "badrow_direction" || key == "peak_sep"; +} + +// Keywords that change where a pixel lands or what it is worth, which this reader does not do. +const std::set REFUSED = {"detector_shift_x", "detector_shift_y"}; + +// Sets one panel property; false if the key is not one. Only where the panel is: its pixel range, axes, +// corner, pixel size, depth offset and camera length. +bool SetPanel(const std::string &filename, Panel &p, const std::string &key, const std::string &value) { if (key == "min_fs") p.min_fs = std::stoll(value); else if (key == "max_fs") p.max_fs = std::stoll(value); else if (key == "min_ss") p.min_ss = std::stoll(value); else if (key == "max_ss") p.max_ss = std::stoll(value); - else if (key == "fs") ParseDirection(value, p.fs_x, p.fs_y); - else if (key == "ss") ParseDirection(value, p.ss_x, p.ss_y); + else if (key == "fs") ParseDirection(filename, key, value, p.fs_x, p.fs_y); + else if (key == "ss") ParseDirection(filename, key, value, p.ss_x, p.ss_y); else if (key == "corner_x") p.corner_x = std::stod(value); else if (key == "corner_y") p.corner_y = std::stod(value); else if (key == "res") p.res = std::stod(value); - else if (key == "clen") p.clen = value; + else if (key == "clen") p.clen = ParseValue(filename, key, value, LENGTH_UNITS, 1e-3); else if (key == "coffset") p.coffset = std::stod(value); - else if (key == "adu_per_eV") p.adu_per_eV = std::stod(value); - else if (key == "adu_per_photon") p.adu_per_photon = std::stod(value); - else if (key == "max_adu") p.max_adu = std::stod(value); - else if (key == "no_index") p.no_index = std::stoi(value) != 0; else return false; return true; @@ -60,17 +119,17 @@ bool SetPanel(Panel &p, const std::string &key, const std::string &value) { } // namespace -Geometry Read(const std::string &filename) { +xfel::Layout Read(const std::string &filename, std::vector &ignored) { std::ifstream in(filename); if (!in) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Cannot open " + filename); - Geometry geom; + xfel::Layout geom; // Global panel keys hold for every panel that appears after them, as in the old syntax; a file in // the current syntax states them first, so this reads both the same. Panel defaults; std::map panel_index; - std::map bad_index; + std::map dims; std::string line; while (std::getline(in, line)) { @@ -82,14 +141,14 @@ Geometry Read(const std::string &filename) { const std::string value = Trim(line.substr(eq + 1)); const auto slash = key.find('/'); + if (REFUSED.contains(key.substr(slash == std::string::npos ? 0 : slash + 1))) + Refuse(filename, key, "is not supported"); + if (slash == std::string::npos) { if (key == "data") geom.data = value; - else if (key == "photon_energy") geom.photon_energy = value; - else if (key == "mask") geom.mask = value; - else if (key == "mask_file") geom.mask_file = value; - else if (key == "mask_good") geom.mask_good = std::stoul(value, nullptr, 0); - else if (key == "mask_bad") geom.mask_bad = std::stoul(value, nullptr, 0); - else SetPanel(defaults, key, value); // anything else (rigid groups, dim, ...) is not needed + else if (key.starts_with("dim")) dims[key] = value; + else if (!SetPanel(filename, defaults, key, value) && !Bookkeeping(key)) + ignored.push_back(key); continue; } @@ -97,29 +156,36 @@ Geometry Read(const std::string &filename) { const std::string prop = key.substr(slash + 1); if (name.starts_with("bad")) { - if (!bad_index.contains(name)) { - bad_index[name] = geom.bad.size(); - geom.bad.push_back({}); - } - auto &b = geom.bad[bad_index[name]]; - if (prop == "panel") b.panel = value; - else if (prop == "min_x") { b.in_lab = true; b.min_x = std::stod(value); } - else if (prop == "max_x") { b.in_lab = true; b.max_x = std::stod(value); } - else if (prop == "min_y") { b.in_lab = true; b.min_y = std::stod(value); } - else if (prop == "max_y") { b.in_lab = true; b.max_y = std::stod(value); } - else if (prop == "min_fs") b.min_fs = std::stoll(value); - else if (prop == "max_fs") b.max_fs = std::stoll(value); - else if (prop == "min_ss") b.min_ss = std::stoll(value); - else if (prop == "max_ss") b.max_ss = std::stoll(value); + ignored.push_back(key); continue; } + if (prop == "data") { + if (value != geom.data) + Refuse(filename, key, "puts a panel in its own dataset"); + continue; + } + if (prop.starts_with("dim")) { + dims[prop] = value; + continue; + } if (!panel_index.contains(name)) { panel_index[name] = geom.panels.size(); geom.panels.push_back(defaults); geom.panels.back().name = name; } - SetPanel(geom.panels[panel_index[name]], prop, value); + if (!SetPanel(filename, geom.panels[panel_index[name]], prop, value) && !Bookkeeping(prop)) + ignored.push_back(key); + } + + // The frame layouts read: a stack of whole frames (%, ss, fs) or one frame (ss, fs). A dimension + // naming a module number lays the panels out in a way the stored frame here cannot express. + for (const auto &[key, value]: dims) { + const bool ok = (key == "dim0" && (value == "%" || (dims.size() == 2 && value == "ss"))) + || (key == "dim1" && (value == "ss" || (dims.size() == 2 && value == "fs"))) + || (key == "dim2" && value == "fs" && dims.size() == 3); + if (!ok) + Refuse(filename, key + " = " + value, "lays the frame out in a way not supported"); } if (geom.panels.empty()) @@ -128,11 +194,30 @@ Geometry Read(const std::string &filename) { if (!(p.res > 0)) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, filename + ": panel " + p.name + " has no pixel size (res)"); - if (p.clen.empty()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - filename + ": panel " + p.name + " has no camera length (clen)"); } return geom; } +DetectorModules ModulesOnImage(const xfel::Layout &geom, double beam_x, double beam_y) { + if (geom.panels.size() < 2) + return {}; + std::vector out; + for (const auto &p: geom.panels) { + DetectorModule m; + m.x0 = p.min_fs; + m.y0 = p.min_ss; + m.width = p.max_fs - p.min_fs + 1; + m.height = p.max_ss - p.min_ss + 1; + // The panel's first pixel has its centre half a pixel along both axes from the corner. + m.ox = beam_x + p.corner_x + 0.5 * (p.fs_x + p.ss_x); + m.oy = beam_y + p.corner_y + 0.5 * (p.fs_y + p.ss_y); + m.fx = p.fs_x; + m.fy = p.fs_y; + m.sx = p.ss_x; + m.sy = p.ss_y; + out.push_back(m); + } + return DetectorModules(out); +} + } // namespace crystfel diff --git a/reader/CrystFELGeometry.h b/reader/CrystFELGeometry.h index 9253f9825..18f1557ac 100644 --- a/reader/CrystFELGeometry.h +++ b/reader/CrystFELGeometry.h @@ -3,55 +3,28 @@ #pragma once -#include -#include #include #include -// A CrystFEL detector geometry file (.geom), as far as it places pixels and says what they are worth: -// the panels with their pixel ranges in the stored frame, fast/slow-scan directions and corners, the -// camera length, the pixel size, the gain and the regions to leave out. Both the current syntax and -// the older one - global keys that only take effect for the panels that follow them, bad regions -// stated in whole-frame coordinates - are read. -// -// Every length is in pixels except clen and coffset, which are in metres. A value CrystFEL takes from -// the image file instead (clen = /path, photon_energy = /path) is kept as that path, for the reader -// to resolve. +#include "DetectorLayout.h" +#include "../common/DetectorModules.h" + +// A CrystFEL detector geometry file (.geom), read for the module geometry only: the panels with their +// pixel ranges in the stored frame, fast/slow-scan directions, corners, pixel size and depth offset +// (coffset), where the frames are (data, dimN), and - optionally - the camera length. Everything else a +// geometry file may hold (photon energy, gain, saturation, masks, bad regions, value flags) is the data +// file's or the command line's to state, or rugnux's to measure: such keywords are listed in `ignored` +// for the log. Both the current syntax and the older one (global keys that take effect for the panels +// after them) are read. A keyword that would move pixels in a way not supported (a tilted panel, a +// frame laid out by a module dimension, panels in different datasets, a detector shift) is refused. namespace crystfel { -struct Panel { - std::string name; - int64_t min_fs = 0, max_fs = 0, min_ss = 0, max_ss = 0; // inclusive, in the stored frame - double fs_x = 1, fs_y = 0, ss_x = 0, ss_y = 1; - double corner_x = 0, corner_y = 0; - double res = 0; // pixels per metre - std::string clen; // metres, or an HDF5 path - double coffset = 0; // metres - std::optional adu_per_eV; - std::optional adu_per_photon; - std::optional max_adu; - bool no_index = false; -}; +// The layout the file describes; keywords read past are added to ignored, for the log. +xfel::Layout Read(const std::string &filename, std::vector &ignored); -// A region to leave out: either a pixel range on one panel (or, in the old syntax, on the whole -// frame), or a box in the laboratory plane, in pixels. -struct BadRegion { - std::string panel; // empty: fs/ss are whole-frame coordinates - bool in_lab = false; - double min_x = 0, max_x = 0, min_y = 0, max_y = 0; - int64_t min_fs = 0, max_fs = 0, min_ss = 0, max_ss = 0; -}; - -struct Geometry { - std::vector panels; - std::vector bad; - std::string data = "/data/data"; - std::string photon_energy; // eV, or an HDF5 path - std::string mask; // HDF5 path, or empty - std::string mask_file; // file holding it, or empty for the image file itself - uint32_t mask_good = 0, mask_bad = 0; -}; - -Geometry Read(const std::string &filename); +// The panels of a geometry file as the modules of an image they are laid on: each panel's pixel range is +// the image's columns (fs) and rows (ss) it covers, and its corner places it relative to the beam, which +// the image has at (beam_x, beam_y). Fewer than two panels: none. +DetectorModules ModulesOnImage(const xfel::Layout &geom, double beam_x, double beam_y); } // namespace crystfel diff --git a/reader/DetectorLayout.h b/reader/DetectorLayout.h new file mode 100644 index 000000000..3b992b516 --- /dev/null +++ b/reader/DetectorLayout.h @@ -0,0 +1,77 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include + +// Where the panels of a stills detector lie and what their pixels are worth, as the XFEL readers +// gather it from a file's own metadata (or, while it is still read, a CrystFEL geometry file): the +// panels with their pixel ranges in the stored frame, fast/slow-scan directions and corners, the camera +// length, the pixel size, the gain and the regions to leave out. The conventions are CrystFEL's - a +// corner is the very corner of the panel's first pixel, x and y in pixels in the detector plane. +// +// Lengths are in pixels, except the camera length, which is held in metres. A value taken from the +// image files (a camera length or photon energy stated per shot) is kept as its path, with the factor +// that turns the stored number into the unit held here. +namespace xfel { + +// A number stated in the layout, or a path to it in the image file, times scale. +struct Value { + std::string path; // empty: value holds the number + double value = 0; + double scale = 1; // what one stored unit at path is worth +}; + +struct Panel { + std::string name; + int64_t min_fs = 0, max_fs = 0, min_ss = 0, max_ss = 0; // inclusive, in the stored frame + double fs_x = 1, fs_y = 0, ss_x = 0, ss_y = 1; + double corner_x = 0, corner_y = 0; + double res = 0; // pixels per metre + std::optional clen; // metres + double coffset = 0; // metres + std::optional adu_per_eV; + std::optional adu_per_photon; + std::optional max_adu; + int64_t mask_edge_pixels = 0; + bool no_index = false; +}; + +// A region to leave out: either a pixel range on one panel (or, in the old syntax, on the whole +// frame), or a box in the laboratory plane, in pixels. +struct BadRegion { + std::string panel; // empty: fs/ss are whole-frame coordinates + bool in_lab = false; + double min_x = 0, max_x = 0, min_y = 0, max_y = 0; + int64_t min_fs = 0, max_fs = 0, min_ss = 0, max_ss = 0; +}; + +// A mask array: a pixel is bad where it has any of bad's bits set, or lacks any of good's - or, with +// nonzero_good, where the array is zero. +struct Mask { + std::string data; // HDF5 path + std::string file; // file holding it, or empty for the image file itself + uint32_t good = 0, bad = 0; + bool nonzero_good = false; + + [[nodiscard]] bool Bad(uint32_t v) const { + return (bad != 0 && (v & bad) != 0) || (good != 0 && (v & good) != good) || (nonzero_good && v == 0); + } +}; + +struct Layout { + std::vector panels; + std::vector bad; + std::vector masks; + std::string data = "/data/data"; + std::optional photon_energy; // eV + std::optional wavelength; // Angstrom; used where no photon_energy is stated + // A stored value below, above or at these marks a pixel unusable in that frame. + std::optional flag_lessthan, flag_morethan, flag_equal; +}; + +} // namespace xfel diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 12ba9a265..33b287c0f 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -6,6 +6,7 @@ #include #include "HDF5MetadataSource.h" +#include "NXmxModules.h" #include "spdlog/fmt/fmt.h" #include "../image_analysis/bragg_integration/CalcISigma.h" #include "../image_analysis/spot_finding/SpotUtils.h" @@ -28,7 +29,9 @@ static Coord McStasToInternal(const std::vector &v) { // separated from the tilt by looking at the module alone - so it is left as it is rather than // approximated. Every real file examined here is exactly discrete. static std::optional ReadModuleOrientation(HDF5Object *file) { - const std::string base = "/entry/instrument/detector/module/"; + // One module for the whole image, or - where the file states its modules one by one - the first. + const std::string base = file->IsDataSet("/entry/instrument/detector/module/fast_pixel_direction") + ? "/entry/instrument/detector/module/" : "/entry/instrument/detector/module_0/"; if (!file->IsDataSet(base + "fast_pixel_direction") || !file->IsDataSet(base + "slow_pixel_direction")) return {}; @@ -1092,6 +1095,9 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen detector.MinCountTime(std::chrono::microseconds(0)); detector.ReadOutTime(std::chrono::nanoseconds(0)); dataset->experiment.Detector(detector); + // The modules, where the file states them one by one (module_0, module_1, ...). + dataset->experiment.Modules(ReadNXmxModules(*master_file, detector.GetPixelSize_mm() / 1000.0, + detector.GetImageOrientation())); // frame_time is the period between frames, count_time the exposure within one. NXmx requires // neither, and a master written outside the DECTRIS toolchain often carries only count_time; diff --git a/reader/JFJochCrystFELReader.cpp b/reader/JFJochCrystFELReader.cpp deleted file mode 100644 index 84d7c4f22..000000000 --- a/reader/JFJochCrystFELReader.cpp +++ /dev/null @@ -1,422 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#include "JFJochCrystFELReader.h" - -#include -#include -#include -#include - -#include - -#include "../common/JFJochException.h" -#include "../common/Logger.h" -#include "../writer/HDF5Objects.h" - -namespace { - -bool IsPath(const std::string &s) { - return !s.empty() && s[0] == '/'; -} - -std::string ReplaceEvent(std::string path, const std::string &group) { - const auto pos = path.find('%'); - if (pos != std::string::npos) - path.replace(pos, 1, group); - return path; -} - -} // namespace - -JFJochCrystFELReader::JFJochCrystFELReader() = default; - -JFJochCrystFELReader::~JFJochCrystFELReader() { - Close(); -} - -// One number for the whole run from a value the geometry takes from the image files: a scalar, a -// per-frame vector, or one scalar per shot group. -double JFJochCrystFELReader::MeanOverFirstFile(const std::string &path) { - auto &file = *open_files_[0]; - std::vector values; - if (path.find('%') != std::string::npos) { - for (const auto &f: frames_) - if (f.file == 0) - values.push_back(HDF5DataSet(file, ReplaceEvent(path, f.group)).ReadScalar()); - } else if (file.GetDimension(path).empty()) - values.push_back(HDF5DataSet(file, path).ReadScalar()); - else - values = file.ReadVector(path); - if (values.empty()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No values in " + path); - return std::accumulate(values.begin(), values.end(), 0.0) / static_cast(values.size()); -} - -void JFJochCrystFELReader::ReadFiles(const std::vector &files, const std::string &geometry_file, - std::optional gain) { - Close(); - geom_ = crystfel::Read(geometry_file); - if (gain) - for (auto &p: geom_.panels) - p.adu_per_photon = *gain; - Open(files, geometry_file); -} - -// A Cheetah CXI file carries what one flat panel needs - the frames at /entry_1/data_1/data, the pixel -// size, the photon energy and a camera length, all in SI units - but not where the beam is. The panel is -// put with its middle on the beam, which is as good as no centre at all: the run measures it. -// -// Only a square frame is taken this way. A Cheetah file from a tiled detector (CSPAD, AGIPD, ePix, -// JUNGFRAU, MPCCD) holds its panels stacked in a slab, which read as one flat panel would place every -// pixel wrongly and say nothing; those slabs are not square, the single-panel CCDs are, and a tiled -// detector needs its geometry file (--crystfel-geometry) anyway. -bool JFJochCrystFELReader::CanReadCXI(const std::string &path) { - try { - if (!std::filesystem::is_regular_file(path)) - return false; - std::unique_lock ul(hdf5_mutex); - HDF5ReadOnlyFile file(path); - if (!file.Exists("/entry_1/data_1/data")) - return false; - const auto dim = file.GetDimension("/entry_1/data_1/data"); - return dim.size() == 3 && dim[1] == dim[2] - && file.Exists("/entry_1/instrument_1/detector_1/x_pixel_size") - && file.Exists("/entry_1/instrument_1/source_1/energy"); - } catch (...) { - return false; - } -} - -void JFJochCrystFELReader::ReadCXI(const std::vector &files, std::optional gain) { - Close(); - if (files.empty()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No CXI file to read"); - { - std::unique_lock ul(hdf5_mutex); - HDF5ReadOnlyFile file(files[0]); - const auto mean = [&](const std::string &name) { - const auto v = file.ReadVector(name); - return std::accumulate(v.begin(), v.end(), 0.0) / static_cast(v.size()); - }; - const auto dim = file.GetDimension("/entry_1/data_1/data"); - constexpr double JOULE_PER_EV = 1.602176634e-19; - const std::string detector = "/entry_1/instrument_1/detector_1/"; - crystfel::Panel p; - p.name = "cxi"; - p.max_ss = static_cast(dim[1]) - 1; - p.max_fs = static_cast(dim[2]) - 1; - p.corner_x = -0.5 * static_cast(dim[2]); - p.corner_y = -0.5 * static_cast(dim[1]); - p.res = 1.0 / mean(detector + "x_pixel_size"); - p.clen = fmt::format("{}", file.Exists(detector + "distance") ? mean(detector + "distance") : 0.1); - if (gain) - p.adu_per_photon = *gain; - geom_ = {}; - geom_.panels.push_back(p); - geom_.data = "/entry_1/data_1/data"; - geom_.photon_energy = fmt::format("{}", mean("/entry_1/instrument_1/source_1/energy") / JOULE_PER_EV); - } - Open(files, "the CXI metadata"); -} - -void JFJochCrystFELReader::Open(const std::vector &files, const std::string &geometry_file) { - if (files.empty()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No file to read"); - files_ = files; - - std::unique_lock ul(hdf5_mutex); - - // Every frame of every file, in file order. - const auto pct = geom_.data.find('%'); - for (size_t f = 0; f < files_.size(); f++) { - open_files_.emplace_back(std::make_unique(files_[f])); - auto &file = *open_files_.back(); - if (pct != std::string::npos) { - const std::string parent = geom_.data.substr(0, pct); - for (const auto &group: file.FindLeafs(parent.empty() ? "/" : parent)) { - const auto dataset = ReplaceEvent(geom_.data, group); - if (file.Exists(dataset)) - frames_.push_back({f, dataset, -1, group}); - } - } else { - const auto dim = file.GetDimension(geom_.data); - if (dim.size() == 3) { - for (hsize_t i = 0; i < dim[0]; i++) - frames_.push_back({f, geom_.data, static_cast(i), ""}); - } else if (dim.size() == 2) - frames_.push_back({f, geom_.data, -1, ""}); - } - } - if (frames_.empty()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "No frames at " + geom_.data + " in the files given"); - - const auto dim = open_files_[frames_[0].file]->GetDimension(frames_[0].dataset); - frame_height_ = dim[dim.size() - 2]; - frame_width_ = dim[dim.size() - 1]; - - const auto &p0 = geom_.panels[0]; - const double clen_m = IsPath(p0.clen) ? MeanOverFirstFile(p0.clen) * 1e-3 : std::stod(p0.clen); - const double distance_mm = (clen_m + p0.coffset) * 1000.0; - if (geom_.photon_energy.empty()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - geometry_file + " states no photon_energy"); - const double photon_energy_eV = IsPath(geom_.photon_energy) ? MeanOverFirstFile(geom_.photon_energy) - : std::stod(geom_.photon_energy); - for (const auto &p: geom_.panels) - if (p.res != p0.res || p.clen != p0.clen || p.coffset != p0.coffset) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - geometry_file + ": panels at different distances or with different " - "pixel sizes cannot be put on one flat detector"); - - std::vector file_mask; - if (!geom_.mask.empty()) { - if (!geom_.mask_file.empty()) { - HDF5ReadOnlyFile mask_file(geom_.mask_file); - file_mask = mask_file.ReadVector(geom_.mask, {0, 0}, {frame_height_, frame_width_}); - } else - file_mask = open_files_[0]->ReadVector(geom_.mask, {0, 0}, {frame_height_, frame_width_}); - } - CheckDirectChunks(); - ul.unlock(); - - // CrystFEL's laboratory frame and the flat image's (column, row, beam) are both right-handed with - // the beam along z, so the point (X, Y) lands at column beam_x + X, row beam_y + Y - and CrystFEL's - // horizontal x stays the image's horizontal, where the polarisation correction expects it. A pixel - // centre is at its index, so beam_x is chosen to put the first panel's pixel centres on whole columns. - auto position = [](const crystfel::Panel &p, double fs, double ss) { - return std::pair{p.corner_x + fs * p.fs_x + ss * p.ss_x, p.corner_y + fs * p.fs_y + ss * p.ss_y}; - }; - double max_x = -1e9, min_x = 1e9, max_y = -1e9, min_y = 1e9; - for (const auto &p: geom_.panels) - for (double fs: {0.0, double(p.max_fs - p.min_fs + 1)}) - for (double ss: {0.0, double(p.max_ss - p.min_ss + 1)}) { - const auto [x, y] = position(p, fs, ss); - max_x = std::max(max_x, x); - min_x = std::min(min_x, x); - max_y = std::max(max_y, y); - min_y = std::min(min_y, y); - } - const auto [x0, y0] = position(p0, 0.5, 0.5); - const double beam_x = std::ceil(x0 - min_x) + 1 - x0; - const double beam_y = std::ceil(y0 - min_y) + 1 - y0; - width_ = static_cast(std::ceil(beam_x + max_x)) + 2; - height_ = static_cast(std::ceil(beam_y + max_y)) + 2; - - // Everything starts as a gap; the panels' pixels are put on top, less the bad regions. - std::vector mask(width_ * height_, 1u << PixelMask::ModuleGapPixelBit); - destination_.assign(frame_width_ * frame_height_, -1); - adu_per_photon_.assign(frame_width_ * frame_height_, 1.0f); - counts_in_photons_ = true; - std::optional saturation; - - for (const auto &p: geom_.panels) { - if (p.max_fs >= static_cast(frame_width_) || p.max_ss >= static_cast(frame_height_)) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Panel " + p.name + " lies outside the stored frame"); - double gain = 1.0; - if (p.adu_per_photon.has_value()) - gain = p.adu_per_photon.value(); - else if (p.adu_per_eV.has_value()) - gain = p.adu_per_eV.value() * photon_energy_eV; - else - counts_in_photons_ = false; - if (p.max_adu.has_value()) - saturation = std::min(saturation.value_or(INT64_MAX), p.max_adu.value() / gain); - if (p.no_index) - continue; - - for (int64_t ss = p.min_ss; ss <= p.max_ss; ss++) { - for (int64_t fs = p.min_fs; fs <= p.max_fs; fs++) { - const auto [x, y] = position(p, fs - p.min_fs + 0.5, ss - p.min_ss + 0.5); - bool bad = false; - for (const auto &b: geom_.bad) { - if (b.in_lab) - bad |= x >= b.min_x && x <= b.max_x && y >= b.min_y && y <= b.max_y; - else if (b.panel.empty() || b.panel == p.name) - bad |= fs >= b.min_fs && fs <= b.max_fs && ss >= b.min_ss && ss <= b.max_ss; - } - if (bad) - continue; - - const size_t src = ss * frame_width_ + fs; - const auto col = std::lround(beam_x + x); - const auto row = std::lround(beam_y + y); - const size_t dst = static_cast(row) * width_ + static_cast(col); - destination_[src] = static_cast(dst); - adu_per_photon_[src] = static_cast(gain); - - const uint32_t m = src < file_mask.size() ? file_mask[src] : 0; - const bool masked = (geom_.mask_bad != 0 && (m & geom_.mask_bad) != 0) - || (geom_.mask_good != 0 && (m & geom_.mask_good) != geom_.mask_good); - mask[dst] = masked ? (1u << PixelMask::ErrorPixelBit) : 0; - } - } - } - - dataset_ = std::make_shared(); - dataset_->experiment = default_experiment; - DetectorSetup detector = DetDECTRIS(static_cast(width_), static_cast(height_), - "CrystFEL geometry", {}); - detector.PixelSize_um(static_cast(1e6 / p0.res)); - // Nothing a parallax or sensor-efficiency correction needs is stated in a geometry file. - detector.SensorThickness_um(0); - // Where the geometry states no limit, the stored integer type still has one - a Cheetah file keeps - // its frames as int16, and a value clipped at 32767 is no measurement either. - if (!saturation.has_value() && !stored_float_ && stored_size_ > 0) { - const double type_max = std::exp2(8.0 * static_cast(stored_size_) - (stored_signed_ ? 1 : 0)) - 1; - saturation = type_max / *std::max_element(adu_per_photon_.begin(), adu_per_photon_.end()); - } - if (saturation.has_value()) - detector.SaturationLimit(static_cast(std::floor(saturation.value()))); - 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(beam_x)); - dataset_->experiment.BeamY_pxl(static_cast(beam_y)); - dataset_->experiment.DetectorDistance_mm(static_cast(distance_mm)); - dataset_->experiment.IncidentEnergy_keV(static_cast(photon_energy_eV / 1000.0)); - dataset_->experiment.ImagesPerTrigger(static_cast(frames_.size())); - dataset_->pixel_mask = std::make_shared(mask); - - Logger("CrystFELReader").Info("{} frames in {} files, {} panels; {} x {} image, beam at ({:.2f}, {:.2f}) px, " - "{:.3f} mm, {:.1f} eV, {}{}", frames_.size(), files_.size(), geom_.panels.size(), - width_, height_, beam_x, beam_y, distance_mm, photon_energy_eV, - counts_in_photons_ ? "converted to photons" : "in detector units (no gain stated)", - direct_chunks_ ? ", chunks decoded in parallel" : ""); - SetStartMessage(dataset_); -} - -void JFJochCrystFELReader::CheckDirectChunks() { - const auto &frame = frames_[0]; - HDF5DataSet dataset(*open_files_[frame.file], frame.dataset); - - const hid_t dcpl = H5Dget_create_plist(dataset.GetID()); - const int rank = frame.index >= 0 ? 3 : 2; - std::vector chunk(rank, 0); - bool ok = H5Pget_layout(dcpl) == H5D_CHUNKED && H5Pget_chunk(dcpl, rank, chunk.data()) == rank - && chunk[rank - 1] == frame_width_ && chunk[rank - 2] == frame_height_ && (rank == 2 || chunk[0] == 1); - bool deflate = false; - shuffled_ = false; - for (int i = 0; ok && i < H5Pget_nfilters(dcpl); i++) { - unsigned int flags; - size_t nelements = 0; - const H5Z_filter_t filter = H5Pget_filter2(dcpl, i, &flags, &nelements, nullptr, 0, nullptr, nullptr); - if (filter == H5Z_FILTER_DEFLATE) - deflate = true; - else if (filter == H5Z_FILTER_SHUFFLE) - shuffled_ = true; - else - ok = false; - } - H5Pclose(dcpl); - - const hid_t type = H5Dget_type(dataset.GetID()); - stored_size_ = H5Tget_size(type); - stored_float_ = H5Tget_class(type) == H5T_FLOAT; - stored_signed_ = H5Tget_sign(type) == H5T_SGN_2; - ok &= (H5Tget_class(type) == H5T_INTEGER && (stored_size_ == 2 || stored_size_ == 4)) - || (stored_float_ && stored_size_ == 4); - H5Tclose(type); - - direct_chunks_ = ok && deflate; -} - -bool JFJochCrystFELReader::CountsInPhotons() const { - return counts_in_photons_; -} - -uint64_t JFJochCrystFELReader::GetNumberOfImages() const { - return frames_.size(); -} - -void JFJochCrystFELReader::Close() { - std::unique_lock ul(hdf5_mutex); - open_files_.clear(); - files_.clear(); - frames_.clear(); - dataset_.reset(); -} - -template -void JFJochCrystFELReader::DecodeInto(int64_t image_number, Buffer &buffer) const { - const auto &frame = frames_.at(image_number); - std::vector stored(frame_width_ * frame_height_); - if (direct_chunks_) { - std::vector compressed; - { - std::unique_lock ul(hdf5_mutex); - HDF5DataSet dataset(*open_files_[frame.file], frame.dataset); - if (frame.index >= 0) - dataset.ReadDirectChunk(compressed, {static_cast(frame.index), 0, 0}); - else - dataset.ReadDirectChunk(compressed, {0, 0}); - } - const size_t n = stored.size(); - std::vector raw(n * stored_size_); - uLongf raw_size = raw.size(); - if (uncompress(raw.data(), &raw_size, compressed.data(), compressed.size()) != Z_OK || raw_size != raw.size()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Cannot inflate " + frame.dataset + " of " + files_[frame.file]); - if (shuffled_) { - std::vector tmp(raw.size()); - for (size_t i = 0; i < n; i++) - for (size_t b = 0; b < stored_size_; b++) - tmp[i * stored_size_ + b] = raw[b * n + i]; - raw.swap(tmp); - } - for (size_t i = 0; i < n; i++) { - const uint8_t *v = raw.data() + i * stored_size_; - if (stored_float_) - stored[i] = *reinterpret_cast(v); - else if (stored_size_ == 2) - stored[i] = stored_signed_ ? *reinterpret_cast(v) : *reinterpret_cast(v); - else - stored[i] = stored_signed_ ? *reinterpret_cast(v) : *reinterpret_cast(v); - } - } else { - std::unique_lock ul(hdf5_mutex); - HDF5DataSet dataset(*open_files_[frame.file], frame.dataset); - if (frame.index >= 0) - dataset.ReadVector(stored, {static_cast(frame.index), 0, 0}, {1, frame_height_, frame_width_}); - else - dataset.ReadVector(stored, {0, 0}, {frame_height_, frame_width_}); - } - - buffer.assign(width_ * height_ * sizeof(int32_t), 0); - auto out = reinterpret_cast(buffer.data()); - for (size_t i = 0; i < stored.size(); i++) - if (destination_[i] >= 0) - out[destination_[i]] = static_cast(std::lround(stored[i] / adu_per_photon_[i])); -} - -bool JFJochCrystFELReader::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; - DecodeInto(image_number, buffer); - message.image = CompressedImage(buffer.data(), buffer.size(), width_, height_, - CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION); - message.number = image_number; - return true; -} - -bool JFJochCrystFELReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) { - if (image_number < 0 || static_cast(image_number) >= frames_.size()) - return false; - DecodeInto(image_number, image.image_buffer); - image.image = CompressedImage(image.image_buffer.data(), image.image_buffer.size(), width_, height_, - CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION); - return true; -} - -std::vector JFJochCrystFELReader::ReadSpots(int64_t) const { - return {}; -} diff --git a/reader/JFJochCrystFELReader.h b/reader/JFJochCrystFELReader.h deleted file mode 100644 index eee018827..000000000 --- a/reader/JFJochCrystFELReader.h +++ /dev/null @@ -1,91 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#pragma once - -#include -#include -#include -#include - -#include "JFJochReader.h" -#include "CrystFELGeometry.h" - -class HDF5ReadOnlyFile; - -// Reads serial (stills) frames from HDF5 files laid out the way CrystFEL reads them, placed on the -// detector by a CrystFEL geometry file: Cheetah CXI files with the frames stacked along the first -// dimension, files holding one frame each, or files with one group per shot (data = /%/data). The -// panels are put on one flat image at their geometry-file positions, rounded to whole pixels, with -// everything between them and every bad region masked. Pixel values are converted to photons with the -// geometry's gain where it states one, and are otherwise left in detector units. -// -// One wavelength and one camera length hold for the whole run: where the geometry takes either from -// the image files, it is the mean over the first file's frames. -// -// A path is one file, or a directory, every .h5/.cxi file of which is read in name order. -class JFJochCrystFELReader : public JFJochReader { - struct Frame { - size_t file; - std::string dataset; - int64_t index; // along the first dimension of a 3D dataset; -1 for a 2D one - std::string group; // the shot group the geometry's % stands for, if any - }; - - crystfel::Geometry geom_; - std::vector files_; - std::vector> open_files_; - std::vector frames_; - std::shared_ptr dataset_; - - size_t frame_width_ = 0, frame_height_ = 0; // the stored frame - size_t width_ = 0, height_ = 0; // the flat image - // For every stored pixel, where it lands on the flat image (-1: nowhere) and what one photon reads. - std::vector destination_; - std::vector adu_per_photon_; - bool counts_in_photons_ = false; - - // Where every frame is stored as one deflated (and optionally byte-shuffled) chunk of a plain - // number type - the Cheetah outputs - the chunk is read raw under the HDF5 lock and decoded - // outside it, so decoding runs on every reading thread instead of one at a time. - bool direct_chunks_ = false; - bool shuffled_ = false; - bool stored_float_ = false, stored_signed_ = false; - size_t stored_size_ = 0; - - void CheckDirectChunks(); - - double MeanOverFirstFile(const std::string &path); - // The files at path, placed by geom_; geometry_file names where geom_ came from, for messages. - void Open(const std::vector &files, const std::string &geometry_file); - - bool LoadImage_i(std::shared_ptr &dataset, - DataMessage &message, - std::vector &buffer, - int64_t image_number, - bool update_dataset) override; - - template - void DecodeInto(int64_t image_number, Buffer &buffer) const; - -public: - JFJochCrystFELReader(); - ~JFJochCrystFELReader() override; - - // gain, where given, is the detector units per photon of every panel, whatever the geometry says. - // files are read in the order given, their frames one after another. - void ReadFiles(const std::vector &files, const std::string &geometry_file, - std::optional gain = std::nullopt); - // Cheetah CXI files, with no geometry file: one flat panel, placed by what the first file states, - // centred on the beam. - static bool CanReadCXI(const std::string &path); - void ReadCXI(const std::vector &files, std::optional gain = std::nullopt); - - // False where the geometry states no gain, so the images are in detector units, not photons. - [[nodiscard]] bool CountsInPhotons() const; - - [[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/JFJochXFELReader.cpp b/reader/JFJochXFELReader.cpp new file mode 100644 index 000000000..abae9e4d8 --- /dev/null +++ b/reader/JFJochXFELReader.cpp @@ -0,0 +1,841 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "JFJochXFELReader.h" +#include "CrystFELGeometry.h" + +#include +#include +#include +#include +#include +#include + +#include + +#include "../common/JFJochException.h" +#include "../common/Logger.h" +#include "../writer/HDF5Objects.h" + +namespace { + +constexpr double JOULE_PER_EV = 1.602176634e-19; +constexpr double PI_D = 3.14159265358979323846; +const std::string CXI_DATA = "/entry_1/data_1/data"; +const std::string CXI_DETECTOR = "/entry_1/instrument_1/detector_1/"; + +// The bits of a Cheetah pixel mask that mark a pixel unusable in every frame: invalid, hot, dead, +// shadowed, bad, missing and noisy (Cheetah's PIXEL_IS_* flags). The peak mask, saturation and the +// resolution limits are left alone. +constexpr uint32_t CHEETAH_BAD_BITS = 1u | 4u | 8u | 16u | 128u | 512u | 1024u; + +std::string ReplaceEvent(std::string path, const std::string &group) { + const auto pos = path.find('%'); + if (pos != std::string::npos) + path.replace(pos, 1, group); + return path; +} + +// A running mean, and the spread, of a value the geometry takes from the image files. +struct Mean { + double sum = 0, min = INFINITY, max = -INFINITY; + size_t n = 0; + void Add(double v) { sum += v; n++; min = std::min(min, v); max = std::max(max, v); } + double Get() const { return sum / static_cast(n); } +}; + +// Adds what one file states at path: one value per shot group (a path with %), a per-frame vector, or +// one scalar for the file. +void Accumulate(HDF5ReadOnlyFile &file, const std::string &path, const std::vector &groups, + Mean &mean) { + if (path.find('%') != std::string::npos) { + for (const auto &g: groups) + mean.Add(HDF5DataSet(file, ReplaceEvent(path, g)).ReadScalar()); + } else if (file.GetDimension(path).empty()) + mean.Add(HDF5DataSet(file, path).ReadScalar()); + else + for (const double v: file.ReadVector(path)) + mean.Add(v); +} + +} // namespace + +JFJochXFELReader::JFJochXFELReader() = default; + +JFJochXFELReader::~JFJochXFELReader() { + Close(); +} + +// The HDF5 filter IDs registered for bitshuffle and for plain LZ4. +constexpr H5Z_filter_t BITSHUFFLE_FILTER_ID = 32008; +constexpr H5Z_filter_t LZ4_FILTER_ID = 32004; + +std::shared_ptr JFJochXFELReader::RawFileAt(size_t file) const { + if (raw_files_.size() != files_.size()) { + raw_files_.assign(files_.size(), nullptr); + user_block_.assign(files_.size(), 0); + } + if (!raw_files_[file]) { + const hid_t fcpl = H5Fget_create_plist(FileAt(file).GetID()); + hsize_t user_block = 0; + if (fcpl >= 0) { + H5Pget_userblock(fcpl, &user_block); + H5Pclose(fcpl); + } + auto raw = std::make_shared(files_[file]); + if (!raw->IsOpen()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Cannot open " + files_[file]); + raw_files_[file] = raw; + user_block_[file] = user_block; + } + return raw_files_[file]; +} + +HDF5ReadOnlyFile &JFJochXFELReader::FileAt(size_t file) const { + auto it = std::find_if(open_files_.begin(), open_files_.end(), [&](const auto &e) { return e.first == file; }); + if (it == open_files_.end()) { + if (open_files_.size() >= MAX_OPEN_FILES) + open_files_.erase(open_files_.begin()); + open_files_.emplace_back(file, std::make_unique(files_[file])); + } else if (it != open_files_.end() - 1) + std::rotate(it, it + 1, open_files_.end()); + return *open_files_.back().second; +} + +// The files given, in the order given; a directory is not read (the caller expands a glob). +static void CheckFiles(const std::vector &files) { + if (files.empty()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No file to read"); +} + +void JFJochXFELReader::ReadFiles(const std::vector &files, const std::string &geometry_file, + std::optional gain) { + Close(); + CheckFiles(files); + std::vector ignored; + geom_ = crystfel::Read(geometry_file, ignored); + source_ = geometry_file; + if (!ignored.empty()) { + std::string list; + for (const auto &key: ignored) + list += (list.empty() ? "" : ", ") + key; + Logger("XFELReader").Info("{}: only the module geometry (and camera length) is taken; ignored: {}", + geometry_file, list); + } + + // What the geometry file no longer states comes from the data files themselves, where they are one of + // the layouts this reader knows: a Cheetah CXI file (photon energy, camera length, the shared mask + // and the per-shot one) or a SACLA Cheetah file (photon energy per shot, camera length, pixel mask). + // Otherwise the photon energy is the command line's to give, and the camera length too where the + // geometry has none. + { + std::unique_lock ul(hdf5_mutex); + HDF5ReadOnlyFile file(files[0]); + if (file.Exists(CXI_DATA)) { + if (file.Exists("/entry_1/instrument_1/source_1/energy")) + geom_.photon_energy = xfel::Value{"/entry_1/instrument_1/source_1/energy", 0, 1.0 / JOULE_PER_EV}; + if (file.Exists(CXI_DETECTOR + "mask_shared")) + geom_.masks.push_back({CXI_DETECTOR + "mask_shared", "", 0, CHEETAH_BAD_BITS}); + if (file.Exists("/entry_1/data_1/mask")) + geom_.masks.push_back({"/entry_1/data_1/mask", "", 0, CHEETAH_BAD_BITS}); + if (!geom_.panels[0].clen && file.Exists(CXI_DETECTOR + "distance")) + for (auto &p: geom_.panels) + p.clen = xfel::Value{CXI_DETECTOR + "distance", 0, 1.0}; + } else if (file.Exists("/metadata/distance_in_mm")) { + geom_.photon_energy = xfel::Value{"/%/photon_energy_ev", 0, 1.0}; + geom_.masks.push_back({"/metadata/pixelmask", "", 0, 0xffffffffu}); + if (!geom_.panels[0].clen) + for (auto &p: geom_.panels) + p.clen = xfel::Value{"/metadata/distance_in_mm", 0, 1e-3}; + } + } + for (auto &p: geom_.panels) { + if (!p.clen) + p.clen = xfel::Value{"", 0.1, 1.0}; // a placeholder: the distance is given on the command line + if (gain) + p.adu_per_photon = *gain; + } + Open(files); + Place(); +} + +// A Cheetah CXI file carries what one flat panel needs - the frames at /entry_1/data_1/data, the pixel +// size, the photon energy and a camera length, all in SI units - but not where the beam is. The panel is +// put with its middle on the beam, which is as good as no centre at all: the run measures it. +// +// Only a square frame is taken this way. A Cheetah file from a tiled detector (CSPAD, AGIPD, ePix, +// JUNGFRAU, MPCCD) holds its panels stacked in a slab, which read as one flat panel would place every +// pixel wrongly and say nothing; those slabs are not square, the single-panel CCDs are. +bool JFJochXFELReader::CanReadCXI(const std::string &path) { + try { + if (!std::filesystem::is_regular_file(path)) + return false; + std::unique_lock ul(hdf5_mutex); + HDF5ReadOnlyFile file(path); + if (!file.Exists(CXI_DATA)) + return false; + const auto dim = file.GetDimension(CXI_DATA); + return dim.size() == 3 && dim[1] == dim[2] + && file.Exists(CXI_DETECTOR + "x_pixel_size") + && file.Exists("/entry_1/instrument_1/source_1/energy"); + } catch (...) { + return false; + } +} + +std::optional JFJochXFELReader::CXIRefusal(const std::string &path) { + try { + if (!std::filesystem::is_regular_file(path) || CanReadCXI(path)) + return std::nullopt; + std::unique_lock ul(hdf5_mutex); + HDF5ReadOnlyFile file(path); + if (!file.Exists(CXI_DATA)) + return std::nullopt; + const auto dim = file.GetDimension(CXI_DATA); + if (dim.size() == 3 && dim[1] != dim[2]) + return fmt::format("{} is a Cheetah CXI file whose {} x {} frame is the panels of a tiled detector " + "stacked in a slab; their layout is not in the file - give it with " + "--crystfel-geometry", path, dim[2], dim[1]); + return fmt::format("{} is a CXI file without the pixel size and photon energy of a single-panel " + "detector; give its geometry with --crystfel-geometry", path); + } catch (...) { + return std::nullopt; + } +} + +void JFJochXFELReader::ReadCXI(const std::vector &files, std::optional gain) { + Close(); + CheckFiles(files); + std::vector dim; + bool has_mask = false; + { + std::unique_lock ul(hdf5_mutex); + HDF5ReadOnlyFile file(files[0]); + dim = file.GetDimension(CXI_DATA); + has_mask = file.Exists(CXI_DETECTOR + "mask_shared"); + const auto pixel = file.ReadVector(CXI_DETECTOR + "x_pixel_size"); + xfel::Panel p; + p.name = "cxi"; + p.max_ss = static_cast(dim[1]) - 1; + p.max_fs = static_cast(dim[2]) - 1; + p.corner_x = -0.5 * static_cast(dim[2]); + p.corner_y = -0.5 * static_cast(dim[1]); + p.res = 1.0 / pixel.at(0); + if (file.Exists(CXI_DETECTOR + "distance")) + p.clen = xfel::Value{CXI_DETECTOR + "distance", 0, 1.0}; + else + p.clen = xfel::Value{"", 0.1, 1.0}; // a placeholder: the distance is given or found + if (gain) + p.adu_per_photon = *gain; + geom_ = {}; + geom_.panels.push_back(p); + geom_.data = CXI_DATA; + geom_.photon_energy = xfel::Value{"/entry_1/instrument_1/source_1/energy", 0, 1.0 / JOULE_PER_EV}; + if (has_mask) + geom_.masks.push_back({CXI_DETECTOR + "mask_shared", "", 0, CHEETAH_BAD_BITS}); + } + source_ = "the CXI metadata"; + Open(files); + Place(); +} + +// A SACLA Cheetah file: one group per shot (tag-N) with the frame and the shot's photon energy, and +// /metadata stating where every sensor of the MPCCD lies (posx/posy in um, the in-plane angle - in +// degrees, whatever its name says - and the pixel size), the camera length and a pixel mask that +// already holds the sensors' edges and the rows behind each amplifier's shield. The sensors are stacked +// one after the other along the frame's slow axis. +bool JFJochXFELReader::CanReadPanelTable(const std::string &path) { + try { + if (!std::filesystem::is_regular_file(path)) + return false; + std::unique_lock ul(hdf5_mutex); + HDF5ReadOnlyFile file(path); + return file.Exists("/metadata/posx_in_um") && file.Exists("/metadata/posy_in_um") + && file.Exists("/metadata/angle_in_rad") && file.Exists("/metadata/pixelsizex_in_um") + && file.Exists("/metadata/distance_in_mm"); + } catch (...) { + return false; + } +} + +void JFJochXFELReader::ReadPanelTable(const std::vector &files, std::optional gain) { + Close(); + CheckFiles(files); + struct Metrology { + std::vector posx, posy, angle, pixel; + double distance_mm = 0; + bool operator==(const Metrology &) const = default; + }; + auto read_metrology = [](HDF5ReadOnlyFile &file) { + Metrology m; + m.posx = file.ReadVector("/metadata/posx_in_um"); + m.posy = file.ReadVector("/metadata/posy_in_um"); + m.angle = file.ReadVector("/metadata/angle_in_rad"); + m.pixel = file.ReadVector("/metadata/pixelsizex_in_um"); + m.distance_mm = HDF5DataSet(file, "/metadata/distance_in_mm").ReadScalar(); + return m; + }; + + Metrology m; + std::string first_group; + { + std::unique_lock ul(hdf5_mutex); + HDF5ReadOnlyFile file(files[0]); + m = read_metrology(file); + // The sensors of one run are placed once: every file has to state the same metrology. + for (size_t f = 1; f < files.size(); f++) { + HDF5ReadOnlyFile other(files[f]); + if (!(read_metrology(other) == m)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + files[f] + " places the sensors or the detector differently from " + + files[0] + "; read the runs separately"); + } + for (const auto &g: file.FindLeafs("/")) + if (file.Exists("/" + g + "/data")) { + first_group = g; + break; + } + if (first_group.empty()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No shot groups in " + files[0]); + } + const size_t nsensor = m.posx.size(); + if (nsensor == 0 || m.posy.size() != nsensor || m.angle.size() != nsensor || m.pixel.size() != nsensor) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + files[0] + ": /metadata does not state the same number of positions, angles and " + "pixel sizes"); + if (!(m.distance_mm > 0)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, files[0] + " states no camera length"); + + std::vector dim; + { + std::unique_lock ul(hdf5_mutex); + HDF5ReadOnlyFile file(files[0]); + dim = file.GetDimension("/" + first_group + "/data"); + } + if (dim.size() != 2 || dim[0] % nsensor != 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + files[0] + ": the frame does not hold " + std::to_string(nsensor) + " sensors stacked"); + const auto rows = static_cast(dim[0] / nsensor); + + geom_ = {}; + geom_.data = "/%/data"; + geom_.photon_energy = xfel::Value{"/%/photon_energy_ev", 0, 1.0}; + geom_.masks.push_back({"/metadata/pixelmask", "", 0, 0xffffffffu}); + // The placement the SACLA pipeline itself writes into its CrystFEL geometry: the corner of the + // sensor's first pixel at (-posx, posy) and its fast axis at the stated angle. + for (size_t i = 0; i < nsensor; i++) { + const double a = m.angle[i] * PI_D / 180.0; + xfel::Panel p; + p.name = "q" + std::to_string(i + 1); + p.min_ss = static_cast(i) * rows; + p.max_ss = p.min_ss + rows - 1; + p.max_fs = static_cast(dim[1]) - 1; + p.corner_x = -m.posx[i] / m.pixel[i]; + p.corner_y = m.posy[i] / m.pixel[i]; + p.fs_x = std::cos(a); + p.fs_y = std::sin(a); + p.ss_x = -std::sin(a); + p.ss_y = std::cos(a); + p.res = 1e6 / m.pixel[i]; + p.clen = xfel::Value{"", m.distance_mm * 1e-3, 1.0}; + if (gain) + p.adu_per_photon = *gain; + geom_.panels.push_back(p); + } + source_ = "the /metadata panel table"; + Open(files); + Place(); +} + +JFJochXFELReader::Layout JFJochXFELReader::ReadLayout(HDF5ReadOnlyFile &file, const Frame &frame) { + Layout l; + HDF5DataSet dataset(file, frame.dataset); + const auto dim = file.GetDimension(frame.dataset); + if (dim.size() < 2 || dim.size() > 3) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + frame.dataset + " is not a stack of frames or one frame"); + l.height = dim[dim.size() - 2]; + l.width = dim[dim.size() - 1]; + + const hid_t dcpl = H5Dget_create_plist(dataset.GetID()); + const int rank = static_cast(dim.size()); + std::vector chunk(rank, 0); + bool ok = H5Pget_layout(dcpl) == H5D_CHUNKED && H5Pget_chunk(dcpl, rank, chunk.data()) == rank + && chunk[rank - 1] == l.width && chunk[rank - 2] == l.height && (rank == 2 || chunk[0] == 1); + bool deflate = false, bitshuffle = false; + for (int i = 0; ok && i < H5Pget_nfilters(dcpl); i++) { + unsigned int flags; + size_t nelements = 0; + const H5Z_filter_t filter = H5Pget_filter2(dcpl, i, &flags, &nelements, nullptr, 0, nullptr, nullptr); + if (filter == H5Z_FILTER_DEFLATE) + deflate = true; + else if (filter == H5Z_FILTER_SHUFFLE) + l.shuffled = true; + else if (filter == BITSHUFFLE_FILTER_ID || filter == LZ4_FILTER_ID) + bitshuffle = true; + else + ok = false; + } + H5Pclose(dcpl); + // Bitshuffle/LZ4 (or zstd, or plain LZ4) is the one filter of its chunk: CompressedImage decodes it. + if (ok && bitshuffle && !deflate && !l.shuffled) + l.algorithm = HDF5Dcpl(dataset).GetCompression(); + + const hid_t type = H5Dget_type(dataset.GetID()); + l.stored_size = H5Tget_size(type); + l.stored_float = H5Tget_class(type) == H5T_FLOAT; + l.stored_signed = H5Tget_sign(type) == H5T_SGN_2; + ok &= ((H5Tget_class(type) == H5T_INTEGER && (l.stored_size == 2 || l.stored_size == 4)) + || (l.stored_float && l.stored_size == 4)) + && H5Tget_order(type) == H5T_ORDER_LE; + H5Tclose(type); + l.deflated = deflate; + l.direct_chunks = ok && (deflate || (bitshuffle && l.algorithm != CompressionAlgorithm::NO_COMPRESSION)); + return l; +} + +void JFJochXFELReader::Open(const std::vector &files) { + files_ = files; + geom_.masks.insert(geom_.masks.end(), extra_masks_.begin(), extra_masks_.end()); + + // Every frame of every file, in file order, with the mean of every value taken from the files. + Mean energy, clen; + const auto &p0 = geom_.panels.at(0); + const auto pct = geom_.data.find('%'); + for (const auto &p: geom_.panels) + if (p.res != p0.res || p.clen->path != p0.clen->path || p.clen->value != p0.clen->value) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + source_ + ": panels with different camera lengths or pixel sizes cannot be put " + "on one flat detector"); + const xfel::Value *energy_value = geom_.photon_energy ? &*geom_.photon_energy + : (geom_.wavelength ? &*geom_.wavelength : nullptr); + + std::unique_lock ul(hdf5_mutex); + for (size_t f = 0; f < files_.size(); f++) { + HDF5ReadOnlyFile file(files_[f]); + const size_t first = frames_.size(); + std::vector groups; + if (pct != std::string::npos) { + const std::string parent = geom_.data.substr(0, pct); + for (const auto &group: file.FindLeafs(parent.empty() ? "/" : parent)) { + const auto dataset = ReplaceEvent(geom_.data, group); + if (file.Exists(dataset)) { + frames_.push_back({f, dataset, -1, group}); + groups.push_back(group); + } + } + } else if (file.Exists(geom_.data)) { + const auto dim = file.GetDimension(geom_.data); + if (dim.size() == 3) { + for (hsize_t i = 0; i < dim[0]; i++) + frames_.push_back({f, geom_.data, static_cast(i), ""}); + } else if (dim.size() == 2) + frames_.push_back({f, geom_.data, -1, ""}); + } + if (frames_.size() == first) + continue; + + // Every file has to store its frames the way the first does: the decoding is set up once. + const Layout l = ReadLayout(file, frames_[first]); + if (first == 0) + layout_ = l; + else if (!(l == layout_)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + files_[f] + " stores its frames differently (size, type or compression) from " + + files_[0] + "; read them separately"); + + if (energy_value != nullptr && !energy_value->path.empty()) + Accumulate(file, energy_value->path, groups, energy); + if (!p0.clen->path.empty()) + Accumulate(file, p0.clen->path, groups, clen); + } + if (frames_.empty()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "No frames at " + geom_.data + " in the files given"); + frame_height_ = layout_.height; + frame_width_ = layout_.width; + + clen_m_ = p0.clen->path.empty() ? p0.clen->value : clen.Get() * p0.clen->scale; + double coffset_sum = 0; + for (const auto &p: geom_.panels) + coffset_sum += p.coffset; + distance_mm_ = (clen_m_ + coffset_sum / static_cast(geom_.panels.size())) * 1000.0; + // Nothing states the photon energy: it is the command line's to give (--wavelength). + photon_energy_eV_ = 0; + if (energy_value != nullptr) { + const double v = energy_value->path.empty() ? energy_value->value : energy.Get() * energy_value->scale; + photon_energy_eV_ = geom_.photon_energy ? v : WVL_1A_IN_KEV * 1000.0 / v; + if (!(photon_energy_eV_ > 0) || !std::isfinite(photon_energy_eV_)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + source_ + ": no usable photon energy in the files given"); + } + if (energy.n > 1) + Logger("XFELReader").Info("Photon energy {:.1f} eV, the mean of {} values from {:.1f} to {:.1f} (the run " + "is read at one wavelength)", photon_energy_eV_, energy.n, + geom_.photon_energy ? energy.min * energy_value->scale + : WVL_1A_IN_KEV * 1000.0 / (energy.max * energy_value->scale), + geom_.photon_energy ? energy.max * energy_value->scale + : WVL_1A_IN_KEV * 1000.0 / (energy.min * energy_value->scale)); +} + +void JFJochXFELReader::Place() { + std::vector file_mask(frame_width_ * frame_height_, 0); + std::vector masked_by_file(frame_width_ * frame_height_, 0); + { + std::unique_lock ul(hdf5_mutex); + frame_masks_.clear(); + for (const auto &m: geom_.masks) { + std::vector v; + if (m.file.empty() && FileAt(0).GetDimension(m.data).size() == 3) { + frame_masks_.push_back(m); // one plane per frame: applied as each frame is read + continue; + } + if (!m.file.empty()) { + HDF5ReadOnlyFile mask_file(m.file); + v = mask_file.ReadVector(m.data, {0, 0}, {frame_height_, frame_width_}); + } else + v = FileAt(0).ReadVector(m.data, {0, 0}, {frame_height_, frame_width_}); + for (size_t i = 0; i < v.size(); i++) + if (m.Bad(v[i])) + masked_by_file[i] = 1; + } + } + + // CrystFEL's laboratory frame and the flat image's (column, row, beam) are both right-handed with + // the beam along z, so the point (X, Y) lies at column beam_x + X, row beam_y + Y of the detector plane + // - and CrystFEL's horizontal x stays the image's horizontal, where the polarisation correction expects + // it. A pixel centre is at its index. A panel set nearer or further than the plane (its own coffset) + // is projected onto it along the ray from the sample, so each pixel keeps its scattering direction. + auto position = [this](const xfel::Panel &p, double fs, double ss) { + const double scale = distance_mm_ / ((clen_m_ + p.coffset) * 1000.0); + return std::pair{scale * (p.corner_x + fs * p.fs_x + ss * p.ss_x), + scale * (p.corner_y + fs * p.fs_y + ss * p.ss_y)}; + }; + // Each panel goes onto the image whole, as a block: its axes turned to the nearest image axis and its + // first pixel put at a whole pixel. Where it really lies - the fraction of a pixel, the small turn - is + // its module's placement (DetectorModules), so no pixel is moved on its own. + struct Block { + int64_t fx = 1, fy = 0, sx = 0, sy = 1; // image step per +1 fs / +1 ss + double x0 = 0, y0 = 0; // plane position of the first pixel's centre, beam at 0 + }; + auto snap = [](double x, double y, int64_t &ix, int64_t &iy) { + if (std::fabs(x) >= std::fabs(y)) { ix = x > 0 ? 1 : -1; iy = 0; } + else { ix = 0; iy = y > 0 ? 1 : -1; } + }; + std::vector blocks; + double min_x = 1e9, min_y = 1e9, max_x = -1e9, max_y = -1e9; + for (const auto &p: geom_.panels) { + Block bl; + snap(p.fs_x, p.fs_y, bl.fx, bl.fy); + snap(p.ss_x, p.ss_y, bl.sx, bl.sy); + if (bl.fx * bl.sy - bl.fy * bl.sx == 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Panel " + p.name + " has its fast and slow axes along the same line"); + std::tie(bl.x0, bl.y0) = position(p, 0.5, 0.5); + const double w = static_cast(p.max_fs - p.min_fs), h = static_cast(p.max_ss - p.min_ss); + for (double fs: {0.0, w}) + for (double ss: {0.0, h}) { + const double x = bl.x0 + fs * static_cast(bl.fx) + ss * static_cast(bl.sx); + const double y = bl.y0 + fs * static_cast(bl.fy) + ss * static_cast(bl.sy); + min_x = std::min(min_x, x); max_x = std::max(max_x, x); + min_y = std::min(min_y, y); max_y = std::max(max_y, y); + } + blocks.push_back(bl); + } + // The first panel's pixel centres on whole columns and rows, two pixels of margin all round. + const double beam_x = std::ceil(blocks[0].x0 - min_x) + 2 - blocks[0].x0; + const double beam_y = std::ceil(blocks[0].y0 - min_y) + 2 - blocks[0].y0; + width_ = static_cast(std::ceil(beam_x + max_x)) + 3; + height_ = static_cast(std::ceil(beam_y + max_y)) + 3; + + // Everything starts as a gap; the panels' pixels are put on top, less the bad regions. + std::vector mask(width_ * height_, 1u << PixelMask::ModuleGapPixelBit); + std::vector taken(width_ * height_, 0); + destination_.assign(frame_width_ * frame_height_, -1); + adu_per_photon_.assign(frame_width_ * frame_height_, 1.0f); + counts_in_photons_ = true; + std::optional saturation; + std::vector modules; + + for (size_t n = 0; n < geom_.panels.size(); n++) { + const auto &p = geom_.panels[n]; + const Block &bl = blocks[n]; + if (p.max_fs >= static_cast(frame_width_) || p.max_ss >= static_cast(frame_height_)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Panel " + p.name + " lies outside the stored frame"); + double gain = 1.0; + if (p.adu_per_photon.has_value()) + gain = p.adu_per_photon.value(); + else if (p.adu_per_eV.has_value()) + gain = p.adu_per_eV.value() * photon_energy_eV_; + else + counts_in_photons_ = false; + if (p.max_adu.has_value()) + saturation = std::min(saturation.value_or(INT64_MAX), p.max_adu.value() / gain); + if (p.no_index) + continue; + + // The image pixel the panel's first pixel lands on, and the block it fills. + const int64_t c0 = std::lround(beam_x + bl.x0), r0 = std::lround(beam_y + bl.y0); + const int64_t w = p.max_fs - p.min_fs, h = p.max_ss - p.min_ss; + const int64_t bx0 = c0 + std::min(0, bl.fx * w) + std::min(0, bl.sx * h); + const int64_t by0 = r0 + std::min(0, bl.fy * w) + std::min(0, bl.sy * h); + + // Its placement: the plane position of image pixel (bx0, by0), and the plane steps of +1 column + // and +1 row - the panel's own axes, scaled by its projection, taken through the block's turn. + const double scale = distance_mm_ / ((clen_m_ + p.coffset) * 1000.0); + DetectorModule m; + m.x0 = bx0; + m.y0 = by0; + m.width = std::abs(bl.fx * w + bl.sx * h) + 1; + m.height = std::abs(bl.fy * w + bl.sy * h) + 1; + // image (dc, dr) = fs (fx, fy) + ss (sx, sy); the inverse of a signed permutation is its transpose. + auto plane_of = [&](double dc, double dr) { + const double fs = dc * static_cast(bl.fx) + dr * static_cast(bl.fy); + const double ss = dc * static_cast(bl.sx) + dr * static_cast(bl.sy); + const auto [x, y] = position(p, fs + 0.5, ss + 0.5); + return std::pair{beam_x + x, beam_y + y}; + }; + std::tie(m.ox, m.oy) = plane_of(static_cast(bx0 - c0), static_cast(by0 - r0)); + m.fx = scale * (static_cast(bl.fx) * p.fs_x + static_cast(bl.sx) * p.ss_x); + m.fy = scale * (static_cast(bl.fx) * p.fs_y + static_cast(bl.sx) * p.ss_y); + m.sx = scale * (static_cast(bl.fy) * p.fs_x + static_cast(bl.sy) * p.ss_x); + m.sy = scale * (static_cast(bl.fy) * p.fs_y + static_cast(bl.sy) * p.ss_y); + modules.push_back(m); + + const int64_t e = p.mask_edge_pixels; + for (int64_t ss = p.min_ss; ss <= p.max_ss; ss++) { + for (int64_t fs = p.min_fs; fs <= p.max_fs; fs++) { + const int64_t col = c0 + (fs - p.min_fs) * bl.fx + (ss - p.min_ss) * bl.sx; + const int64_t row = r0 + (fs - p.min_fs) * bl.fy + (ss - p.min_ss) * bl.sy; + const size_t dst = static_cast(row) * width_ + static_cast(col); + if (taken[dst]) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Panel " + p.name + " overlaps another on the image"); + taken[dst] = 1; + if (fs < p.min_fs + e || fs > p.max_fs - e || ss < p.min_ss + e || ss > p.max_ss - e) + continue; + const auto [x, y] = position(p, fs - p.min_fs + 0.5, ss - p.min_ss + 0.5); + bool bad = false; + for (const auto &b: geom_.bad) { + if (b.in_lab) + bad |= x >= b.min_x && x <= b.max_x && y >= b.min_y && y <= b.max_y; + else if (b.panel.empty() || b.panel == p.name) + bad |= fs >= b.min_fs && fs <= b.max_fs && ss >= b.min_ss && ss <= b.max_ss; + } + if (bad) + continue; + + const size_t src = ss * frame_width_ + fs; + destination_[src] = static_cast(dst); + adu_per_photon_[src] = static_cast(gain); + mask[dst] = masked_by_file[src] ? (1u << PixelMask::ErrorPixelBit) : 0; + } + } + } + + dataset_ = std::make_shared(); + dataset_->experiment = default_experiment; + DetectorSetup detector = DetDECTRIS(static_cast(width_), static_cast(height_), + "placed by " + source_, {}); + detector.PixelSize_um(static_cast(1e6 / geom_.panels[0].res)); + // Nothing a parallax or sensor-efficiency correction needs is stated in these files. + detector.SensorThickness_um(0); + // Where nothing states a limit, the stored integer type still has one - a Cheetah file keeps its + // frames as int16, and a value clipped at 32767 is no measurement either. + if (!saturation.has_value() && !layout_.stored_float && layout_.stored_size > 0) { + const double type_max = std::exp2(8.0 * static_cast(layout_.stored_size) + - (layout_.stored_signed ? 1 : 0)) - 1; + saturation = type_max / *std::max_element(adu_per_photon_.begin(), adu_per_photon_.end()); + } + if (saturation.has_value()) + detector.SaturationLimit(static_cast(std::floor(saturation.value()))); + 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); + // Dark-subtracted counts go below zero; the images are handed out as signed int32. + dataset_->experiment.BitDepthImage(32); + dataset_->experiment.PixelSigned(true); + dataset_->experiment.BeamX_pxl(static_cast(beam_x)); + dataset_->experiment.BeamY_pxl(static_cast(beam_y)); + dataset_->experiment.DetectorDistance_mm(static_cast(distance_mm_)); + if (photon_energy_eV_ > 0) + dataset_->experiment.IncidentEnergy_keV(static_cast(photon_energy_eV_ / 1000.0)); + dataset_->experiment.ImagesPerTrigger(static_cast(frames_.size())); + // Each panel is a module, placed where the layout puts it; one panel square with the image needs none. + const DetectorModules placed(modules); + if (modules.size() >= 2 || placed.Moved()) + dataset_->experiment.Modules(placed); + dataset_->pixel_mask = std::make_shared(mask); + + Logger("XFELReader").Info("{} frames in {} files, {} panels placed by {}; {} x {} image, beam at ({:.2f}, " + "{:.2f}) px, {:.3f} mm, {:.1f} eV, {}{}", frames_.size(), files_.size(), + geom_.panels.size(), source_, width_, height_, beam_x, beam_y, distance_mm_, + photon_energy_eV_, + counts_in_photons_ ? "converted to photons" : "in detector units (no gain known)", + layout_.direct_chunks ? ", chunks decoded in parallel" : ""); + SetStartMessage(dataset_); +} + +void JFJochXFELReader::ApplyGain(float adu_per_photon) { + for (auto &p: geom_.panels) + p.adu_per_photon = adu_per_photon; + Place(); +} + +bool JFJochXFELReader::HasPhotonEnergy() const { + return photon_energy_eV_ > 0; +} + +bool JFJochXFELReader::CountsInPhotons() const { + return counts_in_photons_; +} + +uint64_t JFJochXFELReader::GetNumberOfImages() const { + return frames_.size(); +} + +void JFJochXFELReader::Close() { + std::unique_lock ul(hdf5_mutex); + open_files_.clear(); + raw_files_.clear(); + user_block_.clear(); + files_.clear(); + frames_.clear(); + dataset_.reset(); +} + +template +void JFJochXFELReader::DecodeInto(int64_t image_number, Buffer &buffer) const { + const auto &frame = frames_.at(image_number); + std::vector stored(frame_width_ * frame_height_); + if (layout_.direct_chunks) { + // Under the lock only where the chunk is (a lookup in its index); the read and the decoding not. + std::shared_ptr raw_file; + uint64_t address = 0, chunk_size = 0; + { + std::unique_lock ul(hdf5_mutex); + HDF5DataSet dataset(FileAt(frame.file), frame.dataset); + const hsize_t coord_3d[3] = {static_cast(frame.index), 0, 0}; + const hsize_t coord_2d[2] = {0, 0}; + unsigned filter_mask = 0; + haddr_t addr = HADDR_UNDEF; + hsize_t size = 0; + if (H5Dget_chunk_info_by_coord(dataset.GetID(), frame.index >= 0 ? coord_3d : coord_2d, &filter_mask, + &addr, &size) < 0 || addr == HADDR_UNDEF || size == 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "No stored chunk for " + frame.dataset + " of " + files_[frame.file]); + raw_file = RawFileAt(frame.file); + address = user_block_[frame.file] + addr; + chunk_size = size; + } + std::vector compressed(chunk_size); + raw_file->ReadAt(compressed.data(), compressed.size(), address); + + const size_t n = stored.size(); + const size_t size = layout_.stored_size; + std::vector raw(n * size); + if (layout_.deflated) { + uLongf raw_size = raw.size(); + if (uncompress(raw.data(), &raw_size, compressed.data(), compressed.size()) != Z_OK + || raw_size != raw.size()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Cannot inflate " + frame.dataset + " of " + files_[frame.file]); + } else { + const CompressedImage image(compressed.data(), compressed.size(), frame_width_, frame_height_, + CalcImageMode(size, layout_.stored_float, layout_.stored_signed), + layout_.algorithm); + std::vector decoded; + std::memcpy(raw.data(), image.GetUncompressedPtr(decoded), raw.size()); + } + if (layout_.shuffled) { + std::vector tmp(raw.size()); + for (size_t i = 0; i < n; i++) + for (size_t b = 0; b < size; b++) + tmp[i * size + b] = raw[b * n + i]; + raw.swap(tmp); + } + for (size_t i = 0; i < n; i++) { + const uint8_t *v = raw.data() + i * size; + if (layout_.stored_float) { + float x; + std::memcpy(&x, v, 4); + stored[i] = x; + } else if (size == 2) { + uint16_t x; + std::memcpy(&x, v, 2); + stored[i] = layout_.stored_signed ? static_cast(static_cast(x)) : x; + } else { + uint32_t x; + std::memcpy(&x, v, 4); + stored[i] = layout_.stored_signed ? static_cast(static_cast(x)) + : static_cast(x); + } + } + } else { + std::unique_lock ul(hdf5_mutex); + HDF5DataSet dataset(FileAt(frame.file), frame.dataset); + if (frame.index >= 0) + dataset.ReadVector(stored, {static_cast(frame.index), 0, 0}, {1, frame_height_, frame_width_}); + else + dataset.ReadVector(stored, {0, 0}, {frame_height_, frame_width_}); + } + + // A pixel a per-frame mask marks in this frame is handed out as no measurement. + std::vector masked; + if (!frame_masks_.empty() && frame.index >= 0) { + masked.assign(stored.size(), 0); + std::unique_lock ul(hdf5_mutex); + for (const auto &m: frame_masks_) { + const auto v = FileAt(frame.file).ReadVector(m.data, {static_cast(frame.index), 0, 0}, + {1, frame_height_, frame_width_}); + for (size_t i = 0; i < v.size(); i++) + if (m.Bad(v[i])) + masked[i] = 1; + } + } + + // And so is one whose value the layout flags as no measurement. + const auto flagged = [this](float v) { + return (geom_.flag_lessthan && v < *geom_.flag_lessthan) || (geom_.flag_morethan && v > *geom_.flag_morethan) + || (geom_.flag_equal && v == *geom_.flag_equal); + }; + + buffer.assign(width_ * height_ * sizeof(int32_t), 0); + auto out = reinterpret_cast(buffer.data()); + for (size_t i = 0; i < stored.size(); i++) + if (destination_[i] >= 0) + out[destination_[i]] = (!masked.empty() && masked[i]) || flagged(stored[i]) + ? INT32_MIN + : static_cast(std::lround(stored[i] / adu_per_photon_[i])); +} + +bool JFJochXFELReader::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; + DecodeInto(image_number, buffer); + message.image = CompressedImage(buffer.data(), buffer.size(), width_, height_, + CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION); + message.number = image_number; + return true; +} + +bool JFJochXFELReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) { + if (image_number < 0 || static_cast(image_number) >= frames_.size()) + return false; + DecodeInto(image_number, image.image_buffer); + image.image = CompressedImage(image.image_buffer.data(), image.image_buffer.size(), width_, height_, + CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION); + return true; +} + +std::vector JFJochXFELReader::ReadSpots(int64_t) const { + return {}; +} diff --git a/reader/JFJochXFELReader.h b/reader/JFJochXFELReader.h new file mode 100644 index 000000000..ec79bdfae --- /dev/null +++ b/reader/JFJochXFELReader.h @@ -0,0 +1,146 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include +#include + +#include "JFJochReader.h" +#include "DetectorLayout.h" +#include "HDF5ImageSource.h" + +class HDF5ReadOnlyFile; + +// Reads serial (stills) frames from HDF5 files written by the XFEL pipelines, and places their panels on +// one flat image: +// +// - a SACLA Cheetah file (one group per shot, /tag-N/data, the sensors stacked in the frame): the +// sensors are placed where the file's /metadata puts them, with its camera length and pixel mask, +// and the photon energy of every shot; +// - a Cheetah CXI file of a single-panel detector (the frames stacked at /entry_1/data_1/data): one +// flat panel with the pixel size, photon energy and camera length the file states, centred on the +// beam - the run measures the centre; +// - any of these, or files holding one frame each, laid out by a CrystFEL geometry file instead. +// +// Each panel is put on the flat image whole, as a block, its axes turned to the nearest image axis and its +// first pixel at a whole pixel, with everything between the panels and every masked or bad region masked. +// Where the panel really lies (the fraction of a pixel, a small turn) is its module's placement on the +// detector plane (DetectorModules), which the run maps its spots and predictions through. Pixel values are converted to photons +// with the gain where one is known (the geometry's, or one given or measured later), and are otherwise +// left in detector units. +// +// One wavelength and one camera length hold for the whole run: where the file states one per shot, it +// is the mean over every frame of the run. +// +// A path is one file, or a directory, every .h5/.cxi file of which is read in name order. Only a few +// files are held open at a time, so a run of thousands of files is read as well as one. +class JFJochXFELReader : public JFJochReader { + std::vector extra_masks_; + struct Frame { + size_t file; + std::string dataset; + int64_t index; // along the first dimension of a 3D dataset; -1 for a 2D one + std::string group; // the shot group the geometry's % stands for, if any + }; + + xfel::Layout geom_; + std::string source_; // where geom_ came from, for messages + std::vector files_; + std::vector frames_; + std::shared_ptr dataset_; + + // The most recently used files, held open; the least recent is closed past MAX_OPEN_FILES. + static constexpr size_t MAX_OPEN_FILES = 8; + mutable std::vector>> open_files_; + HDF5ReadOnlyFile &FileAt(size_t file) const; // under hdf5_mutex + std::shared_ptr RawFileAt(size_t file) const; // under hdf5_mutex + + size_t frame_width_ = 0, frame_height_ = 0; // the stored frame + size_t width_ = 0, height_ = 0; // the flat image + // For every stored pixel, where it lands on the flat image (-1: nowhere) and what one photon reads. + std::vector destination_; + std::vector adu_per_photon_; + bool counts_in_photons_ = false; + double clen_m_ = 0, distance_mm_ = 0, photon_energy_eV_ = 0; + // Masks stored with one plane per frame (a Cheetah CXI's per-shot mask), applied as each frame is read. + std::vector frame_masks_; + + // Where every frame is stored as one chunk of a plain little-endian number type - deflated (and + // optionally byte-shuffled), as Cheetah writes it, or bitshuffle/LZ4 (or zstd), as Jungfraujoch does - + // only the chunk's place in the file is asked of HDF5 under its lock; the bytes are read and decoded + // outside it, on every reading thread at once, the way the NXmx reader reads its images + // (HDF5ImageSource: the chunk lookup, its RawFile, and CompressedImage's decoders). + struct Layout { + uint64_t height = 0, width = 0; + bool direct_chunks = false, shuffled = false, deflated = false; + CompressionAlgorithm algorithm = CompressionAlgorithm::NO_COMPRESSION; // where not deflated + bool stored_float = false, stored_signed = false; + size_t stored_size = 0; + bool operator==(const Layout &) const = default; + }; + Layout layout_; + // Positional-read handles on the files, opened as first needed (under hdf5_mutex), with the size of + // each file's user block, which HDF5 chunk addresses count from. + mutable std::vector> raw_files_; + mutable std::vector user_block_; + static Layout ReadLayout(HDF5ReadOnlyFile &file, const Frame &frame); + + // The files, in the order given, every frame in them, and the run's mean of every value the geometry + // takes from the files. + void Open(const std::vector &files); + // Places the panels, given the gains, on the flat image; builds the dataset. + void Place(); + + bool LoadImage_i(std::shared_ptr &dataset, + DataMessage &message, + std::vector &buffer, + int64_t image_number, + bool update_dataset) override; + + template + void DecodeInto(int64_t image_number, Buffer &buffer) const; + +public: + JFJochXFELReader(); + ~JFJochXFELReader() override; + + // gain, where given, is the detector units per photon of every panel, whatever the geometry says. + // files are read in the order given, their frames one after another. + void ReadFiles(const std::vector &files, const std::string &geometry_file, + std::optional gain = std::nullopt); + + // Cheetah CXI files of a single-panel detector: one flat panel, placed by + // what the file states, centred on the beam. + static bool CanReadCXI(const std::string &path); + // Why a Cheetah CXI file is not read without a geometry file (a tiled detector's slab), or nothing + // if the path is not one. + static std::optional CXIRefusal(const std::string &path); + void ReadCXI(const std::vector &files, std::optional gain = std::nullopt); + + // Cheetah HDF5 files with one group per shot and the panels' placement in /metadata (posx/posy in um, the + // in-plane angle, pixel size, camera length) - as SACLA writes them - placed by that table. + static bool CanReadPanelTable(const std::string &path); + void ReadPanelTable(const std::vector &files, std::optional gain = std::nullopt); + + // Mask arrays in other files, in the frames' own layout (--mask-good / --mask-bad), applied with those + // the data files state; set before reading. + void ExtraMasks(const std::vector &masks) { extra_masks_ = masks; } + + // Converts to photons with this gain from here on (detector units per photon, every panel). + void ApplyGain(float adu_per_photon); + + // False where no gain is known, so the images are in detector units, not photons. + [[nodiscard]] bool CountsInPhotons() const; + // False where nothing read states the photon energy (a geometry file over frames of no known layout): + // the caller gives it. + [[nodiscard]] bool HasPhotonEnergy() const; + + [[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/JFJochXTCReader.cpp b/reader/JFJochXTCReader.cpp index 14479da59..772f16b28 100644 --- a/reader/JFJochXTCReader.cpp +++ b/reader/JFJochXTCReader.cpp @@ -3,43 +3,29 @@ #include "JFJochXTCReader.h" +#include #include -#include -#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 {}; +// The pixel size of a Rayonix image, from the model name - an MX340-HS is 340 mm across - and the +// image's width. Only where the width times the binning is a whole readout of 3840, 5760 or 7680 +// pixels, i.e. no region of interest was read out. +std::optional RayonixPixel_um(const std::string &model, uint32_t binning, uint32_t width) { + if (!model.starts_with("MX") || binning == 0) + return std::nullopt; + const auto digits = model.substr(2, model.find_first_not_of("0123456789", 2) - 2); + if (digits.empty()) + return std::nullopt; + const uint64_t native = static_cast(width) * binning; + if (native != 3840 && native != 5760 && native != 7680) + return std::nullopt; + return std::stod(digits) * 1000.0 / width; } } // namespace @@ -48,108 +34,92 @@ 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::vector &streams) { - frames_.clear(); + Close(); + Logger logger("XTCReader"); + std::vector all; + std::vector energies; + std::string model; + uint32_t binning = 0; + size_t damaged = 0; for (const auto &stream: streams) { - auto f = xtc::IndexFrames(stream); - frames_.insert(frames_.end(), f.begin(), f.end()); + auto index = xtc::IndexFile(stream); + all.insert(all.end(), index.frames.begin(), index.frames.end()); + energies.insert(energies.end(), index.photon_energy_eV.begin(), index.photon_energy_eV.end()); + if (!index.rayonix_model.empty()) { + model = index.rayonix_model; + binning = index.rayonix_binning; + } + damaged += index.damaged; + if (index.truncated) + logger.Warning("{} ends inside its last event, which is left out", stream); } - if (frames_.empty()) + + // The diffraction detector is the camera with the largest image. Another camera of the same size + // would make that choice a guess, so it is refused; an event without it is left out. + std::map area; + for (const auto &f: all) + area[f.source] = static_cast(f.width) * f.height; + if (area.empty()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No camera images found in the XTC streams given"); + const auto largest = std::max_element(area.begin(), area.end(), + [](const auto &a, const auto &b) { return a.second < b.second; }); + if (std::count_if(area.begin(), area.end(), [&](const auto &a) { return a.second == largest->second; }) > 1) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "The XTC streams hold images of the same size from more than one camera; which is the " + "diffraction detector is not clear"); + for (const auto &f: all) + if (f.source == largest->first) + frames_.push_back(f); 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"); + if (damaged > 0) + logger.Warning("{} damaged, compressed or unknown-version images left out", damaged); + + const auto &f0 = frames_.front(); + if (!energies.empty()) + photon_energy_eV_ = std::accumulate(energies.begin(), energies.end(), 0.0) / static_cast(energies.size()); + pixel_um_ = RayonixPixel_um(model, binning, f0.width); - const std::string expt = FindExpt(frames_.front().file); dataset_ = std::make_shared(); dataset_->experiment = default_experiment; - has_geometry_ = !expt.empty(); - if (has_geometry_) { - Logger("XTCReader").Info("{} images in XTC, geometry from {}", frames_.size(), expt); - ReadGeometry(expt); - } else { - // XTC states no geometry. One flat panel, square-on, with its middle on the beam - the run - // measures the centre - and the pixel size, distance and wavelength to be given by the caller. - Logger("XTCReader").Info("{} images in XTC, no .expt beside them: one flat panel centred on the beam", - frames_.size()); - const auto &f0 = frames_.front(); - DetectorSetup detector = DetDECTRIS(f0.width, f0.height, "XTC camera", {}); - detector.SensorThickness_um(0); - detector.SaturationLimit(65535); - 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(0.5f * static_cast(f0.width)); - dataset_->experiment.BeamY_pxl(0.5f * static_cast(f0.height)); - } + DetectorSetup detector = DetDECTRIS(f0.width, f0.height, model.empty() ? "XTC camera" : model, {}); + // 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); + detector.MinFrameTime(std::chrono::microseconds(0)); + detector.MinCountTime(std::chrono::microseconds(0)); + detector.ReadOutTime(std::chrono::nanoseconds(0)); + if (pixel_um_) + detector.PixelSize_um(static_cast(*pixel_um_)); + + dataset_->experiment.Detector(detector); + dataset_->experiment.BeamX_pxl(0.5f * static_cast(f0.width)); + dataset_->experiment.BeamY_pxl(0.5f * static_cast(f0.height)); + if (photon_energy_eV_) + dataset_->experiment.IncidentEnergy_keV(static_cast(*photon_energy_eV_ / 1000.0)); + logger.Info("{} images in XTC ({}); {}; {}; one flat panel centred on the beam", frames_.size(), + model.empty() ? "camera" : model, + photon_energy_eV_ ? fmt::format("photon energy {:.1f} eV, the mean of {} beam records", + *photon_energy_eV_, energies.size()) + : std::string("no beam records"), + pixel_um_ ? fmt::format("{:.2f} um pixels ({}x{} binning)", *pixel_um_, binning, binning) + : std::string("pixel size not stated")); + dataset_->error_value = -1; dataset_->experiment.ImagesPerTrigger(static_cast(frames_.size())); - dataset_->pixel_mask = std::make_shared(frames_.front().width, frames_.front().height); + dataset_->pixel_mask = std::make_shared(f0.width, f0.height); + SetStartMessage(dataset_); +} +void JFJochXTCReader::ApplyGain(float adu_per_photon) { + adu_per_photon_ = adu_per_photon; + counts_in_photons_ = true; + dataset_->experiment.Detector().SaturationLimit(static_cast(std::floor(65535.0f / adu_per_photon))); SetStartMessage(dataset_); } @@ -160,6 +130,10 @@ uint64_t JFJochXTCReader::GetNumberOfImages() const { void JFJochXTCReader::Close() { frames_.clear(); dataset_.reset(); + photon_energy_eV_.reset(); + pixel_um_.reset(); + adu_per_photon_ = 1.0f; + counts_in_photons_ = false; } template @@ -169,7 +143,11 @@ CompressedImage JFJochXTCReader::DecodeInto(int64_t image_number, Buffer &buffer 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); + auto out = reinterpret_cast(buffer.data()); + xtc::ReadFrame(f, out, scratch); + if (adu_per_photon_ != 1.0f) + for (size_t i = 0; i < static_cast(f.width) * f.height; i++) + out[i] = static_cast(std::lround(static_cast(out[i]) / adu_per_photon_)); return CompressedImage(buffer.data(), buffer.size(), f.width, f.height, CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION); } diff --git a/reader/JFJochXTCReader.h b/reader/JFJochXTCReader.h index faaa8813e..8da0f99e8 100644 --- a/reader/JFJochXTCReader.h +++ b/reader/JFJochXTCReader.h @@ -3,6 +3,7 @@ #pragma once +#include #include #include @@ -14,19 +15,21 @@ // (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. Where a DIALS experiment list (.expt) lies beside the XTC files, or -// in their parent directory - the geometry a deposition ships for cctbx.xfel - it is read; otherwise -// the panel is placed square-on with its middle on the beam, the run measures the centre, and the -// pixel size, distance and wavelength are given on the command line. +// What the XTC states is used: the photon energy of every event's electron-beam record (the run is +// read at their mean) and, for a Rayonix, the pixel size from its model and binning. The distance is +// not recorded, and the beam centre is not either: the panel is placed square-on with its middle on +// the beam, the run measures the centre, and the distance (and whatever the XTC did not state) is +// given on the command line. // -// 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. +// The images are stills: no goniometer. Each is decoded on demand, the camera's stated pixel offset +// taken off. A damaged or compressed image is left out, and so is an event the file ends inside. class JFJochXTCReader : public JFJochReader { std::vector frames_; std::shared_ptr dataset_; - bool has_geometry_ = false; + std::optional photon_energy_eV_; + std::optional pixel_um_; + float adu_per_photon_ = 1.0f; + bool counts_in_photons_ = false; bool LoadImage_i(std::shared_ptr &dataset, DataMessage &message, @@ -37,8 +40,6 @@ class JFJochXTCReader : public JFJochReader { template CompressedImage DecodeInto(int64_t image_number, Buffer &buffer, std::vector &scratch) const; - void ReadGeometry(const std::string &expt); - public: ~JFJochXTCReader() override = default; @@ -46,9 +47,14 @@ public: // The event streams, read in the order given, their frames one after another. void ReadFiles(const std::vector &streams); - // Whether a DIALS .expt beside the streams gave the geometry; without one the pixel size, the - // distance and the wavelength are the caller's to give. - [[nodiscard]] bool HasGeometry() const { return has_geometry_; } + + // What the XTC stated; anything it did not is the caller's to give. + [[nodiscard]] bool HasWavelength() const { return photon_energy_eV_.has_value(); } + [[nodiscard]] bool HasPixelSize() const { return pixel_um_.has_value(); } + + // Converts to photons with this gain from here on (detector units per photon). + void ApplyGain(float adu_per_photon); + [[nodiscard]] bool CountsInPhotons() const { return counts_in_photons_; } [[nodiscard]] uint64_t GetNumberOfImages() const override; void Close() override; diff --git a/reader/NXmxModules.cpp b/reader/NXmxModules.cpp new file mode 100644 index 000000000..ff07c9be2 --- /dev/null +++ b/reader/NXmxModules.cpp @@ -0,0 +1,68 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "NXmxModules.h" + +#include + +#include "../common/JFJochException.h" + +namespace { + +// A translation dataset as a vector in pixels of the internal frame (x = column, y = row, z = beam). +Coord ReadTranslation_px(HDF5Object &file, const std::string &path, double pixel_size_m, const RotMatrix &delta) { + HDF5DataSet ds(file, path); + std::vector value; + ds.ReadVector(value); + const auto vec = ds.ReadAttrVec("vector"); + if (value.empty() || vec.size() != 3) + throw JFJochException(JFJochExceptionCategory::HDF5, path + ": not a translation"); + double scale = 1.0; + if (ds.AttrExists("units")) { + const std::string units = ds.ReadAttrStr("units"); + if (units == "mm") + scale = 1e-3; + else if (units == "um" || units == "microns") + scale = 1e-6; + } + const double len_px = value[0] * scale / pixel_size_m; + // McStas -> internal is a half turn about z; the discrete image orientation is undone after it. + const Coord mcstas(static_cast(-vec[0]), static_cast(-vec[1]), static_cast(vec[2])); + return delta.transpose() * (static_cast(len_px) * mcstas); +} + +} // namespace + +DetectorModules ReadNXmxModules(HDF5Object &file, double pixel_size_m, const DetectorOrientation &orientation) { + const RotMatrix delta = orientation.Matrix(); + std::vector out; + for (size_t n = 0;; n++) { + const std::string base = "/entry/instrument/detector/module_" + std::to_string(n) + "/"; + if (!file.Exists(base + "data_origin")) + break; + const auto origin = file.ReadVector(base + "data_origin"); + const auto size = file.ReadVector(base + "data_size"); + if (origin.size() != 2 || size.size() != 2) + throw JFJochException(JFJochExceptionCategory::HDF5, base + ": only two-dimensional modules are read"); + const Coord o = ReadTranslation_px(file, base + "module_offset", pixel_size_m, delta); + const Coord f = ReadTranslation_px(file, base + "fast_pixel_direction", pixel_size_m, delta); + const Coord s = ReadTranslation_px(file, base + "slow_pixel_direction", pixel_size_m, delta); + if (std::fabs(o.z) > 1e-3 || std::fabs(f.z) > 1e-6 || std::fabs(s.z) > 1e-6) + throw JFJochException(JFJochExceptionCategory::HDF5, base + ": a module out of the detector plane is not read"); + DetectorModule m; + m.y0 = origin[0]; + m.x0 = origin[1]; + m.height = size[0]; + m.width = size[1]; + m.ox = o.x; + m.oy = o.y; + m.fx = f.x; + m.fy = f.y; + m.sx = s.x; + m.sy = s.y; + out.push_back(m); + } + if (out.size() < 2) + return {}; + return DetectorModules(out); +} diff --git a/reader/NXmxModules.h b/reader/NXmxModules.h new file mode 100644 index 000000000..d68ffd0b8 --- /dev/null +++ b/reader/NXmxModules.h @@ -0,0 +1,14 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include "../common/DetectorModules.h" +#include "../common/DetectorOrientation.h" +#include "../writer/HDF5Objects.h" + +// The modules an NXmx file states, as WriteNXmxModules (writer/NXmxModules.h) writes them: groups +// /entry/instrument/detector/module_0, module_1, ... Their offsets are taken relative to the image's pixel +// (0, 0) and in the detector plane; a module standing out of that plane (an offset or pixel direction +// with a component along the detector normal) is refused. Fewer than two such groups: none. +DetectorModules ReadNXmxModules(HDF5Object &file, double pixel_size_m, const DetectorOrientation &orientation); diff --git a/reader/XTC.cpp b/reader/XTC.cpp index f8b5779aa..8429114e8 100644 --- a/reader/XTC.cpp +++ b/reader/XTC.cpp @@ -4,10 +4,9 @@ #include "XTC.h" #include -#include #include #include -#include +#include #include "../common/JFJochException.h" @@ -18,24 +17,30 @@ 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 ID_EBEAM = 15; +constexpr uint32_t ID_RAYONIX_CONFIG = 73; constexpr uint32_t TRANSITION_L1ACCEPT = 12; constexpr uint32_t TRANSITION_LAST = 12; +constexpr size_t EBEAM_V7_SIZE = 4 + 20 * 8; +constexpr size_t EBEAM_V7_PHOTON_ENERGY = 4 + 17 * 8; -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")); +std::ifstream Open(const std::string &path) { + std::ifstream f(path, std::ios::binary); if (!f) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Cannot open " + path); return f; } +bool ReadAt(std::ifstream &f, uint64_t pos, void *dst, size_t n) { + f.clear(); + f.seekg(static_cast(pos)); + f.read(static_cast(dst), static_cast(n)); + return f.gcount() == static_cast(n); +} + uint32_t Word(const uint8_t *p, size_t i) { uint32_t v; - memcpy(&v, p + 4 * i, 4); + std::memcpy(&v, p + 4 * i, 4); return v; } @@ -50,27 +55,55 @@ 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) { +// Walks the Xtcs in [begin, end) of the file, adding what it knows to index: images and the beam only +// from an event (an L1Accept), the configuration from any transition. +void Walk(std::ifstream &f, const std::string &file, uint64_t begin, uint64_t end, bool event, xtc::Index &index) { 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) + if (!ReadAt(f, pos, h, XTC_HEADER)) return; + const uint32_t damage = Word(h, 0); + const uint64_t source = (static_cast(Word(h, 1)) << 32) | Word(h, 2); const uint32_t type = Word(h, 3) & 0xffff; + const uint32_t version = (Word(h, 3) >> 16) & 0x7fff; + const bool compressed = (Word(h, 3) >> 31) != 0; const uint32_t extent = Word(h, 4); if (extent < XTC_HEADER || pos + extent > end) return; + const uint64_t payload = pos + XTC_HEADER, size = extent - XTC_HEADER; + if (type == ID_XTC) - FindFrame(f, pos + XTC_HEADER, pos + extent, best); - else if (type == ID_FRAME && extent >= XTC_HEADER + FRAME_HEADER) { + Walk(f, file, payload, pos + extent, event, index); + else if (type == ID_FRAME && event) { uint8_t fh[FRAME_HEADER]; - if (fread(fh, 1, FRAME_HEADER, f) == FRAME_HEADER) { + if (size >= FRAME_HEADER && ReadAt(f, payload, fh, 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}; + // Only a whole, undamaged 16-bit image of the version described above is read. + if (depth > 8 && depth <= 16 && npixel > 0) { + if (damage != 0 || compressed || version != 1 || FRAME_HEADER + 2 * npixel > size) + index.damaged++; + else + index.frames.push_back({file, payload + FRAME_HEADER, width, height, Word(fh, 3), source}); + } + } + } else if (type == ID_EBEAM && event && version == 7 && !compressed && damage == 0 && size >= EBEAM_V7_SIZE) { + uint8_t eb[EBEAM_V7_SIZE]; + double energy = 0; + if (ReadAt(f, payload, eb, EBEAM_V7_SIZE) && Word(eb, 0) == 0) { + std::memcpy(&energy, eb + EBEAM_V7_PHOTON_ENERGY, 8); + if (energy > 100.0 && energy < 100000.0) + index.photon_energy_eV.push_back(energy); + } + } else if (type == ID_RAYONIX_CONFIG && version == 2 && !compressed && size >= 2 && size < 4096) { + std::vector cfg(size); + if (ReadAt(f, payload, cfg.data(), size)) { + index.rayonix_binning = cfg[0] == cfg[1] ? cfg[0] : 0; + const std::string text(cfg.begin(), cfg.end()); + const auto mx = text.find("MX"); + if (mx != std::string::npos) + index.rayonix_model = text.substr(mx, text.find('\0', mx) - mx); } } pos += extent; @@ -85,48 +118,45 @@ bool CanRead(const std::string &path) { std::error_code ec; 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")); + std::ifstream f(path, std::ios::binary); uint8_t h[DGRAM_HEADER]; - return f && fread(h, 1, DGRAM_HEADER, f.get()) == DGRAM_HEADER && PlausibleFirstDatagram(h); + return f && ReadAt(f, 0, h, DGRAM_HEADER) && PlausibleFirstDatagram(h); } -std::vector IndexFrames(const std::string &file) { - File f = Open(file); - std::vector out; +Index IndexFile(const std::string &file) { + std::ifstream f = Open(file); + const uint64_t file_size = std::filesystem::file_size(file); + Index index; 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) { + while (ReadAt(f, pos, h, 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); - } + // A partly downloaded or partly written stream ends inside its last datagram. + if (end > file_size) { + index.truncated = true; + break; } + Walk(f, file, payload, end, Transition(h) == TRANSITION_L1ACCEPT, index); pos = end; } - return out; + return index; } 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()) + std::ifstream f = Open(frame.file); + if (!ReadAt(f, frame.offset, scratch.data(), 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; + std::memcpy(&v, scratch.data() + 2 * i, 2); + out[i] = static_cast(v) - static_cast(frame.pixel_offset); } } diff --git a/reader/XTC.h b/reader/XTC.h index 222e38ff1..156db6d44 100644 --- a/reader/XTC.h +++ b/reader/XTC.h @@ -14,32 +14,48 @@ // 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. +// five 32-bit words: damage, the source (two words), the type and version (bit 31: compressed), 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 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. +// - a camera image, Id_Frame (Camera::FrameV1): width, height, bit depth and an offset added to every +// pixel, 32 bits each, then the pixels row by row; +// - the electron beam, Id_EBeam version 7 (Bld::BldDataEBeamV7): a damage mask, then 20 doubles, the +// 18th the photon energy in eV; +// - a Rayonix configuration, Id_RayonixConfig version 2: the binning along fast and slow (a byte each) +// and, further on, the device name ("MX340-HS:..."). +// +// 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 + uint64_t offset = 0; // of the first pixel uint32_t width = 0; uint32_t height = 0; + uint32_t pixel_offset = 0; // added to every pixel by the camera; taken off when read + uint64_t source = 0; // the camera, as the Xtc names it +}; + +// What one stream holds. +struct Index { + std::vector frames; // every undamaged 16-bit image of every event + std::vector photon_energy_eV; // of every event with an undamaged electron-beam record + std::string rayonix_model; // the Rayonix configuration's device name, if any + uint32_t rayonix_binning = 0; + size_t damaged = 0; // images left out: damaged, compressed or of another version + bool truncated = false; // the file ends inside its last datagram, which is left out }; // True if the path is an XTC file - by its first datagram, not only its name. The "smalldata" copies // (.smd.xtc) carry the events without their images and are not read. bool CanRead(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); +Index IndexFile(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); diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 448975f01..e0696c1b3 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -5,6 +5,8 @@ ADD_LIBRARY(JFJochRugnux STATIC Rugnux.cpp RugnuxPipeline.h RugnuxStills.cpp + ModuleDiagnostics.cpp + ModuleDiagnostics.h RugnuxCalibration.cpp RugnuxPreScan.cpp RugnuxPasses.cpp diff --git a/rugnux/ModuleDiagnostics.cpp b/rugnux/ModuleDiagnostics.cpp new file mode 100644 index 000000000..97fb691d6 --- /dev/null +++ b/rugnux/ModuleDiagnostics.cpp @@ -0,0 +1,365 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "ModuleDiagnostics.h" + +#include +#include +#include + +#include + +#include "../common/JFJochMath.h" +#include "../common/Logger.h" +#include "../image_analysis/SensorAbsorption.h" + +namespace { + +constexpr double RAD2DEG = 180.0 / PI; + +struct Fit { + bool ok = false; + std::vector modules; + GeometryRefinerResult r; +}; + +Fit Place(const DiffractionGeometry &geom, const UnitCell &cell, const std::vector &frames, + const std::vector &start, GeometryRefinerSettings st, bool refine) { + st.modules = start; + st.refine_modules = refine; + Fit f; + f.r = RefineGlobalGeometry(geom, cell, frames, st); + f.ok = f.r.ok; + f.modules = f.r.modules; + return f; +} + +// Weighted straight line y = a + b x, with b's standard error and the weighted correlation. +struct Line { double a = 0, b = 0, se_b = 0, corr = 0; }; + +Line FitLine(const std::vector &x, const std::vector &y, const std::vector &w) { + double sw = 0, sx = 0, sy = 0; + for (size_t i = 0; i < x.size(); i++) { sw += w[i]; sx += w[i] * x[i]; sy += w[i] * y[i]; } + Line l; + if (sw <= 0 || x.size() < 3) + return l; + const double mx = sx / sw, my = sy / sw; + double sxx = 0, syy = 0, sxy = 0; + for (size_t i = 0; i < x.size(); i++) { + sxx += w[i] * (x[i] - mx) * (x[i] - mx); + syy += w[i] * (y[i] - my) * (y[i] - my); + sxy += w[i] * (x[i] - mx) * (y[i] - my); + } + if (sxx <= 0) + return l; + l.b = sxy / sxx; + l.a = my - l.b * mx; + l.corr = syy > 0 ? sxy / std::sqrt(sxx * syy) : 0; + double rss = 0; + for (size_t i = 0; i < x.size(); i++) { + const double e = y[i] - l.a - l.b * x[i]; + rss += w[i] * e * e; + } + // Weights normalised to the number of points, so the error is that of an ordinary fit of n points. + const double n = static_cast(x.size()); + l.se_b = std::sqrt(rss / sw * n / (n - 2) / (sxx / sw * n)); + return l; +} + +// Which modules count: those the placement moved (enough spots), excluding the unmoved ones. +std::vector Moved(const GeometryRefinerResult &r, const GeometryRefinerSettings &st) { + std::vector out; + for (size_t m = 0; m < r.modules.size(); m++) + if (r.modules[m].spots >= st.min_module_spots) + out.push_back(m); + return out; +} + +void LogShifts(Logger &log, const std::string &label, const std::vector ¢res, + const std::vector &modules, const std::vector &which) { + for (size_t m: which) + log.Info("ModuleDiag {} module {} dx {:+.3f} dy {:+.3f} turn_deg {:+.4f}", label, m, + modules[m].dx - centres[m].dx, modules[m].dy - centres[m].dy, + (modules[m].rot_rad - centres[m].rot_rad) * RAD2DEG); +} + +// The shift field's global patterns, fitted over the modules (weights: spots): translation, turn about +// the beam, radial scale, shear, and the two tilt-like quadratic patterns. What is left after them is +// what the modules hold on their own. +void Decompose(Logger &log, const std::string &label, double bx, double by, const std::vector &mods, + const std::vector &which, double lever) { + const size_t n = which.size(); + if (n < 4) + return; + // Unknowns: tx, ty, a11, a12, a21, a22, tilt_y, tilt_x (8); two equations per module. + Eigen::MatrixXd A = Eigen::MatrixXd::Zero(2 * n, 8); + Eigen::VectorXd b(2 * n), w(2 * n); + for (size_t k = 0; k < n; k++) { + const auto &m = mods[which[k]]; + const double x = (m.x0 - bx) / lever, y = (m.y0 - by) / lever; + A(2 * k, 0) = 1; A(2 * k, 2) = x; A(2 * k, 3) = y; A(2 * k, 6) = x * x; A(2 * k, 7) = x * y; + A(2 * k + 1, 1) = 1; A(2 * k + 1, 4) = x; A(2 * k + 1, 5) = y; A(2 * k + 1, 6) = x * y; A(2 * k + 1, 7) = y * y; + b(2 * k) = m.dx; + b(2 * k + 1) = m.dy; + w(2 * k) = w(2 * k + 1) = std::sqrt(static_cast(m.spots)); + } + const Eigen::MatrixXd Aw = w.asDiagonal() * A; + const Eigen::VectorXd bw = w.asDiagonal() * b; + // Linear part alone, then with the tilt-like terms: how much the quadratic patterns explain. + const Eigen::VectorXd lin = Aw.leftCols(6).colPivHouseholderQr().solve(bw); + const Eigen::VectorXd full = Aw.colPivHouseholderQr().solve(bw); + auto rms = [&](const Eigen::VectorXd &res) { + double s = 0, sw = 0; + for (Eigen::Index i = 0; i < res.size(); i++) { s += w(i) * w(i) * res(i) * res(i); sw += w(i) * w(i); } + return std::sqrt(s / sw); + }; + const double rms_raw = rms(b); + const double rms_lin = rms(b - A.leftCols(6) * lin); + const double rms_full = rms(b - A * full); + log.Info("ModuleDiag {} decomposition: translation ({:+.3f}, {:+.3f}) px, turn {:+.4f} deg, radial scale " + "{:+.3f} px at the lever, shear ({:+.3f}, {:+.3f}) px; tilt-like ({:+.3f}, {:+.3f}) px at the lever; " + "rms raw {:.3f}, after linear {:.3f}, after linear+tilt {:.3f} px (lever {:.0f} px)", + label, full(0), full(1), 0.5 * (full(4) - full(3)) / lever * RAD2DEG, + 0.5 * (full(2) + full(5)), 0.5 * (full(2) - full(5)), 0.5 * (full(3) + full(4)), + full(6), full(7), rms_raw, rms_lin, rms_full, lever); +} + +void RadialTable(Logger &log, const std::string &label, double bx, double by, const std::vector &mods, + const std::vector &which, const std::vector &mean_d, + const std::vector &rms_before, const std::vector &rms_after, + const std::vector> &sigma) { + std::vector r, rad, tan, mag, turn, w; + for (size_t m: which) { + const auto &s = mods[m]; + const double ox = s.x0 - bx, oy = s.y0 - by, rr = std::hypot(ox, oy); + const double ux = ox / rr, uy = oy / rr; + const double radial = s.dx * ux + s.dy * uy, tangential = -s.dx * uy + s.dy * ux; + r.push_back(rr); + rad.push_back(radial); + tan.push_back(tangential); + mag.push_back(std::hypot(s.dx, s.dy)); + turn.push_back(s.rot_rad * RAD2DEG); + w.push_back(static_cast(s.spots)); + const auto sg = m < sigma.size() ? sigma[m] : std::array{NAN, NAN, NAN}; + log.Info("ModuleDiag {} table module {} r_px {:.0f} mean_d {:.2f} spots {} dx {:+.3f} dy {:+.3f} radial {:+.3f} " + "tangential {:+.3f} turn_deg {:+.4f} sigma_dx {:.3f} sigma_dy {:.3f} sigma_turn_deg {:.4f} " + "rms_before {:.3f} rms_after {:.3f}", label, m, rr, m < mean_d.size() ? mean_d[m] : NAN, s.spots, + s.dx, s.dy, radial, tangential, s.rot_rad * RAD2DEG, sg[0], sg[1], sg[2] * RAD2DEG, + m < rms_before.size() ? rms_before[m] : NAN, m < rms_after.size() ? rms_after[m] : NAN); + } + for (const auto &[name, v]: std::vector *>>{ + {"radial", &rad}, {"tangential", &tan}, {"magnitude", &mag}, {"turn_deg", &turn}}) { + const Line l = FitLine(r, *v, w); + log.Info("ModuleDiag {} trend {} vs r: slope {:+.3e} +- {:.3e} per px ({:+.3f} +- {:.3f} px per 1000 px), " + "corr {:+.2f}", label, name, l.b, l.se_b, 1000 * l.b, 1000 * l.se_b, l.corr); + } +} + +} // namespace + +void RunModuleDiagnostics(const ModuleDiagnosticsInput &in) { + Logger log("ModuleDiag"); + const auto &st = in.settings; + const auto &placed = in.placed; + const std::vector which = Moved(placed, st); + const double bx = placed.beam_x_px, by = placed.beam_y_px; + double lever = 0; + for (size_t m: which) + lever += std::hypot(in.centres[m].x0 - bx, in.centres[m].y0 - by); + lever /= std::max(1, which.size()); + + // Before: the same frames, modules held where given. + const Fit held = Place(in.nominal, in.cell, in.frames, in.centres, st, false); + + // Bootstrap over frames: the placement's spread, per module. + constexpr int BOOT = 20; + std::mt19937 rng(12345); + std::vector> sum(in.centres.size(), {0, 0, 0}), sum2(in.centres.size(), {0, 0, 0}); + int boots = 0; + for (int b = 0; b < BOOT; b++) { + std::vector resample; + std::uniform_int_distribution pick(0, in.frames.size() - 1); + for (size_t i = 0; i < in.frames.size(); i++) + resample.push_back(in.frames[pick(rng)]); + const Fit f = Place(in.nominal, in.cell, resample, in.start, st, true); + if (!f.ok) + continue; + boots++; + for (size_t m = 0; m < f.modules.size(); m++) { + const double v[3] = {f.modules[m].dx, f.modules[m].dy, f.modules[m].rot_rad}; + for (int j = 0; j < 3; j++) { sum[m][j] += v[j]; sum2[m][j] += v[j] * v[j]; } + } + } + std::vector> sigma(in.centres.size(), {NAN, NAN, NAN}); + if (boots >= 2) + for (size_t m = 0; m < sigma.size(); m++) + for (int j = 0; j < 3; j++) { + const double mean = sum[m][j] / boots; + sigma[m][j] = std::sqrt(std::max(0.0, sum2[m][j] / boots - mean * mean) * boots / (boots - 1)); + } + log.Info("ModuleDiag bootstrap: {} of {} resamples placed", boots, BOOT); + + log.Info("ModuleDiag placed: beam ({:.3f}, {:.3f}) distance {:.4f} rot1 {:.5f} rot2 {:.5f} cell {:.3f} {:.3f} {:.3f}", + placed.beam_x_px, placed.beam_y_px, placed.distance_mm, placed.rot1_rad, placed.rot2_rad, + placed.cell.a, placed.cell.b, placed.cell.c); + RadialTable(log, "placed", bx, by, placed.modules, which, placed.module_mean_d_A, held.r.module_capped_rms_px, + placed.module_capped_rms_px, sigma); + Decompose(log, "placed", bx, by, placed.modules, which, lever); + + // Low against high resolution: the placement fitted on the spots coarser than 3 A, and on the finer ones. + for (int half = 0; half < 2; half++) { + std::vector sub = in.frames; + for (auto &f: sub) { + std::vector keep; + for (const auto &s: f.spots) { + const double d = 1.0 / std::max(1e-6f, in.nominal.DetectorToRecip(s.x, s.y).Length()); + if ((half == 0) == (d > 3.0)) + keep.push_back(s); + } + f.spots = keep; + } + const std::string label = half == 0 ? "lowres_d_gt_3A" : "highres_d_lt_3A"; + const Fit f = Place(in.nominal, in.cell, sub, in.start, st, true); + if (!f.ok) { + log.Info("ModuleDiag {}: no placement (too few spots)", label); + continue; + } + log.Info("ModuleDiag {}: beam ({:.3f}, {:.3f}) distance {:.4f} spots {}", label, f.r.beam_x_px, f.r.beam_y_px, + f.r.distance_mm, f.r.spots_used); + LogShifts(log, label, in.centres, f.modules, Moved(f.r, st)); + } + + // A wrong start: the beam 3 px off in x and y, the distance 1 mm long. + { + DiffractionGeometry g = in.nominal; + g.BeamX_pxl(g.GetBeamX_pxl() + 3).BeamY_pxl(g.GetBeamY_pxl() + 3).DetectorDistance_mm(g.GetDetectorDistance_mm() + 1); + const Fit f = Place(g, in.cell, in.frames, in.start, st, true); + if (f.ok) { + log.Info("ModuleDiag perturbed_start: beam ({:.3f}, {:.3f}) distance {:.4f}", f.r.beam_x_px, f.r.beam_y_px, + f.r.distance_mm); + LogShifts(log, "perturbed_start", in.centres, f.modules, Moved(f.r, st)); + } + } + + // The cell free (a weak anchor) instead of pinned. + { + GeometryRefinerSettings s2 = st; + s2.cell_reg_coeff = 0.3; + const Fit f = Place(in.nominal, in.cell, in.frames, in.start, s2, true); + if (f.ok) { + log.Info("ModuleDiag cell_free: beam ({:.3f}, {:.3f}) distance {:.4f} cell {:.3f} {:.3f} {:.3f}", + f.r.beam_x_px, f.r.beam_y_px, f.r.distance_mm, f.r.cell.a, f.r.cell.b, f.r.cell.c); + LogShifts(log, "cell_free", in.centres, f.modules, Moved(f.r, st)); + } + } + + // The tilt refined, and the modules held to no tilt-like pattern. + { + GeometryRefinerSettings s2 = st; + s2.refine_tilt = true; + const Fit f = Place(in.nominal, in.cell, in.frames, in.start, s2, true); + if (f.ok) { + log.Info("ModuleDiag tilt_refined: beam ({:.3f}, {:.3f}) distance {:.4f} rot1_deg {:+.4f} rot2_deg {:+.4f} " + "median_residual {:.3f}", f.r.beam_x_px, f.r.beam_y_px, f.r.distance_mm, f.r.rot1_rad * RAD2DEG, + f.r.rot2_rad * RAD2DEG, f.r.median_residual_px); + const auto w2 = Moved(f.r, st); + RadialTable(log, "tilt_refined", f.r.beam_x_px, f.r.beam_y_px, f.modules, w2, f.r.module_mean_d_A, + held.r.module_capped_rms_px, f.r.module_capped_rms_px, {}); + Decompose(log, "tilt_refined", f.r.beam_x_px, f.r.beam_y_px, f.modules, w2, lever); + } + } + + // Inner against outer modules: the global geometry from each alone, modules held. + { + std::vector r; + for (size_t m: which) + r.push_back(std::hypot(in.centres[m].x0 - bx, in.centres[m].y0 - by)); + std::vector sorted = r; + std::sort(sorted.begin(), sorted.end()); + const double median_r = sorted.empty() ? 0 : sorted[sorted.size() / 2]; + for (int outer = 0; outer < 2; outer++) { + std::vector sub = in.frames; + for (auto &f: sub) { + std::vector keep; + for (const auto &s: f.spots) { + if (s.module < 0) + continue; + const double rr = std::hypot(in.centres[s.module].x0 - bx, in.centres[s.module].y0 - by); + if ((rr >= median_r) == (outer == 1)) + keep.push_back(s); + } + f.spots = keep; + } + const Fit f = Place(in.nominal, in.cell, sub, in.centres, st, false); + log.Info("ModuleDiag {}_modules_only: ok {} beam ({:.3f}, {:.3f}) distance {:.4f} cell {:.3f} {:.3f} {:.3f} " + "spots {}", outer ? "outer" : "inner", f.ok, f.r.beam_x_px, f.r.beam_y_px, f.r.distance_mm, + f.r.cell.a, f.r.cell.b, f.r.cell.c, f.r.spots_used); + } + } + + // Synthetic control: every spot moved onto its prediction at the placed geometry, as found on unmoved + // modules, with Gaussian noise of its module's rms. Placing these must find no module moves; with a + // parallax shift added (the mean conversion depth times tan 2theta, radially), whatever it finds is what + // parallax would put into the modules. + { + const double depth_um = in.sensor_thickness_um > 0 ? in.sensor_thickness_um : 300.0; + const std::string material = in.sensor_material.empty() ? "Si" : in.sensor_material; + const double L = sensor_absorption::AttenuationLength_um(material, in.nominal.GetWavelength_A()); + // Mean depth of an exponential truncated at the sensor's thickness. + const double mean_z_um = L - depth_um * std::exp(-depth_um / L) / (1 - std::exp(-depth_um / L)); + const double mean_z_px = mean_z_um / (in.nominal.GetPixelSize_mm() * 1000.0); + DiffractionGeometry truth = in.nominal; + truth.BeamX_pxl(static_cast(placed.beam_x_px)).BeamY_pxl(static_cast(placed.beam_y_px)) + .DetectorDistance_mm(static_cast(placed.distance_mm)); + for (int parallax = 0; parallax < 2; parallax++) { + std::mt19937 noise_rng(777); + std::vector syn = in.frames; + size_t kept = 0; + for (size_t i = 0; i < syn.size(); i++) { + std::vector keep; + for (size_t k = 0; k < syn[i].spots.size(); k++) { + auto s = syn[i].spots[k]; + if (i >= placed.predicted_px.size() || k >= placed.predicted_px[i].size()) + continue; + const auto p = placed.predicted_px[i][k]; + if (!std::isfinite(p[0]) || s.module < 0) + continue; + // The prediction is where the placed modules would have the spot found; on modules + // with no error it is found where the placed module puts it. + const auto &mv = placed.modules[s.module]; + const double ux = p[0] - mv.x0, uy = p[1] - mv.y0; + const double c = std::cos(mv.rot_rad), sn = std::sin(mv.rot_rad); + double x = mv.x0 + c * ux - sn * uy + mv.dx; + double y = mv.y0 + sn * ux + c * uy + mv.dy; + if (parallax) { + const double ox = x - truth.GetBeamX_pxl(), oy = y - truth.GetBeamY_pxl(); + const double rr = std::hypot(ox, oy); + const double tan2t = rr * truth.GetPixelSize_mm() / truth.GetDetectorDistance_mm(); + if (rr > 0) { + x += mean_z_px * tan2t * ox / rr; + y += mean_z_px * tan2t * oy / rr; + } + } + const double sigma_px = s.module < static_cast(placed.module_capped_rms_px.size()) + ? std::max(0.1, placed.module_capped_rms_px[s.module]) : 0.5; + std::normal_distribution nd(0.0, sigma_px / std::sqrt(2.0)); + s.x = static_cast(x + nd(noise_rng)); + s.y = static_cast(y + nd(noise_rng)); + keep.push_back(s); + } + kept += keep.size(); + syn[i].spots = keep; + } + const std::string label = parallax ? "synthetic_parallax" : "synthetic_no_error"; + const Fit f = Place(truth, placed.cell, syn, in.centres, st, true); + log.Info("ModuleDiag {}: {} spots, mean conversion depth {:.1f} um ({:.3f} px; {} {:.0f} um), ok {}, beam " + "({:.3f}, {:.3f}) distance {:.4f}", label, kept, mean_z_um, mean_z_px, material, depth_um, f.ok, + f.r.beam_x_px, f.r.beam_y_px, f.r.distance_mm); + if (f.ok) { + RadialTable(log, label, f.r.beam_x_px, f.r.beam_y_px, f.modules, Moved(f.r, st), f.r.module_mean_d_A, + {}, f.r.module_capped_rms_px, {}); + Decompose(log, label, f.r.beam_x_px, f.r.beam_y_px, f.modules, Moved(f.r, st), lever); + } + } + } +} diff --git a/rugnux/ModuleDiagnostics.h b/rugnux/ModuleDiagnostics.h new file mode 100644 index 000000000..ec568d7fa --- /dev/null +++ b/rugnux/ModuleDiagnostics.h @@ -0,0 +1,31 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "../common/DiffractionGeometry.h" +#include "../common/UnitCell.h" +#include "../image_analysis/geom_refinement/GeometryRefiner.h" + +// Tests of whether a module placement is the detector's metrology or something else absorbed module by +// module (--refine-modules=diagnose): the placement against the modules' radius and resolution, its global +// patterns, frame-bootstrap errors, fits on low- and high-resolution spots apart, a perturbed start, inner +// against outer modules, a placement of synthetic spots with no module error (and with parallax added), +// the cell set free and the tilt refined. Everything goes to the log, one line per number, under the +// "ModuleDiag" tag, for a script to tabulate. +struct ModuleDiagnosticsInput { + DiffractionGeometry nominal; + UnitCell cell; + std::vector frames; // the placement's frames, each spot's module assigned + std::vector centres; // the modules, unmoved + std::vector start; // where the placement started (after the coarse search) + GeometryRefinerSettings settings; // as the placement ran + GeometryRefinerResult placed; // the placement + std::string sensor_material; + double sensor_thickness_um = 0; +}; + +void RunModuleDiagnostics(const ModuleDiagnosticsInput &in); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 6f1db8cd2..037f225a1 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -477,6 +477,12 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { // determined rather than the long axis that check has to assume. run->JoinBeamCenterCapture(); result.background_beam_centre = run->measured_beam_center_; + result.modules = run->experiment_.GetModules(); + result.modules_placed = !run->module_shifts_.empty(); + // Measured only (a JUNGFRAU's metrology): the run kept the modules where they were given, and the + // placement it found is what is reported. + if (result.modules_placed && !run->config_.apply_modules) + result.modules.Apply(run->module_shifts_); if (result.consensus_cell) header_experiment.SetUnitCell(*result.consensus_cell); result.beam_centre_need_pxl = BeamCenterNeed_pxl(header_experiment); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 7c3282bc7..b0d485b71 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -18,6 +18,7 @@ #include "../common/DiffractionExperiment.h" #include "../common/PixelMask.h" +#include "../common/DetectorModules.h" #include "../common/AzimuthalIntegrationMapping.h" #include "../common/CrystalLattice.h" #include "../common/JFJochException.h" @@ -150,6 +151,19 @@ struct ProcessConfig { // determines the shared beam/distance/cell that per-image refinement cannot reach, and the main pass // re-indexes with it. The value is the number of strong frames fed to the bundle (--refine-geometry). std::optional refine_geometry; + // Within that bundle adjustment, a detector of several modules (the experiment's DetectorModules: the + // file's, the reader's panels, or --module-geometry) has each module placed as a rigid block, where + // moving them predicts the spots better on frames the move was not fitted to (--refine-modules=off: never). + bool refine_modules = true; + // With it, a run of tests on the placement (ModuleDiagnostics.h) goes to the log (--refine-modules=diagnose). + bool module_diagnostics = false; + // Whether the placement found is used from there on, or only measured and reported (and written to + // _modules.h5). Used for a foreign multi-panel detector (SACLA, a geometry file's panels); a + // JUNGFRAU's metrology is measured only. + bool apply_modules = true; + // The file's detector distance is nominal (SACLA's /metadata): with a reference cell it is refined + // from the cell, as where the frames have shown a stated distance to be wrong. + bool stills_distance_nominal = false; // Beam-stop shadow detection (--detect-beam-stop). When set, a pre-pass projects this many // frames, finds where the beam stop and its holder shadow the detector and marks them in the @@ -407,6 +421,11 @@ struct ProcessResult { // beam-centre pre-scan that runs before the first one. processing_time_s is the last pass alone // and so understates a two-pass run by more than a factor of two. double total_time_s = 0.0; + // Where the detector's modules ended up (the placement the run started from, moved where it placed + // them); empty for a detector of one module. + DetectorModules modules; + // The run moved them (RefineStillsGeometry). + bool modules_placed = false; // The end message of this pass, as the _process.h5 records it: the lattice, refined geometry, // reindexing and the run means - what the broker puts into its own when it runs rugnux live. EndMessage end_message; @@ -925,6 +944,11 @@ class Rugnux { // fed to the second pass's Bragg prediction. Empty outside the rotation two-pass. std::vector prepass_mosaicity_; + // The steps RefineStillsGeometry moved the modules by (empty: not moved). + std::vector module_shifts_; + // Placed once, by the first RefineStillsGeometry. + bool modules_placed_ = false; + // Two-pass geometry pre-pass: the detector geometry the first pass ends on - beam_x, beam_y in px, // distance in mm, then rot1 and rot2 in rad - applied by Run() before the second pass. Empty when // the pre-pass committed no change (the second pass then reproduces the first). diff --git a/rugnux/RugnuxStills.cpp b/rugnux/RugnuxStills.cpp index a7a5c3ced..1bc05a7b9 100644 --- a/rugnux/RugnuxStills.cpp +++ b/rugnux/RugnuxStills.cpp @@ -11,6 +11,7 @@ #include "../image_analysis/hot_pixels/HotPixels.h" #include "../image_analysis/SensorAbsorption.h" #include "DiagnosticOutput.h" +#include "ModuleDiagnostics.h" #include #include @@ -21,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -85,6 +87,18 @@ using namespace rugnux_internal; +namespace { + +// The module of each spot, by where it lies on the detector plane (spots are on the plane already). +void AssignModules(std::vector &frames, const DetectorModules &modules) { + for (auto &f: frames) + for (auto *list: {&f.spots, &f.unindexed}) + for (auto &s: *list) + s.module = static_cast(modules.ModuleOfPlanePoint(s.x, s.y)); +} + +} // namespace + void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_process, RugnuxObserver *observer) { Logger logger("Rugnux"); @@ -131,7 +145,11 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ // least 8000), so on data that indexes well it meant re-indexing the ENTIRE run to keep 200 frames - // 99% of this pass was sampling. Four times the bundle leaves the "strongest N" selection a real // pool to choose from. - const int strong_target = refine_frames * 4; + // Placing the modules of a modular detector needs spots on the outer ones, which only the strong + // high-resolution frames reach, and the placement takes every strong frame sampled (see `placing`): + // there the sample goes five times further. + const bool placing_modules = config_.refine_modules && !modules_placed_ && experiment_.GetModules().size() >= 2; + const int strong_target = refine_frames * (placing_modules ? 20 : 4); // The sample is cut into a FIXED number of interleaved stripes - stripe s is sample positions // s, s+STRIPES, s+2*STRIPES, ... - and each stripe stops once it has contributed its share of the @@ -196,10 +214,13 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ GeomRefineFrame f; f.lattice = *msg.indexing_lattice; f.ordinal = ordinal; - for (const auto &s : msg.spots) + for (const auto &s : msg.spots) { if (s.indexed && s.lattice == 0) f.spots.push_back(GeomRefineSpot{s.x, s.y, static_cast(s.h), static_cast(s.k), static_cast(s.l)}); + else + f.unindexed.push_back(GeomRefineSpot{s.x, s.y}); + } if (f.spots.size() >= 6) { if (static_cast(f.spots.size()) >= MIN_STRONG_SPOTS) stripe_strong++; @@ -234,6 +255,9 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ if (a.spots.size() != b.spots.size()) return a.spots.size() > b.spots.size(); return a.ordinal < b.ordinal; }); + // Placing modules needs spots on every module, and the outer ones see few: that fit takes every strong + // frame the sampling found, not only the bundle's share. + std::vector placing = strong; if (static_cast(strong.size()) > refine_frames) strong.resize(refine_frames); @@ -252,18 +276,121 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ // the user's cell only the starting point. Only where the frames have already shown the file's // distance to be wrong (DecideStillsGeometry, which is what -C buys on a header with a bad distance) // is the cell the better of the two, and then it is anchored and the distance refined from it. - gr_settings.refine_distance = stills_distance_from_cell_; - if (!stills_distance_from_cell_) + const bool distance_from_cell = stills_distance_from_cell_ + || (config_.stills_distance_nominal && !config_.detector_distance_given); + gr_settings.refine_distance = distance_from_cell; + if (!distance_from_cell) gr_settings.cell_reg_coeff = 0.0; - const GeometryRefinerResult r = RefineGlobalGeometry( - experiment_.GetDiffractionGeometry(), *cell, strong, gr_settings); + const DiffractionGeometry nominal = experiment_.GetDiffractionGeometry(); + GeometryRefinerResult r = RefineGlobalGeometry(nominal, *cell, strong, gr_settings); if (!r.ok) { logger.Warning("Geometry refinement did not run (too few strong frames / did not converge); " "keeping the input geometry"); return; } + // A detector of several modules: each placed as a rigid block, where that predicts the spots better on + // frames the placement was not fitted to. The strong frames are dealt alternately into two halves; each + // half is fitted with the modules held and with them moved, and each fit is scored on the other half + // with its geometry frozen. The modules move if the moved fits predict the held-out spots better. + // The modules are those the run was given (the file's, the reader's panels or --module-geometry), never + // guessed from the image; they are placed once, by the first RefineStillsGeometry. + const DetectorModules &given = experiment_.GetModules(); + const std::vector centres = config_.refine_modules && !modules_placed_ && given.size() >= 2 + ? given.Centres() : std::vector{}; + modules_placed_ = true; + if (!centres.empty()) { + AssignModules(placing, given); + AssignModules(strong, given); + + // Where a module is off by more than the spots' spacing allows the refinement to pull, its spots + // went unindexed: each module is first searched for over +-SCAN_RANGE_PX on the bundle's frames, + // and the placement starts from the shift that puts the most of its spots on the lattices. + constexpr double SCAN_RANGE_PX = 25.0; + std::vector start = centres; + const auto scan = ScanModuleShifts(nominal, strong, centres, SCAN_RANGE_PX); + for (size_t m = 0; m < start.size(); m++) + if (scan[m].on_lattice_at_best > scan[m].on_lattice_at_zero) { + start[m].dx = scan[m].dx; + start[m].dy = scan[m].dy; + if (std::hypot(scan[m].dx, scan[m].dy) >= 1.0) + logger.Info("Module {} at ({:.0f}, {:.0f}): {} spots on the lattices as placed, {} shifted by " + "({:+.1f}, {:+.1f}) px", m, start[m].x0, start[m].y0, scan[m].on_lattice_at_zero, + scan[m].on_lattice_at_best, scan[m].dx, scan[m].dy); + } + // Once the modules are placed, spots that went unindexed may sit on the lattice: the placement + // takes every spot and lets the fit sort them. + for (auto &f: placing) { + f.spots.insert(f.spots.end(), f.unindexed.begin(), f.unindexed.end()); + f.unindexed.clear(); + } + + std::vector half[2]; + for (size_t i = 0; i < placing.size(); i++) + half[i % 2].push_back(placing[i]); + + // Held-out residual of the geometry fitted on half `fit`, scored on every spot of the other half. + auto held_out = [&](const std::vector &start, bool refine) { + double sum = 0; + for (int fit = 0; fit < 2; fit++) { + GeometryRefinerSettings st = gr_settings; + st.modules = start; + st.refine_modules = refine; + const auto f = RefineGlobalGeometry(nominal, *cell, half[fit], st); + if (!f.ok) + return std::numeric_limits::infinity(); + DiffractionGeometry g = nominal; + g.BeamX_pxl(f.beam_x_px).BeamY_pxl(f.beam_y_px).DetectorDistance_mm(f.distance_mm); + GeometryRefinerSettings ev = gr_settings; + ev.refine_beam_distance = false; + ev.modules = f.modules; + const auto e = RefineGlobalGeometry(g, f.cell, half[1 - fit], ev); + if (!e.ok) + return std::numeric_limits::infinity(); + sum += e.capped_rms_px; + } + return sum / 2.0; + }; + const double held = held_out(centres, false); + const double moved = held_out(start, true); + logger.Info("Module placement: {} modules; held-out spot residual (capped rms) {:.3f} px as given, " + "{:.3f} px with each placed", centres.size(), held, moved); + + if (moved < held) { + GeometryRefinerSettings st = gr_settings; + st.modules = start; + st.refine_modules = true; + const auto rm = RefineGlobalGeometry(nominal, *cell, placing, st); + if (rm.ok) { + for (size_t m = 0; m < rm.modules.size(); m++) { + const auto &ms = rm.modules[m]; + logger.Info("Module {} at ({:.0f}, {:.0f}): shift ({:+.2f}, {:+.2f}) px, turn {:+.3f} deg, {} spots", + m, ms.x0, ms.y0, ms.dx, ms.dy, ms.rot_rad * 180.0 / PI, ms.spots); + } + // The move is a step from where the modules were given; the experiment now has them there, + // so every spot found and every prediction integrated from here on uses it. The image is + // not touched. + if (config_.module_diagnostics) { + ModuleDiagnosticsInput diag{nominal, *cell, placing, centres, start, st, rm, + experiment_.GetDetectorSetup().GetSensorMaterial(), + experiment_.GetDetectorSetup().GetSensorThickness_um()}; + RunModuleDiagnostics(diag); + } + module_shifts_ = rm.modules; + if (config_.apply_modules) { + DetectorModules placed = given; + placed.Apply(rm.modules); + experiment_.Modules(placed); + r = rm; + logger.Info("Module placement: modules placed (median residual {:.3f} px)", rm.median_residual_px); + } else + logger.Info("Module placement: measured and reported only - the run keeps the modules where " + "the file has them (median residual {:.3f} px with them placed)", rm.median_residual_px); + } + } + } + logger.Info("Geometry refinement: beam ({:.2f}, {:.2f}) -> ({:.2f}, {:.2f}) px, distance {:.4f} -> {:.4f} mm", experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), r.beam_x_px, r.beam_y_px, experiment_.GetDetectorDistance_mm(), r.distance_mm); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 0e88d9a5d..70dda57f5 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -29,11 +29,15 @@ #include "../reader/JFJochCBFReader.h" #include "../reader/JFJochMarCCDReader.h" #include "../reader/JFJochSMVReader.h" -#include "../reader/JFJochCrystFELReader.h" +#include "../reader/JFJochXFELReader.h" #include "../reader/JFJochXTCReader.h" +#include "../reader/CrystFELGeometry.h" +#include "../reader/NXmxModules.h" #include "../reader/JFJochStartMessageReader.h" +#include "../image_analysis/bragg_integration/BackgroundGain.h" #include "BrokerFeed.h" #include "../writer/HDF5NXmx.h" +#include "../writer/NXmxModules.h" #include #include #include @@ -157,11 +161,13 @@ void print_usage() { std::cout << std::endl; std::cout << " Detector geometry" << std::endl; - std::cout << " --crystfel-geometry Read the HDF5/CXI files as CrystFEL does with this geometry file: frames from its data path, panels placed at its corners, camera length, pixel size, gain, photon energy and bad regions from it. One wavelength and one camera length for the run; panel positions are taken as given" << std::endl; - std::cout << " (Cheetah CXI files with a square frame are read without one: one flat panel, with the pixel size, photon energy and camera length the file states and the beam centre measured; --detector-distance overrides a camera length the file has wrong)" << std::endl; + std::cout << " (Serial XFEL stills are read without a geometry file: SACLA Cheetah files are placed by their own /metadata; Cheetah CXI files of a single-panel detector (square frame) are one flat panel with the pixel size, photon energy and camera length they state; LCLS-I XTC streams of a Rayonix take the photon energy and pixel size from the XTC and need --detector-distance. The beam centre is measured, the run is read at the mean photon energy, and SACLA or CXI files that state no gain are converted to photons with the gain their background shows)" << std::endl; + std::cout << " --crystfel-geometry Read the HDF5/CXI files with their panels placed by this CrystFEL geometry file: only the module geometry (panel ranges, fast/slow directions, corners, pixel size, coffset, the frames' data path) and, if stated, the camera length are taken from it; photon energy, masks and camera length otherwise come from the data files (Cheetah CXI, SACLA) or the command line, gain and saturation are measured. Tilted panels, dimN layouts, per-panel data and detector shifts are refused" << std::endl; std::cout << std::endl; std::cout << " Detector mask" << std::endl; + std::cout << " --mask-good : HDF5/CXI stills: a mask array in another HDF5 file, in the frames' own layout (as stored, before the panels are placed); a pixel is usable only where it is non-zero. May be given more than once" << std::endl; + std::cout << " --mask-bad : As --mask-good, but a pixel is unusable where the array is non-zero" << std::endl; std::cout << " --detect-beam-stop[=N|off] Find the beam stop and its holder in a projection of N images and add them to the pixel mask (bit 9), so nothing shadowed by them is integrated. ON by default (60 images); =off disables. Reflections behind the stop are attenuated but not flagged, so they are integrated low with a plausible sigma and no existing rejection catches them" << std::endl; std::cout << " --beam-center-check[=off] Measure the beam centre from the isotropy of the scattered background on EVERY run, report how far the file's value is from it, and index a second first pass with it to see whether the two centres give the same lattice. The fit reads the projection --detect-beam-stop already builds, so it costs no extra frames. On a run that indexes, the measured centre is adopted where the file's centre indexes nothing and the measured one indexes a majority, and where the two centres return cells related by an integer volume factor, in either direction, and the measured centre's lattice carries materially more of the pooled spots against its own chance level - the axis harmonic a centre error along the spindle produces, which nothing downstream repairs. Where the indexing cannot tell the two centres apart - the same cell with a different metric symmetry, or the same lattice found at both only by reading less of each frame - the first pass is run at both and the run keeps the one that merges better (on the same lattice: the one that measures more reflections at I/sigma >= 2). ON by default; =off disables On stills, where there is no first pass, the file's centre, the measured one and the point-symmetry capture's are each indexed on a sample of frames and the run takes the one that puts materially more spots on the lattice of their frame; with a reference cell (-C) the detector distance is searched the same way, over 60-160% of the file's, before the stills geometry refinement" << std::endl; std::cout << std::endl; @@ -197,6 +203,8 @@ void print_usage() { std::cout << " --fft-min-unit-cell Shortest cell axis the FFT search accepts, in A (default: 10). A candidate with a shorter axis is discarded, but a de-novo rotation run also tries a second first-pass hypothesis with the floor at 5 A, so a small-molecule cell is indexed without this flag; -C lowers it on its own to cover the cell given" << std::endl; std::cout << " -r, --refine Geometry refinement algorithm (none|orientation|beam_and_lattice|flex); flex tries all three per image and keeps whichever indexes the most spots (alias: multi)" << std::endl; std::cout << " --refine-geometry[=N|off] Stills: extra first pass that bundle-adjusts the shared beam/distance/cell from N strongly-indexed frames (default: 200) then re-indexes (lifts weak-stills indexing). Default ON for stills when a reference cell is given (-C / reference MTZ); =off disables" << std::endl; + std::cout << " --refine-modules[=on|off|diagnose] Stills, in that bundle adjustment: a detector of several modules (the input's own - NXdetector_module groups, SACLA sensors, geometry-file panels - or --module-geometry) has each placed as a rigid block - shift and in-plane turn - where that predicts the spots better on frames held out of the fit (default: on). Used from there on for the panels of a foreign detector (SACLA sensors, geometry-file panels); for a JUNGFRAU measured and reported only. The placement is written to _modules.h5" << std::endl; + std::cout << " --module-geometry Where the detector's modules are: an NXmx file with one NXdetector_module group per module (module_0, module_1, ...; what a run writes as _modules.h5), or a CrystFEL geometry file whose panels are laid on the image (fs = column, ss = row). Replaces the input's own modules" << std::endl; std::cout << " --index-ice-rings[=on|off] Index on the spots flagged as sitting on an ice ring too, instead of setting them aside (default: off; no effect without --detect-ice-rings)" << std::endl; std::cout << std::endl; @@ -267,7 +275,7 @@ void print_usage() { std::cout << " --detector-distance Detector distance (mm)" << std::endl; std::cout << " --wavelength Wavelength (A)" << std::endl; std::cout << " --pixel-size Pixel size (um)" << std::endl; - std::cout << " --gain Detector units per photon, for a file that states none (CXI and --crystfel-geometry input; the frames are converted to photons with it)" << std::endl; + std::cout << " --gain Detector units per photon, for a file that states none (SACLA, CXI, XTC and --crystfel-geometry input; the frames are converted to photons with it). Default: for SACLA and CXI input, measured on the background of a few frames; XTC frames stay in detector units" << std::endl; std::cout << " --saturation Pixel value from which a pixel is overloaded, in the units the frames are read in (photons where they are converted). Replaces the file's limit and the one stills read off the frames" << std::endl; std::cout << " --rot1 PONI rotation 1 (rad)" << std::endl; std::cout << " --rot2 PONI rotation 2 (rad)" << std::endl; @@ -366,6 +374,10 @@ enum { OPT_PIXEL_SIZE, OPT_GAIN, OPT_SATURATION, + OPT_MASK_GOOD, + OPT_MASK_BAD, + OPT_REFINE_MODULES, + OPT_MODULE_GEOMETRY, OPT_ROT1, OPT_ROT2, OPT_ROT3, @@ -434,6 +446,8 @@ static option long_options[] = { {"pixel-size", required_argument, nullptr, OPT_PIXEL_SIZE}, {"gain", required_argument, nullptr, OPT_GAIN}, {"saturation", required_argument, nullptr, OPT_SATURATION}, + {"mask-good", required_argument, nullptr, OPT_MASK_GOOD}, + {"mask-bad", required_argument, nullptr, OPT_MASK_BAD}, {"rot1", required_argument, nullptr, OPT_ROT1}, {"rot2", required_argument, nullptr, OPT_ROT2}, {"rot3", required_argument, nullptr, OPT_ROT3}, @@ -446,6 +460,8 @@ static option long_options[] = { {"prepass-fraction", required_argument, nullptr, OPT_PREPASS_FRACTION}, {"rotation-scale", required_argument, nullptr, OPT_ROTATION_SCALE}, {"refine-geometry", optional_argument, nullptr, OPT_REFINE_GEOMETRY}, + {"refine-modules", optional_argument, nullptr, OPT_REFINE_MODULES}, + {"module-geometry", required_argument, nullptr, OPT_MODULE_GEOMETRY}, {"detect-beam-stop", optional_argument, nullptr, OPT_DETECT_BEAM_STOP}, {"beam-center-check", optional_argument, nullptr, OPT_BEAM_CENTER_CHECK}, @@ -697,6 +713,31 @@ void configure_offline_output(DiffractionExperiment &experiment, const std::stri } namespace { + +// The gain of a detector whose file states none - ADU per photon, from the background's variance per +// unit mean (BackgroundGain.h) on a few frames spread over the run - so a foreign stills file can be +// converted to photons as it is read. +std::optional MeasureGain(JFJochReader &reader) { + const auto dataset = reader.GetDataset(); + const size_t n = reader.GetNumberOfImages(); + const size_t frames = std::min(n, background_gain::GAIN_MAX_FRAMES); + std::vector levels; + JFJochReaderRawImage raw; + std::vector buffer; + std::vector pixels; + for (size_t i = 0; i < frames; i++) { + if (!reader.ReadRawImage(static_cast(i * n / frames), raw) + || raw.image.GetMode() != CompressedImageMode::Int32) + continue; + const size_t w = raw.image.GetWidth(), h = raw.image.GetHeight(); + pixels.resize(w * h); + std::memcpy(pixels.data(), raw.image.GetUncompressedPtr(buffer), w * h * sizeof(int32_t)); + levels.push_back(background_gain::MeasureFrame(pixels, w, h, dataset->pixel_mask->GetMask(), + dataset->experiment.GetSaturationLimit())); + } + return background_gain::Gain(levels); +} + std::atomic g_active_process{nullptr}; // --simulate-live: a stored sweep read as the collection it was while being written - image i is @@ -996,6 +1037,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { // Geometry overrides (default: keep the value stored in the input file) std::optional beam_x, beam_y, detector_distance_mm, wavelength_A, rot1_rad, rot2_rad, rot3_rad, polarization_factor; std::optional pixel_size_um, gain_per_photon, saturation_value; + std::vector extra_masks; // --mask-good / --mask-bad bool detector_mirror_y = false; int64_t detector_quarter_turns = 0; std::optional smooth_g_deg_arg; // --smooth-g[=deg]; default 5 deg for rot3d, 0 (off) otherwise @@ -1034,6 +1076,9 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { bool beam_center_check = true; // --beam-center-check[=off]; on by default std::optional refine_geometry; // --refine-geometry[=N]: stills global geometry-refinement pass bool refine_geometry_disabled = false; // --refine-geometry=off: opt out of the stills default-on + bool refine_modules = true; // --refine-modules=off + bool module_diagnostics = false; // --refine-modules=diagnose + std::string module_geometry_file; // --module-geometry std::optional bandwidth_fwhm; // relative FWHM of dlambda/lambda @@ -1153,6 +1198,13 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { : RUGNUX_REFINE_GEOMETRY_FRAMES; break; } + case OPT_REFINE_MODULES: + refine_modules = !optarg || std::string(optarg) != "off"; + module_diagnostics = optarg && std::string(optarg) == "diagnose"; + break; + case OPT_MODULE_GEOMETRY: + module_geometry_file = optarg; + break; case OPT_FORCE_ROTATION_LATTICE: { if (rotation_indexing) { logger.Error("Rotation indexing already enabled"); @@ -1593,6 +1645,26 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { case OPT_PIXEL_SIZE: pixel_size_um = parse_float_arg(optarg, "--pixel-size", logger); break; case OPT_GAIN: gain_per_photon = parse_float_arg(optarg, "--gain", logger); break; case OPT_SATURATION: saturation_value = parse_float_arg(optarg, "--saturation", logger); break; + case OPT_MASK_GOOD: + case OPT_MASK_BAD: { + // :, split at the last colon (the dataset is an HDF5 path, starting with /). + const std::string arg = optarg; + const auto colon = arg.rfind(':'); + if (colon == std::string::npos || colon + 1 >= arg.size() || arg[colon + 1] != '/') { + logger.Error("{} takes :, e.g. mask.h5:/mask", + opt == OPT_MASK_GOOD ? "--mask-good" : "--mask-bad"); + exit(EXIT_FAILURE); + } + xfel::Mask m; + m.file = arg.substr(0, colon); + m.data = arg.substr(colon + 1); + if (opt == OPT_MASK_GOOD) + m.nonzero_good = true; + else + m.bad = 0xffffffffu; + extra_masks.push_back(m); + break; + } case OPT_WAVELENGTH: { // Guard > 0: wavelength is used as a divisor (WVL_1A_IN_KEV / wavelength) below, and a 0 // would produce a non-finite incident energy that throws unguarded and aborts the process. @@ -1784,18 +1856,17 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { SimulatedLive cbf_reader; SimulatedLive marccd_reader; SimulatedLive smv_reader; - SimulatedLive crystfel_reader; + SimulatedLive xfel_reader; SimulatedLive xtc_reader; JFJochStartMessageReader live_reader; std::unique_ptr multi_reader; std::unique_ptr broker_feed; JFJochReader *reader_ptr = nullptr; - // Which format this is. Asked INSIDE the try: these open and read the file - a directory whose - // 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. + // Which format this is. Asked INSIDE the try: these open and read the file, 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; - bool input_is_cxi = false; + bool input_is_cxi = false, input_is_xtc = false, input_is_panel_table = false; try { // Every input is the format the first one is. const auto all_inputs = [&input_files](auto can_read) { @@ -1806,8 +1877,11 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { "kind - Jungfraujoch/NXmx masters, SACLA HDF5, CXI or XTC"); exit(EXIT_FAILURE); }; - if (input_files.size() > 1 && crystfel_geometry.empty() && !JFJochCrystFELReader::CanReadCXI(input_file) - && !JFJochXTCReader::CanRead(input_file)) { + const auto is_foreign = [](const std::string &f) { + return JFJochXFELReader::CanReadPanelTable(f) || JFJochXFELReader::CanReadCXI(f) || JFJochXTCReader::CanRead(f); + }; + input_is_crystfel = !crystfel_geometry.empty(); + if (input_files.size() > 1 && !input_is_crystfel && !is_foreign(input_file)) { // One frame of a CBF, marCCD or SMV sweep names the whole sweep. if (JFJochCBFReader::CanRead(input_file) || JFJochMarCCDReader::CanRead(input_file) || JFJochSMVReader::CanRead(input_file)) { @@ -1815,9 +1889,8 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { "frames", input_file); exit(EXIT_FAILURE); } - if (!all_inputs([](const std::string &f) { - return !JFJochCrystFELReader::CanReadCXI(f) && !JFJochXTCReader::CanRead(f) - && !JFJochCBFReader::CanRead(f) && !JFJochMarCCDReader::CanRead(f) + if (!all_inputs([&is_foreign](const std::string &f) { + return !is_foreign(f) && !JFJochCBFReader::CanRead(f) && !JFJochMarCCDReader::CanRead(f) && !JFJochSMVReader::CanRead(f); })) refuse_mixed(); std::vector> inputs; @@ -1838,30 +1911,52 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { multi_reader->FirstImage(k) + static_cast(multi_reader->ImagesInInput(k)) - 1); reader_ptr = multi_reader.get(); } else { - input_is_crystfel = !crystfel_geometry.empty(); - if (gain_per_photon && !input_is_crystfel && !JFJochCrystFELReader::CanReadCXI(input_file)) { - logger.Error("--gain is for a CXI file or --crystfel-geometry input; {} is read in the units it states", - input_file); + input_is_panel_table = !input_is_crystfel && JFJochXFELReader::CanReadPanelTable(input_file); + input_is_cxi = !input_is_crystfel && !input_is_panel_table && JFJochXFELReader::CanReadCXI(input_file); + if ((input_is_panel_table && !all_inputs(JFJochXFELReader::CanReadPanelTable)) + || (input_is_cxi && !all_inputs(JFJochXFELReader::CanReadCXI))) + refuse_mixed(); + if (!extra_masks.empty() && !input_is_crystfel && !input_is_panel_table && !input_is_cxi) { + logger.Error("--mask-good and --mask-bad are for HDF5/CXI stills (SACLA, CXI or --crystfel-geometry " + "input); {} is not one", input_file); exit(EXIT_FAILURE); } - input_is_cxi = !input_is_crystfel && JFJochCrystFELReader::CanReadCXI(input_file); - input_is_cbf = !input_is_crystfel && !input_is_cxi && 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_cxi && !input_is_cbf && !input_is_marccd && !input_is_smv - && JFJochXTCReader::CanRead(input_file); - if ((input_is_cxi && !all_inputs(JFJochCrystFELReader::CanReadCXI)) - || (input_is_xtc && !all_inputs(JFJochXTCReader::CanRead))) - refuse_mixed(); + xfel_reader.ExtraMasks(extra_masks); if (input_is_crystfel) { - crystfel_reader.ReadFiles(input_files, crystfel_geometry, gain_per_photon); - reader_ptr = &crystfel_reader; + xfel_reader.ReadFiles(input_files, crystfel_geometry, gain_per_photon); + reader_ptr = &xfel_reader; + if (!xfel_reader.HasPhotonEnergy() && !wavelength_A) { + logger.Error("Nothing read states the photon energy of {}; give it with --wavelength", input_file); + exit(EXIT_FAILURE); + } logger.Info("Read {} frames from {} files with the CrystFEL geometry {}", - crystfel_reader.GetNumberOfImages(), input_files.size(), crystfel_geometry); + xfel_reader.GetNumberOfImages(), input_files.size(), crystfel_geometry); + } else if (input_is_panel_table) { + xfel_reader.ReadPanelTable(input_files, gain_per_photon); + reader_ptr = &xfel_reader; + logger.Info("Read {} frames from {} HDF5 files placed by their /metadata panel table", xfel_reader.GetNumberOfImages(), input_files.size()); } else if (input_is_cxi) { - crystfel_reader.ReadCXI(input_files, gain_per_photon); - reader_ptr = &crystfel_reader; - logger.Info("Read {} frames from {} CXI files", crystfel_reader.GetNumberOfImages(), input_files.size()); + xfel_reader.ReadCXI(input_files, gain_per_photon); + reader_ptr = &xfel_reader; + logger.Info("Read {} frames from {} CXI files", xfel_reader.GetNumberOfImages(), input_files.size()); + } else if (const auto refusal = JFJochXFELReader::CXIRefusal(input_file)) { + logger.Error("{}", *refusal); + exit(EXIT_FAILURE); + } + input_is_cbf = !reader_ptr && JFJochCBFReader::CanRead(input_file); + input_is_marccd = !reader_ptr && !input_is_cbf && JFJochMarCCDReader::CanRead(input_file); + input_is_smv = !reader_ptr && !input_is_cbf && !input_is_marccd && JFJochSMVReader::CanRead(input_file); + input_is_xtc = !reader_ptr && !input_is_cbf && !input_is_marccd && !input_is_smv + && JFJochXTCReader::CanRead(input_file); + if (input_is_xtc && !all_inputs(JFJochXTCReader::CanRead)) + refuse_mixed(); + if (gain_per_photon && !reader_ptr && !input_is_xtc) { + logger.Error("--gain is for SACLA, CXI, XTC or --crystfel-geometry input; {} is read in the units it " + "states", input_file); + exit(EXIT_FAILURE); + } + if (reader_ptr) { + // read above } else if (input_is_cbf) { cbf_reader.ReadFiles(input_file); reader_ptr = &cbf_reader; @@ -1879,11 +1974,24 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { smv_reader.GetNumberOfImages(), input_file); } else if (input_is_xtc) { xtc_reader.ReadFiles(input_files); - if (!xtc_reader.HasGeometry() && (!pixel_size_um || !detector_distance_mm || !wavelength_A)) { - logger.Error("XTC states no geometry and no DIALS .expt lies beside it: give --pixel-size, " - "--detector-distance and --wavelength (the beam centre is measured)"); + std::vector missing; + if (!detector_distance_mm) + missing.emplace_back("--detector-distance"); + if (!xtc_reader.HasPixelSize() && !pixel_size_um) + missing.emplace_back("--pixel-size"); + if (!xtc_reader.HasWavelength() && !wavelength_A) + missing.emplace_back("--wavelength"); + if (!missing.empty()) { + std::string options; + for (const auto &m: missing) + options += (options.empty() ? "" : ", ") + m; + logger.Error("The XTC does not state the detector distance{}{}: give {} (the beam centre is " + "measured)", xtc_reader.HasPixelSize() ? "" : ", the pixel size", + xtc_reader.HasWavelength() ? "" : ", the wavelength", options); exit(EXIT_FAILURE); } + if (gain_per_photon) + xtc_reader.ApplyGain(*gain_per_photon); reader_ptr = &xtc_reader; logger.Info("Read {} XTC still images from {} streams", xtc_reader.GetNumberOfImages(), input_files.size()); @@ -1932,6 +2040,20 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { } JFJochReader &reader = *reader_ptr; + // A SACLA or Cheetah CXI file that states no gain is converted to photons with the one its background + // shows, unless --gain gave it. Not an XTC Rayonix: there the conversion cost R-free (one LCLS set, 0.28 -> + // 0.30) where it gained it on the PAL Rayonix, so XTC frames stay in detector units, with the + // background clip widened by the gain the pre-scan measures, as for any CCD. + if (!gain_per_photon && (input_is_panel_table || input_is_cxi || input_is_crystfel) && !xfel_reader.CountsInPhotons()) { + if (const auto g = MeasureGain(reader)) { + logger.Info("Gain: {:.3f} detector units per photon, measured on the background; the frames are " + "converted to photons with it", *g); + xfel_reader.ApplyGain(*g); + } else + logger.Info("Gain: the file states none and the background does not show it; the frames stay in " + "detector units"); + } + const auto dataset = reader.GetDataset(); if (!dataset) { logger.Error("No experiment dataset found in the input file"); @@ -2055,7 +2177,8 @@ 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 || input_is_cxi || input_is_xtc) { + if (input_is_cbf || input_is_marccd || input_is_smv || input_is_crystfel || input_is_panel_table || input_is_cxi + || 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); @@ -2646,6 +2769,42 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { if (detector_mirror_y || detector_quarter_turns != 0) experiment.Detector().ImageOrientation( DetectorOrientation(detector_mirror_y, detector_quarter_turns)); + // Where the modules are, if given: after the geometry above, which a geometry file's panels are laid + // relative to. + if (!module_geometry_file.empty()) { + try { + DetectorModules modules; + if (module_geometry_file.ends_with(".geom")) { + std::vector ignored; + modules = crystfel::ModulesOnImage(crystfel::Read(module_geometry_file, ignored), + experiment.GetBeamX_pxl(), experiment.GetBeamY_pxl()); + } else { + std::unique_lock ul(hdf5_mutex); + HDF5ReadOnlyFile file(module_geometry_file); + modules = ReadNXmxModules(file, experiment.GetPixelSize_mm() / 1000.0, + experiment.GetDiffractionGeometry().GetOrientation()); + } + if (modules.empty()) { + logger.Error("--module-geometry {} describes fewer than two modules", module_geometry_file); + exit(EXIT_FAILURE); + } + for (const auto &m: modules.Get()) + if (m.x0 < 0 || m.y0 < 0 || m.x0 + m.width > experiment.GetXPixelsNum() + || m.y0 + m.height > experiment.GetYPixelsNum()) { + logger.Error("--module-geometry {}: a module lies outside the {} x {} image", module_geometry_file, + experiment.GetXPixelsNum(), experiment.GetYPixelsNum()); + exit(EXIT_FAILURE); + } + experiment.Modules(modules); + logger.Info("--module-geometry: {} modules from {}", modules.size(), module_geometry_file); + } catch (const std::exception &e) { + logger.Error("Cannot read --module-geometry {}: {}", module_geometry_file, e.what()); + exit(EXIT_FAILURE); + } + } + if (!experiment.GetModules().empty()) + logger.Info("Detector: {} modules{}", experiment.GetModules().size(), + experiment.GetModules().Moved() ? ", placed off the image's own grid" : ""); // --polarization is applied after configure_offline_output below, which sets the rugnux default. // Calibrating from the run-summed profile needs the profile to be binned in azimuth - see @@ -3289,11 +3448,11 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { config.beam_center_check = beam_center_check; config.detector_distance_given = detector_distance_mm.has_value(); config.detect_saturation = !saturation_value - && (input_is_cbf || input_is_marccd || input_is_smv - || input_is_crystfel || input_is_cxi || input_is_xtc); + && (input_is_cbf || input_is_marccd || input_is_smv || input_is_crystfel + || input_is_panel_table || input_is_cxi || input_is_xtc); config.adaptive_integration_radius = adaptive_integration_radius; - config.counts_in_adu = input_is_smv || input_is_marccd || input_is_xtc - || ((input_is_crystfel || input_is_cxi) && !crystfel_reader.CountsInPhotons()); + config.counts_in_adu = input_is_smv || input_is_marccd || (input_is_xtc && !xtc_reader.CountsInPhotons()) + || ((input_is_crystfel || input_is_panel_table || input_is_cxi) && !xfel_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(); @@ -3302,6 +3461,10 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { config.prepass_fraction = prepass_fraction; config.forced_rotation_lattice = forced_rotation_lattice; config.refine_geometry = refine_geometry; + config.refine_modules = refine_modules; + config.apply_modules = input_is_panel_table || input_is_cxi || input_is_crystfel; + config.stills_distance_nominal = input_is_panel_table; + config.module_diagnostics = module_diagnostics; config.run_scaling = run_scaling; config.scaling_iter = scaling_iter; config.reference_data = reference_data; @@ -3381,6 +3544,15 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { } g_active_process = nullptr; + // Where the modules are, for a later run's --module-geometry (and any NXmx reader). + if (!output_prefix.empty() && result.modules_placed) { + try { + WriteModuleGeometryFile(output_prefix + "_modules.h5", result.modules, experiment.GetDiffractionGeometry()); + } catch (const std::exception &e) { + logger.Warning("Cannot write {}_modules.h5: {}", output_prefix, e.what()); + } + } + // Every file the run made, from what it returned. if (!output_prefix.empty()) { try { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d5970ae84..da48248f1 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -11,9 +11,12 @@ ADD_EXECUTABLE(jfjoch_portable_test EXCLUDE_FROM_ALL ZSTDCompressorTest.cpp JFJochReaderTest.cpp CrystFELGeometryTest.cpp + XFELReaderTest.cpp SaturationCeilingTest.cpp + DetectorModulesTest.cpp CellPopulationTest.cpp MultiReaderTest.cpp + ModuleRefinementTest.cpp SweepLayoutTest.cpp CBORTest.cpp ImageSpotFinderCPUTest.cpp @@ -105,6 +108,11 @@ ADD_EXECUTABLE(jfjoch_test JFJochReaderTest.cpp MultiReaderTest.cpp CrystFELGeometryTest.cpp + XFELReaderTest.cpp + SaturationCeilingTest.cpp + DetectorModulesTest.cpp + CellPopulationTest.cpp + ModuleRefinementTest.cpp LiveCollectionTest.cpp RugnuxLiveStatusTest.cpp SweepLayoutTest.cpp diff --git a/tests/CrystFELGeometryTest.cpp b/tests/CrystFELGeometryTest.cpp index 9b5ef3f81..baa2a7235 100644 --- a/tests/CrystFELGeometryTest.cpp +++ b/tests/CrystFELGeometryTest.cpp @@ -2,6 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include +#include #include #include @@ -19,19 +20,19 @@ std::string WriteGeom(const std::string &name, const std::string &text) { TEST_CASE("CrystFELGeometry_OldSyntax", "[CrystFELGeometry][portable]") { // One panel and whole-frame bad strips, global keys only: the old syntax. - const auto geom = crystfel::Read("../../tests/crystfel/jf4m.geom"); + std::vector ignored; + const auto geom = crystfel::Read("../../tests/crystfel/jf4m.geom", ignored); REQUIRE(geom.panels.size() == 1); CHECK(geom.panels[0].res == Catch::Approx(13333.3)); - CHECK(geom.panels[0].clen == "0.075"); - CHECK(geom.panels[0].adu_per_photon.value() == Catch::Approx(1.0)); - CHECK(geom.photon_energy == "12400"); + CHECK(geom.panels[0].clen->path.empty()); + CHECK(geom.panels[0].clen->value == Catch::Approx(0.075)); CHECK(geom.data == "/entry/data/data"); - REQUIRE(!geom.bad.empty()); - CHECK(geom.bad[0].panel.empty()); - CHECK(!geom.bad[0].in_lab); - CHECK(geom.bad[0].min_fs == 1030); - CHECK(geom.bad[0].max_fs == 1037); - CHECK(geom.bad[0].max_ss == 2163); + // Only the module geometry is read: gain, photon energy, value flags and bad regions are ignored. + CHECK(!geom.panels[0].adu_per_photon.has_value()); + CHECK(!geom.photon_energy.has_value()); + CHECK(!geom.flag_lessthan.has_value()); + CHECK(geom.bad.empty()); + CHECK(std::find(ignored.begin(), ignored.end(), "photon_energy") != ignored.end()); } TEST_CASE("CrystFELGeometry_Panels", "[CrystFELGeometry][portable]") { @@ -51,8 +52,8 @@ q1/min_ss = 0 q1/max_ss = 1023 q1/fs = +0.000390x -1.000000y q1/ss = +1.000000x +0.000390y -q1/corner_x = 32.9283 -q1/corner_y = 1033.68 +q1/corner_x = 25.5 +q1/corner_y = 1040.25 q1/coffset = 0.000067 q2/min_ss = 1024 q2/max_ss = 2047 @@ -70,18 +71,20 @@ badbox/max_x = 10 badbox/min_y = -5 badbox/max_y = 5 )"); - const auto geom = crystfel::Read(path); + std::vector ignored; + const auto geom = crystfel::Read(path, ignored); REQUIRE(geom.panels.size() == 2); const auto &q1 = geom.panels[0]; CHECK(q1.name == "q1"); - CHECK(q1.clen == "/LCLS/detector_1/EncoderValue"); + CHECK(q1.clen->path == "/LCLS/detector_1/EncoderValue"); + CHECK(q1.clen->scale == Catch::Approx(1e-3)); CHECK(q1.coffset == Catch::Approx(0.000067)); - CHECK(q1.adu_per_eV.value() == Catch::Approx(0.001)); + CHECK(!q1.adu_per_eV.has_value()); CHECK(q1.fs_x == Catch::Approx(0.00039)); CHECK(q1.fs_y == Catch::Approx(-1.0)); CHECK(q1.ss_x == Catch::Approx(1.0)); CHECK(q1.ss_y == Catch::Approx(0.00039)); - CHECK(q1.corner_y == Catch::Approx(1033.68)); + CHECK(q1.corner_y == Catch::Approx(1040.25)); CHECK(q1.max_ss == 1023); const auto &q2 = geom.panels[1]; @@ -89,17 +92,53 @@ badbox/max_y = 5 CHECK(q2.fs_y == Catch::Approx(0.0)); CHECK(q2.ss_y == Catch::Approx(1.0)); CHECK(q2.min_ss == 1024); - CHECK(q2.adu_per_photon.value() == Catch::Approx(9.5)); + CHECK(!q2.adu_per_photon.has_value()); CHECK(q2.res == Catch::Approx(20000)); CHECK(geom.data == "/%/data"); - CHECK(geom.mask_file == "mask.h5"); - CHECK(geom.mask_bad == 1); - REQUIRE(geom.bad.size() == 2); - CHECK(geom.bad[0].panel == "q1"); - CHECK(geom.bad[0].max_fs == 4); - CHECK(geom.bad[1].in_lab); - CHECK(geom.bad[1].min_x == Catch::Approx(-10)); - CHECK(geom.bad[1].max_y == Catch::Approx(5)); + CHECK(geom.masks.empty()); + CHECK(geom.bad.empty()); + for (const std::string key: {"adu_per_eV", "photon_energy", "mask", "mask_file", "mask_good", "mask_bad", + "q2/adu_per_photon", "badq1/min_fs", "badbox/max_y"}) + CHECK(std::find(ignored.begin(), ignored.end(), key) != ignored.end()); + std::filesystem::remove(path); +} + +TEST_CASE("CrystFELGeometry_UnitsAndRefusals", "[CrystFELGeometry][portable]") { + std::vector ignored; + const auto path = WriteGeom("jfjoch_crystfel_geometry_units.geom", R"( +clen = 70.5 mm +photon_energy = 9.3 keV +res = 10000 +mask0_data = /entry/mask +mask0_badbits = 0x6 +p0/min_fs = 0 +p0/max_fs = 99 +p0/min_ss = 0 +p0/max_ss = 99 +p0/mask_edge_pixels = 2 +p0/no_index = false +some_new_keyword = 1 +)"); + const auto geom = crystfel::Read(path, ignored); + CHECK(geom.panels[0].clen->value == Catch::Approx(0.0705)); + CHECK(!geom.photon_energy.has_value()); + CHECK(geom.masks.empty()); + CHECK(std::find(ignored.begin(), ignored.end(), "some_new_keyword") != ignored.end()); + CHECK(std::find(ignored.begin(), ignored.end(), "p0/mask_edge_pixels") != ignored.end()); + + // A camera length is optional. + const auto no_clen = WriteGeom("jfjoch_crystfel_geometry_noclen.geom", + "res = 10000\np0/min_fs = 0\np0/max_fs = 9\np0/min_ss = 0\np0/max_ss = 9\n"); + CHECK(!crystfel::Read(no_clen, ignored).panels[0].clen.has_value()); + std::filesystem::remove(no_clen); + + for (const std::string bad: {"p0/fs = +1.0x +0.1z", "dim1 = 0", "p0/data = /other", "detector_shift_x = /x"}) { + const auto p = WriteGeom("jfjoch_crystfel_geometry_refused.geom", + "clen = 0.1\nres = 10000\nphoton_energy = 9000\np0/min_fs = 0\np0/max_fs = 9\n" + "p0/min_ss = 0\np0/max_ss = 9\n" + bad + "\n"); + CHECK_THROWS(crystfel::Read(p, ignored)); + std::filesystem::remove(p); + } std::filesystem::remove(path); } diff --git a/tests/DetectorModulesTest.cpp b/tests/DetectorModulesTest.cpp new file mode 100644 index 000000000..9a3478a29 --- /dev/null +++ b/tests/DetectorModulesTest.cpp @@ -0,0 +1,112 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include +#include + +#include "../common/DetectorGeometryModular.h" +#include "../common/DetectorModules.h" +#include "../common/DiffractionExperiment.h" +#include "../common/JFJochMath.h" +#include "../reader/NXmxModules.h" +#include "../writer/NXmxModules.h" + +// Two JUNGFRAU modules side by side, 8 pixels apart: each a 1030 x 514 box of the assembled image, placed +// as imaged, so every pixel maps onto itself. +TEST_CASE("DetectorModules_JungfrauAsImaged", "[DetectorModules]") { + const DetectorGeometryModular geom(2, 2, 8, 36, false); + const DetectorModules modules = ModulesAsImaged(geom); + REQUIRE(modules.size() == 2); + CHECK(modules[0].x0 == 0); + CHECK(modules[1].x0 == 1038); + for (const auto &m: modules.Get()) { + CHECK(m.width == 1030); + CHECK(m.height == 514); + } + CHECK_FALSE(modules.Moved()); + + CHECK(modules.ModuleOfImagePixel(1033, 100) == -1); // in the gap + CHECK(modules.ModuleOfImagePixel(1040, 100) == 1); + const auto [px, py] = modules.ImageToPlane(1040.25, 100.5); + CHECK(px == Catch::Approx(1040.25)); + CHECK(py == Catch::Approx(100.5)); + + // Not modular: no modules. + CHECK(ModulesAsImaged(DetectorGeometryModular(1)).empty()); +} + +// A module moved by a rigid step: the image -> plane map moves its pixels, and plane -> image undoes it. +TEST_CASE("DetectorModules_Apply", "[DetectorModules]") { + DetectorModules modules({DetectorModule::AsImaged(0, 0, 100, 50), DetectorModule::AsImaged(108, 0, 100, 50)}); + const auto centres = modules.Centres(); + std::vector step = centres; + step[1].dx = 1.5; + step[1].dy = -0.25; + step[1].rot_rad = 0.2 * PI / 180.0; + modules.Apply(step); + CHECK(modules.Moved()); + + // Module 0 is where it was. + const auto [ax, ay] = modules.ImageToPlane(10, 10); + CHECK(ax == Catch::Approx(10)); + CHECK(ay == Catch::Approx(10)); + + // Module 1's centre moved by the shift, a corner by the shift and the turn. + const auto [cx, cy] = modules[1].PlaneCentre(); + CHECK(cx == Catch::Approx(centres[1].x0 + 1.5)); + CHECK(cy == Catch::Approx(centres[1].y0 - 0.25)); + const double c = std::cos(step[1].rot_rad), s = std::sin(step[1].rot_rad); + const double ux = 108 - centres[1].x0, uy = 0 - centres[1].y0; + const auto [qx, qy] = modules.ImageToPlane(108, 0); + CHECK(qx == Catch::Approx(centres[1].x0 + c * ux - s * uy + 1.5)); + CHECK(qy == Catch::Approx(centres[1].y0 + s * ux + c * uy - 0.25)); + + // And back: a plane point inside module 1's footprint goes to the image pixel that put it there. + const auto [ix, iy] = modules.PlaneToImage(qx, qy); + CHECK(ix == Catch::Approx(108).margin(1e-9)); + CHECK(iy == Catch::Approx(0).margin(1e-9)); + const auto [jx, jy] = modules.PlaneToImage(modules.ImageToPlane(150.3, 20.7).first, + modules.ImageToPlane(150.3, 20.7).second); + CHECK(jx == Catch::Approx(150.3)); + CHECK(jy == Catch::Approx(20.7)); +} + +// Written as NXdetector_module groups and read back, the modules come back as they were - here two +// JUNGFRAU modules, one of them moved off the image's grid, on a detector stored mirrored. +TEST_CASE("DetectorModules_NXmxRoundTrip", "[DetectorModules]") { + DetectorModules modules = ModulesAsImaged(DetectorGeometryModular(2, 2, 8, 36, false)); + auto step = modules.Centres(); + step[0].dx = -0.7; + step[0].dy = 0.35; + step[0].rot_rad = -0.05 * PI / 180.0; + modules.Apply(step); + + DiffractionGeometry geom; + geom.BeamX_pxl(1100).BeamY_pxl(250).DetectorDistance_mm(120).PixelSize_mm(0.075) + .Orientation(DetectorOrientation(true, 1)); + + const std::string path = "detector_modules_roundtrip.h5"; + WriteModuleGeometryFile(path, modules, geom); + DetectorModules back; + { + HDF5ReadOnlyFile file(path); + back = ReadNXmxModules(file, 0.075e-3, geom.GetOrientation()); + } + std::filesystem::remove(path); + + REQUIRE(back.size() == 2); + for (size_t i = 0; i < 2; i++) { + const auto &a = modules[i], &b = back[i]; + CHECK(b.x0 == a.x0); + CHECK(b.y0 == a.y0); + CHECK(b.width == a.width); + CHECK(b.height == a.height); + CHECK(b.ox == Catch::Approx(a.ox).margin(1e-3)); + CHECK(b.oy == Catch::Approx(a.oy).margin(1e-3)); + CHECK(b.fx == Catch::Approx(a.fx).margin(1e-6)); + CHECK(b.fy == Catch::Approx(a.fy).margin(1e-6)); + CHECK(b.sx == Catch::Approx(a.sx).margin(1e-6)); + CHECK(b.sy == Catch::Approx(a.sy).margin(1e-6)); + } +} diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index b8705be01..1c83d63f7 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -4836,29 +4836,44 @@ TEST_CASE("JFJochReader_DECTRISMultichannel", "[HDF5][Full]") { } // LCLS-I XTC: datagrams of a 40-byte header and a tree of Xtcs; a camera image is a Camera::FrameV1 -// leaf. The fixture writes a Configure and two L1Accepts. The first carries a second, smaller camera -// beside the diffraction detector, which the reader must pass over; the geometry comes from a -// one-panel DIALS experiment list beside the stream. +// leaf. The fixture writes a Configure carrying a Rayonix configuration, and L1Accepts with the +// electron beam's photon energy. The first event carries a second, smaller camera beside the +// diffraction detector, which the reader must pass over; one event's image is damaged, and the stream +// ends inside a last event, as a partly downloaded file does. namespace { void PutWords(std::string &s, std::initializer_list words) { for (uint32_t w: words) s.append(reinterpret_cast(&w), 4); } - std::string XtcLeaf(uint32_t type, const std::string &payload) { + std::string XtcLeaf(uint32_t type, const std::string &payload, uint32_t source = 0, uint32_t damage = 0) { std::string s; - PutWords(s, {0, 0, 0, type, static_cast(20 + payload.size())}); + PutWords(s, {damage, 0, source, type, static_cast(20 + payload.size())}); return s + payload; } - std::string XtcFrame(uint32_t width, uint32_t height, uint16_t base) { + std::string XtcFrame(uint32_t width, uint32_t height, uint16_t base, uint32_t source, uint32_t damage = 0) { std::string p; - PutWords(p, {width, height, 16, 0}); + PutWords(p, {width, height, 16, 3}); // an offset of 3 counts, taken off when read for (uint32_t i = 0; i < width * height; i++) { - const uint16_t v = base + i; + const uint16_t v = base + i % 100; p.append(reinterpret_cast(&v), 2); } - return XtcLeaf(2 | (1u << 16), p); + return XtcLeaf(2 | (1u << 16), p, source, damage); + } + + std::string XtcEBeam(double photon_energy_eV) { + std::string p(4 + 20 * 8, '\0'); + std::memcpy(p.data() + 4 + 17 * 8, &photon_energy_eV, 8); + return XtcLeaf(15 | (7u << 16), p); + } + + std::string XtcRayonixConfig(const std::string &model, uint8_t binning) { + std::string p(60, '\0'); + p[0] = static_cast(binning); + p[1] = static_cast(binning); + p.replace(20, model.size(), model); + return XtcLeaf(73 | (2u << 16), p); } std::string XtcDatagram(uint32_t transition, const std::string &payload) { @@ -4869,55 +4884,50 @@ namespace { } } -TEST_CASE("JFJochXTCReader_Geometry", "[portable]") { - const uint32_t nx = 24, ny = 16; +TEST_CASE("JFJochXTCReader", "[portable]") { + // 1920 wide at 4x4 binning: a whole MX340-HS readout, so 340 mm / 1920 pixels. + const uint32_t nx = 1920, ny = 4; + std::filesystem::remove_all("xtc_test"); std::filesystem::create_directories("xtc_test"); { std::ofstream f("xtc_test/run-s00-c00.xtc", std::ios::binary); - f << XtcDatagram(4, ""); - f << XtcDatagram(12, XtcLeaf(1, XtcLeaf(13, std::string(32, '\0')) + XtcFrame(4, 4, 7)) - + XtcLeaf(1, XtcFrame(nx, ny, 100))); - f << XtcDatagram(12, XtcLeaf(1, XtcFrame(nx, ny, 200))); - } - { - // The beam meets the panel 12 pixels along fast and 8 along slow from the first pixel's - // corner, 150 mm away; 0.1 mm pixels. - std::ofstream f("xtc_test/geometry.expt"); - f << R"({"beam": [{"direction": [0, 0, 1], "wavelength": 1.3, "polarization_normal": [0, 1, 0], - "polarization_fraction": 0.999}], - "detector": [{"panels": [{"fast_axis": [1, 0, 0], "slow_axis": [0, -1, 0], - "origin": [-1.2, 0.8, -150], "pixel_size": [0.1, 0.1], - "image_size": [24, 16]}]}]})"; + f << XtcDatagram(4, XtcLeaf(1, XtcRayonixConfig("MX340-HS:125", 4))); + f << XtcDatagram(12, XtcLeaf(1, XtcLeaf(13, std::string(32, '\0')) + XtcFrame(4, 4, 7, 2) + XtcEBeam(9500)) + + XtcLeaf(1, XtcFrame(nx, ny, 100, 1))); + f << XtcDatagram(12, XtcLeaf(1, XtcFrame(nx, ny, 300, 1, 0x4000) + XtcEBeam(9520))); // damaged + f << XtcDatagram(12, XtcLeaf(1, XtcFrame(nx, ny, 200, 1) + XtcEBeam(9510))); + const std::string cut = XtcDatagram(12, XtcLeaf(1, XtcFrame(nx, ny, 400, 1))); + f << cut.substr(0, cut.size() / 2); } REQUIRE(JFJochXTCReader::CanRead("xtc_test/run-s00-c00.xtc")); CHECK_FALSE(JFJochXTCReader::CanRead("xtc_test")); // a directory is not read - CHECK_FALSE(JFJochXTCReader::CanRead("xtc_test/geometry.expt")); - JFJochXTCReader reader; - REQUIRE_NOTHROW(reader.ReadFiles({"xtc_test/run-s00-c00.xtc"})); - CHECK(reader.GetNumberOfImages() == 2); - const auto x = reader.GetDataset()->experiment; - CHECK(x.GetXPixelsNum() == nx); - CHECK(x.GetYPixelsNum() == ny); - CHECK(x.GetDetectorDistance_mm() == Catch::Approx(150.0)); - CHECK(x.GetPixelSize_mm() == Catch::Approx(0.1)); - CHECK(x.GetWavelength_A() == Catch::Approx(1.3).epsilon(1e-4)); - CHECK(x.GetBeamX_pxl() == Catch::Approx(12.0)); - CHECK(x.GetBeamY_pxl() == Catch::Approx(8.0)); - CHECK_FALSE(x.GetGoniometer().has_value()); // stills - REQUIRE(x.GetPolarizationFactor().has_value()); - CHECK(x.GetPolarizationFactor().value() == Catch::Approx(0.998)); - - for (int i = 0; i < 2; i++) { - auto raw = reader.GetRawImage(i); - REQUIRE(raw); - CHECK(raw->image.GetWidth() == nx); - CHECK(raw->image.GetHeight() == ny); - const auto *pixels = reinterpret_cast(raw->image_buffer.data()); - CHECK(pixels[0] == 100 * (i + 1)); - CHECK(pixels[nx * ny - 1] == 100 * (i + 1) + nx * ny - 1); + { + // The XTC's own beam records and Rayonix configuration. + JFJochXTCReader reader; + REQUIRE_NOTHROW(reader.ReadFiles({"xtc_test/run-s00-c00.xtc"})); + CHECK(reader.GetNumberOfImages() == 2); + CHECK(reader.HasWavelength()); + CHECK(reader.HasPixelSize()); + const auto x = reader.GetDataset()->experiment; + CHECK(x.GetXPixelsNum() == nx); + CHECK(x.GetPixelSize_mm() == Catch::Approx(340.0 / 1920.0)); + // The damaged image's beam record still counts: three records, 9500, 9520 and 9510 eV. + CHECK(x.GetIncidentEnergy_keV() == Catch::Approx(9.510).epsilon(1e-5)); + CHECK(x.GetBeamX_pxl() == Catch::Approx(nx / 2.0)); + for (int i = 0; i < 2; i++) { + auto raw = reader.GetRawImage(i); + REQUIRE(raw); + const auto *pixels = reinterpret_cast(raw->image_buffer.data()); + CHECK(pixels[0] == 100 * (i + 1) - 3); + CHECK(pixels[99] == 100 * (i + 1) + 99 - 3); + } + reader.ApplyGain(2.0f); + CHECK(reader.CountsInPhotons()); + auto raw = reader.GetRawImage(0); + CHECK(reinterpret_cast(raw->image_buffer.data())[1] == 49); // (100 + 1 - 3) / 2 } - reader.Close(); + std::filesystem::remove_all("xtc_test"); } diff --git a/tests/ModuleRefinementTest.cpp b/tests/ModuleRefinementTest.cpp new file mode 100644 index 000000000..5b4886055 --- /dev/null +++ b/tests/ModuleRefinementTest.cpp @@ -0,0 +1,94 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include +#include + +#include "../common/JFJochMath.h" +#include "../image_analysis/geom_refinement/GeometryRefiner.h" +#include "../image_analysis/bragg_prediction/BraggPrediction.h" + +// A detector of four modules (the quadrants around the beam), each displaced by a known rigid move: the +// bundle adjustment with module placement must find the moves, and the fit without it must be worse. +TEST_CASE("RefineGlobalGeometry_ModulePlacement") { + DiffractionExperiment exp; + exp.IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(150); + const DiffractionGeometry geom = exp.GetDiffractionGeometry(); + + // Moves that hold none of what the gauge takes away: zero in sum, in turn and outwards. + std::vector truth(4); + const double centre[4][2] = {{500, 500}, {1500, 500}, {500, 1500}, {1500, 1500}}; + const double move[4][3] = {{2.0, -2.0, 0.2}, {2.0, 2.0, -0.2}, {-2.0, -2.0, 0.1}, {-2.0, 2.0, -0.1}}; + for (int m = 0; m < 4; m++) { + truth[m].x0 = centre[m][0]; + truth[m].y0 = centre[m][1]; + truth[m].dx = move[m][0]; + truth[m].dy = move[m][1]; + truth[m].rot_rad = move[m][2] * PI / 180.0; + } + auto module_of = [](double x, double y) { return (x < 1000 ? 0 : 1) + (y < 1000 ? 0 : 2); }; + + const CrystalLattice latt(50, 60, 70, 90, 90, 90); + BraggPredictionSettings settings{.high_res_A = 2.0, .ewald_dist_cutoff = 0.001}; + std::mt19937 rng(7); + std::uniform_real_distribution u(-1.0f, 1.0f); + + std::vector frames; + for (int f = 0; f < 60; f++) { + const Coord axis = Coord(u(rng), u(rng), u(rng)).Normalize(); + const CrystalLattice rotated = latt.Multiply(RotMatrix(static_cast(PI) * u(rng), axis)); + BraggPrediction prediction; + const int n = prediction.Calc(exp, rotated, settings); + GeomRefineFrame frame; + frame.lattice = rotated; + frame.ordinal = f; + for (int i = 0; i < n; i++) { + const auto &r = prediction.GetReflections().at(i); + const int m = module_of(r.predicted_x, r.predicted_y); + // where the image shows a spot whose true place is the prediction: the inverse of the move + const auto &t = truth[m]; + const double c = std::cos(t.rot_rad), s = std::sin(t.rot_rad); + const double vx = r.predicted_x - t.x0 - t.dx, vy = r.predicted_y - t.y0 - t.dy; + GeomRefineSpot spot; + spot.x = static_cast(t.x0 + c * vx + s * vy); + spot.y = static_cast(t.y0 - s * vx + c * vy); + spot.h = r.h; + spot.k = r.k; + spot.l = r.l; + spot.module = static_cast(m); + frame.spots.push_back(spot); + } + if (frame.spots.size() >= 20) + frames.push_back(frame); + } + REQUIRE(frames.size() > 30); + + GeometryRefinerSettings st; + st.crystal_system = gemmi::CrystalSystem::Orthorhombic; + st.num_threads = 2; + for (const auto &t: truth) { + ModuleShift start; + start.x0 = t.x0; + start.y0 = t.y0; + st.modules.push_back(start); + } + + const auto held = RefineGlobalGeometry(geom, latt.GetUnitCell(), frames, st); + REQUIRE(held.ok); + + st.refine_modules = true; + const auto placed = RefineGlobalGeometry(geom, latt.GetUnitCell(), frames, st); + REQUIRE(placed.ok); + CHECK(placed.median_residual_px < 0.1); + CHECK(held.median_residual_px > 5 * placed.median_residual_px); + + for (int m = 0; m < 4; m++) { + CHECK(placed.modules[m].dx == Catch::Approx(truth[m].dx).margin(0.1)); + CHECK(placed.modules[m].dy == Catch::Approx(truth[m].dy).margin(0.1)); + CHECK(placed.modules[m].rot_rad * 180.0 / PI == Catch::Approx(move[m][2]).margin(0.01)); + } + CHECK(placed.beam_x_px == Catch::Approx(1000).margin(0.1)); + CHECK(placed.beam_y_px == Catch::Approx(1000).margin(0.1)); + CHECK(placed.distance_mm == Catch::Approx(150).margin(0.05)); +} diff --git a/tests/XFELReaderTest.cpp b/tests/XFELReaderTest.cpp new file mode 100644 index 000000000..8bd0caa49 --- /dev/null +++ b/tests/XFELReaderTest.cpp @@ -0,0 +1,184 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include +#include +#include +#include + +#include "../reader/JFJochXFELReader.h" + +// Synthetic SACLA Cheetah and Cheetah CXI files, written with the HDF5 C API the way those pipelines +// write them: whole frames as single deflated (and byte-shuffled) chunks, so the direct-chunk path is +// what is read. +namespace { + +void WriteFrames(hid_t loc, const std::string &name, hid_t type, const void *data, std::vector dims) { + const hid_t space = H5Screate_simple(static_cast(dims.size()), dims.data(), nullptr); + std::vector chunk = dims; + if (chunk.size() == 3) + chunk[0] = 1; + const hid_t dcpl = H5Pcreate(H5P_DATASET_CREATE); + H5Pset_chunk(dcpl, static_cast(chunk.size()), chunk.data()); + H5Pset_shuffle(dcpl); + H5Pset_deflate(dcpl, 4); + const hid_t ds = H5Dcreate2(loc, name.c_str(), type, space, H5P_DEFAULT, dcpl, H5P_DEFAULT); + H5Dwrite(ds, type, H5S_ALL, H5S_ALL, H5P_DEFAULT, data); + H5Dclose(ds); + H5Pclose(dcpl); + H5Sclose(space); +} + +void WriteVector(hid_t loc, const std::string &name, const std::vector &v) { + const hsize_t n = v.size(); + const hid_t space = H5Screate_simple(1, &n, nullptr); + const hid_t ds = H5Dcreate2(loc, name.c_str(), H5T_NATIVE_DOUBLE, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + H5Dwrite(ds, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, v.data()); + H5Dclose(ds); + H5Sclose(space); +} + +void WriteScalar(hid_t loc, const std::string &name, double v) { + const hid_t space = H5Screate(H5S_SCALAR); + const hid_t ds = H5Dcreate2(loc, name.c_str(), H5T_NATIVE_DOUBLE, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + H5Dwrite(ds, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &v); + H5Dclose(ds); + H5Sclose(space); +} + +int32_t PixelAt(JFJochXFELReader &reader, int64_t image, size_t col, size_t row) { + auto raw = reader.GetRawImage(image); + REQUIRE(raw); + return reinterpret_cast(raw->image_buffer.data())[row * raw->image.GetWidth() + col]; +} + +} // namespace + +TEST_CASE("JFJochXFELReader_SACLA", "[portable]") { + // Two sensors of 8 rows x 4 columns stacked in a 16 x 4 frame, 50 um pixels. The first lies square + // with its first corner on the beam; the second is turned by 90 degrees with its corner 20 pixels + // along x (posx = -1000 um), so its fast axis runs along y and its slow axis along -x. + const std::string name = "xfel_reader_sacla.h5"; + { + const hid_t f = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + const hid_t meta = H5Gcreate2(f, "metadata", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + WriteVector(meta, "posx_in_um", {0, -1000}); + WriteVector(meta, "posy_in_um", {0, 0}); + WriteVector(meta, "posz_in_um", {0, 0}); + WriteVector(meta, "angle_in_rad", {0, 90}); // degrees, whatever the name says + WriteVector(meta, "pixelsizex_in_um", {50, 50}); + WriteVector(meta, "pixelsizey_in_um", {50, 50}); + WriteScalar(meta, "distance_in_mm", 60.0); + std::vector pixelmask(16 * 4, 0); + pixelmask[1] = 1; // first sensor, row 0, column 1 + WriteFrames(meta, "pixelmask", H5T_NATIVE_UINT16, pixelmask.data(), {16, 4}); + H5Gclose(meta); + for (int shot = 0; shot < 2; shot++) { + const std::string group = "tag-" + std::to_string(100 + shot); + const hid_t g = H5Gcreate2(f, group.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + std::vector frame(16 * 4, 10); + frame[0] = static_cast(20 + shot); // first sensor (fs 0, ss 0) + frame[8 * 4] = -5; // second sensor (fs 0, ss 0): below zero + frame[8 * 4 + 1] = 33; // second sensor (fs 1, ss 0) + frame[9 * 4] = 44; // second sensor (fs 0, ss 1) + WriteFrames(g, "data", H5T_NATIVE_INT16, frame.data(), {16, 4}); + WriteScalar(g, "photon_energy_ev", 10000.0 + 10.0 * shot); + H5Gclose(g); + } + H5Fclose(f); + } + + REQUIRE(JFJochXFELReader::CanReadPanelTable(name)); + CHECK_FALSE(JFJochXFELReader::CanReadCXI(name)); + JFJochXFELReader reader; + REQUIRE_NOTHROW(reader.ReadPanelTable({name})); + CHECK(reader.GetNumberOfImages() == 2); + CHECK_FALSE(reader.CountsInPhotons()); + const auto x = reader.GetDataset()->experiment; + CHECK(x.GetIncidentEnergy_keV() == Catch::Approx(10.005)); + CHECK(x.GetDetectorDistance_mm() == Catch::Approx(60.0)); + CHECK(x.GetPixelSize_mm() == Catch::Approx(0.05)); + // The first sensor's first pixel centre is at (0.5, 0.5) from the beam, which the flat image puts on + // a whole column: the beam at 1.5, the pixel at column and row 2. + CHECK(x.GetBeamX_pxl() == Catch::Approx(1.5)); + CHECK(x.GetBeamY_pxl() == Catch::Approx(1.5)); + CHECK(PixelAt(reader, 0, 2, 2) == 20); + CHECK(PixelAt(reader, 1, 2, 2) == 21); + // Second sensor: (fs 0, ss 0) centre at x = 20 - 0.5, y = 0.5; fs runs along +y, ss along -x. + CHECK(PixelAt(reader, 0, 21, 2) == -5); + CHECK(PixelAt(reader, 0, 21, 3) == 33); + CHECK(PixelAt(reader, 0, 20, 2) == 44); + // The file's mask holds the first sensor's (fs 1, ss 0), at column 3, row 2. + const auto &mask = reader.GetDataset()->pixel_mask->GetMask(); + CHECK(mask[2 * x.GetXPixelsNum() + 3] != 0); + CHECK(mask[2 * x.GetXPixelsNum() + 2] == 0); + + reader.ApplyGain(2.0f); + CHECK(reader.CountsInPhotons()); + CHECK(PixelAt(reader, 0, 2, 2) == 10); + reader.Close(); + std::filesystem::remove(name); +} + +TEST_CASE("JFJochXFELReader_CXI", "[portable]") { + // A single-panel Cheetah CXI file: three 6 x 6 frames of signed 32-bit counts, 75 um pixels, + // 9 keV stated in joules, a mask_shared with one bad pixel. + const std::string name = "xfel_reader_test.cxi"; + { + const hid_t f = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + for (const char *g: {"entry_1", "entry_1/data_1", "entry_1/instrument_1", "entry_1/instrument_1/detector_1", + "entry_1/instrument_1/source_1"}) + H5Gclose(H5Gcreate2(f, g, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)); + std::vector frames(3 * 36, 5); + frames[0] = -7; // frame 0, pixel (0, 0): dark-subtracted below zero + frames[36 + 7] = 70000; // frame 1, pixel (1, 1) + WriteFrames(f, "entry_1/data_1/data", H5T_NATIVE_INT32, frames.data(), {3, 6, 6}); + const std::string det = "entry_1/instrument_1/detector_1/"; + WriteVector(f, det + "x_pixel_size", {75e-6, 75e-6, 75e-6}); + WriteVector(f, det + "distance", {0.1, 0.1, 0.1}); + WriteVector(f, "entry_1/instrument_1/source_1/energy", std::vector(3, 9000.0 * 1.602176634e-19)); + std::vector mask(36, 0); + mask[35] = 8; // dead + WriteFrames(f, det + "mask_shared", H5T_NATIVE_UINT16, mask.data(), {6, 6}); + H5Fclose(f); + } + + REQUIRE(JFJochXFELReader::CanReadCXI(name)); + CHECK_FALSE(JFJochXFELReader::CanReadPanelTable(name)); + CHECK_FALSE(JFJochXFELReader::CXIRefusal(name).has_value()); + JFJochXFELReader reader; + REQUIRE_NOTHROW(reader.ReadCXI({name})); + CHECK(reader.GetNumberOfImages() == 3); + const auto x = reader.GetDataset()->experiment; + CHECK(x.GetIncidentEnergy_keV() == Catch::Approx(9.0)); + CHECK(x.GetPixelSize_mm() == Catch::Approx(0.075)); + CHECK(x.GetDetectorDistance_mm() == Catch::Approx(100.0)); + // Corner at (-3, -3): the pixel (fs, ss) lands at column and row 2 + fs, 2 + ss. + CHECK(PixelAt(reader, 0, 2, 2) == -7); + CHECK(PixelAt(reader, 1, 3, 3) == 70000); + CHECK(PixelAt(reader, 2, 4, 4) == 5); + const auto &m = reader.GetDataset()->pixel_mask->GetMask(); + CHECK(m[7 * x.GetXPixelsNum() + 7] != 0); + CHECK(m[6 * x.GetXPixelsNum() + 6] == 0); + reader.Close(); + std::filesystem::remove(name); +} + +TEST_CASE("JFJochXFELReader_CXI_TiledRefused", "[portable]") { + // A Cheetah CXI slab of a tiled detector is not square, and is not read as one flat panel. + const std::string name = "xfel_reader_tiled.cxi"; + { + const hid_t f = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + for (const char *g: {"entry_1", "entry_1/data_1"}) + H5Gclose(H5Gcreate2(f, g, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)); + std::vector frames(2 * 8 * 4, 1); + WriteFrames(f, "entry_1/data_1/data", H5T_NATIVE_INT16, frames.data(), {2, 8, 4}); + H5Fclose(f); + } + CHECK_FALSE(JFJochXFELReader::CanReadCXI(name)); + const auto refusal = JFJochXFELReader::CXIRefusal(name); + REQUIRE(refusal.has_value()); + CHECK(refusal->find("tiled detector") != std::string::npos); + std::filesystem::remove(name); +} diff --git a/writer/CMakeLists.txt b/writer/CMakeLists.txt index ad4a0a3fa..003a5eaab 100644 --- a/writer/CMakeLists.txt +++ b/writer/CMakeLists.txt @@ -10,6 +10,7 @@ TARGET_LINK_LIBRARIES(JFJochHDF5Wrappers Compression hdf5-static) ADD_LIBRARY(JFJochWriter STATIC HDF5DataFile.h HDF5DataFile.cpp HDF5NXmx.cpp HDF5NXmx.h + NXmxModules.cpp NXmxModules.h FileWriter.cpp FileWriter.h HDF5DataFilePlugin.cpp HDF5DataFilePlugin.h diff --git a/writer/NXmxModules.cpp b/writer/NXmxModules.cpp new file mode 100644 index 000000000..dd822e860 --- /dev/null +++ b/writer/NXmxModules.cpp @@ -0,0 +1,70 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "NXmxModules.h" + +#include + +namespace { + +// Internal frame (x = column, y = row downward, z = beam) -> McStas, a half turn about z. +std::vector ToMcStas(const Coord &v) { + return {0.0 - v.x, 0.0 - v.y, 0.0 + v.z}; +} + +// A translation of `length` metres along v (pixels, internal frame), as NXmx writes one: the length as +// the value, the direction as a unit vector - and a zero translation with a unit vector all the same. +void SaveTranslation(HDF5Group &group, const std::string &name, const RotMatrix &delta, double dx_px, double dy_px, + double pixel_size_m, const std::string &depends_on) { + const double len_px = std::hypot(dx_px, dy_px); + const Coord dir = len_px > 0 ? Coord(static_cast(dx_px / len_px), static_cast(dy_px / len_px), 0) + : Coord(0, 0, 1); + SaveScalar(group, name, len_px * pixel_size_m)-> + Transformation("m", depends_on, "", "", "translation", ToMcStas(delta * dir), {0, 0, 0}, ""); +} + +} // namespace + +void WriteNXmxModules(HDF5Object &file, const DetectorModules &modules, double pixel_size_m, + const DetectorOrientation &orientation, const std::string &depends_on) { + const RotMatrix delta = orientation.Matrix(); + for (size_t n = 0; n < modules.size(); n++) { + const DetectorModule &m = modules[n]; + const std::string path = "/entry/instrument/detector/module_" + std::to_string(n); + HDF5Group group(file, path); + group.NXClass("NXdetector_module"); + group.SaveVector("data_origin", std::vector{static_cast(m.y0), static_cast(m.x0)}); + group.SaveVector("data_size", std::vector{static_cast(m.height), + static_cast(m.width)}); + SaveTranslation(group, "module_offset", delta, m.ox, m.oy, pixel_size_m, depends_on); + SaveTranslation(group, "fast_pixel_direction", delta, m.fx, m.fy, pixel_size_m, path + "/module_offset"); + SaveTranslation(group, "slow_pixel_direction", delta, m.sx, m.sy, pixel_size_m, path + "/module_offset"); + } +} + +void WriteModuleGeometryFile(const std::string &path, const DetectorModules &modules, + const DiffractionGeometry &geom) { + HDF5File file(path); + HDF5Group(file, "/entry").NXClass("NXentry"); + HDF5Group(file, "/entry/instrument").NXClass("NXinstrument"); + HDF5Group(file, "/entry/instrument/detector").NXClass("NXdetector"); + HDF5Group transformations(file, "/entry/instrument/detector/transformations"); + transformations.NXClass("NXtransformations"); + + // Sample -> image pixel (0, 0), square to the beam: what the modules' offsets are measured from. The + // tilt, if any, is the master file's to state; this file carries the modules only. + const double pixel_m = geom.GetPixelSize_mm() / 1000.0; + const RotMatrix delta = geom.GetOrientation().Matrix(); + const auto v = ToMcStas(delta * Coord(-geom.GetBeamX_pxl() * static_cast(pixel_m), + -geom.GetBeamY_pxl() * static_cast(pixel_m), + geom.GetDetectorDistance_mm() / 1000.0f)); + const double len = std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + SaveScalar(transformations, "translation", len)-> + Transformation("m", ".", "detector", "detector_arm", "translation", + std::vector{v[0] / len, v[1] / len, v[2] / len}); + SaveScalar(file, "/entry/instrument/detector/depends_on", + std::string("/entry/instrument/detector/transformations/translation")); + + WriteNXmxModules(file, modules, pixel_m, geom.GetOrientation(), + "/entry/instrument/detector/transformations/translation"); +} diff --git a/writer/NXmxModules.h b/writer/NXmxModules.h new file mode 100644 index 000000000..f2790855e --- /dev/null +++ b/writer/NXmxModules.h @@ -0,0 +1,30 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "../common/DetectorModules.h" +#include "../common/DetectorOrientation.h" +#include "../common/DiffractionGeometry.h" +#include "HDF5Objects.h" + +// The detector's modules as NXmx describes them: one NXdetector_module group per module, +// /entry/instrument/detector/module_, after the NeXus base class NXdetector_module +// (https://manual.nexusformat.org/classes/base_classes/NXdetector_module.html) as NXmx uses it +// (https://manual.nexusformat.org/classes/applications/NXmx.html): +// - data_origin, data_size: where the module is in the image array, slowest index first ([row, column]); +// - module_offset: translation from the image's pixel (0, 0) to the module's first pixel, depends_on +// the detector's own translation (the sample -> pixel (0, 0) vector the master writes); +// - fast_pixel_direction, slow_pixel_direction: translation of one pixel along the module's columns and +// rows, depends_on module_offset. +// Vectors are in the McStas frame, lengths in metres. A pixel (row i, column j) of a module is then at +// translation + module_offset + j fast + i slow - the convention dxtbx reads. +void WriteNXmxModules(HDF5Object &file, const DetectorModules &modules, double pixel_size_m, + const DetectorOrientation &orientation, const std::string &depends_on); + +// A file holding only the module geometry (and the detector translation it hangs from), for +// rugnux --module-geometry. +void WriteModuleGeometryFile(const std::string &path, const DetectorModules &modules, + const DiffractionGeometry &geom);