diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index db0f96f86..95c79ba8b 100644 --- a/rugnux/Rugnux.cpp +++ b/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(&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 indexer_pool; std::unique_ptr 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(indexing_settings, IndexerConstruction::OnFirstUse); indexer = std::make_unique(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(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(validation.size()), + alt.score, static_cast(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(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; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 1e863dc5c..ff607d58f 100644 --- a/rugnux/Rugnux.h +++ b/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. diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 5019e0a05..e5a9bff6f 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -18,6 +18,7 @@ #include #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(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);