rugnux: accept a CBF sweep, decide the rotation axis sign, and search as far as a given cell needs

Three changes to the rotation path, all in the same first-pass block.

Input. Rugnux now holds a JFJochReader rather than the concrete HDF5 reader, so a sweep can come
from a directory of miniCBF as well as from a master file; the CLI picks by looking at the file.
Only the _process.h5 image links need the concrete reader and they ask for it by dynamic_cast -
there is nothing for a virtual dataset to point at in a directory of CBFs, so that file writes its
images instead of linking to them. --mode scale refuses a CBF sweep, there being no stored
reflections to re-scale.

Search bound. A reference cell longer than the FFT bound cannot be represented by the search meant
to find it, so the run reported "found no lattice" - strictly worse than not giving the cell at
all. Raise the bound to cover a given cell (margin 1.3, as DIALS uses when a cell is known), and
raise it in the long-axis rescue too, where it costs nothing: that pass already coarsens the
resolution to 3.5 A, which shrinks the transform by 1.75x, so spending it back on reach leaves the
rescue at 1.14x the standard pass. The rescue's constrained re-index gets a pool wide enough to
hold what the coarse pass recovered, which it did not before. ROTATION ONLY - stills fire the FFT
once per image across every worker and routinely run with a known cell, so widening the transform
there would cost the whole serial run for nothing. Clamped rather than throwing: the rescue can
recover an implausible axis, and a bad candidate must not take the run down with it.

Axis sign. The rotation-axis sign is a convention the input often cannot settle - a miniCBF header
names the axis but gives no direction, and an NXmx vector is only meaningful together with the
detector orientation, which nothing here reads. It is decidable from the data, though: the wrong
sign does not index at all (measured on one sweep, changing nothing else: every frame indexed one
way, none the other). So after a poor first pass, try the other sign and keep whichever indexes
more validation frames - the same count the scheme choice already uses, so no new metric and no new
threshold. It flips the AXIS, not the angles, because prediction reads the axis too; the decision
then holds for the rest of the run. Runs before the long-axis rescue, since with the sign wrong
every candidate lattice is wrong. Only after a poor pass, so a correctly-signed file costs nothing,
and spot finding is not repeated.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

rugnux/Rugnux.cpp

@@ -290,7 +290,7 @@ namespace {
 
 }
 
-Rugnux::Rugnux(JFJochHDF5Reader &reader, DiffractionExperiment experiment,
+Rugnux::Rugnux(JFJochReader &reader, DiffractionExperiment experiment,
                              PixelMask pixel_mask, ProcessConfig config)
     : reader_(reader), experiment_(std::move(experiment)),
       pixel_mask_(std::move(pixel_mask)), config_(std::move(config)),
@@ -354,6 +354,38 @@ static void LogXDSGeometry(const DiffractionExperiment &experiment) {
                     goniometer->GetAxis().x, goniometer->GetAxis().y, goniometer->GetAxis().z);
 }
 
