Files
leonarski_fandClaude Opus 5.5 ebce906025 Serial XFEL stills inputs: readers from the files themselves, module geometry and its refinement
Combines the stills reader and module work (stills-readers, stills-modgeom, stills-euxfel) onto
rc178 with the owner's decisions for this release:

- Inputs are files, never a directory; several files of one format are one dataset (the multi-input
  CLI of the previous commit), so every reader takes an explicit list: JFJochXFELReader::ReadFiles,
  ReadCXI, ReadPanelTable and JFJochXTCReader::ReadFiles take std::vector<std::string>; ListFiles and
  xtc::CollectStreams are gone, CanRead* accept a regular file only.
- One HDF5 frame reader (JFJochCrystFELReader -> JFJochXFELReader) parameterised by what the file holds:
  per-shot groups or a 3D stack, a /metadata panel table (posx/posy/angle, pixel size, nominal camera
  length, pixel mask: as SACLA writes), a Cheetah CXI single panel (square frame; a tiled slab is
  refused with a message). At most 8 files open; every file checked against the first. Signed 32-bit
  decode, per-frame 3D masks, panels at their own coffset (from the EuXFEL work).
- Frames stored one chunk each are read the way the NXmx reader reads its images: only the chunk's
  place is asked of HDF5 under its lock (H5Dget_chunk_info_by_coord), the bytes are read through
  HDF5ImageSource::RawFile (positional, no lock) and decoded outside it - deflate(+shuffle) by zlib,
  bitshuffle/LZ4 (and zstd, plain LZ4) by CompressedImage's decoders. Bitshuffle/LZ4 used to go through
  HDF5's filter under the global lock, one frame at a time.
- LCLS-I XTC (Rayonix): wavelength from the beam record, pixel size from the Rayonix model and
  binning, only --detector-distance needed; damaged/compressed frames skipped, a truncated last event
  dropped; std::ifstream, no POSIX calls: builds with MSVC. No DIALS .expt support at all.
- CrystFEL .geom: module geometry only (panel ranges, fs/ss, corners, pixel size, coffset, data path)
  and optionally the camera length; everything else is ignored with a note in the log.
