From 01e4f81e4430e5b34d6e3ca19daf5b88f82797c9 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 8 Jul 2026 14:09:15 +0200 Subject: [PATCH] Consolidate offline tools into rugnux; merge on by default Fold jfjoch_azint and rugnux_scale into the single rugnux CLI, and flip scaling/merging from opt-in to on-by-default. - rugnux --azint-only replaces jfjoch_azint (reuses ProcessMode:: AzimuthalIntegration; adds the correction toggles and geometry overrides the old tool carried). - rugnux --scale replaces rugnux_scale (re-scale/merge stored reflections, same simpler scale/merge as before: no space-group search or rot3d defaults). Its workflow is folded into the CLI verbatim. - Merging is now on by default for rotation and stills; --no-merge disables it, replacing -M/--scale-merge. Also fix the batch ReadReflections reader so it works on rugnux's own NXmxIntegrated _process.h5: read per-image reflections/mosaicity/lattice from the master (global index) first, falling back to the source-data locator for legacy/VDS datasets. Without this, both --scale and the former rugnux_scale read zero reflections from an integrated snapshot. RugnuxCommandLine now emits `rugnux --azint-only` and `--no-merge`; docs and tests updated to match. Co-Authored-By: Claude Fable 5 --- CLAUDE.md | 6 +- CMakeLists.txt | 2 +- docs/CHANGELOG.md | 1 + docs/RUGNUX.md | 36 ++-- docs/TOOLS.md | 49 +---- reader/HDF5MetadataSource.cpp | 19 +- rugnux/Rugnux.h | 6 +- rugnux/RugnuxCommandLine.cpp | 16 +- rugnux/RugnuxCommandLine.h | 2 +- tests/RugnuxTest.cpp | 7 +- tools/CMakeLists.txt | 12 +- tools/jfjoch_azint.cpp | 296 -------------------------- tools/rugnux_cli.cpp | 327 +++++++++++++++++++++++++--- tools/rugnux_scale.cpp | 390 ---------------------------------- 14 files changed, 359 insertions(+), 810 deletions(-) delete mode 100644 tools/jfjoch_azint.cpp delete mode 100644 tools/rugnux_scale.cpp diff --git a/CLAUDE.md b/CLAUDE.md index a8754e28..c331f080 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -123,8 +123,10 @@ generated API model and internal types. - `jfjoch_broker` — online, real-time (FPGA + GPU). - `jfjoch_viewer` — interactive Qt desktop (`viewer/`), results not persisted. - `rugnux` (`tools/rugnux_cli.cpp`, built on the `Rugnux` library in `rugnux/`) — offline batch over - a stored HDF5; writes `_process.h5` and `.mtz`/`.cif`/`.hkl`. `rugnux_scale` re-scales/merges - already-integrated data. (rugnux = the data-processing half of the system; see `docs/NAMING.md`.) + a stored HDF5; writes `_process.h5` and `.mtz`/`.cif`/`.hkl`. Merging is on by default (`--no-merge` + to disable); `--azint-only` runs only azimuthal integration and `--scale` re-scales/merges the + already-integrated reflections in a `_process.h5`. (rugnux = the data-processing half of the system; + see `docs/NAMING.md`.) **`image_analysis/` pipeline** (subdirs): `spot_finding`, `indexing` (`ffbidx`/`fft` GPU, `fftw` CPU), `lattice_search`, `geom_refinement`, `pixel_refinement`, `bragg_prediction`, diff --git a/CMakeLists.txt b/CMakeLists.txt index cd78188b..b4fc549a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -279,7 +279,7 @@ IF (JFJOCH_VIEWER_ONLY) ADD_SUBDIRECTORY(reader) ADD_SUBDIRECTORY(rugnux) ADD_SUBDIRECTORY(viewer) - ADD_SUBDIRECTORY(tools) # builds only the portable analysis tools (rugnux/rugnux_scale/azint/extract_hkl) + ADD_SUBDIRECTORY(tools) # builds only the portable analysis tools (rugnux/extract_hkl) ELSE() ADD_SUBDIRECTORY(jungfrau) ADD_SUBDIRECTORY(compression) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index b85f78d2..ecb7f121 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.156 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. +* rugnux: Consolidate the offline analysis tools into one `rugnux` binary - the former `jfjoch_azint` is now `rugnux --azint-only` and `rugnux_scale` is now `rugnux --scale`. Scaling and merging are on by default (`--no-merge` disables them), replacing the previous opt-in `-M, --scale-merge`. * jfjoch_process: Major rotation (rot3d) data processing overhaul - robust profile-fit integration, Cauchy-loss scaling with optional absorption surface, de-novo indexing and space-group/centering determination fixes, and merging statistics + ISa in the mmCIF output. * jfjoch_process: Bragg integration now runs on the GPU in the offline/non-FPGA workflow (one box-sum + profile-fit engine, GPU when available, CPU otherwise); the FPGA workflow integrates on the CPU directly from the assembled image. The previous standalone integrators are removed. * jfjoch_process: Deterministic Bragg prediction - when more reflections are predicted than fit the output, they are ranked by distance to the Ewald sphere before truncation, so repeated runs produce identical reflections. diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 7dddc771..494702fa 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -42,9 +42,10 @@ the first pass. - `_process.h5` — NXmx-compliant HDF5 with derived metadata (spots, indexing, integration, azimuthal integration, per-image statistics). See - [HDF5 / NeXus data format](HDF5.md) for the layout. -- When merging (`-M`, or whenever a `--reference-mtz` is supplied), the merged reflections are - written as `.cif` (mmCIF — the default), or `.mtz` / `.hkl` depending on + [HDF5 / NeXus data format](HDF5.md) for the layout. Written by default only when **not** merging + (i.e. under `--no-merge`); add `--write-process-h5` to also write it when merging. +- Merging is **on by default** (`--no-merge` disables it). The merged reflections are written as + `.cif` (mmCIF — the default), or `.mtz` / `.hkl` depending on `--scaling-output`. Both the mmCIF and the MTZ carry the **refined unit cell** (from rotation indexing) and the **space group determined from systematic absences** (constrained to the indexed lattice symmetry). No-reference scaling additionally emits per-iteration `_iterN_scale.dat`. @@ -52,12 +53,13 @@ the first pass. Merged statistics (⟨I/σ⟩, CC1/2, completeness, …), the error model and timing are printed to the console. -## Re-scaling and re-merging (`rugnux_scale`) +## Re-scaling and re-merging (`rugnux --scale`) -The companion tool `rugnux_scale` re-scales and merges the *already-integrated* reflections stored -in one or more `_process.h5` files, without re-running spot finding or integration. Use it to -re-merge quickly with a different space group, resolution limit, anomalous setting or reference MTZ, -or to combine several processed runs into one set of merged intensities. +The `--scale` mode re-scales and merges the *already-integrated* reflections stored in a +`_process.h5` file, without re-running spot finding or integration. Use it to re-merge quickly with a +different space group, resolution limit, anomalous setting or reference MTZ. It reuses the same +`-o/-N/-s/-e/-S/-A/-B/-z/--scaling-*` options as the full run, and (unlike the full pipeline) does +not run a space-group search, so pass `-S` for the correct symmetry. ## Quick start @@ -73,9 +75,10 @@ rugnux rotation_master.h5 \ Because the dataset carries a rotation goniometer axis, it is processed as **rotation data by default**: two-pass rotation indexing (index the sweep once, then process every frame against that -lattice) with the **`rot3d`** partiality model (rotation partials combined into 3D fulls). Rotation -data is **scaled and merged automatically** (as if `-M` were given); the unit cell is taken from the -rotation indexer and the space group is determined from systematic absences, and both are written +lattice) with the **`rot3d`** partiality model (rotation partials combined into 3D fulls). Scaling +and merging run **by default** (for both rotation and stills; `--no-merge` turns them off); the unit +cell is taken from the rotation indexer and the space group is determined from systematic absences, +and both are written into the merged `.cif`. Run **fully de novo** (no `-C`/`-S`) for the best result — supplying a cell or space group up front @@ -95,7 +98,7 @@ rugnux serial_master.h5 \ -o lyso_serial -N 32 \ -X ffbidx -C 79,79,38,90,90,90 -S 96 \ --spot-sigma 4 \ - -M -z reference.mtz \ + -z reference.mtz \ --scaling-high-resolution 1.8 ``` @@ -117,6 +120,13 @@ General: | `-t, --stride ` | Process every *n*-th image (default: 1) | | `-v, --verbose` | Verbose output | +Modes (default: full analysis — spot finding, indexing, integration and merging): + +| Option | Description | +| --- | --- | +| `--azint-only` | Only run azimuthal integration (no spot finding/indexing); writes `_process.h5` | +| `--scale` | Only re-scale/merge the already-integrated reflections in the input `_process.h5` (no re-integration) | + Spot finding: | Option | Description | @@ -163,7 +173,7 @@ Scaling and merging: | Option | Description | | --- | --- | -| `-M, --scale-merge` | Scale and merge (automatic for rotation data; only needed to force scaling on stills) | +| `--no-merge` | Skip scaling and merging (on by default); write only the per-image `_process.h5` | | `-A, --anomalous` | Anomalous mode (keep Friedel pairs separate) | | `-B, --refine-bfactor` | Refine a per-image B-factor | | `--scale-fulls` / `--no-scale-fulls` | rot3d: refit a per-frame scale on the combined fulls (XDS order, Unity model); on by default for rotation data, off for stills | diff --git a/docs/TOOLS.md b/docs/TOOLS.md index 345ed34e..d9f74dcf 100644 --- a/docs/TOOLS.md +++ b/docs/TOOLS.md @@ -9,12 +9,10 @@ number of command-line tools. Each prints its own usage when run with `-h` or wi ### rugnux Offline CLI tool that runs the full crystallographic analysis pipeline (spot finding, indexing, integration, scaling/merging) on a stored HDF5 dataset, producing a `_process.h5` file and, when -merging, reflection files. See [rugnux](RUGNUX.md). - -### rugnux_scale -Re-scales and merges the already-integrated reflections from one or more `_process.h5` files -(no re-integration). Useful to re-merge with a different space group, partiality, resolution limit -or reference MTZ, or to combine several runs. See [rugnux](RUGNUX.md). +merging, reflection files. Merging is on by default (`--no-merge` disables it). Two extra modes +narrow the work: `--azint-only` runs only azimuthal integration (no spot finding/indexing), and +`--scale` re-scales/merges the already-integrated reflections in a `_process.h5` without +re-integrating. See [rugnux](RUGNUX.md). ### jfjoch_extract_hkl Extracts reflections (HKL list) from a Jungfraujoch master file; can sum the same HKL across @@ -66,41 +64,4 @@ writing. Tests single-threaded HDF5 writer performance. ### jfjoch_simplon_test -Minimal test client for a DECTRIS SIMPLON detector API. - -### jfjoch_azint - -Runs CPU azimuthal integration on an existing Jungfraujoch HDF5 file using `N` threads and writes the -result to `_process.h5`. Example: -``` -./jfjoch_azint -o output -N 8 input.h5 -``` - -The integration settings and geometry default to the values stored in the input file; any of them can be -overridden on the command line: -``` -jfjoch_azint {} - -o, --output-prefix Output file prefix (default: output) - -N, --threads Number of threads (default: 1) - -s, --start-image Start image number (default: 0) - -e, --end-image End image number (default: all) - -t, --stride Image stride (default: 1) - -v, --verbose Verbose output - - Azimuthal integration: - --min-q Minimum Q for integration (1/A) - --max-q Maximum Q for integration (1/A) - --q-spacing Q bin spacing (1/A) - --azimuthal-bins Number of azimuthal bins (default: 1) - --polarization-correction Enable/disable polarization correction - --solid-angle-correction Enable/disable solid angle correction - - Geometry overrides: - --beam-x Beam center X (pixel) - --beam-y Beam center Y (pixel) - --detector-distance Detector distance (mm) - --wavelength Wavelength (A) - --rot1 PONI rotation 1 (rad) - --rot2 PONI rotation 2 (rad) - --polarization Polarization factor -``` \ No newline at end of file +Minimal test client for a DECTRIS SIMPLON detector API. \ No newline at end of file diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 811ab091..f5b8af3a 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -1102,14 +1102,21 @@ std::vector HDF5MetadataSource::ReadReflections(size_t start // Generic (non-image-specific) detector geometry from experiment setup. outcome.geom = cached_geom; - // Per-image reflections and MX metadata live in the same file as the image pixels, - // at the source-local index (the locator keeps the data-file handle cached). - const auto loc = ResolveMeta(static_cast(img)); - HDF5Object *meta_file = loc.file.get(); - const size_t meta_image_id = loc.local_index; + // Per-image reflections and MX metadata are stored in this master at the global index for a + // self-contained integrated _process.h5 snapshot, or co-located with the pixels in the source + // data file at the source-local index for a legacy/VDS dataset. Prefer the master (so an + // integrated snapshot reads without its linked source data present); fall back to the source. + HDF5ReadOnlyFile *meta_file = master_file.get(); + size_t meta_image_id = img; + std::string refl_group = fmt::format("/entry/reflections/image_{:06d}", img); + if (!master_file->Exists(refl_group)) { + const auto loc = ResolveMeta(static_cast(img)); + meta_file = loc.file.get(); + meta_image_id = loc.local_index; + refl_group = fmt::format("/entry/reflections/image_{:06d}", meta_image_id); + } // ── reflections ────────────────────────────────────────────────────── - const std::string refl_group = fmt::format("/entry/reflections/image_{:06d}", meta_image_id); ReadReflectionsFromGroup(*meta_file, refl_group, outcome.reflections); // ── per-image mosaicity ─────────────────────────────────────────────── diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 1ba2d77e..69e9a059 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -20,14 +20,14 @@ class JFJochHDF5Reader; -// Offline reprocessing of a stored Jungfraujoch HDF5 dataset, shared by rugnux, -// jfjoch_azint and (later) the viewer. The full processing workflow lives here, not in the CLIs: +// Offline reprocessing of a stored Jungfraujoch HDF5 dataset, shared by rugnux (including its +// --azint-only mode) and (later) the viewer. The full processing workflow lives here, not in the CLIs: // setup, an optional two-pass rotation-indexing pre-pass, a parallel per-image loop (std::thread), // an optional scaling/merging post-pass, and the _process.h5 output. The detector geometry and all // algorithm settings are configured on the DiffractionExperiment by the caller; ProcessConfig only // carries run control. Cancellable from any thread (e.g. SIGINT or a GUI button) via Cancel(). enum class ProcessMode { - AzimuthalIntegration, // preprocess + azimuthal integration only (jfjoch_azint) + AzimuthalIntegration, // preprocess + azimuthal integration only (rugnux --azint-only) FullAnalysis // spot finding + indexing + refinement + integration (rugnux) }; diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 0c7aa6bd..315bc4cd 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -53,7 +53,9 @@ std::string RugnuxCommandLine(const ProcessConfig &config, const std::string &input_file) { std::vector args; const bool azint = (config.mode == ProcessMode::AzimuthalIntegration); - args.emplace_back(azint ? "jfjoch_azint" : "rugnux"); + args.emplace_back("rugnux"); + if (azint) + args.emplace_back("--azint-only"); auto add = [&](const std::string &flag, const std::string &val) { args.push_back(flag); @@ -72,10 +74,10 @@ std::string RugnuxCommandLine(const ProcessConfig &config, if (azint) { const auto a = experiment.GetAzimuthalIntegrationSettings(); - add("--min-q", num(a.GetLowQ_recipA())); - add("--max-q", num(a.GetHighQ_recipA())); - add("--q-spacing", num(a.GetQSpacing_recipA())); - add("--azimuthal-bins", std::to_string(a.GetAzimuthalBinCount())); + add("--azim-min-q", num(a.GetLowQ_recipA())); + add("--azim-max-q", num(a.GetHighQ_recipA())); + add("--azim-q-spacing", num(a.GetQSpacing_recipA())); + add("--azim-phi-bins", std::to_string(a.GetAzimuthalBinCount())); add("--polarization-correction", a.IsPolarizationCorrection() ? "on" : "off"); add("--solid-angle-correction", a.IsSolidAngleCorrection() ? "on" : "off"); } else { @@ -116,13 +118,15 @@ std::string RugnuxCommandLine(const ProcessConfig &config, args.emplace_back("--force-still"); } + // Merging is on by default; emit --no-merge only when it was turned off. if (config.run_scaling) { - args.emplace_back("-M"); const auto sc = experiment.GetScalingSettings(); if (!sc.GetMergeFriedel()) args.emplace_back("-A"); if (sc.GetRefineB()) args.emplace_back("-B"); + } else { + args.emplace_back("--no-merge"); } } diff --git a/rugnux/RugnuxCommandLine.h b/rugnux/RugnuxCommandLine.h index 197198f6..b8362f70 100644 --- a/rugnux/RugnuxCommandLine.h +++ b/rugnux/RugnuxCommandLine.h @@ -9,7 +9,7 @@ class DiffractionExperiment; -// Reconstruct an equivalent rugnux / jfjoch_azint command line for a configured run, so a +// Reconstruct an equivalent rugnux command line (including --azint-only) for a configured run, so a // job set up in the GUI can be handed off to a cluster. Covers the settings that matter for the // run, not every obscure flag; geometry is taken from the input file, so geometry overrides are // not emitted. diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index f5bf5620..553810fb 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -174,8 +174,9 @@ TEST_CASE("RugnuxCommandLine_AzInt", "[process]") { config.output_prefix = "az"; const std::string cmd = RugnuxCommandLine(config, x, "in.h5"); - CHECK(cmd.rfind("jfjoch_azint", 0) == 0); - CHECK(cmd.find("--azimuthal-bins 4") != std::string::npos); - CHECK(cmd.find("--min-q") != std::string::npos); + CHECK(cmd.rfind("rugnux", 0) == 0); + CHECK(cmd.find("--azint-only") != std::string::npos); + CHECK(cmd.find("--azim-phi-bins 4") != std::string::npos); + CHECK(cmd.find("--azim-min-q") != std::string::npos); CHECK(cmd.find("in.h5") != std::string::npos); } diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 8984eb0c..4194ea3a 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -8,18 +8,12 @@ ADD_EXECUTABLE(jfjoch_extract_hkl jfjoch_extract_hkl.cpp TARGET_LINK_LIBRARIES(jfjoch_extract_hkl JFJochReader) INSTALL(TARGETS jfjoch_extract_hkl RUNTIME COMPONENT viewer) +# rugnux is the single offline analysis CLI: full pipeline by default, plus --azint-only +# (azimuthal integration) and --scale (re-scale/merge stored reflections) modes. ADD_EXECUTABLE(rugnux rugnux_cli.cpp) TARGET_LINK_LIBRARIES(rugnux Rugnux JFJochReader JFJochImageAnalysis JFJochWriter) INSTALL(TARGETS rugnux RUNTIME COMPONENT viewer) -ADD_EXECUTABLE(rugnux_scale rugnux_scale.cpp) -TARGET_LINK_LIBRARIES(rugnux_scale JFJochReader JFJochImageAnalysis JFJochWriter) -INSTALL(TARGETS rugnux_scale RUNTIME COMPONENT viewer) - -ADD_EXECUTABLE(jfjoch_azint jfjoch_azint.cpp) -TARGET_LINK_LIBRARIES(jfjoch_azint Rugnux JFJochReader JFJochImageAnalysis JFJochWriter) -INSTALL(TARGETS jfjoch_azint RUNTIME COMPONENT viewer) - # In-place re-compress /entry/data/data of a _data file from bitshuffle/LZ4 to bitshuffle/zstd. ADD_EXECUTABLE(jfjoch_recompress jfjoch_recompress.cpp) TARGET_LINK_LIBRARIES(jfjoch_recompress JFJochHDF5Wrappers) @@ -29,8 +23,6 @@ INSTALL(TARGETS jfjoch_recompress RUNTIME COMPONENT viewer) IF (WIN32) TARGET_LINK_LIBRARIES(jfjoch_extract_hkl wingetopt) TARGET_LINK_LIBRARIES(rugnux wingetopt) - TARGET_LINK_LIBRARIES(rugnux_scale wingetopt) - TARGET_LINK_LIBRARIES(jfjoch_azint wingetopt) ENDIF() # Online / hardware tools (broker, FPGA, receiver, detector). They link the non-portable diff --git a/tools/jfjoch_azint.cpp b/tools/jfjoch_azint.cpp deleted file mode 100644 index c525a6ac..00000000 --- a/tools/jfjoch_azint.cpp +++ /dev/null @@ -1,296 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#include -#include -#include -#include -#include -#include -#include - -#include "../reader/JFJochHDF5Reader.h" -#include "../common/Logger.h" -#include "../common/Definitions.h" -#include "../common/DiffractionExperiment.h" -#include "../common/PixelMask.h" -#include "../common/print_license.h" -#include "../rugnux/Rugnux.h" - -void print_usage() { - std::cout << "Usage ./jfjoch_azint {} " << std::endl; - std::cout << "Runs CPU azimuthal integration on a Jungfraujoch HDF5 file and writes _process.h5" << std::endl; - std::cout << "Options:" << std::endl; - std::cout << " -o, --output-prefix Output file prefix (default: output)" << std::endl; - std::cout << " -N, --threads Number of threads (default: 1)" << std::endl; - std::cout << " -s, --start-image Start image number (default: 0)" << std::endl; - std::cout << " -e, --end-image End image number (default: all)" << std::endl; - std::cout << " -t, --stride Image stride (default: 1)" << std::endl; - std::cout << " -v, --verbose Verbose output" << std::endl; - std::cout << std::endl; - - std::cout << " Azimuthal integration (defaults taken from the input file)" << std::endl; - std::cout << " --min-q Minimum Q for integration (1/A)" << std::endl; - std::cout << " --max-q Maximum Q for integration (1/A)" << std::endl; - std::cout << " --q-spacing Q bin spacing (1/A)" << std::endl; - std::cout << " --azimuthal-bins Number of azimuthal bins (default: 1)" << std::endl; - std::cout << " --polarization-correction Enable/disable polarization correction" << std::endl; - std::cout << " --solid-angle-correction Enable/disable solid angle correction" << std::endl; - std::cout << std::endl; - - std::cout << " Geometry overrides (defaults taken from the input file)" << std::endl; - std::cout << " --beam-x Beam center X (pixel)" << std::endl; - std::cout << " --beam-y Beam center Y (pixel)" << std::endl; - std::cout << " --detector-distance Detector distance (mm)" << std::endl; - std::cout << " --wavelength Wavelength (A)" << std::endl; - std::cout << " --rot1 PONI rotation 1 (rad)" << std::endl; - std::cout << " --rot2 PONI rotation 2 (rad)" << std::endl; - std::cout << " --polarization Polarization factor" << std::endl; -} - -enum { - OPT_MIN_Q = 1000, - OPT_MAX_Q, - OPT_Q_SPACING, - OPT_AZIMUTHAL_BINS, - OPT_POLARIZATION_CORRECTION, - OPT_SOLID_ANGLE_CORRECTION, - OPT_BEAM_X, - OPT_BEAM_Y, - OPT_DETECTOR_DISTANCE, - OPT_WAVELENGTH, - OPT_ROT1, - OPT_ROT2, - OPT_POLARIZATION -}; - -static option long_options[] = { - {"verbose", no_argument, nullptr, 'v'}, - {"output-prefix", required_argument, nullptr, 'o'}, - {"threads", required_argument, nullptr, 'N'}, - {"start-image", required_argument, nullptr, 's'}, - {"end-image", required_argument, nullptr, 'e'}, - {"stride", required_argument, nullptr, 't'}, - - {"min-q", required_argument, nullptr, OPT_MIN_Q}, - {"max-q", required_argument, nullptr, OPT_MAX_Q}, - {"q-spacing", required_argument, nullptr, OPT_Q_SPACING}, - {"azimuthal-bins", required_argument, nullptr, OPT_AZIMUTHAL_BINS}, - {"polarization-correction", required_argument, nullptr, OPT_POLARIZATION_CORRECTION}, - {"solid-angle-correction", required_argument, nullptr, OPT_SOLID_ANGLE_CORRECTION}, - - {"beam-x", required_argument, nullptr, OPT_BEAM_X}, - {"beam-y", required_argument, nullptr, OPT_BEAM_Y}, - {"detector-distance", required_argument, nullptr, OPT_DETECTOR_DISTANCE}, - {"wavelength", required_argument, nullptr, OPT_WAVELENGTH}, - {"rot1", required_argument, nullptr, OPT_ROT1}, - {"rot2", required_argument, nullptr, OPT_ROT2}, - {"polarization", required_argument, nullptr, OPT_POLARIZATION}, - {nullptr, 0, nullptr, 0} -}; - -bool parse_on_off(const char *arg, bool &out) { - std::string s = arg ? arg : ""; - std::transform(s.begin(), s.end(), s.begin(), - [](unsigned char c) { return static_cast(std::tolower(c)); }); - if (s == "on" || s == "1" || s == "true" || s == "yes") { - out = true; - return true; - } - if (s == "off" || s == "0" || s == "false" || s == "no") { - out = false; - return true; - } - return false; -} - -namespace { - std::atomic g_active_process{nullptr}; - void handle_sigint(int) { - if (auto *p = g_active_process.load()) - p->Cancel(); - } -} - -int main(int argc, char **argv) { - for (int i = 0; i < argc; i++) - std::cout << argv[i] << " "; - std::cout << std::endl << std::endl; - - RegisterHDF5Filter(); - print_license("jfjoch_azint"); - - Logger logger("jfjoch_azint"); - - std::string output_prefix = "output"; - int nthreads = 1; - int start_image = 0; - int end_image = -1; // -1 indicates process until end - int image_stride = 1; - bool verbose = false; - - // Azimuthal integration overrides (default: keep value from input file) - std::optional min_q; - std::optional max_q; - std::optional q_spacing; - std::optional azimuthal_bins; - std::optional polarization_correction; - std::optional solid_angle_correction; - - // Geometry overrides (default: keep value from input file) - std::optional beam_x; - std::optional beam_y; - std::optional detector_distance_mm; - std::optional wavelength_A; - std::optional rot1_rad; - std::optional rot2_rad; - std::optional polarization_factor; - - if (argc == 1) { - print_usage(); - exit(EXIT_FAILURE); - } - - int opt; - int option_index = 0; - const char *short_opts = "vo:N:s:e:t:"; - - while ((opt = getopt_long(argc, argv, short_opts, long_options, &option_index)) != -1) { - switch (opt) { - case 'o': output_prefix = optarg; break; - case 'v': verbose = true; break; - case 'N': nthreads = atoi(optarg); break; - case 's': start_image = atoi(optarg); break; - case 'e': end_image = atoi(optarg); break; - case 't': image_stride = atoi(optarg); break; - - case OPT_MIN_Q: min_q = atof(optarg); break; - case OPT_MAX_Q: max_q = atof(optarg); break; - case OPT_Q_SPACING: q_spacing = atof(optarg); break; - case OPT_AZIMUTHAL_BINS: azimuthal_bins = atoi(optarg); break; - case OPT_POLARIZATION_CORRECTION: { - bool value; - if (!parse_on_off(optarg, value)) { - logger.Error("Invalid polarization correction value (expected on|off): {}", optarg); - exit(EXIT_FAILURE); - } - polarization_correction = value; - break; - } - case OPT_SOLID_ANGLE_CORRECTION: { - bool value; - if (!parse_on_off(optarg, value)) { - logger.Error("Invalid solid angle correction value (expected on|off): {}", optarg); - exit(EXIT_FAILURE); - } - solid_angle_correction = value; - break; - } - - case OPT_BEAM_X: beam_x = atof(optarg); break; - case OPT_BEAM_Y: beam_y = atof(optarg); break; - case OPT_DETECTOR_DISTANCE: detector_distance_mm = atof(optarg); break; - case OPT_WAVELENGTH: wavelength_A = atof(optarg); break; - case OPT_ROT1: rot1_rad = atof(optarg); break; - case OPT_ROT2: rot2_rad = atof(optarg); break; - case OPT_POLARIZATION: polarization_factor = atof(optarg); break; - - default: - print_usage(); - exit(EXIT_FAILURE); - } - } - - if (optind != argc - 1) { - logger.Error("Input file not specified"); - print_usage(); - exit(EXIT_FAILURE); - } - - const std::string input_file = argv[optind]; - logger.Verbose(verbose); - - if (image_stride <= 0) { - logger.Error("Image stride must be positive"); - exit(EXIT_FAILURE); - } - - // 1. Read input file - JFJochHDF5Reader reader; - try { - reader.ReadFile(input_file); - } catch (const std::exception &e) { - logger.Error("Error reading input file: {}", e.what()); - exit(EXIT_FAILURE); - } - - const auto dataset = reader.GetDataset(); - if (!dataset) { - logger.Error("No experiment dataset found in the input file"); - exit(EXIT_FAILURE); - } - logger.Info("Loaded dataset from {}", input_file); - - // 2. Build experiment: defaults from the input file, overridden by command line. Output and - // runtime invariants are set inside Rugnux; here we only configure the geometry and - // azimuthal-integration settings. - DiffractionExperiment experiment(dataset->experiment); - if (beam_x.has_value()) experiment.BeamX_pxl(beam_x.value()); - if (beam_y.has_value()) experiment.BeamY_pxl(beam_y.value()); - if (detector_distance_mm.has_value()) experiment.DetectorDistance_mm(detector_distance_mm.value()); - if (wavelength_A.has_value()) experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / wavelength_A.value()); - if (rot1_rad.has_value()) experiment.PoniRot1_rad(rot1_rad.value()); - if (rot2_rad.has_value()) experiment.PoniRot2_rad(rot2_rad.value()); - if (polarization_factor.has_value()) experiment.PolarizationFactor(polarization_factor.value()); - - AzimuthalIntegrationSettings azint_settings = experiment.GetAzimuthalIntegrationSettings(); - if (min_q.has_value() || max_q.has_value()) - azint_settings.QRange_recipA(min_q.value_or(azint_settings.GetLowQ_recipA()), - max_q.value_or(azint_settings.GetHighQ_recipA())); - if (q_spacing.has_value()) azint_settings.QSpacing_recipA(q_spacing.value()); - if (azimuthal_bins.has_value()) azint_settings.AzimuthalBinCount(azimuthal_bins.value()); - if (polarization_correction.has_value()) azint_settings.PolarizationCorrection(polarization_correction.value()); - if (solid_angle_correction.has_value()) azint_settings.SolidAngleCorrection(solid_angle_correction.value()); - experiment.ImportAzimuthalIntegrationSettings(azint_settings); - - logger.Info("Geometry: beam ({:.2f}, {:.2f}) pxl, distance {:.2f} mm, wavelength {:.5f} A", - experiment.GetBeamX_pxl(), experiment.GetBeamY_pxl(), - experiment.GetDetectorDistance_mm(), experiment.GetWavelength_A()); - logger.Info("Azimuthal integration: Q range [{:.4f}, {:.4f}] 1/A, spacing {:.4f} 1/A, {} Q bins x {} azimuthal bins", - azint_settings.GetLowQ_recipA(), azint_settings.GetHighQ_recipA(), - azint_settings.GetQSpacing_recipA(), azint_settings.GetQBinCount(), - azint_settings.GetAzimuthalBinCount()); - logger.Info("Corrections: polarization {}, solid angle {}", - azint_settings.IsPolarizationCorrection() ? "on" : "off", - azint_settings.IsSolidAngleCorrection() ? "on" : "off"); - - // 3. Run the shared azimuthal-integration workflow. - ProcessConfig config; - config.mode = ProcessMode::AzimuthalIntegration; - config.start_image = start_image; - config.end_image = end_image; - config.stride = image_stride; - config.nthreads = nthreads; - config.output_prefix = output_prefix; - - Rugnux process(reader, experiment, dataset->pixel_mask, config); - - g_active_process = &process; - std::signal(SIGINT, handle_sigint); - - ProcessResult result; - try { - result = process.Run(); - } catch (const std::exception &e) { - logger.Error("Processing failed: {}", e.what()); - exit(EXIT_FAILURE); - } - g_active_process = nullptr; - - // 4. Report statistics - std::cout << fmt::format("Processing time: {:.2f} s", result.processing_time_s) << std::endl; - std::cout << fmt::format("Frame rate: {:.2f} Hz", result.frame_rate_hz) << std::endl; - std::cout << fmt::format("Total throughput: {:.2f} MB/s", result.throughput_MBs) << std::endl; - if (result.cancelled) - logger.Warning("Processing was cancelled after {} images", result.images_processed); -} diff --git a/tools/rugnux_cli.cpp b/tools/rugnux_cli.cpp index bcddb708..2f7448f3 100644 --- a/tools/rugnux_cli.cpp +++ b/tools/rugnux_cli.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -12,10 +13,17 @@ #include "../reader/JFJochHDF5Reader.h" #include "../common/Logger.h" +#include "../common/Definitions.h" #include "../common/DiffractionExperiment.h" #include "../common/PixelMask.h" #include "../common/print_license.h" #include "../image_analysis/LoadFCalcFromMtz.h" +#include "../image_analysis/UpdateReflectionResolution.h" +#include "../image_analysis/WriteReflections.h" +#include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/ScaleOnTheFly.h" +#include "../image_analysis/scale_merge/RotationScaleMerge.h" +#include "../image_analysis/scale_merge/TwinningAnalysis.h" #include "../rugnux/Rugnux.h" // Default rot3d per-frame scale-G smoothing range (XDS DELPHI-like), in degrees of rotation. @@ -36,6 +44,11 @@ void print_usage() { std::cout << " -v, --verbose Verbose output" << std::endl; std::cout << std::endl; + std::cout << " Modes (default: full analysis - spot finding, indexing, integration and merging)" << std::endl; + std::cout << " --azint-only Only run azimuthal integration (no spot finding/indexing); writes _process.h5" << std::endl; + std::cout << " --scale Only re-scale/merge the already-integrated reflections in (no re-integration)" << std::endl; + std::cout << std::endl; + std::cout << " Spot finding" << std::endl; std::cout << " --spot-sigma Noise sigma level for spot finding (default: 3.0)" << std::endl; std::cout << " --spot-threshold Photon count threshold for spot finding (default: 10)" << std::endl; @@ -57,9 +70,9 @@ void print_usage() { std::cout << " -r, --refine Geometry refinement algorithm (none|orientation|beam_and_lattice)" << std::endl; std::cout << std::endl; - std::cout << " Scaling and merging" << std::endl; - std::cout << " -M, --scale-merge Scale and merge (refine mosaicity); automatic for rotation data" << std::endl; - std::cout << " --scale-fulls After -P rot3d combine, refit a per-frame scale on the fulls (XDS order, Unity model); implies -M. Default ON for rot3d" << std::endl; + std::cout << " Scaling and merging (on by default)" << std::endl; + std::cout << " --no-merge Skip scaling and merging; write only the per-image _process.h5" << std::endl; + std::cout << " --scale-fulls rot3d: after the 3D combine, refit a per-frame scale on the fulls (XDS order, Unity model). Default ON for rot3d" << std::endl; std::cout << " --no-scale-fulls Disable the rot3d scale-fulls refit (it is on by default for rot3d)" << std::endl; std::cout << " --write-process-h5 Also write the (large) _process.h5 when merging (default: only .mtz/.cif when merging)" << std::endl; std::cout << " --smooth-g[=deg] rot3d: smooth per-frame scale G over a deg-degree rotation range (XDS DELPHI-like) before the combine (default: 5 for rot3d; 0 = off)" << std::endl; @@ -86,6 +99,18 @@ void print_usage() { std::cout << " --azim-min-q Azimuthal-integration minimum Q (1/A)" << std::endl; std::cout << " --azim-max-q Azimuthal-integration maximum Q (1/A)" << std::endl; std::cout << " --azim-phi-bins Number of azimuthal (phi) bins (default: 1)" << std::endl; + std::cout << " --polarization-correction Enable/disable azimuthal polarization correction" << std::endl; + std::cout << " --solid-angle-correction Enable/disable azimuthal solid angle correction" << std::endl; + std::cout << std::endl; + + std::cout << " Geometry overrides (defaults taken from the input file)" << std::endl; + std::cout << " --beam-x Beam center X (pixel)" << std::endl; + std::cout << " --beam-y Beam center Y (pixel)" << std::endl; + std::cout << " --detector-distance Detector distance (mm)" << std::endl; + std::cout << " --wavelength Wavelength (A)" << std::endl; + std::cout << " --rot1 PONI rotation 1 (rad)" << std::endl; + std::cout << " --rot2 PONI rotation 2 (rad)" << std::endl; + std::cout << " --polarization Polarization factor" << std::endl; } enum { @@ -118,7 +143,19 @@ enum { OPT_FORCE_STILL, OPT_AZIM_MIN_Q, OPT_AZIM_MAX_Q, - OPT_AZIM_PHI_BINS + OPT_AZIM_PHI_BINS, + OPT_AZINT_ONLY, + OPT_SCALE, + OPT_NO_MERGE, + OPT_POLARIZATION_CORRECTION, + OPT_SOLID_ANGLE_CORRECTION, + OPT_BEAM_X, + OPT_BEAM_Y, + OPT_DETECTOR_DISTANCE, + OPT_WAVELENGTH, + OPT_ROT1, + OPT_ROT2, + OPT_POLARIZATION }; static option long_options[] = { @@ -136,7 +173,9 @@ static option long_options[] = { {"space-group", required_argument, nullptr, 'S'}, {"anomalous", no_argument, nullptr, 'A'}, {"refine-bfactor", no_argument, nullptr, 'B'}, - {"scale-merge", no_argument, nullptr, 'M'}, + {"azint-only", no_argument, nullptr, OPT_AZINT_ONLY}, + {"scale", no_argument, nullptr, OPT_SCALE}, + {"no-merge", no_argument, nullptr, OPT_NO_MERGE}, {"scale-fulls", no_argument, nullptr, OPT_SCALE_FULLS}, {"no-scale-fulls", no_argument, nullptr, OPT_NO_SCALE_FULLS}, {"write-process-h5", no_argument, nullptr, OPT_WRITE_PROCESS_H5}, @@ -150,6 +189,15 @@ static option long_options[] = { {"azim-min-q", required_argument, nullptr, OPT_AZIM_MIN_Q}, {"azim-max-q", required_argument, nullptr, OPT_AZIM_MAX_Q}, {"azim-phi-bins", required_argument, nullptr, OPT_AZIM_PHI_BINS}, + {"polarization-correction", required_argument, nullptr, OPT_POLARIZATION_CORRECTION}, + {"solid-angle-correction", required_argument, nullptr, OPT_SOLID_ANGLE_CORRECTION}, + {"beam-x", required_argument, nullptr, OPT_BEAM_X}, + {"beam-y", required_argument, nullptr, OPT_BEAM_Y}, + {"detector-distance", required_argument, nullptr, OPT_DETECTOR_DISTANCE}, + {"wavelength", required_argument, nullptr, OPT_WAVELENGTH}, + {"rot1", required_argument, nullptr, OPT_ROT1}, + {"rot2", required_argument, nullptr, OPT_ROT2}, + {"polarization", required_argument, nullptr, OPT_POLARIZATION}, {"redo-rotation-spots", no_argument, nullptr, OPT_REDO_ROTATION_SPOTS}, {"force-rotation-lattice", required_argument, nullptr, OPT_FORCE_ROTATION_LATTICE}, @@ -192,6 +240,21 @@ bool parse_float_strict(const std::string &t, float &out) { } }; +bool parse_on_off(const char *arg, bool &out) { + std::string s = arg ? arg : ""; + std::transform(s.begin(), s.end(), s.begin(), + [](unsigned char c) { return static_cast(std::tolower(c)); }); + if (s == "on" || s == "1" || s == "true" || s == "yes") { + out = true; + return true; + } + if (s == "off" || s == "0" || s == "false" || s == "no") { + out = false; + return true; + } + return false; +} + std::optional parse_unit_cell_arg(const char *arg) { if (!arg) return std::nullopt; @@ -301,18 +364,25 @@ int main(int argc, char **argv) { int image_stride = 1; bool verbose = false; + bool azint_only = false; // --azint-only: azimuthal integration only (no spot finding/indexing) + bool scale_only = false; // --scale: re-scale/merge stored reflections only (no re-integration) bool rotation_indexing = false; bool force_still = false; // --force-still: process a rotation dataset as stills (indexing + scaling) bool two_pass_rotation = true; bool reuse_rotation_spots = true; int rotation_indexing_image_count = 100; std::optional rotation_indexing_range; - bool run_scaling = false; + bool run_scaling = true; // merge is on by default; --no-merge turns it off std::optional scale_fulls_arg; // --scale-fulls / --no-scale-fulls; default on for rot3d bool write_process_h5_flag = false; // --write-process-h5; also write _process.h5 when merging std::optional detect_ice_rings; // --detect-ice-rings[=on|off]; unset => use the dataset (file) value std::optional min_q, max_q, q_spacing; // azimuthal integration range / -q spacing (1/A) std::optional azimuthal_bins; // --azimuthal-bins + std::optional polarization_correction; // --polarization-correction (azimuthal integration) + std::optional solid_angle_correction; // --solid-angle-correction (azimuthal integration) + + // 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, polarization_factor; std::optional smooth_g_deg_arg; // --smooth-g[=deg]; default 5 deg for rot3d, 0 (off) otherwise bool anomalous_mode = false; std::optional space_group_number; @@ -352,7 +422,7 @@ int main(int argc, char **argv) { int opt; int option_index = 0; - const char *short_opts = "vo:N:s:e:t:R::X:C:z:FABS:Mr:q:"; + const char *short_opts = "vo:N:s:e:t:R::X:C:z:FABS:r:q:"; while ((opt = getopt_long(argc, argv, short_opts, long_options, &option_index)) != -1) { switch (opt) { @@ -513,11 +583,16 @@ int main(int argc, char **argv) { max_spot_count_override = atoll(optarg); logger.Info("Max spot count overridden to {}", max_spot_count_override.value()); break; - case 'M': - run_scaling = true; + case OPT_AZINT_ONLY: + azint_only = true; + break; + case OPT_SCALE: + scale_only = true; + break; + case OPT_NO_MERGE: + run_scaling = false; break; case OPT_SCALE_FULLS: - run_scaling = true; scale_fulls_arg = true; break; case OPT_NO_SCALE_FULLS: @@ -537,8 +612,6 @@ int main(int argc, char **argv) { case OPT_WRITE_PROCESS_H5: write_process_h5_flag = true; break; - run_scaling = true; - break; case OPT_SMOOTH_G: smooth_g_deg_arg = optarg ? std::stod(optarg) : SMOOTH_G_DEFAULT_DEG; break; @@ -599,6 +672,31 @@ int main(int argc, char **argv) { case OPT_AZIM_PHI_BINS: azimuthal_bins = atoi(optarg); break; + case OPT_POLARIZATION_CORRECTION: { + bool value; + if (!parse_on_off(optarg, value)) { + logger.Error("Invalid polarization correction value (expected on|off): {}", optarg); + exit(EXIT_FAILURE); + } + polarization_correction = value; + break; + } + case OPT_SOLID_ANGLE_CORRECTION: { + bool value; + if (!parse_on_off(optarg, value)) { + logger.Error("Invalid solid angle correction value (expected on|off): {}", optarg); + exit(EXIT_FAILURE); + } + solid_angle_correction = value; + break; + } + case OPT_BEAM_X: beam_x = atof(optarg); break; + case OPT_BEAM_Y: beam_y = atof(optarg); break; + case OPT_DETECTOR_DISTANCE: detector_distance_mm = atof(optarg); break; + case OPT_WAVELENGTH: wavelength_A = atof(optarg); break; + case OPT_ROT1: rot1_rad = atof(optarg); break; + case OPT_ROT2: rot2_rad = atof(optarg); break; + case OPT_POLARIZATION: polarization_factor = atof(optarg); break; case OPT_SCALING_ITERATIONS: scaling_iter = atoi(optarg); if (scaling_iter <= 0) { @@ -629,6 +727,11 @@ int main(int argc, char **argv) { input_file = argv[optind]; logger.Verbose(verbose); + if (azint_only && scale_only) { + logger.Error("--azint-only and --scale are mutually exclusive"); + exit(EXIT_FAILURE); + } + // Validate space group number early const gemmi::SpaceGroup *space_group = nullptr; if (space_group_number.has_value()) { @@ -695,6 +798,115 @@ int main(int argc, char **argv) { } } + // --scale: re-scale and merge the already-integrated reflections stored in the input file, + // without re-running spot finding or integration (folded in from the former rugnux_scale tool). + if (scale_only) { + const auto total_images = static_cast(reader.GetNumberOfImages()); + const int last_image = (end_image < 0 || end_image >= total_images) ? total_images - 1 : end_image; + auto reflections = reader.ReadReflections(start_image, last_image); + + DiffractionExperiment experiment(dataset->experiment); + experiment.BitDepthImage(32).Compression(CompressionAlgorithm::BSHUF_LZ4); + experiment.FilePrefix(output_prefix); + experiment.Mode(DetectorMode::Standard); + experiment.PixelSigned(true); + experiment.OverwriteExistingFiles(true); + experiment.PolarizationFactor(0.99); + experiment.SetFileWriterFormat(FileWriterFormat::NXmxLegacy); + // Keep the space group stored in the input file (written by the full pipeline) unless -S overrides. + if (space_group_number.has_value()) + experiment.SpaceGroupNumber(space_group_number); + experiment.NumTriggers(1); + + // A rotation (goniometer) dataset uses RotationScaleMerge unless --force-still asks for stills scaling. + IndexingSettings indexing_settings; + indexing_settings.RotationIndexing(experiment.GetGoniometer().has_value() && !force_still); + experiment.ImportIndexingSettings(indexing_settings); + + ScalingSettings scaling_settings; + if (d_min_scale_merge) + scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value()); + scaling_settings.MergeFriedel(!anomalous_mode); + scaling_settings.RefineB(refine_bfactor); + scaling_settings.MinPartiality(min_partiality); + scaling_settings.MinCCForImage(min_image_cc / 100.0); // --min-image-cc is percent; the setting is a fraction + if (intensity_format) + scaling_settings.FileFormat(intensity_format.value()); + experiment.ImportScalingSettings(scaling_settings); + + if (!experiment.GetUnitCell()) { + logger.Error("Experiment unit cell not found, cannot update reflection resolution"); + exit(EXIT_FAILURE); + } + auto refl_stats = UpdateReflectionResolution(experiment.GetUnitCell().value(), reflections); + logger.Info("Read {} reflections from {} images", refl_stats.n_reflections, refl_stats.n_images); + experiment.ImagesPerTrigger(refl_stats.n_images); + + const auto scale_start = std::chrono::steady_clock::now(); + std::vector merged_reflections; + MergeStatistics merged_statistics; + double error_model_isa = 0.0; + + // Rotation (rot3d): the dedicated RotationScaleMerge does the whole self-scale -> 3D combine -> + // merge. It does not support external-reference scaling, B-factor or wedge refinement. Everything + // else (stills, reference scaling) uses ScaleOnTheFly + MergeOnTheFly. + const bool is_rotation = experiment.IsRotationIndexing(); + if (is_rotation) { + if (!reference_data.empty() || experiment.GetScalingSettings().GetRefineB() + || experiment.GetRefineRotationWedgeInScaling() + || experiment.GetScalingSettings().GetRotationWedgeForScaling().has_value()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Rotation scaling/merging (RotationScaleMerge) does not support reference " + "scaling, B-factor refinement or wedge refinement"); + RotationScaleMerge rsm(experiment, reflections, experiment.GetUnitCell(), + scaling_iter, 0.0f, nthreads, logger); + rsm.Ingest(); + auto r = rsm.Run(false); + merged_reflections = std::move(r.merged); + merged_statistics = std::move(r.statistics); + error_model_isa = r.isa; + } else { + for (int i = 0; i < scaling_iter; i++) { + if (reference_data.empty()) + ScaleOnTheFly(experiment, MergeAll(experiment, reflections)).Scale(reflections, nthreads); + else + ScaleOnTheFly(experiment, reference_data).Scale(reflections, nthreads); + } + MergeOnTheFly merge_engine(experiment); + merge_engine.ReferenceCell(experiment.GetUnitCell()); + for (size_t i = 0; i < reflections.size(); ++i) + merge_engine.AddImage(reflections[i], static_cast(i)); + merged_reflections = merge_engine.ExportReflections(); + merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data); + error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0; + } + + logger.Info("Scale + merge completed in {:.2f} s ({} unique reflections)", + std::chrono::duration(std::chrono::steady_clock::now() - scale_start).count(), + merged_reflections.size()); + + std::cout << merged_statistics; + + // Space-group determination lives in the full rugnux pipeline; --scale only consumes a space + // group (from the file or -S) and merges in it. + const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value(); + if (!fixed_space_group) + logger.Warning("No space group in the input file or on the command line - merged in P1. " + "Re-run rugnux (which determines and stores the space group) or pass " + "-S to scale and merge in the correct symmetry."); + + const auto twin_sg_number = experiment.GetSpaceGroupNumber(); + const gemmi::SpaceGroup *twin_sg = twin_sg_number + ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr; + std::cout << std::endl << TwinningAnalysisToText(AnalyzeTwinning(merged_reflections, twin_sg)) << std::endl; + + if (!output_prefix.empty()) + WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics, + error_model_isa > 0 ? fmt::format("{:.2f}", error_model_isa) : "?", + output_prefix); + return 0; + } + uint64_t total_images_in_file = reader.GetNumberOfImages(); if (end_image < 0 || end_image > total_images_in_file) end_image = total_images_in_file; @@ -722,6 +934,72 @@ int main(int argc, char **argv) { // 2. Setup Experiment & Components DiffractionExperiment experiment(dataset->experiment); + + // Geometry overrides (default: keep the value stored in the input file). Applied before the + // azimuthal-integration settings are derived, which depend on the geometry. + if (beam_x) experiment.BeamX_pxl(beam_x.value()); + if (beam_y) experiment.BeamY_pxl(beam_y.value()); + if (detector_distance_mm) experiment.DetectorDistance_mm(detector_distance_mm.value()); + if (wavelength_A) experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / wavelength_A.value()); + if (rot1_rad) experiment.PoniRot1_rad(rot1_rad.value()); + if (rot2_rad) experiment.PoniRot2_rad(rot2_rad.value()); + if (polarization_factor) experiment.PolarizationFactor(polarization_factor.value()); + + // Azimuthal integration (default q-spacing 0.01 1/A, from AzimuthalIntegrationSettings): the profile + // resolves the narrow ice rings for the ice-ring score. Shared by --azint-only and full analysis. + // -q / --azim-* / correction flags override; defaults come from the input file. + { + AzimuthalIntegrationSettings azint_settings = experiment.GetAzimuthalIntegrationSettings(); + if (min_q || max_q) + azint_settings.QRange_recipA(min_q.value_or(azint_settings.GetLowQ_recipA()), + max_q.value_or(azint_settings.GetHighQ_recipA())); + if (q_spacing) + azint_settings.QSpacing_recipA(q_spacing.value()); + if (azimuthal_bins) + azint_settings.AzimuthalBinCount(azimuthal_bins.value()); + if (polarization_correction) + azint_settings.PolarizationCorrection(polarization_correction.value()); + if (solid_angle_correction) + azint_settings.SolidAngleCorrection(solid_angle_correction.value()); + experiment.ImportAzimuthalIntegrationSettings(azint_settings); + logger.Info("Azimuthal integration: Q [{:.4f}, {:.4f}] 1/A, spacing {:.4f}, {} Q x {} azimuthal bins", + azint_settings.GetLowQ_recipA(), azint_settings.GetHighQ_recipA(), + azint_settings.GetQSpacing_recipA(), azint_settings.GetQBinCount(), + azint_settings.GetAzimuthalBinCount()); + } + + // --azint-only: azimuthal integration only (no spot finding / indexing / scaling). Rugnux reads + // the geometry and azimuthal-integration settings configured above off the experiment. + if (azint_only) { + ProcessConfig config; + config.mode = ProcessMode::AzimuthalIntegration; + config.start_image = start_image; + config.end_image = end_image; + config.stride = image_stride; + config.nthreads = nthreads; + config.output_prefix = output_prefix; + + Rugnux process(reader, experiment, dataset->pixel_mask, config); + g_active_process = &process; + std::signal(SIGINT, handle_sigint); + + ProcessResult result; + try { + result = process.Run(); + } catch (const std::exception &e) { + logger.Error("Processing failed: {}", e.what()); + exit(EXIT_FAILURE); + } + g_active_process = nullptr; + + std::cout << fmt::format("Processing time: {:.2f} s", result.processing_time_s) << std::endl; + std::cout << fmt::format("Frame rate: {:.2f} Hz", result.frame_rate_hz) << std::endl; + std::cout << fmt::format("Total throughput: {:.2f} MB/s", result.throughput_MBs) << std::endl; + if (result.cancelled) + logger.Warning("Processing was cancelled after {} images", result.images_processed); + return 0; + } + experiment.BitDepthImage(32).Compression(CompressionAlgorithm::BSHUF_LZ4); experiment.FilePrefix(output_prefix); experiment.Mode(DetectorMode::Standard); // Ensure full image analysis @@ -770,10 +1048,8 @@ int main(int argc, char **argv) { "indexing). Use --force-still to treat it as stills."); } - // Rotation data is scaled and merged automatically (as if -M were given); stills still require - // an explicit -M. --force-still keeps rotation_indexing false, so stills stay opt-in. - if (rotation_indexing) - run_scaling = true; + // Scaling and merging are on by default (run_scaling initialised true); --no-merge turns them off + // for both rotation and stills, in which case only the per-image _process.h5 is written. // Configure Indexing IndexingSettings indexing_settings; @@ -842,25 +1118,6 @@ int main(int argc, char **argv) { : "profile (empirical)"); } - // Azimuthal integration (default q-spacing 0.01 1/A, from AzimuthalIntegrationSettings): the profile - // resolves the narrow ice rings for the ice-ring score. -q / --azim-* override. Applied before - // Rugnux builds the azint mapping from the experiment. - { - AzimuthalIntegrationSettings azint_settings = experiment.GetAzimuthalIntegrationSettings(); - if (min_q || max_q) - azint_settings.QRange_recipA(min_q.value_or(azint_settings.GetLowQ_recipA()), - max_q.value_or(azint_settings.GetHighQ_recipA())); - if (q_spacing) - azint_settings.QSpacing_recipA(q_spacing.value()); - if (azimuthal_bins) - azint_settings.AzimuthalBinCount(azimuthal_bins.value()); - experiment.ImportAzimuthalIntegrationSettings(azint_settings); - logger.Info("Azimuthal integration: Q [{:.4f}, {:.4f}] 1/A, spacing {:.4f}, {} Q x {} azimuthal bins", - azint_settings.GetLowQ_recipA(), azint_settings.GetHighQ_recipA(), - azint_settings.GetQSpacing_recipA(), azint_settings.GetQBinCount(), - azint_settings.GetAzimuthalBinCount()); - } - SpotFindingSettings spot_settings; spot_settings.enable = true; spot_settings.indexing = true; diff --git a/tools/rugnux_scale.cpp b/tools/rugnux_scale.cpp deleted file mode 100644 index 67a2a4ab..00000000 --- a/tools/rugnux_scale.cpp +++ /dev/null @@ -1,390 +0,0 @@ -// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#include -#include -#include -#include -#include -#include -#include - -#include "../reader/JFJochHDF5Reader.h" -#include "../common/Logger.h" -#include "../common/DiffractionExperiment.h" -#include "../common/time_utc.h" -#include "../common/print_license.h" -#include "../image_analysis/MXAnalysisWithoutFPGA.h" -#include "../writer/FileWriter.h" -#include "../image_analysis/IndexAndRefine.h" -#include "../common/JFJochReceiverPlots.h" -#include "../compression/JFJochCompressor.h" -#include "../image_analysis/LoadFCalcFromMtz.h" -#include "../image_analysis/scale_merge/Merge.h" -#include "../image_analysis/scale_merge/TwinningAnalysis.h" -#include "../image_analysis/scale_merge/ScaleOnTheFly.h" -#include "../image_analysis/scale_merge/RotationScaleMerge.h" -#include "../image_analysis/WriteReflections.h" -#include "../image_analysis/UpdateReflectionResolution.h" - -void print_usage() { - std::cout << "Usage ./rugnux_scale {} " << std::endl; - std::cout << "Options:" << std::endl; - std::cout << " -o, --output-prefix Output file prefix (default: output)" << std::endl; - std::cout << " -N, --threads Number of threads (default: 1)" << std::endl; - std::cout << " -s, --start-image Start image number (default: 0)" << std::endl; - std::cout << " -e, --end-image End image number (default: all)" << std::endl; - std::cout << " -v, --verbose Verbose output" << std::endl; - std::cout << "" << std::endl; - std::cout << " Scaling and merging" << std::endl; - std::cout << " -S, --space-group Space group number" << std::endl; - std::cout << " --scaling-high-resolution High resolution limit for scaling/merging (default: 0.0; no limit)" - << std::endl; - std::cout << " --force-still Scale as independent stills (ScaleOnTheFly) even for a rotation dataset " - "(default: rotation data uses RotationScaleMerge)" << std::endl; - std::cout << " -A, --anomalous Anomalous mode (don't merge Friedel pairs)" << std::endl; - std::cout << " -B, --refine-bfactor Refine per image B-factor" << std::endl; - std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << - std::endl; - std::cout << " --min-image-cc Per-image CC limit in percent (default: no limit)" << std::endl; - std::cout << " --scaling-iterations Number of scaling iterations with no reference data (default: 3)" - << std::endl; - std::cout << " --scaling-output Output format for scaling results mtz|cif|txt (default: cif)" << std::endl; - std::cout << " -z, --reference-mtz Reference MTZ file" << std::endl; - std::cout << " --reference-column