+// The FFT recovers a basis vector only up to fft_max_unit_cell_A: FFTIndexer sizes its projected
+// histogram from that value and the transform's last usable bin IS that length. So a reference cell
+// longer than the bound cannot be represented by the search meant to find it - the run then reports
+// "found no lattice", which is strictly worse than not giving the cell at all (measured on a 545 A
+// axis: a partial triclinic solution without -C, nothing with it). Raise the bound to cover the cell,
+// with the margin DIALS uses when a cell is known (max_cell = 1.3 * longest axis, see
+// dials/algorithms/indexing/indexer.py find_max_cell).
+//
+// ROTATION ONLY. Stills fire the FFT once per image across every worker and routinely run with a
+// known cell, so widening the transform there would cost the whole serial run for nothing; the
+// caller applies this inside its rotation branch.
+static void RaiseFFTBoundForKnownCell(IndexingSettings &settings, const DiffractionExperiment &x) {
+    const auto cell = x.GetUnitCell();
+    if (!cell.has_value() || !cell->is_finite())
+        return;
+
+    constexpr float KNOWN_CELL_MARGIN = 1.3f;
+    const float longest = std::max({cell->a, cell->b, cell->c});
+    // Clamp rather than let the setter throw. A cell that needs more than the search can ever reach
+    // is a real possibility - the long-axis rescue calls this with whatever its coarse pass recovered,
+    // and that has been seen to be an implausible harmonic - and a bad candidate must not take the
+    // whole run down. Searching as far as the indexer goes is the best available answer; if the cell
+    // really is that long the pass simply fails to index, which the caller already handles.
+    const float needed = std::min(KNOWN_CELL_MARGIN * longest,
+                                  IndexingSettings::fft_max_unit_cell_limit_A);
+    if (needed > settings.GetFFT_MaxUnitCell_A()) {
+        settings.FFT_MaxUnitCell_A(needed);
+        Logger("Rugnux").Info("FFT search bound raised to {:.0f} A to cover the given cell "
+                              "(longest axis {:.1f} A)", needed, longest);
+    }
+}
+
 void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, RugnuxObserver *observer) {
     Logger logger("Rugnux");
 
@@ -1390,29 +1422,37 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
     start_message.file_format = FileWriterFormat::NXmxIntegrated;
     start_message.write_master_file = true;
     start_message.write_images = false;
-    // With the stride: image i of the process file is source image start_image + i * stride, so the
-    // linked pictures line up with the per-image analysis written beside them.
-    start_message.hdf5_source_data =
-        reader_.GetHDF5DataSource(start_image, images_to_process, config_.stride);
+    // The process file LINKS to the original images rather than writing its own, which only an HDF5
+    // input can offer: there is nothing for a virtual dataset to point at in a directory of CBFs. So
+    // ask for the links only when the source is HDF5, and otherwise write the images into the file -
+    // the per-image analysis is the point of it either way, and without pictures nothing downstream
+    // could show the frame a reflection came from.
+    if (auto *hdf5 = dynamic_cast<JFJochHDF5Reader *>(&reader_)) {
+        // With the stride: image i of the process file is source image start_image + i * stride, so
+        // the linked pictures line up with the per-image analysis written beside them.
+        start_message.hdf5_source_data =
+            hdf5->GetHDF5DataSource(start_image, images_to_process, config_.stride);
 
-    // This file links to the ORIGINAL images instead of writing its own (write_images = false), so
-    // everything describing the pixels must describe those files. experiment_ cannot: it is pinned
-    // to signed 32-bit, the container the reader hands images out in, which would type the virtual
-    // dataset int32 over unsigned 16- or 32-bit sources and silently convert every value on read.
-    // Only consumers going through the virtual view ever saw it - our own reader resolves the
-    // mapping and opens the source itself.
-    const auto stored = reader_.GetStoredPixelFormat();
-    start_message.bit_depth_image = stored.bit_depth;
-    // Unsigned only, for the reason FillMessage gives: DIALS remaps the top two codes whenever the
-    // field is present, and on a signed image those land inside the trusted range.
-    start_message.bit_depth_readout.reset();
-    if (!stored.is_signed)
-        start_message.bit_depth_readout = stored.bit_depth;
-    start_message.pixel_signed = stored.is_signed;
-    start_message.error_value = stored.is_signed
-                                ? -(int64_t(1) << (stored.bit_depth - 1))          // INTx_MIN
-                                : (int64_t(1) << stored.bit_depth) - 1;            // UINTx_MAX
-    start_message.underload_value = stored.is_signed ? start_message.error_value.value() + 1 : 0;
+        // Everything describing the pixels must describe those linked files. experiment_ cannot: it
+        // is pinned to signed 32-bit, the container the reader hands images out in, which would type
+        // the virtual dataset int32 over unsigned 16- or 32-bit sources and silently convert every
+        // value on read. Only consumers going through the virtual view ever saw it - our own reader
+        // resolves the mapping and opens the source itself.
+        const auto stored = hdf5->GetStoredPixelFormat();
+        start_message.bit_depth_image = stored.bit_depth;
+        // Unsigned only, for the reason FillMessage gives: DIALS remaps the top two codes whenever
+        // the field is present, and on a signed image those land inside the trusted range.
+        start_message.bit_depth_readout.reset();
+        if (!stored.is_signed)
+            start_message.bit_depth_readout = stored.bit_depth;
+        start_message.pixel_signed = stored.is_signed;
+        start_message.error_value = stored.is_signed
+                                    ? -(int64_t(1) << (stored.bit_depth - 1))      // INTx_MIN
+                                    : (int64_t(1) << stored.bit_depth) - 1;        // UINTx_MAX
+        start_message.underload_value = stored.is_signed ? start_message.error_value.value() + 1 : 0;
+    } else {
+        start_message.write_images = true;   // no file to link to; the process file carries them
+    }
 
     // The full-detector mask is 4 B/px and the ROI map 2 B/px - tens of megabytes each - and nothing
     // but the process file reads them. Copy them where a writer is actually built rather than here:
@@ -1448,6 +1488,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
     // Full-analysis shared engines.
     std::unique_ptr<IndexerThreadPool> indexer_pool;
     std::unique_ptr<IndexAndRefine> indexer;
+    // The FFT search bound indexer_pool was built with. The long-axis rescue needs it to tell whether a
+    // cell it recovered still fits, since a pool's transform is sized once, at construction.
+    float indexing_bound_A = 0.0f;
     if (per_image_analysis) {
         auto indexing_settings = experiment_.GetIndexingSettings();
         // Two-pass rotation fires the FFT indexer only twice - the two first-pass schemes, run in
@@ -1462,7 +1505,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
             // Clamped to what IndexingSettings accepts - a big host would otherwise ask for more than
             // the setting's limit and throw before a single image is read.
             indexing_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS));