- --mask-good / --mask-bad <file>:<dataset>: a mask array in another HDF5 file, in the frames' stored
  layout (what a .geom's mask_file used to give).
- Gain measured on the background for panel-table and CXI input that states none (--gain overrides).
- Detector modules from explicit geometry only (DetectorModules: the reader's panels, an NXmx file with
  one NXdetector_module per module, or --module-geometry); never from the image. Spots and predictions
  are mapped through the modules, the image is not re-assembled.
- Module refinement in the stills bundle adjustment (Ceres; per-module shift and in-plane turn, gauge
  on mean shift/turn/scale and on the two tilt-like patterns when the tilt is refined; adopted only
  where it predicts held-out frames better). Applied for a foreign multi-panel detector (panel table,
  .geom panels); for a JUNGFRAU measured and reported only. Written to <prefix>_modules.h5 (NXmx
  NXdetector_module groups), read back with --module-geometry.
- A panel-table file's camera length is nominal: with a reference cell it is refined from the cell.
- Broker and writer output unchanged: masters keep the single NXdetector_module group, the start
  message carries no module list.

Tests: [CrystFELGeometry], JFJochXFELReader_*, JFJochXTCReader, [JFJochMultiReader], [CellPopulation],
SaturationCeiling_*, [DetectorModules], RefineGlobalGeometry_ModulePlacement, [reindex], Sublattice.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
2026-10-20 23:30:48 +02:00

134 lines
7.4 KiB
C++

// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
#include <mutex>
#include "../common/JFJochMessages.h"
#include "../common/DiffractionExperiment.h"
#include "../common/AzimuthalIntegrationMapping.h"
#include "../common/PixelMask.h"
#include "../common/AzimuthalIntegrationProfile.h"
#include "bragg_prediction/BraggPrediction.h"
#include "bragg_integration/BraggIntegrationEngine.h"
#include "bragg_integration/BraggIntegrationEngineCPU.h"
#include "spot_finding/ImageSpotFinder.h"
#include "spot_finding/AdaptiveSpotFinderCPU.h"
#include "indexing/IndexerThreadPool.h"
#include "azint/AzIntEngine.h"
#include "roi/ROIIntegration.h"
#include "IndexAndRefine.h"
#include "image_preprocessing/ImagePreprocessor.h"
#include "image_preprocessing/ImagePreprocessorBuffer.h"
#include "image_preprocessing/ImagePreprocessorCPU.h"
class CudaStream;
class AdaptiveSpotFinderGPU;
// MXAnalysisWithoutFPGA is not thread safe - it has to owned by a single thread
class MXAnalysisWithoutFPGA {
const DiffractionExperiment &experiment;
const AzimuthalIntegrationMapping &integration;
std::vector<uint8_t> decompression_buffer;
std::unique_ptr<ImagePreprocessor> preprocessor;
// The preprocessor, where it is the CPU one.
ImagePreprocessorCPU *preprocessor_cpu = nullptr;
size_t npixels;
size_t xpixels;
// Built on first use: the fused adaptive finder produces the azimuthal profile as a by-product,
// so on the rugnux path this engine is constructed and then never run.
std::unique_ptr<AzIntEngine> azint;
AzIntEngine &AzInt();
std::unique_ptr<ROIIntegration> roi;
// Built on first use. Which finder an image takes arrives with its SpotFindingSettings, and
// with adaptive detection on - the default everywhere but the broker - this one is never asked
// for; on the GPU it is ~14 MB and 15 device allocations per worker.
std::unique_ptr<ImageSpotFinder> spotFinder;
ImageSpotFinder &FixedThresholdFinder();
// Self-calibrating finder, used when spot settings request adaptive detection. Kept alongside the
// default finder because the choice arrives with the per-image settings, not at construction. It is
// an AdaptiveSpotFinderCPU by default; on the GPU path, when the fused engine is enabled (rugnux
// offline only), it is instead an AdaptiveSpotFinderGPU that also computes the azimuthal profile,
// aliased through fused_adaptive so Analyze() can take that profile and skip the separate azint pass.
std::unique_ptr<ImageSpotFinder> adaptiveSpotFinder;
AdaptiveSpotFinderGPU *fused_adaptive = nullptr;
// The CPU finder, where it gives the profile too (the azimuthal integration being on the CPU).
AdaptiveSpotFinderCPU *fused_adaptive_cpu = nullptr;
const bool enable_fused_adaptive_gpu;
IndexAndRefine &indexer;
std::unique_ptr<BraggPrediction> prediction;
std::unique_ptr<BraggIntegrationEngine> 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<Reflection> IntegrateOnImage(const std::vector<Reflection> &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),
// on the images for which count_crowding is on.
std::unique_ptr<BraggIntegrationEngineCPU> crowding_engine;
bool count_crowding = false;
std::unique_ptr<ImagePreprocessorBuffer> preprocessor_buffer;
const PixelMask &mask;
// Decompress the image into decompression_buffer (or read it straight from the message, when it is
// not compressed) and return where it landed.
const uint8_t *Decompress(const CompressedImage &image);
// The CPU preprocessing. A bitshuffled image is decoded a block (~32 KB) at a time and each block
// is preprocessed while it is in cache, then - with ring_pass - put through the adaptive finder's
// plain ring pass as well, so the image is never written out whole before it is preprocessed and
// the preprocessed pixels are read back from cache, not from memory. The blocks come in pixel
// order, so every sum is taken in the same order as by separate passes.
ImageStatistics PreprocessCPU(const CompressedImage &image, bool ring_pass);
// Pixels outside the resolution limits, bit-packed. Built by the integration mapping, which is
// shared by every worker's engine and hands out the same mask to all of them.
std::shared_ptr<const std::vector<uint32_t>> mask_resolution;
// The limits mask_resolution was built for. Kept as the OPTIONAL the caller passed, so an unset
// high-resolution limit compares equal to itself and the mask is not rebuilt on every image.
std::optional<float> mask_high_res;
std::optional<float> mask_low_res;
void UpdateMaskResolution(const SpotFindingSettings& settings);
#ifdef JFJOCH_USE_CUDA
std::shared_ptr<CudaStream> stream; // kept so RebuildROI() can recreate the GPU ROI engine
#endif
public:
// enable_fused_adaptive_gpu turns on the fused GPU azint+adaptive spot finder (only takes effect on
// the GPU path with adaptive detection). The rugnux offline path and the interactive viewer enable
// it by default, as does the online receiver. It only changes performance - the fused engine
// reproduces the CPU finder's spots. Note it also decides whether the preprocessed image is copied
// back to the host each frame: that copy exists only for a CPU engine to read, and with the flag on
// no CPU engine is built, so the copy is skipped.
MXAnalysisWithoutFPGA(const DiffractionExperiment &experiment, const AzimuthalIntegrationMapping &integration,
const PixelMask &mask, IndexAndRefine &indexer, bool enable_fused_adaptive_gpu = false);
void Analyze(DataMessage &output, AzimuthalIntegrationProfile &profile, const SpotFindingSettings &spot_finding_settings);
// Surgical ROI-only paths used when a full re-analysis is not wanted: rebuild the
// ROI engine after the ROI set changes, recompute ROIs after preprocessing a new
// image (reanalyze off), or just rerun ROIs on the current preprocessed image (an
// interactive ROI move). A full Analyze() already computes ROIs, so needs nothing.
void RebuildROI();
void AnalyzeROIOnly(DataMessage &output);
void RunROIOnly(DataMessage &output);
// What this worker's Bragg integrator counted (BraggIntegrationCounts). Each worker builds its own
// analysis, so a caller that wants the run's totals sums this over the workers it started.
[[nodiscard]] BraggIntegrationCounts BraggCounts() const;
// Also count, on the images Analyze is given while CountCrowding(true), how crowded the background
// rings of the integrated reflections would be at `settings` rather than the experiment's own
// (BraggIntegrationEngineCPU::CountNeighbourCrowding), read back with CrowdingCounts.
void CountCrowdingAt(const BraggIntegrationSettings &settings);
void CountCrowding(bool on) { count_crowding = on; }
[[nodiscard]] BraggIntegrationCounts CrowdingCounts() const;
};