From 7c23b4befee2ec4d6e0196b6c9911811b02c3fe6 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 21 Oct 2026 12:42:08 +0200 Subject: [PATCH] CrystFEL geometry: read the pixel masks and bad regions instead of ignoring them The .geom parser now takes the masking information the format defines: - mask arrays: mask/mask_file/mask_good/mask_bad (mask 0) and maskN_data/maskN_file/maskN_goodbits/ maskN_badbits for N = 0..7; global keys hold for the panels after them, a panel may override. A mask that is the same on every panel is one whole-frame mask, otherwise one per panel, applied on that panel's pixels only. - a relative mask_file is looked for in the working directory (as CrystFEL does), then beside the .geom. - bad regions (bad*): a lab-plane box (min_x..max_y) or a fs/ss range on a panel (or, old syntax, on the whole frame); mask_edge_pixels and no_index. JFJochXFELReader already had the static pixel mask and the per-frame mask path for these; the masks go through them unchanged. A mask whose path holds % (one per shot group) is applied frame by frame like a 3D one. Where the geometry names mask arrays, the data file's own Cheetah masks (CXI mask_shared and per-shot mask with Cheetah's generic bad bits, /metadata/pixelmask) are not added: the geometry says which bits matter (e.g. an AGIPD mask_bad that includes bit 32). --mask-good / --mask-bad still add their masks to these. Tests: [CrystFELGeometry] (new keys, per-panel override, relative mask file), JFJochXFELReader_* (new GeometryMasks case on a synthetic CXI + .geom), [JFJochMultiReader], JFJochXTCReader*. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- docs/RUGNUX_ADVANCED.md | 2 +- docs/RUGNUX_FORMATS.md | 19 ++++-- reader/CrystFELGeometry.cpp | 100 +++++++++++++++++++++++++-- reader/CrystFELGeometry.h | 19 +++--- reader/DetectorLayout.h | 5 +- reader/JFJochXFELReader.cpp | 57 +++++++++++----- reader/JFJochXFELReader.h | 2 + rugnux/rugnux_cli.cpp | 2 +- tests/CrystFELGeometryTest.cpp | 103 +++++++++++++++++++++++++--- tests/XFELReaderTest.cpp | 120 +++++++++++++++++++++++++++++++++ 10 files changed, 380 insertions(+), 49 deletions(-) diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index cf26ddde2..f4f0c1ad7 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -495,7 +495,7 @@ Geometry overrides (defaults are taken from the input file; override them to rep | `--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)) | +| `--crystfel-geometry ` | Place the panels of an HDF5/CXI input by a CrystFEL geometry file instead of the file's own metadata: module geometry, camera length, mask arrays and bad regions (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 the frames' stored layout, in `` or in each data file; usable only where non-zero (or with all ``) / unusable where non-zero (or with any ``). Repeatable ([Serial XFEL inputs](RUGNUX_FORMATS.md#serial-xfel-inputs)) | diff --git a/docs/RUGNUX_FORMATS.md b/docs/RUGNUX_FORMATS.md index 3b87e75b6..945cfcb1f 100644 --- a/docs/RUGNUX_FORMATS.md +++ b/docs/RUGNUX_FORMATS.md @@ -102,10 +102,17 @@ Damaged or compressed XTC images are left out, and so is a last event that a par 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 +`--crystfel-geometry ` places the panels of an HDF5/CXI input by a CrystFEL geometry file. Taken +from it are the module geometry (panel pixel ranges, fast/slow directions, corners, pixel size, +per-panel `coffset`, the frames' data path), optionally the camera length, and the pixels to leave out: +the mask arrays (`mask`, `mask_file`, `mask_good`, `mask_bad`, or up to eight as `maskN_data`, +`maskN_file`, `maskN_goodbits`, `maskN_badbits`; a panel may state its own), the bad regions (`bad*`, +a box in the laboratory plane or a pixel range on a panel), `mask_edge_pixels` and `no_index`. A mask +dataset with one plane per frame, or with `%` in its path, is applied frame by frame. A relative +`mask_file` is looked for in the working directory, as CrystFEL does, and then beside the geometry file. +Where the geometry names mask arrays, the data file's own (the CXI `mask_shared` and per-shot mask, a +panel table's `/metadata/pixelmask`) are not taken: the geometry says which bits matter. Everything else +it may hold (gain, saturation, 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 @@ -113,8 +120,8 @@ given with `--mask-good [:][@]` (a pixel is usable only whe or has all of ``) or `--mask-bad [:][@]` (unusable where non-zero, or where it has any of ``, e.g. `@0xfeff`), in the frames' own layout, as stored. Without `` the dataset is read from each data file - a per-frame mask stored beside the frames. 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). A gain the background shows at or below one detector unit +and add to the masks the data files or the geometry file state (the CXI `mask_shared` and per-shot mask, +the `/metadata/pixelmask` of a panel-table file, a `--crystfel-geometry` file's masks). A gain the background shows at or below one detector unit per photon is not applied: such frames stay in detector units. ### Module geometry diff --git a/reader/CrystFELGeometry.cpp b/reader/CrystFELGeometry.cpp index f66e39e22..2f61c68e7 100644 --- a/reader/CrystFELGeometry.cpp +++ b/reader/CrystFELGeometry.cpp @@ -3,6 +3,9 @@ #include "CrystFELGeometry.h" +#include +#include +#include #include #include #include @@ -98,8 +101,8 @@ bool Bookkeeping(const std::string &key) { // 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. +// Sets one panel property; false if the key is not one. Where the panel is - its pixel range, axes, +// corner, pixel size, depth offset and camera length - and which of its pixels are left out. 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); @@ -112,11 +115,58 @@ bool SetPanel(const std::string &filename, Panel &p, const std::string &key, con else if (key == "res") p.res = std::stod(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 == "mask_edge_pixels") p.mask_edge_pixels = std::stoll(value); + else if (key == "no_index") p.no_index = value == "1" || value == "true"; else return false; return true; } +// The up to eight masks a panel has (maskN_data, maskN_file, maskN_goodbits, maskN_badbits, N = 0..7; +// mask, mask_file, mask_good and mask_bad are mask 0's). Sets the one the key names; false if the key +// is not a mask key. +using Masks = std::array; +bool SetMask(Masks &masks, const std::string &key, const std::string &value) { + static const std::map OLD = { + {"mask", "data"}, {"mask_file", "file"}, {"mask_good", "goodbits"}, {"mask_bad", "badbits"}}; + static const std::regex CURRENT(R"(mask([0-7])_(data|file|goodbits|badbits))"); + size_t n = 0; + std::string field; + std::smatch match; + if (OLD.contains(key)) + field = OLD.at(key); + else if (std::regex_match(key, match, CURRENT)) { + n = std::stoul(match[1].str()); + field = match[2].str(); + } else + return false; + if (field == "data") masks[n].data = value; + else if (field == "file") masks[n].file = value; + else if (field == "goodbits") masks[n].good = static_cast(std::stoul(value, nullptr, 0)); + else masks[n].bad = static_cast(std::stoul(value, nullptr, 0)); + return true; +} + +// Sets one property of a bad region: a box in the laboratory plane (min_x .. max_y, in pixels) or a +// pixel range in the stored frame (min_fs .. max_ss, inclusive) on one panel - or, in the old syntax, +// with no panel, on the whole frame. False if the key is not one. +bool SetBad(xfel::BadRegion &b, const std::string &key, const std::string &value) { + if (key == "min_x") b.min_x = std::stod(value); + else if (key == "max_x") b.max_x = std::stod(value); + else if (key == "min_y") b.min_y = std::stod(value); + else if (key == "max_y") b.max_y = std::stod(value); + else if (key == "min_fs") b.min_fs = std::stoll(value); + else if (key == "max_fs") b.max_fs = std::stoll(value); + else if (key == "min_ss") b.min_ss = std::stoll(value); + else if (key == "max_ss") b.max_ss = std::stoll(value); + else if (key == "panel") b.panel = value; + else + return false; + if (key.ends_with("_x") || key.ends_with("_y")) + b.in_lab = true; + return true; +} + } // namespace xfel::Layout Read(const std::string &filename, std::vector &ignored) { @@ -128,7 +178,9 @@ xfel::Layout Read(const std::string &filename, std::vector &ignored // 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; + Masks default_masks; + std::vector panel_masks; // each panel's, beside geom.panels + std::map panel_index, bad_index; std::map dims; std::string line; @@ -147,7 +199,8 @@ xfel::Layout Read(const std::string &filename, std::vector &ignored if (slash == std::string::npos) { if (key == "data") geom.data = value; else if (key.starts_with("dim")) dims[key] = value; - else if (!SetPanel(filename, defaults, key, value) && !Bookkeeping(key)) + else if (!SetPanel(filename, defaults, key, value) && !SetMask(default_masks, key, value) + && !Bookkeeping(key)) ignored.push_back(key); continue; } @@ -156,7 +209,12 @@ xfel::Layout Read(const std::string &filename, std::vector &ignored const std::string prop = key.substr(slash + 1); if (name.starts_with("bad")) { - ignored.push_back(key); + if (!bad_index.contains(name)) { + bad_index[name] = geom.bad.size(); + geom.bad.emplace_back(); + } + if (!SetBad(geom.bad[bad_index[name]], prop, value)) + ignored.push_back(key); continue; } @@ -173,8 +231,11 @@ xfel::Layout Read(const std::string &filename, std::vector &ignored panel_index[name] = geom.panels.size(); geom.panels.push_back(defaults); geom.panels.back().name = name; + panel_masks.push_back(default_masks); } - if (!SetPanel(filename, geom.panels[panel_index[name]], prop, value) && !Bookkeeping(prop)) + const size_t n = panel_index[name]; + if (!SetPanel(filename, geom.panels[n], prop, value) && !SetMask(panel_masks[n], prop, value) + && !Bookkeeping(prop)) ignored.push_back(key); } @@ -195,6 +256,33 @@ xfel::Layout Read(const std::string &filename, std::vector &ignored throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, filename + ": panel " + p.name + " has no pixel size (res)"); } + + // Each mask once for the whole frame where every panel has it alike (the usual case: stated once, + // before the panels), otherwise once for each panel that has it, on that panel's pixels only. + for (size_t n = 0; n < default_masks.size(); n++) { + const auto &first = panel_masks[0][n]; + const bool alike = std::all_of(panel_masks.begin(), panel_masks.end(), + [&](const Masks &m) { return m[n] == first; }); + if (alike) { + if (!first.data.empty()) + geom.masks.push_back(first); + continue; + } + for (size_t i = 0; i < geom.panels.size(); i++) { + if (panel_masks[i][n].data.empty()) + continue; + geom.masks.push_back(panel_masks[i][n]); + geom.masks.back().panel = geom.panels[i].name; + } + } + // A mask file is named as it is, or relative to the working directory, as CrystFEL reads it; a + // relative one not found there is looked for beside the geometry file. + for (auto &m: geom.masks) { + const std::filesystem::path file(m.file); + const auto beside = std::filesystem::path(filename).parent_path() / file; + if (!m.file.empty() && file.is_relative() && !std::filesystem::exists(file) && std::filesystem::exists(beside)) + m.file = beside.string(); + } return geom; } diff --git a/reader/CrystFELGeometry.h b/reader/CrystFELGeometry.h index 18f1557ac..9e951a286 100644 --- a/reader/CrystFELGeometry.h +++ b/reader/CrystFELGeometry.h @@ -9,14 +9,17 @@ #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. +// A CrystFEL detector geometry file (.geom): 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), optionally the camera length - and which pixels to leave out: the mask arrays (mask, mask_file, +// mask_good, mask_bad, or maskN_data, maskN_file, maskN_goodbits, maskN_badbits for up to eight, per +// panel where a panel states its own), the bad regions (bad*, in the laboratory plane or as a pixel +// range on a panel), mask_edge_pixels and no_index. Photon energy, gain, saturation and value flags are +// 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 { // The layout the file describes; keywords read past are added to ignored, for the log. diff --git a/reader/DetectorLayout.h b/reader/DetectorLayout.h index 7f594eb4c..90cd56c2b 100644 --- a/reader/DetectorLayout.h +++ b/reader/DetectorLayout.h @@ -53,10 +53,13 @@ struct BadRegion { // 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 data; // HDF5 path; a % in it stands for the shot group, as in the data path std::string file; // file holding it, or empty for the image file itself uint32_t good = 0, bad = 0; bool nonzero_good = false; + std::string panel; // empty: every panel; otherwise only this panel's pixels + + bool operator==(const Mask &) const = default; [[nodiscard]] bool Bad(uint32_t v) const { return (bad != 0 && (v & bad) != 0) || (good != 0 && (v & good) != good) || (nonzero_good && v == 0); diff --git a/reader/JFJochXFELReader.cpp b/reader/JFJochXFELReader.cpp index 8fb6ad658..b3bd3eb59 100644 --- a/reader/JFJochXFELReader.cpp +++ b/reader/JFJochXFELReader.cpp @@ -137,31 +137,36 @@ void JFJochXFELReader::ReadFiles(const std::vector &files, const st 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: {}", + Logger("XFELReader").Info("{}: the module geometry, camera length and masks are taken; ignored: {}", geometry_file, list); } + Logger("XFELReader").Info("{}: {} mask array(s), {} bad region(s)", geometry_file, geom_.masks.size(), + geom_.bad.size()); - // What the geometry file no longer states comes from the data files themselves, where they are one of + // What the geometry file does not state 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. + // The file's own masks are taken only where the geometry names none: a geometry that does names the + // mask arrays with the bits that matter for this detector. 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]); + const bool file_masks = geom_.masks.empty(); 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")) + if (file_masks && 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")) + if (file_masks && 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 (file_masks) + 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}; @@ -511,8 +516,9 @@ void JFJochXFELReader::Place() { 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 + if (m.file.empty() && (m.data.find('%') != std::string::npos + || FileAt(0).GetDimension(m.data).size() == 3)) { + frame_masks_.push_back(m); // one per shot or one plane per frame: applied as each is read continue; } if (!m.file.empty()) { @@ -520,9 +526,7 @@ void JFJochXFELReader::Place() { 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; + MarkBad(m, v, masked_by_file); } } @@ -701,6 +705,22 @@ void JFJochXFELReader::Place() { SetStartMessage(dataset_); } +void JFJochXFELReader::MarkBad(const xfel::Mask &m, const std::vector &v, std::vector &bad) const { + int64_t min_fs = 0, max_fs = static_cast(frame_width_) - 1; + int64_t min_ss = 0, max_ss = static_cast(frame_height_) - 1; + for (const auto &p: geom_.panels) + if (!m.panel.empty() && p.name == m.panel) { + min_fs = p.min_fs; max_fs = p.max_fs; + min_ss = p.min_ss; max_ss = p.max_ss; + } + for (int64_t ss = min_ss; ss <= max_ss; ss++) + for (int64_t fs = min_fs; fs <= max_fs; fs++) { + const size_t i = static_cast(ss) * frame_width_ + static_cast(fs); + if (m.Bad(v[i])) + bad[i] = 1; + } +} + void JFJochXFELReader::ApplyGain(float adu_per_photon) { for (auto &p: geom_.panels) p.adu_per_photon = adu_per_photon; @@ -823,15 +843,16 @@ CompressedImage JFJochXFELReader::DecodeInto(int64_t image_number, Buffer &buffe // 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) { + if (!frame_masks_.empty()) { masked.assign(frame_width_ * frame_height_, 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; + const auto path = ReplaceEvent(m.data, frame.group); + const auto v = frame.index >= 0 + ? FileAt(frame.file).ReadVector(path, {static_cast(frame.index), 0, 0}, + {1, frame_height_, frame_width_}) + : FileAt(frame.file).ReadVector(path, {0, 0}, {frame_height_, frame_width_}); + MarkBad(m, v, masked); } } diff --git a/reader/JFJochXFELReader.h b/reader/JFJochXFELReader.h index 84fb09923..bf135d031 100644 --- a/reader/JFJochXFELReader.h +++ b/reader/JFJochXFELReader.h @@ -75,6 +75,8 @@ class JFJochXFELReader : public JFJochReader { 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_; + // Marks in bad every stored pixel the mask m, read as v, leaves out (on its panel, where it has one). + void MarkBad(const xfel::Mask &m, const std::vector &v, std::vector &bad) const; // 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 - diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 32629a6ee..b2e641876 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -162,7 +162,7 @@ void print_usage() { std::cout << " Detector geometry" << 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 << " --crystfel-geometry Read the HDF5/CXI files with their panels placed by this CrystFEL geometry file: the module geometry (panel ranges, fast/slow directions, corners, pixel size, coffset, the frames' data path), if stated the camera length, and the masks (mask arrays, bad regions, mask_edge_pixels, no_index) are taken from it; photon energy, camera length and - where the geometry names no mask array - masks 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; diff --git a/tests/CrystFELGeometryTest.cpp b/tests/CrystFELGeometryTest.cpp index adfc2a419..b229119b5 100644 --- a/tests/CrystFELGeometryTest.cpp +++ b/tests/CrystFELGeometryTest.cpp @@ -27,11 +27,18 @@ TEST_CASE("CrystFELGeometry_OldSyntax", "[CrystFELGeometry][portable]") { CHECK(geom.panels[0].clen->path.empty()); CHECK(geom.panels[0].clen->value == Catch::Approx(0.075)); CHECK(geom.data == "/entry/data/data"); - // Only the module geometry is read: gain, photon energy, value flags and bad regions are ignored. + // Gain, photon energy and value flags are ignored; the bad strips, with no panel, are on the whole frame. CHECK(!geom.panels[0].adu_per_photon.has_value()); CHECK(!geom.photon_energy.has_value()); - CHECK(geom.bad.empty()); CHECK(std::find(ignored.begin(), ignored.end(), "photon_energy") != ignored.end()); + CHECK(std::find(ignored.begin(), ignored.end(), "flag_lessthan") != ignored.end()); + REQUIRE(geom.bad.size() == 4); + 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); + CHECK(geom.bad[3].min_ss == 1614); } TEST_CASE("CrystFELGeometry_Panels", "[CrystFELGeometry][portable]") { @@ -95,11 +102,25 @@ badbox/max_y = 5 CHECK(q2.res == Catch::Approx(20000)); CHECK(geom.data == "/%/data"); - 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"}) + // The old mask keywords are mask 0, the same for both panels: one mask for the whole frame. + REQUIRE(geom.masks.size() == 1); + CHECK(geom.masks[0].data == "/data/data"); + CHECK(geom.masks[0].file == "mask.h5"); + CHECK(geom.masks[0].good == 0); + CHECK(geom.masks[0].bad == 1); + CHECK(geom.masks[0].panel.empty()); + REQUIRE(geom.bad.size() == 2); + CHECK(geom.bad[0].panel == "q1"); + CHECK(!geom.bad[0].in_lab); + CHECK(geom.bad[0].max_fs == 4); + CHECK(geom.bad[0].max_ss == 1023); + CHECK(geom.bad[1].in_lab); + CHECK(geom.bad[1].min_x == Catch::Approx(-10)); + CHECK(geom.bad[1].max_y == Catch::Approx(5)); + for (const std::string key: {"adu_per_eV", "photon_energy", "q2/adu_per_photon"}) CHECK(std::find(ignored.begin(), ignored.end(), key) != ignored.end()); + for (const std::string key: {"mask", "mask_file", "mask_good", "mask_bad", "badq1/min_fs", "badbox/max_y"}) + CHECK(std::find(ignored.begin(), ignored.end(), key) == ignored.end()); std::filesystem::remove(path); } @@ -122,9 +143,13 @@ 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()); + REQUIRE(geom.masks.size() == 1); + CHECK(geom.masks[0].data == "/entry/mask"); + CHECK(geom.masks[0].file.empty()); + CHECK(geom.masks[0].bad == 6); + CHECK(geom.panels[0].mask_edge_pixels == 2); + CHECK(!geom.panels[0].no_index); 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", @@ -141,3 +166,65 @@ some_new_keyword = 1 } std::filesystem::remove(path); } + +TEST_CASE("CrystFELGeometry_Masks", "[CrystFELGeometry][portable]") { + // Two masks: mask 0 for every panel, with its bad bits overridden on panel a; mask 1 in a file named + // relative to the geometry file's directory. A bad region on panel b, and panel c left out. + const auto dir = std::filesystem::temp_directory_path() / "jfjoch_crystfel_geometry_masks"; + std::filesystem::create_directories(dir); + std::ofstream(dir / "pixels.h5") << ""; + const auto path = (dir / "masks.geom").string(); + std::ofstream(path) << R"( +res = 10000 +mask0_data = /entry/mask +mask0_goodbits = 0x1 +mask1_data = /%/pixel_mask +mask1_file = pixels.h5 +mask1_badbits = 0xfeff +a/min_fs = 0 +a/max_fs = 9 +a/min_ss = 0 +a/max_ss = 9 +a/mask0_badbits = 4 +b/min_fs = 10 +b/max_fs = 19 +b/min_ss = 0 +b/max_ss = 9 +c/min_fs = 20 +c/max_fs = 29 +c/min_ss = 0 +c/max_ss = 9 +c/no_index = true +bad_b/min_fs = 12 +bad_b/max_fs = 13 +bad_b/min_ss = 0 +bad_b/max_ss = 9 +bad_b/panel = b +)"; + std::vector ignored; + const auto geom = crystfel::Read(path, ignored); + CHECK(ignored.empty()); + REQUIRE(geom.panels.size() == 3); + CHECK(geom.panels[2].no_index); + CHECK(!geom.panels[1].no_index); + + // Mask 0 differs between panels: one per panel. Mask 1 is the same everywhere: one for the frame. + REQUIRE(geom.masks.size() == 4); + for (size_t i = 0; i < 3; i++) { + CHECK(geom.masks[i].data == "/entry/mask"); + CHECK(geom.masks[i].good == 1); + CHECK(geom.masks[i].panel == geom.panels[i].name); + } + CHECK(geom.masks[0].bad == 4); + CHECK(geom.masks[1].bad == 0); + CHECK(geom.masks[3].data == "/%/pixel_mask"); + CHECK(geom.masks[3].bad == 0xfeff); + CHECK(geom.masks[3].panel.empty()); + CHECK(std::filesystem::path(geom.masks[3].file) == dir / "pixels.h5"); + + REQUIRE(geom.bad.size() == 1); + CHECK(geom.bad[0].panel == "b"); + CHECK(geom.bad[0].min_fs == 12); + CHECK(geom.bad[0].max_fs == 13); + std::filesystem::remove_all(dir); +} diff --git a/tests/XFELReaderTest.cpp b/tests/XFELReaderTest.cpp index bd2661ede..2d73446b7 100644 --- a/tests/XFELReaderTest.cpp +++ b/tests/XFELReaderTest.cpp @@ -2,6 +2,8 @@ // SPDX-License-Identifier: GPL-3.0-only #include +#include +#include #include #include #include @@ -245,6 +247,124 @@ TEST_CASE("JFJochXFELReader_CXI_Int16", "[portable]") { std::filesystem::remove(name); } +TEST_CASE("JFJochXFELReader_GeometryMasks", "[portable]") { + // A Cheetah CXI file of two 3-column x 4-row panels side by side in a 4 x 6 frame, read through a + // CrystFEL geometry file that names the masks: the per-frame mask with its own bad bits, a mask in a + // file beside the geometry file, a box in the laboratory plane and a pixel range on panel p1. The + // file's own mask_shared is not taken (the geometry names the masks); --mask-bad still adds its own. + const auto dir = std::filesystem::path("xfel_reader_geom_masks"); + std::filesystem::create_directories(dir); + const std::string name = (dir / "frames.cxi").string(); + { + 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(2 * 24); + for (size_t i = 0; i < frames.size(); i++) + frames[i] = static_cast(i % 24) + 100; // every pixel its own value, in both frames + WriteFrames(f, "entry_1/data_1/data", H5T_NATIVE_INT32, frames.data(), {2, 4, 6}); + std::vector mask(2 * 24, 0); + mask[1 * 6 + 1] = 0x20; // frame 0 only, (fs 1, ss 1): a bit Cheetah's own set does not include + WriteFrames(f, "entry_1/data_1/mask", H5T_NATIVE_UINT16, mask.data(), {2, 4, 6}); + const std::string det = "entry_1/instrument_1/detector_1/"; + std::vector shared(24, 0); + shared[3 * 6 + 5] = 8; // (fs 5, ss 3) dead to Cheetah + WriteFrames(f, det + "mask_shared", H5T_NATIVE_UINT16, shared.data(), {4, 6}); + WriteVector(f, "entry_1/instrument_1/source_1/energy", std::vector(2, 9000.0 * 1.602176634e-19)); + H5Fclose(f); + } + { + // Good where bit 0 is set; bit 1 is what --mask-bad below looks at. + const hid_t f = H5Fcreate((dir / "pixels.h5").string().c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + std::vector pixels(24, 1); + pixels[0 * 6 + 4] = 0; // (fs 4, ss 0) + pixels[2 * 6 + 3] = 3; // (fs 3, ss 2) + WriteFrames(f, "mask", H5T_NATIVE_UINT16, pixels.data(), {4, 6}); + H5Fclose(f); + } + const std::string geom = (dir / "frames.geom").string(); + std::ofstream(geom) << R"( +clen = 0.1 +res = 10000 +dim0 = % +dim1 = ss +dim2 = fs +data = /entry_1/data_1/data +mask = /entry_1/data_1/mask +mask_good = 0x0 +mask_bad = 0xffff +mask1_data = /mask +mask1_file = pixels.h5 +mask1_goodbits = 0x1 +p0/min_fs = 0 +p0/max_fs = 2 +p0/min_ss = 0 +p0/max_ss = 3 +p0/fs = x +p0/ss = y +p0/corner_x = -3 +p0/corner_y = -2 +p1/min_fs = 3 +p1/max_fs = 5 +p1/min_ss = 0 +p1/max_ss = 3 +p1/fs = x +p1/ss = y +p1/corner_x = 0 +p1/corner_y = -2 +badbox/min_x = -3 +badbox/max_x = -2 +badbox/min_y = -2 +badbox/max_y = -1 +bad_p1/min_fs = 5 +bad_p1/max_fs = 5 +bad_p1/min_ss = 0 +bad_p1/max_ss = 0 +bad_p1/panel = p1 +)"; + + JFJochXFELReader reader; + reader.ExtraMasks({{"/mask", (dir / "pixels.h5").string(), 0, 2}}); + REQUIRE_NOTHROW(reader.ReadFiles({name}, geom)); + REQUIRE(reader.GetNumberOfImages() == 2); + + // Where each stored pixel lands, found by its value in frame 1 (which the per-frame mask leaves alone). + auto raw = reader.GetRawImage(1); + REQUIRE(raw); + const auto *frame1 = reinterpret_cast(raw->image_buffer.data()); + const size_t n = raw->image.GetWidth() * raw->image.GetHeight(); + auto where = [&](size_t fs, size_t ss) { + const auto it = std::find(frame1, frame1 + n, static_cast(ss * 6 + fs) + 100); + REQUIRE(it != frame1 + n); + return static_cast(it - frame1); + }; + const auto &m = reader.GetDataset()->pixel_mask->GetMask(); + CHECK(m[where(2, 2)] == 0); + CHECK(m[where(4, 0)] != 0); // mask 1, beside the geometry file: lacks the good bit + CHECK(m[where(3, 2)] != 0); // --mask-bad: bit 1 + CHECK(m[where(5, 3)] == 0); // the file's mask_shared: not taken + CHECK(m[where(5, 1)] == 0); + CHECK(m[where(0, 0)] != 0); // the lab box: p0's first pixel (centre at x -2.5, y -1.5) and no other + CHECK(m[where(1, 0)] == 0); + CHECK(m[where(0, 1)] == 0); + CHECK(m[where(5, 0)] != 0); // bad_p1 + size_t usable = 0; + for (size_t ss = 0; ss < 4; ss++) + for (size_t fs = 0; fs < 6; fs++) + usable += m[where(fs, ss)] == 0; + CHECK(usable == 24 - 4); + + // The per-frame mask: (1, 1) is no measurement in frame 0 only. + auto raw0 = reader.GetRawImage(0); + REQUIRE(raw0); + const auto *frame0 = reinterpret_cast(raw0->image_buffer.data()); + CHECK(frame0[where(1, 1)] == INT32_MIN); + CHECK(frame0[where(2, 1)] == 1 * 6 + 2 + 100); + reader.Close(); + std::filesystem::remove_all(dir); +} + namespace { // An LCLS-I XTC stream of the shape JFJochXTCReader reads: a Configure datagram, then one event (an