+            RaiseFFTBoundForKnownCell(indexing_settings, experiment_);
         }
+        indexing_bound_A = indexing_settings.GetFFT_MaxUnitCell_A();
         indexer_pool = std::make_unique<IndexerThreadPool>(indexing_settings, IndexerConstruction::OnFirstUse);
         indexer = std::make_unique<IndexAndRefine>(experiment_, indexer_pool.get());
         // With no per-image file to write, nothing reads the message's reflection list.
@@ -1858,6 +1903,40 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
 
         FirstPass best = pick_best(*indexer_pool, *indexer);
 
+        // Rotation-axis SIGN rescue. The sign is a convention the input cannot settle: a PILATUS
+        // miniCBF header carries only "Oscillation_axis omega" and no direction at all, and facilities
+        // disagree about the sign of the NXmx vector (measured across ten of them: most want the
+        // stored [-1,0,0], APS 19-ID and SPring-8 BL41XU want the opposite). Nothing downstream can
+        // repair it either - the default spindle fit refines the axis DIRECTION by fractions of a
+        // degree and cannot flip it, that being a discrete change.
+        //
+        // It is, however, decidable from the data: the wrong sign does not index at all. Measured on
+        // one dataset by changing nothing but the sign: 100% of frames indexed one way, 0 of 60 the
+        // other. So try the other sign and keep whichever indexes more frames - the same
+        // validation-frame count the scheme choice above already uses, so no new metric and no new
+        // threshold. Runs before the long-axis rescue: with the sign wrong every candidate lattice is
+        // wrong, and the coarse pass below would spend itself on a mis-signed cloud.
+        //
+        // Only after a poor pass, so a correctly-signed file costs nothing; and spot finding is not
+        // repeated, the cache being keyed by image and the spots not depending on the axis.
+        if (!cancelled_ && best.score < 0.5 * static_cast<double>(validation.size())) {
+            if (const auto gon = experiment_.GetGoniometer(); gon.has_value() && gon->IsScanning()) {
+                GoniometerAxis flipped = *gon;
+                flipped.Axis(-gon->GetAxis());
+                experiment_.Goniometer(flipped);
+                const FirstPass alt = pick_best(*indexer_pool, *indexer);
+                if (alt.result.has_value() && alt.score > best.score) {
+                    logger.Info("Rotation axis sign: the file's axis indexes {}/{} validation frames, "
+                                "the opposite sign {}/{} - adopting the opposite sign for this run",
+                                best.score, static_cast<int>(validation.size()),
+                                alt.score, static_cast<int>(validation.size()));
+                    best = alt;
+                } else {
+                    experiment_.Goniometer(*gon);   // the file was right; put it back
+                }
+            }
+        }
+
         // Long-axis rescue. When the de-novo cell indexes few validation frames, a long, finely-spaced
         // axis was likely lost: the unconstrained FFT either collapsed it to a short sub-multiple or let
         // a denser supercell over-fit the accumulated cloud (a small global-orientation error throws the
@@ -1875,6 +1954,16 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
             && best.score < 0.5 * static_cast<double>(validation.size())) {
             auto coarse_settings = experiment_.GetIndexingSettings();
             coarse_settings.FFT_HighResolution_A(3.5f);   // low-order reflections only -> robust long axis
+            // Widen the search bound too, and note this costs nothing: the histogram length is
+            // 4*pi*max_length/high_resolution, so coarsening the resolution to 3.5 A has already
+            // SHRUNK the transform by 1.75x - spending that back on reach leaves the coarse pass at
+            // 1.14x the standard one. Without it the rescue cannot recover what it exists to recover:
+            // a cell whose long axis is past the default 500 A bound is unrepresentable, so the
+            // "recovered" metric would still be a sub-multiple. Only reached after a poor pass, so a
+            // well-indexing crystal never pays. Stills never reach here.
+            constexpr float RESCUE_MAX_UNIT_CELL_A = 1000.0f;
+            if (coarse_settings.GetFFT_MaxUnitCell_A() < RESCUE_MAX_UNIT_CELL_A)
+                coarse_settings.FFT_MaxUnitCell_A(RESCUE_MAX_UNIT_CELL_A);
             coarse_settings.IndexingThreads(2);
             coarse_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS));
             IndexerThreadPool coarse_pool(coarse_settings, IndexerConstruction::OnFirstUse);
@@ -1906,7 +1995,22 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
                 logger.Info("Long-axis rescue: coarse pass recovered a {:.0f} A axis (was {:.0f} A); "
                             "re-indexing constrained by that cell", max_axis(*coarse_ref), max_axis(*best.result));
                 experiment_.SetUnitCell(coarse_ref->lattice.GetUnitCell());
-                const FirstPass constrained = pick_best(*indexer_pool, *indexer);
+                // The re-index must be able to REPRESENT the cell the coarse pass just recovered. The
+                // standard pool was sized at the default bound, so re-using it here would search for a
+                // long axis it cannot hold and hand back the sub-multiple the rescue exists to replace.
+                // count_indexed forces the lattice and never touches the pool, so the existing
+                // IndexAndRefine is reused; only the FFT needs the wider one.
+                auto wide_settings = experiment_.GetIndexingSettings();
+                wide_settings.IndexingThreads(2);
+                wide_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS));
+                RaiseFFTBoundForKnownCell(wide_settings, experiment_);
+                FirstPass constrained;
+                if (wide_settings.GetFFT_MaxUnitCell_A() > indexing_bound_A) {
+                    IndexerThreadPool wide_pool(wide_settings, IndexerConstruction::OnFirstUse);
+                    constrained = pick_best(wide_pool, *indexer);
+                } else {
+                    constrained = pick_best(*indexer_pool, *indexer);
+                }
                 experiment_.SetUnitCell(std::nullopt);   // leave the space-group / cell determination de-novo
                 if (constrained.result.has_value() && constrained.score > best.score)
                     best = constrained;

rugnux/Rugnux.h

@@ -24,6 +24,7 @@
 #include "../image_analysis/rotation_indexer/RotationIndexer.h"   // RotationIndexerResult
 #include "RugnuxCalibration.h"                                   // CalibrationMethod, CalibrationResult
 
+class JFJochReader;
 class JFJochHDF5Reader;
 
 // Offline reprocessing of a stored Jungfraujoch HDF5 dataset, shared by rugnux (all of its --mode
@@ -263,7 +264,10 @@ public:
 };
 
 class Rugnux {
-    JFJochHDF5Reader &reader_;
+    // The base interface, so a sweep can come from any reader: an HDF5 master, or a
+    // directory of PILATUS miniCBF. Only the _process.h5 image links need the concrete
+    // HDF5 reader, and that block asks for it by dynamic_cast.
+    JFJochReader &reader_;
     DiffractionExperiment experiment_;
     PixelMask pixel_mask_;
     ProcessConfig config_;
@@ -369,7 +373,7 @@ class Rugnux {
     ProcessResult RunAllPasses(RugnuxObserver *observer);
 
 public:
-    Rugnux(JFJochHDF5Reader &reader, DiffractionExperiment experiment,
+    Rugnux(JFJochReader &reader, DiffractionExperiment experiment,
                   PixelMask pixel_mask, ProcessConfig config);
 
     // Runs the configured workflow to completion or until Cancel(). Throws on setup failure.

rugnux/rugnux_cli.cpp

@@ -18,6 +18,7 @@
 #include <vector>
 
 #include "../reader/JFJochHDF5Reader.h"
+#include "../reader/JFJochCBFReader.h"
 #include "../common/Logger.h"
 #include "../common/Definitions.h"
 #include "../common/DiffractionExperiment.h"
@@ -1279,14 +1280,27 @@ static int RunRugnux(int argc, char **argv) {
     // and the beam-stop pre-scan instead of landing inside the first pass.
     auto cufft_warmup = std::async(std::launch::async, WarmUpCuFFT);
 
-    // 1. Read Input File
-    JFJochHDF5Reader reader;
+    // 1. Read Input File. Either an HDF5 master, or one PILATUS miniCBF frame - naming any frame of a
+    // sweep reads the whole sweep, which is how these are archived and how a user thinks of them.
+    JFJochHDF5Reader hdf5_reader;
+    JFJochCBFReader cbf_reader;
+    JFJochReader *reader_ptr = nullptr;
+    const bool input_is_cbf = JFJochCBFReader::CanRead(input_file);
     try {
-        reader.ReadFile(input_file);
+        if (input_is_cbf) {
+            cbf_reader.ReadFiles(input_file);
+            reader_ptr = &cbf_reader;
+            logger.Info("Read {} CBF images of the sweep containing {}",
+                        cbf_reader.GetNumberOfImages(), input_file);
+        } else {
+            hdf5_reader.ReadFile(input_file);
+            reader_ptr = &hdf5_reader;
+        }
     } catch (const std::exception &e) {
         logger.Error("Error reading input file: {}", e.what());
         exit(EXIT_FAILURE);
     }
+    JFJochReader &reader = *reader_ptr;
 
     const auto dataset = reader.GetDataset();
     if (!dataset) {
@@ -1362,9 +1376,15 @@ static int RunRugnux(int argc, char **argv) {
     // --mode 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 (mode == RugnuxMode::Scale) {
+        // Re-scaling reads reflections a previous run integrated, which only a _process.h5 holds.
+        if (input_is_cbf) {
+            logger.Error("--mode scale needs the integrated reflections in a _process.h5; "
+                         "a CBF sweep holds none. Run --mode mx on it first.");
+            exit(EXIT_FAILURE);
+        }
         const auto total_images = static_cast<int>(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);
+        auto reflections = hdf5_reader.ReadReflections(start_image, last_image);
 
         DiffractionExperiment experiment(dataset->experiment);
         configure_offline_output(experiment, output_prefix);