mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 72 Wiki Activity
Files
996cd20106d92f5dd7b5a712bbce4b98ae302eda
Jungfraujoch/rugnux/Rugnux.cpp
T
jungfrauandClaude Opus 5 996cd20106 Make the beam-stop mask O(pixels) and parallel 
GetMask() was 3.28 s of the 4.78 s pre-scan on a 16M-pixel detector, all on
one thread. Four changes, none of which alters the mask:

dilate() was a multi-source BFS. On a full rectangle with no obstacles the
8-connected graph distance IS the Chebyshev distance - a path stepping towards
the target never has to leave the frame - so the result is a dilation by the
(2r+1) square clipped to the frame, which separates into a pass along x and a
pass along y. That is O(1) per pixel whatever r is, with no queue and no
4-bytes-per-pixel distance array (72 MB, allocated and filled five times per
call). The erode() case is the one that hurt: it dilates the COMPLEMENT, so on
a detector whose shadow is under 1% of the pixels it seeded the BFS from
essentially every pixel.

fill_holes() floods the background from the border. It now floods the bounding
box of the region grown by one: everything outside that box is background and
the box's own ring is background, so the whole outside is one border-connected
component and a background pixel inside the box is border-connected exactly
when it reaches the ring.

The three baseline iterations re-binned every pixel by radius and re-took a
median each time. The iteration only ever excludes pixels whose background is
below a cut, and dividing by a positive baseline is monotone, so a ring's
excluded pixels are exactly its lowest ones and the next median is an order
statistic of the same, unchanging ring. The rings are binned and sorted once;
each iteration then picks a rank and counts a prefix. Nine full-image passes
become one.

box_sum's vertical pass walked one column at a time, striding a whole row per
step and missing on every access; it now carries a strip of columns together.
Each row's and each column's running sum keeps its terms in its order, so the
floating-point rounding is unchanged - only the traversal differs. The pooled
COUNT is a count of at most 25 pixels, so it is an exact integer box sum now
rather than a floating-point one; the background itself stays in double,
because its running sum adds and subtracts across a whole row and in float the
two roundings would not cancel.

The per-pixel passes then run on all threads, and GetMask takes a thread count.

Measured on a 16M-pixel rotation dataset: GetMask 3.28 s -> 0.99 s, whole
pre-scan 4.78 s -> 2.37 s, whole run 1m10s -> 1m03s. The mask is unchanged on
both a 16M and a 2M-pixel dataset (139126 and 22143 shadow pixels), as are the
space group, the merged reflection count and the merging statistics.

Also corrected the comment on erode(): the dilation cannot seed outside the
frame, so outside behaves as foreground, not as complement.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-23 12:50:12 +02:00

2632 lines
166 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "Rugnux.h"
#include "ModelValidation.h"
#include <algorithm>
#include <atomic>
#include <chrono>
#include <cmath>
#include <functional>
#include <future>
#include <mutex>
#include <numeric>
#include <set>
#include <sstream>
#include "../reader/JFJochHDF5Reader.h"
#include "../common/JFJochMath.h"
#include "../common/Logger.h"
#include "../common/AzimuthalIntegrationMapping.h"
#include "../common/AzimuthalIntegrationProfile.h"
#include "../common/CUDAWrapper.h"
#include "../common/time_utc.h"
#include "../writer/FileWriter.h"
#include "../image_analysis/MXAnalysisWithoutFPGA.h"
#include "../image_analysis/beam_stop/ShadowFinder.h"
#include "../image_analysis/IndexAndRefine.h"
#include "../image_analysis/geom_refinement/BeamCenterFromBackground.h"
#include "../image_analysis/geom_refinement/BeamCenterFromSpots.h"
#include "../image_analysis/geom_refinement/GeometryRefiner.h"
#include "../image_analysis/indexing/IndexerThreadPool.h"
#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h"
#include "../image_analysis/azint/AzIntEngineCPU.h"
#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h"
#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h"
#include "../image_analysis/scale_merge/Merge.h"
#include "../image_analysis/scale_merge/RfreeFlags.h"
#include "../image_analysis/scale_merge/RotationScaleMerge.h"
#include "../image_analysis/scale_merge/ResolutionCutoff.h"
#include "../image_analysis/scale_merge/ReindexAmbiguity.h"
#include "../image_analysis/scale_merge/ScalingResult.h"
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
#include "../image_analysis/geom_refinement/PostRefine.h"
#include "../image_analysis/lattice_search/LatticeSearch.h"
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
#include "../image_analysis/scale_merge/HKLKey.h"
#include "../image_analysis/scale_merge/ScaleOnTheFly.h"
#include "../image_analysis/scale_merge/StillsPartialityRefine.h"
#include "../image_analysis/WriteReflections.h"
#include "../image_analysis/bragg_integration/CalcISigma.h"
#include "../common/Definitions.h"
#include "../common/CorrelationCoefficient.h"
#include <array>
#include <map>
namespace {
// How precisely the beam centre has to be placed before the estimate is used at all. An estimate
// that fails it is treated as no estimate: the spot symmetry falls through to the background and
// the background to the header. A fit's own scatter will not serve as the gate - the ways either
// estimator goes badly wrong all leave a small one - so the spot symmetry reports instead how
// far its answer moves when the search is started elsewhere, which does separate them. On the 38
// rotation regression crystals that stays below 0.35 px against errors of at most 0.74; over a
// ring of injected geometries it rises past 1.0 px on 96 % of the answers that are 2 px or more
// out, and on none of the 38 as they are.
constexpr float MAX_BEAM_CENTER_SIGMA_PXL = 1.0f;
// Images the projection is built from when the beam centre is wanted but the beam-stop pre-pass
// is off; with it on, that pre-pass's image count is used and one projection serves both.
constexpr int BEAM_CENTER_PROJECTION_IMAGES = 60;
// Spots per pre-scan frame kept for the beam centre, the strongest first.
constexpr size_t BEAM_CENTER_SPOTS_PER_IMAGE = 500;
// Frames read for the spot symmetry alone when the first sample could not measure it. Both
// estimators are counting statistics over PAIRS of frames - the Friedel one over frames half a
// turn apart, the second crossing over any two frames a reflection is recorded on - so a sweep
// whose frames carry few spots can only be measured by reading more of them. On the two
// rotation regression crystals that decline, this is what turns the decline into an answer
// 0.04 px from the post-refined centre; the sample the beam stop is built from is left alone.
constexpr int BEAM_CENTER_SPARSE_IMAGES = 400;
// Frames the spot symmetry reads, as a multiple of the beam stop's. Its two estimators are
// counts of matched PAIRS and the pairs grow faster than the frames do, so this is the one knob
// that buys independence from WHICH frames were drawn - and that dependence is the estimator's
// largest remaining failure. Asking the same code for one frame more or less re-draws the whole
// sample: over five draws of 38 sweeps, 60 frames leaves three answers more than 1 px out, none
// of which reports an unusual sigma, and 120 leaves none in 186. It is also what makes the end
// margin below free: at 60 frames excluding the ends moved three centres, at 120 it moves none.
constexpr int BEAM_CENTER_IMAGE_FACTOR = 2;
// Images left out at each end of the sweep. Shutter synchronisation disturbs the first and last
// frame or two of a run, and a partly exposed frame is dark over the whole detector - which is
// exactly the kind of error a projection cannot average away and a background comparison reads
// as signal. The margin is capped at a tenth of the sweep so a short run still has a sample.
constexpr int PRESCAN_END_MARGIN_IMAGES = 5;
// Workers reading the pre-scan sample. Each owns a shard of the beam-stop projection so no two
// threads touch the same accumulator, and a shard costs 20 bytes per pixel - 362 MB on a 16M
// detector - so this is capped well below the worker count of the run proper. The accumulation
// is memory-bound rather than compute-bound, so a handful of workers already saturates it.
constexpr size_t PRESCAN_MAX_WORKERS = 8;
// Pick up to requested_images ordinals spread evenly across [0, images_to_process) for the
// first pass of two-pass rotation indexing.
std::vector<int> select_equally_spaced_image_ordinals(int images_to_process, int requested_images) {
std::vector<int> ret;
if (images_to_process <= 0 || requested_images <= 0)
return ret;
const int n = std::min(images_to_process, requested_images);
if (n == 1) {
ret.push_back(0);
return ret;
}
std::set<int> unique_ordinals;
for (int i = 0; i < n; i++)
unique_ordinals.insert(static_cast<int>(
std::llround(static_cast<double>(i) * static_cast<double>(images_to_process - 1) /
static_cast<double>(n - 1))));
ret.assign(unique_ordinals.begin(), unique_ordinals.end());
return ret;
}
// A sample arranged as pairs of images half a turn apart: half the ordinals spread over the
// part of the sweep that has a partner 180 deg later, and those partners. That is what the
// Friedel mirror needs. It is NOT the plain selection rearranged - even on a full turn the two
// sets are all but disjoint, sharing 2 of 60 ordinals on a 3600 image sweep - so the caller
// keeps them apart and gives each consumer the one it needs. Empty when the sweep is shorter
// than half a turn, which is the caller's signal that the spot symmetry cannot be measured here.
std::vector<int> select_half_turn_paired_ordinals(int images_to_process, int requested_images,
float degrees_per_ordinal) {
if (!(degrees_per_ordinal > 0.0f))
return {};
const int half_turn = static_cast<int>(std::lround(180.0 / degrees_per_ordinal));
if (half_turn < 1 || images_to_process - half_turn < 1)
return {};
std::set<int> ordinals;
for (const int ordinal : select_equally_spaced_image_ordinals(images_to_process - half_turn,
requested_images / 2)) {
ordinals.insert(ordinal);
ordinals.insert(ordinal + half_turn);
}
return {ordinals.begin(), ordinals.end()};
}
// Apply a goniometer rotation SCALE: the stage turned k times the angle the file records. What is
// wrong is the SWEEP, not where it began, so the per-frame increment and the per-frame oscillation
// width both scale by k while the starting angle is left alone - the stage reached that position
// before the sweep started. A constant offset in phi is in any case exactly degenerate with a
// rotation of the crystal orientation about the spindle, which indexing refines away, so anchoring
// the stretched sweep at image 0 costs nothing and keeps the first frame at the angle the file gives
// it. Every angle in the pipeline comes from this object, so scaling it is the whole correction.
GoniometerAxis ScaleRotation(const GoniometerAxis &g, float k) {
GoniometerAxis scaled(g.GetName(), g.GetStart_deg(), g.GetIncrement_deg() * k,
g.GetAxis(), g.GetHelicalStep());
if (const auto wedge = g.GetScreeningWedge())
scaled.ScreeningWedge(*wedge * k);
return scaled;
}
}
Rugnux::Rugnux(JFJochHDF5Reader &reader, DiffractionExperiment experiment,
PixelMask pixel_mask, ProcessConfig config)
: reader_(reader), experiment_(std::move(experiment)),
pixel_mask_(std::move(pixel_mask)), config_(std::move(config)) {
// Bit 9 describes where THIS run found the beam stop, so a mask read back from a file that
// already carries one starts clear; the user mask (bit 8) is left as it was loaded.
pixel_mask_.ClearBeamStopMask(experiment_);
// A manually asserted stage calibration applies to everything, before anything reads an angle.
if (config_.rotation_scale.has_value())
if (const auto g = experiment_.GetGoniometer()) {
Logger("Rugnux").Info("Goniometer rotation scale {:.5f} applied from the command line",
*config_.rotation_scale);
experiment_.Goniometer(ScaleRotation(*g, *config_.rotation_scale));
}
}
// Report the detector geometry the way XDS states it, so a run can be compared line by line against
// the IDXREF.LP / CORRECT.LP of XDS on the same data. XDS is never given this geometry - the durin
// plugin hands it image data only (plugin_get_header returns dimensions, bytes per pixel, pixel size
// and frame count, nothing else) and XDS refines its own from XDS.INP - which is precisely why
// having ours in the same convention is worth printing.
//
// The two laboratory frames coincide: x along increasing detector column, y along increasing row,
// z along the beam. XDS places a pixel at
// x_lab(i,j) = (i-ORGX)*QX*X_axis + (j-ORGY)*QY*Y_axis + DETECTOR_DISTANCE*(X_axis x Y_axis)
// and DiffractionGeometry::LabCoord is that same expression with X_axis = poni_rot*(1,0,0) and
// Y_axis = poni_rot*(0,1,0). So the tilt goes out as the two axis vectors rather than as angles,
// which is also how XDS reports it after refinement. Taking the axes as differences of LabCoord
// keeps them in step with whatever the geometry currently is, tilt included.
//
// ORGX/ORGY are 1-based - XDS counts pixels from 1 - and are the foot of the perpendicular from the
// crystal rather than the direct beam, which is what our beam centre is too
// (docs/DETECTOR_GEOMETRY.md), so they map across without correction.
//
// ROTATION_AXIS is printed as stored. Its direction is right, the frames being shared; its SIGN has
// not been cross-checked against an XDS refinement, so treat a flip there as unconfirmed rather than
// as a discrepancy in the data.
static void LogXDSGeometry(const DiffractionExperiment &experiment) {
Logger logger("Rugnux");
const auto geom = experiment.GetDiffractionGeometry();
const float beam_x = geom.GetBeamX_pxl();
const float beam_y = geom.GetBeamY_pxl();
const float pixel_mm = geom.GetPixelSize_mm();
const Coord origin = geom.LabCoord(beam_x, beam_y);
const Coord x_axis = ((geom.LabCoord(beam_x + 1.0f, beam_y) - origin) / pixel_mm).Normalize();
const Coord y_axis = ((geom.LabCoord(beam_x, beam_y + 1.0f) - origin) / pixel_mm).Normalize();
logger.Info("XDS convention: ORGX= {:.2f} ORGY= {:.2f} DETECTOR_DISTANCE= {:.4f}",
beam_x + 1.0f, beam_y + 1.0f, geom.GetDetectorDistance_mm());
logger.Info("XDS convention: DIRECTION_OF_DETECTOR_X-AXIS= {:.6f} {:.6f} {:.6f}",
x_axis.x, x_axis.y, x_axis.z);
logger.Info("XDS convention: DIRECTION_OF_DETECTOR_Y-AXIS= {:.6f} {:.6f} {:.6f}",
y_axis.x, y_axis.y, y_axis.z);
logger.Info("XDS convention: INCIDENT_BEAM_DIRECTION= 0 0 1 X-RAY_WAVELENGTH= {:.6f} QX= QY= {:.6f}",
geom.GetWavelength_A(), pixel_mm);
if (const auto goniometer = experiment.GetGoniometer())
logger.Info("XDS convention: ROTATION_AXIS= {:.6f} {:.6f} {:.6f}",
goniometer->GetAxis().x, goniometer->GetAxis().y, goniometer->GetAxis().z);
}
void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, RugnuxObserver *observer) {
Logger logger("Rugnux");
// Two-pass rotation runs this twice. The shadow does not move, and re-detecting would find
// nothing (its pixels are masked by now) and so clear the mask the first pass established;
// the beam centre must not be re-read either, because by then the geometry is post-refined.
const auto &current = pixel_mask_.GetMask();
const bool have_shadow = std::any_of(current.begin(), current.end(),
[](uint32_t v) { return (v & (1u << PixelMask::BeamStopPixelBit)) != 0; });
const bool want_shadow = config_.detect_beam_stop.has_value() && !have_shadow;
const bool want_beam_center = config_.estimate_beam_center && !beam_center_placed_;
if (!want_shadow && !want_beam_center)
return;
// The two consumers want different frames and each gets its own set. The shadow is built from
// the even spread over the sweep it has always been built from, so that asking for a beam
// centre cannot move the mask; the spot symmetry needs frames paired half a turn apart, where
// the sweep reaches that far. The two sets barely overlap, so the loop below walks their union
// and every frame is read once.
//
// Both samples are taken from the middle of the sweep, both ends left out: an even spread over
// [0, images_to_process) always picks the very first and the very last image, and those are the
// ones a shutter out of step spoils - the last carries under 90 % of the sample's median counts
// on 12 of 38 regression sweeps and is essentially unexposed on 5.
const auto goniometer = experiment_.GetGoniometer();
const int margin = std::min(PRESCAN_END_MARGIN_IMAGES, images_to_process / 10);
auto shadow_sample = select_equally_spaced_image_ordinals(images_to_process - 2 * margin,
frame_count);
auto spot_sample = want_beam_center && goniometer
? select_half_turn_paired_ordinals(
images_to_process - 2 * margin,
BEAM_CENTER_IMAGE_FACTOR * frame_count,
// Degrees per ORDINAL: the shifted goniometer's increment already
// is that, and rotation data cannot be strided anyway.
std::abs(goniometer->GetIncrement_deg()))
: std::vector<int>();
for (int &ordinal : shadow_sample)
ordinal += margin;
for (int &ordinal : spot_sample)
ordinal += margin;
const bool want_spot_symmetry = !spot_sample.empty();
const std::set<int> shadow_set(shadow_sample.begin(), shadow_sample.end());
const std::set<int> spot_set(spot_sample.begin(), spot_sample.end());
std::set<int> sample = shadow_set;
sample.insert(spot_set.begin(), spot_set.end());
if (sample.empty())
return;
// Say which of the two this is reading for. It runs before the first image of the run proper is
// processed, so without this the job sits at its initial status for the whole pre-scan - which is
// tens of seconds once the beam centre asks for its own frames on top of the projection.
if (observer)
observer->OnPhase(want_shadow ? (want_beam_center ? "Beam stop and beam centre" : "Beam stop")
: "Beam centre");
ShadowFinder finder(experiment_, pixel_mask_);
// The spot symmetry is read in the same pass as the projection rather than in a pre-pass of its
// own. The image is preprocessed a second time for it: the shadow is accumulated from raw counts
// and the spot finder works on the converted image.
std::unique_ptr<AzimuthalIntegrationMapping> prescan_mapping;
if (want_spot_symmetry)
prescan_mapping = std::make_unique<AzimuthalIntegrationMapping>(experiment_, pixel_mask_);
// Everything reading one frame needs that cannot be shared between threads: its own scratch for
// each of the two decompressions, its own preprocessed image and its own spot finder.
struct PreScanWorker {
std::vector<uint8_t> shadow_buffer;
std::vector<uint8_t> decompression_buffer;
std::unique_ptr<ImagePreprocessorCPU> preprocessor;
std::unique_ptr<ImagePreprocessorBuffer> preprocessed;
std::unique_ptr<ImageSpotFinder> spot_finder;
};
const auto make_worker = [&] {
PreScanWorker w;
if (want_spot_symmetry) {
w.preprocessor = std::make_unique<ImagePreprocessorCPU>(experiment_, pixel_mask_);
w.preprocessed = std::make_unique<ImagePreprocessorBuffer>(experiment_.GetPixelsNum());
if (config_.spot_finding.adaptive_threshold)
w.spot_finder = std::make_unique<AdaptiveSpotFinderCPU>(*prescan_mapping);
else
w.spot_finder = std::make_unique<ImageSpotFinderCPU>(experiment_.GetXPixelsNumConv(),
experiment_.GetYPixelsNumConv());
}
return w;
};
// The spots of one image. Taken as a step of its own because the sparse-sweep pass below reads
// images the shadow does not. The frame they belong to is stamped by the caller, which is what
// keeps the pool independent of the order the workers finish in.
const auto find_spots = [&](PreScanWorker &w, CompressedImage &image, int image_idx,
std::vector<BeamCenterSpot> &out) {
try {
w.preprocessor->Analyze(*w.preprocessed,
image.GetUncompressedPtr(w.decompression_buffer), image.GetMode());
} catch (const std::exception &e) {
logger.Warning("Pre-scan: failed to preprocess image {}: {}", image_idx, e.what());
return;
}
auto spots = w.spot_finder->Run(*w.preprocessed, config_.spot_finding);
// The strongest of a crowded frame: the symmetry is over-determined either way, and the
// matching is quadratic in the spots of one frame.
if (spots.size() > BEAM_CENTER_SPOTS_PER_IMAGE) {
std::partial_sort(spots.begin(), spots.begin() + BEAM_CENTER_SPOTS_PER_IMAGE, spots.end(),
[](const DiffractionSpot &a, const DiffractionSpot &b) {
return a.Count() > b.Count(); });
spots.resize(BEAM_CENTER_SPOTS_PER_IMAGE);
}
for (const auto &spot : spots) {
const Coord centroid = spot.RawCoord();
out.push_back({centroid.x, centroid.y, static_cast<float>(spot.Count()), 0});
}
};
std::vector<float> frame_angle_deg;
std::vector<BeamCenterSpot> beam_center_spots;
// Read the sample on several workers. The reader serialises on the HDF5 lock, but the
// decompression, the projection and the spot finding - which is all of the cost on a large
// detector - run in parallel. Each worker accumulates into a shard of its own, so nothing is
// locked while an image is added, and the per-frame results are stitched together in sample
// order below so the beam centre sees the same input however the workers interleaved.
{
const std::vector<int> ordinals(sample.begin(), sample.end());
std::vector<std::vector<BeamCenterSpot>> spots_of(ordinals.size());
// A frame joins the pool only if it could be read, as it did when this was a serial loop.
std::vector<char> spot_read(ordinals.size(), 0);
const size_t nworkers = std::min<size_t>(std::max<size_t>(config_.nthreads, 1),
std::min(PRESCAN_MAX_WORKERS, ordinals.size()));
finder.SetShardCount(nworkers);
std::atomic<size_t> next{0};
std::vector<std::future<void>> futures;
futures.reserve(nworkers);
for (size_t t = 0; t < nworkers; t++)
futures.emplace_back(std::async(std::launch::async, [&, t] {
PreScanWorker w = make_worker();
for (size_t i = next.fetch_add(1); i < ordinals.size(); i = next.fetch_add(1)) {
const int ordinal = ordinals[i];
const int image_idx = start_image + ordinal * config_.stride;
std::shared_ptr<JFJochReaderRawImage> img;
try {
img = reader_.GetRawImage(image_idx);
} catch (const std::exception &e) {
logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what());
continue;
}
if (!img) continue;
DataMessage msg{};
msg.image = img->image;
msg.number = ordinal;
msg.original_number = image_idx;
if (shadow_set.contains(ordinal))
finder.AddImage(msg, w.shadow_buffer, t);
if (spot_set.contains(ordinal)) {
spot_read[i] = 1;
find_spots(w, msg.image, image_idx, spots_of[i]);
}
}
}));
for (auto &f : futures) f.get();
// Frame numbering follows the sample order, exactly as the serial read did.
for (size_t i = 0; i < ordinals.size(); i++) {
if (!spot_read[i]) continue;
const int frame = static_cast<int>(frame_angle_deg.size());
// The goniometer has been shifted so that a LOCAL ordinal maps to the angle of the
// original image it came from (see RunPipeline), so the angle is asked for by ordinal,
// as everywhere else.
frame_angle_deg.push_back(goniometer->GetAngle_deg(static_cast<float>(ordinals[i])));
for (auto &spot : spots_of[i]) {
spot.frame = frame;
beam_center_spots.push_back(spot);
}
}
}
if (finder.GetFrameCount() == 0) {
logger.Warning("Pre-scan: no image could be read. Skipping.");
return;
}
if (want_shadow) {
const auto shadow = finder.GetMask(config_.nthreads);
const auto shadowed = std::count(shadow.begin(), shadow.end(), 1u);
pixel_mask_.LoadBeamStopMask(experiment_, shadow);
logger.Info("Beam stop shadow: {} pixels ({:.2f}% of the detector) found in {} images",
shadowed, 100.0 * static_cast<double>(shadowed) / static_cast<double>(shadow.size()),
finder.GetFrameCount());
}
// The beam centre comes after, so that the shadow is in the mask by the time it is placed and
// the holder arm - the largest azimuthal asymmetry on the detector - is out of the way of the
// background the centre is read from.
if (!want_beam_center)
return;
beam_center_placed_ = true;
// The spot symmetry first where the sweep supports it: it is exact geometry rather than a
// property of the scattering, and it measures the sets the background cannot - a flat radial
// profile at long wavelength leaves that fit nothing to bite on. The background answers
// wherever the sweep does not reach half a turn, which is every still and every screening
// wedge, so the two cover between them what neither covers alone.
std::string source = "spot symmetry";
std::optional<BeamCenterEstimate> estimate;
SpindleEstimate spindle;
if (want_spot_symmetry)
estimate = FindBeamCenterFromSpotSymmetry(experiment_, frame_angle_deg, beam_center_spots,
config_.fit_spindle ? &spindle : nullptr);
// Where it does not come out - no answer at all, or one that moves when the search is started
// elsewhere - read more of the sweep and ask again. Both estimators are counts of matched spot
// PAIRS, the Friedel one over frames half a turn apart and the second crossing over any two
// frames one reflection is recorded on, so a sweep whose frames carry few spots, or whose pairs
// match rarely, is measurable only from more of them. More data first and a weaker estimator
// second: this runs before the fall-through, and only what survives it falls through. The extra
// frames go to the spot finder alone, so the beam-stop projection keeps the images it was
// validated on.
if (want_spot_symmetry && (!estimate || estimate->sigma_pxl > MAX_BEAM_CENTER_SIGMA_PXL)) {
// What has already been read FOR THE SPOTS, which is not the whole sample: the frames the
// shadow was built from carry none of them, and testing against the union would silently
// skip the ones this pass exists to read.
std::vector<int> extra;
for (int ordinal : select_half_turn_paired_ordinals(
images_to_process - 2 * margin, BEAM_CENTER_SPARSE_IMAGES,
std::abs(goniometer->GetIncrement_deg()))) {
ordinal += margin;
if (!spot_set.contains(ordinal))
extra.push_back(ordinal);
}
if (!extra.empty()) {
logger.Info("Beam centre: the spot symmetry does not come out on {} frames "
"({} spots); reading {} more frames",
frame_angle_deg.size(), beam_center_spots.size(), extra.size());
// Its own phase: this reads several times the frames the pass above did, and is the one
// step of the pre-scan a user is likely to sit through wondering what is happening.
if (observer)
observer->OnPhase("Beam centre (reading more frames)");
// Read on workers as above - this pass is several times the size of the first one - and
// append in `extra` order so the pool does not depend on how they interleaved.
std::vector<std::vector<BeamCenterSpot>> extra_spots(extra.size());
std::vector<char> extra_read(extra.size(), 0);
const size_t nworkers = std::min<size_t>(std::max<size_t>(config_.nthreads, 1),
std::min(PRESCAN_MAX_WORKERS, extra.size()));
std::atomic<size_t> next{0};
std::vector<std::future<void>> futures;
futures.reserve(nworkers);
for (size_t t = 0; t < nworkers; t++)
futures.emplace_back(std::async(std::launch::async, [&] {
PreScanWorker w = make_worker();
for (size_t i = next.fetch_add(1); i < extra.size(); i = next.fetch_add(1)) {
const int image_idx = start_image + extra[i] * config_.stride;
std::shared_ptr<JFJochReaderRawImage> img;
try {
img = reader_.GetRawImage(image_idx);
} catch (const std::exception &e) {
logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what());
continue;
}
if (img) {
extra_read[i] = 1;
find_spots(w, img->image, image_idx, extra_spots[i]);
}
}
}));
for (auto &f : futures) f.get();
for (size_t i = 0; i < extra.size(); i++) {
if (!extra_read[i]) continue;
const int frame = static_cast<int>(frame_angle_deg.size());
frame_angle_deg.push_back(goniometer->GetAngle_deg(static_cast<float>(extra[i])));
for (auto &spot : extra_spots[i]) {
spot.frame = frame;
beam_center_spots.push_back(spot);
}
}
estimate = FindBeamCenterFromSpotSymmetry(experiment_, frame_angle_deg,
beam_center_spots,
config_.fit_spindle ? &spindle : nullptr);
}
}
if (config_.fit_spindle && std::isfinite(spindle.azimuth_rad))
logger.Info("Spindle: azimuth about the beam {:.3f} mrad, tip towards it {:.3f} mrad "
"(Friedel vote {:.0f} -> {:.0f})",
1e3 * spindle.azimuth_rad, 1e3 * spindle.tip_rad,
spindle.vote_excess_nominal, spindle.vote_excess);
if (!estimate || estimate->sigma_pxl > MAX_BEAM_CENTER_SIGMA_PXL) {
source = "background";
estimate = FindBeamCenterFromBackground(experiment_, pixel_mask_, finder.GetMeanProjection());
}
if (!estimate) {
logger.Info("Beam centre: not measurable, keeping ({:.2f},{:.2f})",
experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl());
return;
}
const float moved = std::hypot(estimate->beam_x_pxl - experiment_.GetBeamX_pxl(),
estimate->beam_y_pxl - experiment_.GetBeamY_pxl());
const bool commit = estimate->sigma_pxl <= MAX_BEAM_CENTER_SIGMA_PXL;
logger.Info("Beam centre from {}: ({:.2f},{:.2f}) -> ({:.2f},{:.2f}), moved {:.2f} px, "
"sigma {:.2f} px => {}",
source, experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(),
estimate->beam_x_pxl, estimate->beam_y_pxl, moved, estimate->sigma_pxl,
commit ? "COMMIT" : "reject (kept header)");
if (commit)
experiment_.BeamX_pxl(estimate->beam_x_pxl).BeamY_pxl(estimate->beam_y_pxl);
}
void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_process,
RugnuxObserver *observer) {
Logger logger("Rugnux");
// Rotation has its own two-pass; this is stills only.
if (experiment_.IsRotationIndexing()) {
logger.Warning("--refine-geometry is stills-only; rotation data uses two-pass indexing. Ignoring.");
return;
}
// The bundle anchors a known cell; without one there is nothing to refine against.
const auto cell = experiment_.GetUnitCell();
if (!cell.has_value()) {
logger.Warning("--refine-geometry needs a reference unit cell (-C / -S / reference MTZ). Skipping.");
return;
}
const int refine_frames = std::max(1, config_.refine_geometry.value());
if (observer)
observer->OnPhase("Geometry refinement (first pass)");
const auto dataset = reader_.GetDataset();
// Index a spread sample large enough to yield ~refine_frames strong frames even at a low hit rate.
const int sample_budget = std::min(images_to_process, std::max(refine_frames * 50, 8000));
const std::vector<int> sample = select_equally_spaced_image_ordinals(images_to_process, sample_budget);
// First-pass engines, built from the current (nominal) geometry: an index-only pass (no
// integration) to obtain each frame's spots + assigned HKL + orientation.
AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_);
IndexerThreadPool pool(experiment_.GetIndexingSettings(), IndexerConstruction::OnFirstUse);
IndexAndRefine indexer(experiment_, &pool, /*retain_outcomes=*/false);
auto pass_settings = config_.spot_finding;
pass_settings.enable = true;
pass_settings.indexing = true;
pass_settings.quick_integration = false;
// A frame is worth bundling only if it indexed this many spots; the bundle then takes the
// refine_frames strongest of those.
constexpr int MIN_STRONG_SPOTS = 20;
// Stop sampling once there are comfortably more strong frames than the bundle will use, instead of
// draining the whole budget. The budget is sized for a low-hit-rate dataset (refine_frames * 50, at
// least 8000), so on data that indexes well it meant re-indexing the ENTIRE run to keep 200 frames -
// 99% of this pass was sampling. Four times the bundle leaves the "strongest N" selection a real
// pool to choose from.
const int strong_target = refine_frames * 4;
// The sample is cut into a FIXED number of interleaved stripes - stripe s is sample positions
// s, s+STRIPES, s+2*STRIPES, ... - and each stripe stops once it has contributed its share of the
// strong frames. Two things follow, both of which the previous shared-cursor-with-a-shared-stop
// arrangement got wrong. Every stripe is spread over the WHOLE run, so stopping early no longer
// means fitting the geometry to the beginning of it. And which frames get examined depends only
// on the data: it is not a race between the workers, and it does not change with -N, because a
// stripe is processed identically whichever worker happens to claim it.
constexpr int STRIPES = 32;
const int stripe_strong_target = std::max(1, (strong_target + STRIPES - 1) / STRIPES);
std::vector<GeomRefineFrame> frames;
std::mutex frames_mutex;
std::atomic<int> next_stripe = 0;
std::atomic<int> examined = 0;
auto worker = [&]() {
pin_gpu(); // round-robin per worker thread; must precede engine construction
MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, indexer,
/*enable_fused_adaptive_gpu=*/true);
AzimuthalIntegrationProfile profile(mapping);
for (int stripe = next_stripe.fetch_add(1); stripe < STRIPES && !cancelled_;
stripe = next_stripe.fetch_add(1)) {
int stripe_strong = 0;
for (int idx = stripe; idx < static_cast<int>(sample.size()) && !cancelled_
&& stripe_strong < stripe_strong_target; idx += STRIPES) {
examined.fetch_add(1, std::memory_order_relaxed);
const int ordinal = sample[idx];
const int image_idx = start_image + ordinal * config_.stride;
std::shared_ptr<JFJochReaderRawImage> img;
try {
img = reader_.GetRawImage(image_idx);
} catch (const std::exception &e) {
logger.Warning("Geometry refinement: failed to load image {}: {}", image_idx, e.what());
continue;
}
if (!img) continue;
DataMessage msg{};
msg.image = img->image;
msg.number = ordinal;
msg.original_number = image_idx;
if (dataset->efficiency.size() > image_idx)
msg.image_collection_efficiency = dataset->efficiency[image_idx];
try {
analysis.Analyze(msg, profile, pass_settings);
} catch (const std::exception &e) {
continue;
}
if (!msg.indexing_result.value_or(false) || !msg.indexing_lattice.has_value())
continue;
GeomRefineFrame f;
f.lattice = *msg.indexing_lattice;
f.ordinal = ordinal;
for (const auto &s : msg.spots)
if (s.indexed && s.lattice == 0)
f.spots.push_back(GeomRefineSpot{s.x, s.y,
static_cast<int32_t>(s.h), static_cast<int32_t>(s.k), static_cast<int32_t>(s.l)});
if (f.spots.size() >= 6) {
if (static_cast<int>(f.spots.size()) >= MIN_STRONG_SPOTS)
stripe_strong++;
const std::unique_lock ul(frames_mutex);
frames.push_back(std::move(f));
}
}
}
};
std::vector<std::future<void>> futures;
futures.reserve(config_.nthreads);
for (int i = 0; i < config_.nthreads; ++i)
futures.push_back(std::async(std::launch::async, worker));
for (auto &fut : futures)
fut.get();
if (cancelled_)
return;
// Keep the strongest frames (most indexed spots) for the bundle - the true cell indexes many
// spots per frame, and orientation diversity comes for free from independent serial stills.
std::vector<GeomRefineFrame> strong;
for (auto &f : frames)
if (static_cast<int>(f.spots.size()) >= MIN_STRONG_SPOTS)
strong.push_back(std::move(f));
// frames is in the order the workers happened to finish, so the spot count alone does not order
// it: without the ordinal tie-break, equally strong frames would swap places between runs and a
// different bundle would be refined.
std::sort(strong.begin(), strong.end(),
[](const GeomRefineFrame &a, const GeomRefineFrame &b) {
if (a.spots.size() != b.spots.size()) return a.spots.size() > b.spots.size();
return a.ordinal < b.ordinal;
});
if (static_cast<int>(strong.size()) > refine_frames)
strong.resize(refine_frames);
logger.Info("Geometry refinement: indexed {} of {} sampled frames (budget {}), bundling {} strong "
"frames (>= {} spots)",
frames.size(), std::min<int>(examined.load(), static_cast<int>(sample.size())),
sample.size(), strong.size(), MIN_STRONG_SPOTS);
GeometryRefinerSettings gr_settings;
gr_settings.crystal_system = experiment_.GetCrystalSystem();
gr_settings.num_threads = config_.nthreads;
const GeometryRefinerResult r = RefineGlobalGeometry(
experiment_.GetDiffractionGeometry(), *cell, strong, gr_settings);
if (!r.ok) {
logger.Warning("Geometry refinement did not run (too few strong frames / did not converge); "
"keeping the input geometry");
return;
}
logger.Info("Geometry refinement: beam ({:.2f}, {:.2f}) -> ({:.2f}, {:.2f}) px, distance {:.4f} -> {:.4f} mm",
experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), r.beam_x_px, r.beam_y_px,
experiment_.GetDetectorDistance_mm(), r.distance_mm);
logger.Info("Geometry refinement: cell a/b/c {:.3f}/{:.3f}/{:.3f} -> {:.3f}/{:.3f}/{:.3f} A "
"(median residual {:.3f} px, {} frames, {} spots)",
cell->a, cell->b, cell->c, r.cell.a, r.cell.b, r.cell.c,
r.median_residual_px, r.frames_used, r.spots_used);
experiment_.BeamX_pxl(r.beam_x_px).BeamY_pxl(r.beam_y_px).DetectorDistance_mm(r.distance_mm);
experiment_.SetUnitCell(r.cell);
// The beam centre now comes from indexed spots, which the background cannot better.
beam_center_placed_ = true;
}
ProcessResult Rugnux::Run(RugnuxObserver *observer) {
// Rotation two-pass geometry post-refinement: the first pass integrates at the header geometry and
// post-refines the detector geometry (distance + beam from the observed spot positions, cell scale + axis
// from phi_obs); the second pass re-indexes and re-integrates with the refined geometry. BOTH passes are
// written out - the second (refined-geometry) pass is the CANONICAL result and takes the plain "<prefix>_*"
// name, while the first (header-geometry) pass is kept as "<prefix>_01_*" for comparison - so the two can be
// compared and the ice-ring flags / geometry stored in each master reflect the geometry that pass used.
if (config_.mode == ProcessMode::FullAnalysis && config_.rotation_postrefine_geometry
&& experiment_.IsRotationIndexing()) {
Logger logger("Rugnux");
const std::string base_prefix = config_.output_prefix;
const auto gonio_snapshot = experiment_.GetGoniometer();
prepass_detector_geometry_.reset();
prepass_rotation_scale_.reset();
prepass_result_.reset();
force_rotation_result_.reset();
prepass_merge_sg_.reset();
prepass_promoted_point_group_ = false;
logger.Info("Rotation two-pass geometry post-refinement: first pass (header geometry) -> {}_01_*",
base_prefix);
// Keep an empty prefix empty (the "compute stats, persist nothing" mode): the "_01" suffix would
// make it non-empty and RunPipeline would then write stray _01_* files despite no output wanted.
config_.output_prefix = base_prefix.empty() ? base_prefix : base_prefix + "_01";
auto pass1 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/true);
pass1.pass_number = 1;
pass1.pass_count = 2;
if (cancelled_) { config_.output_prefix = base_prefix; return pass1; }
if (gonio_snapshot) experiment_.Goniometer(*gonio_snapshot); // undo the pre-pass goniometer shift
// Apply the post-refined detector geometry for the second pass, keeping the header geometry so
// the run can go back to it if the refined pass turns out worse (see the quality guard below).
const std::array<float, 3> header_geometry{experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(),
experiment_.GetDetectorDistance_mm()};
if (prepass_detector_geometry_) {
const auto &g = *prepass_detector_geometry_;
experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]);
}
// ... and the goniometer rotation scale, on the same measure-then-re-integrate footing: the angles
// in the file are the commanded ones, so a stage that ran fast is a geometry error like any other.
if (prepass_rotation_scale_ && gonio_snapshot) {
experiment_.Goniometer(ScaleRotation(*gonio_snapshot, *prepass_rotation_scale_));
// The pre-pass mosaicity is a width in degrees fitted against the angles the second pass has
// just stopped using, and the override can only ever raise the second pass's own estimate (it
// takes the larger of the two). Carrying it over would hold the second pass at the rocking
// width the uncorrected angles produced - the correction half-applied. Drop it and let the
// second pass fit its own.
prepass_mosaicity_.clear();
logger.Info("Two-pass: goniometer rotation scale {:.5f} will drive the second integration pass "
"(oscillation {:.4f} -> {:.4f} deg per image)", *prepass_rotation_scale_,
gonio_snapshot->GetIncrement_deg(),
gonio_snapshot->GetIncrement_deg() * *prepass_rotation_scale_);
}
// Space group: the second pass RE-INDEXES DE NOVO (clear the group here) so the indexer's pseudo-
// symmetry safeguards recover the true cell - reusing pass-1's group in the indexer forces build_sr
// onto a wrong / doubled cell (a huge oblique cell collapses; a pseudo-centred cell doubles). Pass-1's
// group is instead reinstated for the MERGE ONLY (RunPipeline, before RotationScaleMerge), so the
// symmetry stays fixed with no re-search / no flip.
//
// This holds for a group the centred-lattice test determined too, even though its conventional
// setting is not the one the de-novo primitive frame comes back in. That is precisely what the
// primitive->conventional reindex below does, and when the reindex declines - pass 2's metric is
// not the group's - the centring check right after it stops the pass being adopted at all. Letting
// such a group be re-searched instead put the second pass's determination back in play at a
// slightly different geometry, which on a lattice whose centring is pseudo-symmetric to a couple of
// tenths of a percent is a coin toss: measured on a C-centred monoclinic crystal, pass 1 confirmed
// the centring and pass 2 dropped it, and the run merged in P1.
experiment_.SpaceGroupNumber(std::nullopt);
prepass_merge_sg_ = pass1.space_group_number;
prepass_promoted_point_group_ = pass1.twinning.laue_class_was_chosen_by_promotion;
logger.Info("Rotation two-pass geometry post-refinement: second pass (refined geometry, canonical) -> {}_*",
base_prefix);
config_.output_prefix = base_prefix; // the refined pass is the canonical result (no _02 suffix)
auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
// The post-refinement is measured by pass 1 and consumed by pass 2, so carry it onto whichever
// result is returned - it is a result of the run, not of the pass that happened to fit it.
pass2.post_refine = pass1.post_refine;
pass2.pass_number = 2;
pass2.pass_count = 2;
// Supercell-collapse guard for the second pass. The de-novo lattice search can be BISTABLE on a
// periodically-modulated dataset: the tiny post-refine geometry move tips it out of the true-cell
// basin onto a spurious SUPERCELL (a large multiple of pass-1's cell volume), and the space group
// reused from pass-1 is then stamped onto that wrong-shaped cell - catastrophic. Pass 1 (header
// geometry, un-perturbed) found the correct lattice, so re-run the second pass with pass-1's FULL
// indexing result forced (its refined lattice + orientation + rotation axis + geometry - not just the
// cell, which loses the per-frame orientation and integrates worse): the re-run then matches de-novo
// quality but at the refined geometry. Only fires when pass 2 is markedly LARGER than pass 1 - a
// genuine de-novo demotion to a smaller primitive (the pseudo-symmetry case the second pass re-indexes
// de novo to catch) has pass2 <= pass1 and is left untouched.
if (!cancelled_ && prepass_result_ && pass1.consensus_cell && pass2.consensus_cell) {
const double v1 = gemmi::UnitCell(*pass1.consensus_cell).volume;
const double v2 = gemmi::UnitCell(*pass2.consensus_cell).volume;
if (v1 > 0.0 && v2 > 1.5 * v1) {
logger.Info("Two-pass: second-pass cell volume {:.0f} A^3 is {:.2f}x pass-1 ({:.0f} A^3) - a "
"spurious supercell; re-running the second pass with pass-1's result forced",
v2, v2 / v1, v1);
force_rotation_result_ = *prepass_result_;
// Pass 1's result carries pass 1's goniometer, and the per-image path takes its angles
// from there - forcing it whole would put the uncorrected angles back and silently undo
// the rotation scale for the re-run.
if (prepass_rotation_scale_ && force_rotation_result_->axis)
force_rotation_result_->axis =
ScaleRotation(*force_rotation_result_->axis, *prepass_rotation_scale_);
pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
pass2.post_refine = pass1.post_refine;
pass2.pass_number = 3;
pass2.pass_count = 3;
pass2.pass_decision = fmt::format(
"post-refined geometry adopted, re-run with pass-1's lattice forced: the de-novo "
"cell volume {:.0f} A^3 was {:.2f}x pass-1's ({:.0f} A^3), a spurious supercell",
v2, v2 / v1, v1);
force_rotation_result_.reset();
}
}
// The space group was DECIDED in pass 1 and only reused by pass 2, which therefore has no search
// of its own to report. Carry pass 1's over, or the evidence behind the decision (operator CCs,
// absence counts, the refused higher symmetry) is never shown on the default rotation path.
if (!pass2.space_group_search.has_value())
pass2.space_group_search = pass1.space_group_search;
// Pass 2 is normally the better answer, which is why it is the canonical output - but it is not
// guaranteed to be, and until now it was adopted whatever it produced. Two ways it can be wrong:
// it merges more unique reflections than the cell it settled on can hold (completeness above
// 100% is arithmetically impossible and means the cell is wrong), or its CC1/2 collapses
// relative to pass 1. Measured on a large-cell crystal: pass 1 195538 unique at 92.6% and
// CC1/2 0.98, pass 2 134667 at "117%" and CC1/2 0.62. Both bounds are set where only a failure
// reaches them, so a normal run - where pass 2 is a little better - keeps pass 2.
if (!cancelled_ && pass1.has_merge_statistics && pass2.has_merge_statistics) {
const auto completeness = [](const ProcessResult &r) {
const auto &o = r.merge_statistics.overall;
return o.possible_unique_reflections > 0
? 100.0 * o.unique_reflections / o.possible_unique_reflections : 0.0;
};
const double compl2 = completeness(pass2);
const double cc1 = pass1.merge_statistics.overall.cc_half;
const double cc2 = pass2.merge_statistics.overall.cc_half;
constexpr double MAX_CREDIBLE_COMPLETENESS = 100.5; // rounding headroom, nothing more
constexpr double MAX_CC_HALF_LOSS = 0.05;
if (pass2.lattice_conflicts_with_prepass_sg
|| compl2 > MAX_CREDIBLE_COMPLETENESS || cc2 < cc1 - MAX_CC_HALF_LOSS) {
logger.Warning("Two-pass: the refined pass is worse than the header-geometry pass "
"(completeness {:.1f}% vs {:.1f}%, CC1/2 {:.3f} vs {:.3f}) - going back to "
"the header geometry. The refined geometry did not help this crystal.",
compl2, completeness(pass1), cc2, cc1);
// The refined pass has already written the canonical files, so re-run at the header
// geometry to replace them - the same remedy the supercell collapse above uses, and it
// only costs a pass on a crystal that was going to be wrong otherwise.
experiment_.BeamX_pxl(header_geometry[0]).BeamY_pxl(header_geometry[1])
.DetectorDistance_mm(header_geometry[2]);
if (prepass_rotation_scale_ && gonio_snapshot)
experiment_.Goniometer(*gonio_snapshot); // and the header rotation angles with it
experiment_.SpaceGroupNumber(std::nullopt);
config_.output_prefix = base_prefix;
auto redo = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
if (!redo.space_group_search.has_value())
redo.space_group_search = pass1.space_group_search;
redo.post_refine = pass1.post_refine;
redo.pass_number = pass2.pass_count + 1;
redo.pass_count = pass2.pass_count + 1;
redo.pass_decision = fmt::format(
"header geometry re-adopted: the post-refined pass was worse "
"(completeness {:.1f}% vs {:.1f}%, CC1/2 {:.3f} vs {:.3f})",
compl2, completeness(pass1), cc2, cc1);
return redo;
}
if (pass2.pass_decision.empty())
pass2.pass_decision = fmt::format(
"post-refined geometry adopted (completeness {:.1f}% vs {:.1f}%, "
"CC1/2 {:.3f} vs {:.3f})", compl2, completeness(pass1), cc2, cc1);
}
if (pass2.pass_decision.empty())
pass2.pass_decision = (prepass_detector_geometry_ || prepass_rotation_scale_)
? "post-refined geometry adopted"
: "the post-refinement committed no geometry change, so this pass reproduces the first";
return pass2;
}
return RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
}
namespace {
// Angle between the crystal's symmetry axes and the spindle. A rotation sweep never records the
// reflections whose reciprocal vector lies close to the spindle - the blind cusp - and normally
// symmetry fills it in from an equivalent elsewhere in reciprocal space. It cannot when a symmetry
// axis IS the spindle: the cusp then maps onto itself, every reflection in it is equivalent only to
// other reflections in it, and the hole stays empty however long the sweep runs. Returns the
// smallest angle in degrees between any proper rotation axis of the space group and the goniometer
// axis, with the order of the axis that achieves it.
std::optional<std::pair<double, int>> ClosestSymmetryAxisToSpindle(const gemmi::SpaceGroup &sg,
const CrystalLattice &lattice,
const Coord &spindle) {
const double spindle_len = spindle.Length();
if (spindle_len < 1e-9)
return std::nullopt;
std::optional<std::pair<double, int>> best;
for (const auto &op : sg.operations().derive_symmorphic().sym_ops) {
if (op.rot == gemmi::Op::identity().rot)
continue;
// Order of the rotation, then project onto its invariant direction by summing its powers:
// (1/n) sum_k W^k annihilates everything except the axis.
const gemmi::Op w{op.rot, {0, 0, 0}, op.notation};
int order = 0;
gemmi::Op cur = gemmi::Op::identity();
std::array<std::array<double, 3>, 3> sum{};
for (int k = 0; k < 6; ++k) {
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
sum[i][j] += static_cast<double>(cur.rot[i][j]) / gemmi::Op::DEN;
cur = cur.combine(w).wrap();
++order;
if (cur.rot == gemmi::Op::identity().rot)
break;
}
if (order < 2)
continue;
// Any non-degenerate column of the projector is the axis in fractional direct coordinates.
for (int c = 0; c < 3; ++c) {
const Coord axis = lattice.Vec0() * static_cast<float>(sum[0][c])
+ lattice.Vec1() * static_cast<float>(sum[1][c])
+ lattice.Vec2() * static_cast<float>(sum[2][c]);
const double len = axis.Length();
if (len < 1e-6 * lattice.Vec0().Length())
continue;
const double cosang = std::abs((axis * spindle) / (len * spindle_len));
const double ang = std::acos(std::min(1.0, cosang)) * 180.0 / PI;
if (!best.has_value() || ang < best->first)
best = std::make_pair(ang, order);
break;
}
}
return best;
}
}
ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, bool geometry_prepass) {
Logger logger("Rugnux");
ProcessResult result;
const auto dataset = reader_.GetDataset();
if (!dataset)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"No experiment dataset found in the input file");
if (config_.stride <= 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Image stride must be positive");
const auto total_images_in_file = static_cast<int>(reader_.GetNumberOfImages());
int end_image = config_.end_image;
if (end_image < 0 || end_image > total_images_in_file)
end_image = total_images_in_file;
const int start_image = config_.start_image;
// Round UP: the worker loop below walks ordinals until start + ordinal * stride reaches end_image,
// so a stride that does not divide the range evenly still processes a final partial step. Rounding
// down here would size the per-image arrays one short of the ordinals that index them.
const int images_to_process = (end_image - start_image + config_.stride - 1) / config_.stride;
if (images_to_process <= 0) {
logger.Warning("No images to process (start {}, end {}, stride {}, total {})",
start_image, end_image, config_.stride, total_images_in_file);
return result;
}
const bool full = (config_.mode == ProcessMode::FullAnalysis);
// Powder calibration by spots finds spots and nothing else, so it runs the same per-image engine as
// the full analysis (spot finding lives inside MXAnalysisWithoutFPGA) with indexing switched off in
// its settings; by rings it only needs the azimuthal profile, which is the azint worker exactly.
const bool calibration = (config_.mode == ProcessMode::Calibration);
const bool calibration_spots = calibration && config_.calibration_method == CalibrationMethod::Spots;
const bool per_image_analysis = full || calibration_spots;
const bool write_files = write_output && !config_.output_prefix.empty();
// Output/runtime invariants. Algorithm settings (indexing, scaling, integration, polarization,
// space group, unit cell, ...) are configured on experiment_ by the caller.
experiment_.BitDepthImage(32).PixelSigned(true);
experiment_.Mode(DetectorMode::Standard);
experiment_.OverwriteExistingFiles(true);
// Offline processing: the output prefix is a local path the operator chose (rugnux -o or the
// viewer's "next to the input file"), so use the trusted setter that permits an absolute path.
experiment_.FilePrefixTrusted(config_.output_prefix.empty() ? "output" : config_.output_prefix);
experiment_.SetFileWriterFormat(FileWriterFormat::NXmxLegacy);
experiment_.ImagesPerTrigger(images_to_process);
experiment_.NumTriggers(1);
if (full)
experiment_.Compression(CompressionAlgorithm::BSHUF_LZ4);
// The pipeline indexes images 0..N-1 within this run; if we process a sub-range/strided
// selection, shift the goniometer so local index i maps to the angle of original image
// start+i*stride (keeping the per-image rotation wedge), otherwise rotation angles would be
// wrong for any start_image != 0.
if (const auto g = experiment_.GetGoniometer();
g.has_value() && (start_image != 0 || config_.stride != 1)) {
const float incr = g->GetIncrement_deg();
GoniometerAxis shifted(g->GetName(),
g->GetStart_deg() + incr * static_cast<float>(start_image),
incr * static_cast<float>(config_.stride),
g->GetAxis(), g->GetHelicalStep());
shifted.ScreeningWedge(g->GetScreeningWedge().value_or(incr));
experiment_.Goniometer(shifted);
}
// Stills global geometry-refinement first pass: bundle-adjust the shared beam/distance/cell from a
// sample of strongly-indexed frames and apply it to experiment_, so everything below (azimuthal
// mapping, the main index+integrate pass, scaling) uses the refined geometry. Rotation is skipped
// inside (it has its own two-pass).
if (full && config_.refine_geometry.has_value())
RefineStillsGeometry(start_image, end_image, images_to_process, observer);
// After any geometry refinement, so the shadow is found about the beam centre actually used,
// and before the azimuthal mapping and the output mask below, which both read pixel_mask_.
if (config_.detect_beam_stop.has_value() || config_.estimate_beam_center)
PreScan(start_image, images_to_process,
config_.detect_beam_stop.value_or(BEAM_CENTER_PROJECTION_IMAGES), observer);
AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_);
JFJochReceiverPlots plots;
plots.Setup(experiment_, mapping);
// Output file (NXmxIntegrated master that links back to the original images).
StartMessage start_message;
experiment_.FillMessage(start_message);
start_message.arm_date = dataset->arm_date;
start_message.az_int_bin_to_q = mapping.GetBinToQ();
start_message.az_int_bin_to_two_theta = mapping.GetBinToTwoTheta();
start_message.az_int_q_bin_count = mapping.GetQBinCount();
start_message.az_int_phi_bin_count = mapping.GetAzimuthalBinCount();
if (mapping.GetAzimuthalBinCount() > 1)
start_message.az_int_bin_to_phi = mapping.GetBinToPhi();
start_message.pixel_mask["default"] = pixel_mask_.GetMask(experiment_);
if (full) {
start_message.rois = experiment_.ROI().ExportMetadata();
if (!experiment_.ROI().empty())
start_message.roi_map = experiment_.ExportROIMap();
start_message.max_spot_count = experiment_.GetMaxSpotCount();
}
start_message.master_suffix = "process";
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);
// 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;
std::unique_ptr<FileWriter> writer;
if (write_files && config_.write_process_h5)
writer = std::make_unique<FileWriter>(start_message, /*check_overwrite_at_start=*/true,
/*trusted_path=*/true);
const char *mode_name = full ? "full analysis"
: calibration ? "powder calibration" : "azimuthal integration";
logger.Info("Processing {} images (range {}-{}, stride {}) using {} threads [{}]",
images_to_process, start_image, end_image, config_.stride, config_.nthreads, mode_name);
LogXDSGeometry(experiment_);
if (observer)
observer->OnPhase(full ? "Full analysis" : calibration ? "Powder calibration" : "Azimuthal integration");
// Full-analysis shared engines.
std::unique_ptr<IndexerThreadPool> indexer_pool;
std::unique_ptr<IndexAndRefine> indexer;
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
// parallel (the main per-image loop then uses the forced lattice, no FFT). Size the pool to
// those two so a rotation run doesn't pay to initialise (cuFFT plan + GPU alloc) indexers it
// never uses. Stills fire the FFT once per image across all workers, so keep the default there.
if (config_.rotation_indexing && config_.two_pass_rotation && experiment_.IsRotationIndexing()) {
indexing_settings.IndexingThreads(2);
// Two indexer threads leave the rest of the machine idle for the whole first pass, which is
// a third of the run. Split each scheme's candidate-cell refinement over -N/2 threads so the
// two together use -N. The split is by candidate cell, so the result does not depend on it.
// 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));
}
indexer_pool = std::make_unique<IndexerThreadPool>(indexing_settings, IndexerConstruction::OnFirstUse);
indexer = std::make_unique<IndexAndRefine>(experiment_, indexer_pool.get());
if (!config_.reference_data.empty())
indexer->ReferenceIntensities(config_.reference_data);
// Second pass of the rotation two-pass: predict with the smoothed mosaicity fitted in the pre-pass.
if (!geometry_prepass && !prepass_mosaicity_.empty())
indexer->SetPredictionMosaicityOverride(prepass_mosaicity_);
}
const auto start_time = std::chrono::steady_clock::now();
// First pass of two-pass rotation indexing (full analysis only).
if (full && force_rotation_result_.has_value()) {
// Supercell-collapse fallback: force pass-1's WHOLE indexing result (lattice + refined orientation /
// search metadata / axis), not just its lattice - a lattice-only force loses that metadata and
// integrates markedly worse.
indexer->ForceRotationIndexerResult(*force_rotation_result_);
logger.Info("Rotation indexer result forced from pass 1 (lattice + geometry) - skipping first pass");
} else if (full && config_.forced_rotation_lattice.has_value()) {
indexer->ForceRotationIndexerLattice(*config_.forced_rotation_lattice);
logger.Info("Rotation indexer lattice forced externally - skipping first pass");
} else if (full && config_.rotation_indexing && config_.two_pass_rotation) {
if (observer)
observer->OnPhase("Rotation indexing (first pass)");
// Mid-exposure goniometer angle for an image ordinal (matches the per-image path).
auto rot_angle = [&](int ordinal) -> std::optional<float> {
if (const auto g = experiment_.GetGoniometer())
return g->GetAngle_deg(static_cast<float>(ordinal)) + g->GetWedge_deg() / 2.0f;
return std::nullopt;
};
// Spots for a first-pass image ordinal, cached (a frame is reused across schemes and the
// validation set; re-finding is expensive on the --redo-rotation-spots path). Accessed only
// from single-threaded sections (the scheme feed and the validation loop), so no locking.
std::map<int, std::vector<SpotToSave>> spot_cache;
// One analysis engine for the whole first pass. Constructing it allocates a CUDA stream and a
// full set of GPU buffers (preprocessing, spot finding, azimuthal integration, Bragg
// integration) - far more work than analysing a frame - so building it per image would pay for
// hundreds of them, serially. This section is single-threaded, so one engine is enough.
auto analysis = std::make_unique<MXAnalysisWithoutFPGA>(experiment_, mapping, pixel_mask_, *indexer,
/*enable_fused_adaptive_gpu=*/true);
auto profile = std::make_unique<AzimuthalIntegrationProfile>(mapping);
auto get_spots = [&](int ordinal) -> const std::vector<SpotToSave> & {
auto it = spot_cache.find(ordinal);
if (it != spot_cache.end())
return it->second;
const int image_idx = start_image + ordinal * config_.stride;
std::vector<SpotToSave> spots;
try {
if (auto img = reader_.GetRawImage(image_idx)) {
DataMessage m{};
m.number = ordinal;
m.original_number = image_idx;
auto first_pass = config_.spot_finding;
first_pass.indexing = false;
first_pass.quick_integration = false;
m.image = img->image;
if (dataset->efficiency.size() > image_idx)
m.image_collection_efficiency = dataset->efficiency[image_idx];
analysis->Analyze(m, *profile, first_pass);
spots = std::move(m.spots);
}
} catch (const std::exception &e) {
logger.Warning("First-pass spot read failed for image {}: {}", image_idx, e.what());
}
return spot_cache.emplace(ordinal, std::move(spots)).first->second;
};
// How many validation frames a candidate global lattice actually indexes per image, scored on
// the *real* per-image path (force the candidate, then run the same indexing+refinement the
// main loop uses). A cell that merely fits the accumulated first-pass cloud (a spurious
// sub-lattice from an under-sampled scheme) indexes few frames, while the true cell indexes
// many - this is the discriminator that lets the best scheme win.
const std::vector<int> validation = select_equally_spaced_image_ordinals(
images_to_process, std::min(images_to_process, 60));
auto validation_settings = config_.spot_finding;
validation_settings.indexing = true;
validation_settings.quick_integration = false;
auto count_indexed = [&](IndexAndRefine &idx, const RotationIndexerResult &r) -> int {
idx.ForceRotationIndexerResult(r);
int count = 0;
for (const int ordinal : validation) {
DataMessage m{};
m.number = ordinal;
m.spots = get_spots(ordinal);
if (idx.IndexFrameOnly(m, validation_settings))
count++;
}
return count;
};
// Feed one first-pass scheme (a set of image ordinals) into its own rotation indexer, ready to
// be indexed. Spots are pulled from the single-threaded cache here; the FFT + refinement runs
// separately (RunIndexing) so the two schemes' indexing can overlap.
auto feed_scheme = [&](IndexerThreadPool &pool, const std::vector<int> &ordinals) {
auto ri = std::make_unique<RotationIndexer>(experiment_, pool);
for (const int ordinal : ordinals) {
if (cancelled_ || ri->AccumulationFull())
break;
ri->ProcessImage(ordinal, get_spots(ordinal), rot_angle(ordinal));
}
return ri;
};
// Sampling schemes, tried in order of preference. We index each and keep the one that indexes
// the most frames; the schemes cover complementary failure modes and no single one is best for
// every crystal:
// - "spread": images spread over the whole rotation - full angular range, coarse stride; the
// long-standing default, best for typical crystals.
// - "wedge": consecutive frames from the clean start - native stride keeps each reflection's
// rocking curve continuous across frames, so the FFT can trace a long axis whose fine
// reciprocal spacing the coarse spread cannot resolve (e.g. an F-cubic cell viewed down
// [111]); bounded by the accumulation cap to the least radiation-damaged early wedge.
std::vector<int> all_consecutive(images_to_process);
std::iota(all_consecutive.begin(), all_consecutive.end(), 0);
const std::vector<int> spread_ordinals = select_equally_spaced_image_ordinals(
images_to_process, config_.rotation_indexing_image_count);
const std::vector<std::pair<std::string, std::vector<int>>> schemes = {
{"spread", spread_ordinals},
{"wedge", all_consecutive},
};
struct FirstPass { std::optional<RotationIndexerResult> result; int score = -1; std::string name; double vol = 0.0; };
// Feed both schemes (single-threaded spot cache), run their RunIndexing() concurrently, then keep
// the lattice that indexes the most validation frames on the real per-image path.
auto pick_best = [&](IndexerThreadPool &pool, IndexAndRefine &idx) -> FirstPass {
std::vector<std::unique_ptr<RotationIndexer>> ris;
ris.reserve(schemes.size());
for (const auto &[name, ordinals] : schemes) {
if (cancelled_)
break;
ris.push_back(feed_scheme(pool, ordinals));
}
std::vector<std::future<void>> index_futures;
index_futures.reserve(ris.size());
for (auto &ri : ris) {
RotationIndexer *rp = ri.get();
index_futures.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); }));
}
for (auto &f : index_futures)
f.get();
FirstPass bp;
for (size_t i = 0; i < ris.size(); i++) {
if (cancelled_)
break;
auto found = ris[i]->GetLattice();
if (!found.has_value())
continue;
RotationIndexerResult result = std::move(*found);
int score = count_indexed(idx, result);
const std::string &name = schemes[i].first;
logger.Info("First-pass scheme '{}': indexes {}/{} validation frames", name, score,
static_cast<int>(validation.size()));
// Metric symmetry that costs nearly every frame is not this crystal's symmetry. The
// Bravais class is decided on the UNREFINED FFT candidate against a fixed angular
// tolerance, so a lattice that is pseudo-symmetric to a few tenths of a degree gets
// promoted a class too far; the constraint then snaps a real angle to the ideal one,
// which throws nearly every reflection of every frame out of tolerance (measured: 2/60
// frames for a monoclinic lattice promoted to C-orthorhombic, 39/60 for its primitive
// cell). The indexer's own guard compares the two on the accumulated first-pass cloud,
// where they differ by less than a factor 2 - far too little to act on.
// The bar is a clear majority, not a margin: the unconstrained refinement holds NO cell
// parameter fixed, so it can only index at least as many frames as the constrained one,
// and on genuine symmetry it does index a few more (measured: 47->54 frames on a real
// I-centred orthorhombic crystal, 49->58 on an F-cubic one - a 10% margin demotes both).
// Only a constrained cell that fails outright while its unconstrained cell works is
// evidence of a false promotion, so demand exactly that. It is the same "fails to index
// half the frames" test the long-axis rescue below uses, and it leaves genuine symmetry
// alone; the space group is in any case decided from the intensities later.
const int majority = static_cast<int>(validation.size()) / 2;
if (result.unconstrained && score < majority) {
RotationIndexerResult alt = *result.unconstrained;
const int alt_score = count_indexed(idx, alt);
if (alt_score > majority) {
logger.Info("Scheme '{}': {}-centred {} indexes {}/{} frames but its unconstrained "
"cell indexes {}/{} - the metric symmetry is a false promotion, dropping it",
name, result.search_result.centering,
gemmi::crystal_system_str(result.search_result.system), score,
static_cast<int>(validation.size()), alt_score,
static_cast<int>(validation.size()));
result = std::move(alt);
score = alt_score;
}
}
// Compare PRIMITIVE volumes: two schemes can find the same lattice in different
// settings, and a centred setting's cell is an exact integer multiple of its primitive
// one - a rhombohedral lattice in hexagonal axes is exactly 3x its primitive
// rhombohedral cell. Comparing the centred volumes makes the supercell test below fire
// on that pair and "demote" a perfectly good setting to a threefold-smaller merge, which
// is enough to change the space group the search then picks.
const double vol = std::abs(
result.lattice.ToPrimitive(result.search_result.centering).CalcVolume());
// A later scheme wins if it indexes clearly more frames (>10%).
const bool clearly_more = static_cast<float>(score) > bp.score * 1.1f + 0.5f;
// Integer-supercell tie-break. The validation-frame count saturates - a spurious axis
// multiple (2x/3x...) indexes every frame its true cell does, so both schemes reach the
// same frame total and the count alone cannot tell them apart; the default then keeps
// whichever ran first. When the two schemes tie on frames but their cell volumes are
// related by an integer factor >=2, the larger cell is that spurious supercell and the
// smaller is the true reduced cell: take it, regardless of scheme order. A near-integer
// ratio separates a real axis multiplication from a centering coincidence. Volumes are
// primitive (see above), so a pure setting difference is a ratio of 1 and never fires.
bool integer_subcell = false;
if (bp.result.has_value() && !clearly_more && vol > 1.0 && bp.vol > 1.0) {
const bool tied = static_cast<float>(score) >= bp.score * 0.9f - 0.5f;
const double ratio = (vol < bp.vol) ? bp.vol / vol : vol / bp.vol;
const double nearest = std::round(ratio);
const bool integer_multiple = nearest >= 2.0 && std::abs(ratio - nearest) < 0.15;
integer_subcell = tied && integer_multiple && vol < bp.vol;
}
if (!bp.result.has_value() || clearly_more || integer_subcell) {
if (integer_subcell)
logger.Info("Scheme '{}' cell (vol {:.0f}) is a {:.0f}x sub-cell of '{}' (vol {:.0f}) at "
"equal frame count - adopting the smaller true cell", name, vol,
std::round(bp.vol / vol), bp.name, bp.vol);
bp.score = score;
bp.vol = vol;
bp.name = name;
bp.result = std::move(result);
}
}
return bp;
};
FirstPass best = pick_best(*indexer_pool, *indexer);
// 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
// many high-order reflections off along the fine axis, so the true cell scores worst on the raw
// cloud). Recover the true metric with a COARSE-resolution pass - only low-order reflections,
// where the fine axis stays robust - then RE-INDEX at full resolution constrained by that cell as
// a reference (the -C path): the reference filter drops the collapsed/supercell candidates and
// refines an accurate global lattice. Only runs after a poor standard pass, so well-indexing
// crystals never pay for it and their result is untouched.
auto max_axis = [](const RotationIndexerResult &r) {
const auto uc = r.lattice.GetUnitCell();
return std::max({uc.a, uc.b, uc.c});
};
if (!cancelled_ && best.result.has_value()
&& 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
coarse_settings.IndexingThreads(2);
coarse_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS));
IndexerThreadPool coarse_pool(coarse_settings, IndexerConstruction::OnFirstUse);
// Coarse first pass: keep the recovered cell with the LONGEST axis directly. Its full-
// resolution per-frame validation would be low (the coarse cell is metrically right but
// imprecise), so do NOT score it here - it is only the reference for the constrained re-index.
std::vector<std::unique_ptr<RotationIndexer>> cris;
cris.reserve(schemes.size());
for (const auto &[name, ordinals] : schemes) {
if (cancelled_) break;
cris.push_back(feed_scheme(coarse_pool, ordinals));
}
std::vector<std::future<void>> cfut;
cfut.reserve(cris.size());
for (auto &ri : cris) {
RotationIndexer *rp = ri.get();
cfut.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); }));
}
for (auto &f : cfut) f.get();
std::optional<RotationIndexerResult> coarse_ref;
for (auto &ri : cris) {
auto r = ri->GetLattice();
if (r.has_value() && (!coarse_ref.has_value() || max_axis(*r) > max_axis(*coarse_ref)))
coarse_ref = r;
}
// A recovered axis at least 30% longer than the standard pass found signals a rescued metric.
if (!cancelled_ && coarse_ref.has_value()
&& max_axis(*coarse_ref) > 1.3 * max_axis(*best.result)) {
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);
experiment_.SetUnitCell(std::nullopt); // leave the space-group / cell determination de-novo
if (constrained.result.has_value() && constrained.score > best.score)
best = constrained;
}
}
if (!cancelled_) {
// The two ways the first pass can end without a usable lattice. They are one situation from
// the user's side - no crystal was found in these frames - so both name the cause and then
// give the same remedies.
constexpr const char *next_steps =
" Nothing was integrated. Check the beam centre (--beam-x / --beam-y), raise "
"--max-spots if the frames are well populated with spots, and use --mode azint to "
"write the per-image analysis for inspection.";
if (!best.result.has_value())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
std::string("Two-pass rotation indexing found no lattice.") + next_steps);
// A lattice that indexes almost none of the validation frames is not this crystal's lattice.
// Forcing it anyway integrates every image against it; the few partials that survive yield no
// fulls in the 3D combine, and the run then dies inside the merge complaining about resolution -
// far from the cause. Stop here instead, while the reason is still visible. The floor is one
// frame in six, several times below the weakest crystal that still merges. It is checked AFTER
// the long-axis rescue, so a metric that rescue recovers is never rejected on its pre-rescue
// score.
if (best.score < static_cast<int>(validation.size()) / 6) {
// Name the cell and Bravais class that was rejected. The commonest cause is a metric
// symmetry promoted one class too far - the constrained cell then misses every
// reflection by the small angle the constraint snapped away - and without the cell in
// the message there is nothing to see that from.
const auto &c = best.result->search_result.conventional.GetUnitCell();
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
fmt::format("Two-pass rotation indexing found only a lattice that "
"indexes {}/{} validation frames - too few for it to be "
"this crystal's lattice. It was {}-centred {}, "
"{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}.{}",
best.score, static_cast<int>(validation.size()),
best.result->search_result.centering,
gemmi::crystal_system_str(best.result->search_result.system),
c.a, c.b, c.c, c.alpha, c.beta, c.gamma, next_steps));
}
indexer->ForceRotationIndexerResult(*best.result);
logger.Info("Two-pass rotation indexing found lattice (scheme '{}': {}/{} validation frames)",
best.name, best.score, static_cast<int>(validation.size()));
// Second pass: compare the de-novo lattice with pass 1's HERE, before integrating every image
// with it. A markedly larger cell is the bistable supercell collapse, and a centring the group
// carried over from pass 1 cannot describe is the same disagreement seen from the symmetry side.
// Both were already rejected - the first in Run(), the second before the merge - but only after
// the whole pass had run, so each cost an entire extra pass on a dataset that ended up on pass-1's
// lattice anyway. Take pass-1's result now instead: the same remedy, one pass earlier. A triclinic
// de-novo cell is left alone, since that is the demotion the merge reindexes into the group's
// conventional setting.
if (!geometry_prepass && prepass_result_.has_value()) {
const auto *reuse_sg = prepass_merge_sg_.has_value()
? gemmi::find_spacegroup_by_number(static_cast<int>(*prepass_merge_sg_)) : nullptr;
// Primitive volumes, like the scheme comparison above: a centred conventional cell is an
// exact integer multiple of its primitive one, so two settings of the SAME lattice differ
// by that factor and comparing the conventional cells reads a mere change of setting as a
// supercell - which then forces pass 1's setting and, with it, its lower symmetry.
const double v1 = std::abs(prepass_result_->lattice
.ToPrimitive(prepass_result_->search_result.centering).CalcVolume());
const double v2 = std::abs(best.result->lattice
.ToPrimitive(best.result->search_result.centering).CalcVolume());
const bool supercell = v1 > 1.0 && v2 > 1.5 * v1;
const bool wrong_centering = reuse_sg != nullptr
&& best.result->search_result.system != gemmi::CrystalSystem::Triclinic
&& best.result->search_result.centering != reuse_sg->centring_type();
if (supercell || wrong_centering) {
logger.Info("Two-pass: that lattice ({}-centred, volume {:.0f} A^3) disagrees with pass 1 "
"({}-centred, {:.0f} A^3) - integrating with pass-1's lattice instead",
best.result->search_result.centering, v2,
prepass_result_->search_result.centering, v1);
indexer->ForceRotationIndexerResult(*prepass_result_);
}
}
}
}
// Main per-image loop, spread over N worker threads pulling from a shared counter. HDF5 reads
// are serialized by the global hdf5_mutex; the analysis runs in parallel.
std::atomic<int> next_ordinal = 0;
std::atomic<int> finished_count = 0;
// Images are written at their own ordinal, so a frame that fails leaves a HOLE rather than
// shifting everything after it up. The extent of the written datasets is therefore the highest
// ordinal reached, not the number that succeeded - counting successes makes the file one image
// short for every failure, and the image it loses is the LAST one, not the one that failed.
std::atomic<int> max_written_ordinal = -1;
auto note_written = [&max_written_ordinal](int ordinal) {
int prev = max_written_ordinal.load(std::memory_order_relaxed);
while (prev < ordinal
&& !max_written_ordinal.compare_exchange_weak(prev, ordinal, std::memory_order_relaxed)) {
}
};
std::atomic<uint64_t> total_uncompressed_bytes = 0;
// Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and
// the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of
// a fixed budget instead of everything it found - a powder ring is over-determined either way.
constexpr size_t CALIBRATION_MAX_SPOTS = 20000;
const size_t calibration_spots_per_image =
std::max<size_t>(1, CALIBRATION_MAX_SPOTS / static_cast<size_t>(images_to_process));
std::mutex calibration_spot_m;
std::vector<SpotToSave> calibration_spot_list;
auto azint_worker = [&]() {
std::vector<uint8_t> decompression_buffer;
ImagePreprocessorCPU preprocessor(experiment_, pixel_mask_);
ImagePreprocessorBuffer buffer(experiment_.GetPixelsNum());
AzIntEngineCPU azint(mapping);
AzimuthalIntegrationProfile profile(mapping);
while (!cancelled_) {
const int ordinal = next_ordinal.fetch_add(1);
const int image_idx = start_image + ordinal * config_.stride;
if (image_idx >= end_image) break;
std::shared_ptr<JFJochReaderRawImage> img;
try {
img = reader_.GetRawImage(image_idx);
} catch (const std::exception &e) {
logger.Error("Failed to load image {}: {}", image_idx, e.what());
continue;
}
if (!img) continue;
DataMessage msg{};
msg.image = img->image;
msg.number = ordinal;
msg.original_number = image_idx;
if (dataset->efficiency.size() > image_idx)
msg.image_collection_efficiency = dataset->efficiency[image_idx];
total_uncompressed_bytes += msg.image.GetUncompressedSize();
const auto t0 = std::chrono::steady_clock::now();
try {
const uint8_t *image_ptr = msg.image.GetUncompressedPtr(decompression_buffer);
preprocessor.Analyze(buffer, image_ptr, msg.image.GetMode());
azint.Run(buffer, profile);
} catch (const std::exception &e) {
logger.Error("Error integrating image {}: {}", image_idx, e.what());
continue;
}
msg.azint_time_s = std::chrono::duration<float>(std::chrono::steady_clock::now() - t0).count();
msg.processing_time_s = msg.azint_time_s;
msg.az_int_profile = profile.GetResult();
msg.az_int_profile_count = profile.GetPixelCount();
msg.az_int_profile_std = profile.GetStd();
msg.bkg_estimate = profile.GetBkgEstimate(mapping.Settings());
msg.ice_ring_score = profile.GetIceRingScore(mapping.Settings(),
config_.spot_finding.ice_ring_width_Q_recipA);
msg.run_number = experiment_.GetRunNumber();
msg.run_name = experiment_.GetRunName();
plots.Add(msg, profile);
if (writer) writer->Write(msg);
note_written(ordinal);
if (observer) observer->OnImageProcessed(msg);
const int done = finished_count.fetch_add(1) + 1;
if (observer) observer->OnProgress(done, images_to_process);
}
};
auto full_worker = [&]() {
pin_gpu(); // round-robin per worker thread; must precede engine construction
MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, *indexer,
/*enable_fused_adaptive_gpu=*/true);
AzimuthalIntegrationProfile profile(mapping);
while (!cancelled_) {
const int ordinal = next_ordinal.fetch_add(1);
const int image_idx = start_image + ordinal * config_.stride;
if (image_idx >= end_image) break;
std::shared_ptr<JFJochReaderRawImage> img;
try {
img = reader_.GetRawImage(image_idx);
} catch (const std::exception &e) {
logger.Error("Failed to load image {}: {}", image_idx, e.what());
continue;
}
if (!img) continue;
DataMessage msg{};
msg.image = img->image;
msg.number = ordinal;
msg.original_number = image_idx;
if (dataset->efficiency.size() > image_idx)
msg.image_collection_efficiency = dataset->efficiency[image_idx];
total_uncompressed_bytes += msg.image.GetUncompressedSize();
const auto t0 = std::chrono::steady_clock::now();
try {
analysis.Analyze(msg, profile, config_.spot_finding);
} catch (const std::exception &e) {
logger.Error("Error analyzing image {}: {}", image_idx, e.what());
continue;
}
msg.processing_time_s = std::chrono::duration<float>(std::chrono::steady_clock::now() - t0).count();
msg.run_number = experiment_.GetRunNumber();
msg.run_name = experiment_.GetRunName();
if (calibration_spots) {
std::lock_guard lock(calibration_spot_m);
for (size_t i = 0; i < msg.spots.size() && i < calibration_spots_per_image; ++i)
calibration_spot_list.push_back(msg.spots[i]);
}
plots.Add(msg, profile);
if (writer) writer->Write(msg);
note_written(ordinal);
if (observer) observer->OnImageProcessed(msg);
const int done = finished_count.fetch_add(1) + 1;
if (observer) observer->OnProgress(done, images_to_process);
}
};
if (observer)
observer->OnPhase("Processing images");
std::function<void()> worker = per_image_analysis ? std::function<void()>(full_worker)
: std::function<void()>(azint_worker);
std::vector<std::future<void> > futures;
futures.reserve(config_.nthreads);
const auto image_loop_start = std::chrono::steady_clock::now();
for (int i = 0; i < config_.nthreads; ++i)
futures.push_back(std::async(std::launch::async, worker));
for (auto &f: futures)
f.get();
// Wall time of the per-image loop alone. The per-stage means below cover only what runs inside it,
// so without this there is nothing to compare them against and no way to see how much of a run is
// spent outside it - on the first-pass indexing and on scaling/merging.
result.image_loop_time_s = std::chrono::duration<double>(
std::chrono::steady_clock::now() - image_loop_start).count();
result.cancelled = cancelled_;
result.images_processed = finished_count.load();
// Every image failing is a total failure, not a run that produced nothing: it used to be reported
// only as per-image log lines while the process still exited 0 with no output file.
if (!cancelled_ && images_to_process > 0 && result.images_processed == 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"No image could be analyzed - see the per-image errors above");
result.mean_processing_time = plots.GetMeanProcessingTime();
result.indexing_rate = plots.GetIndexingRate();
// Powder calibration: the whole run has now been accumulated - the summed (q x azimuth) profile in
// plots, or the pooled spot list - so fit the detector geometry to it. The caller reports it and
// writes the PONI file; nothing here is stored in the _process.h5.
if (calibration && !cancelled_) {
if (observer)
observer->OnPhase("Powder calibration");
if (calibration_spots) {
logger.Info("Powder calibration from {} pooled spots", calibration_spot_list.size());
result.calibration = CalibrateFromSpots(calibration_spot_list,
experiment_.GetDiffractionGeometry(),
config_.calibrant_ring_q);
} else {
const auto *summed = plots.GetAzIntProfileObject();
if (!summed)
throw JFJochException(JFJochExceptionCategory::CalibrationError,
"No azimuthal profile was accumulated");
result.calibration = CalibrateFromProfile(summed->GetResult(), mapping,
experiment_.GetDiffractionGeometry(),
config_.calibrant_ring_q);
}
}
// End message (also written to the file).
EndMessage end_msg;
// The EXTENT of the per-image datasets (see note_written), not the success count: those two
// differ by exactly the number of frames that failed, and the writer sizes its virtual datasets
// from this.
end_msg.max_image_number = static_cast<uint64_t>(max_written_ordinal.load() + 1);
end_msg.images_collected_count = result.images_processed;
end_msg.images_sent_to_write_count = result.images_processed;
end_msg.end_date = time_UTC(std::chrono::system_clock::now());
end_msg.run_number = experiment_.GetRunNumber();
end_msg.run_name = experiment_.GetRunName();
end_msg.bkg_estimate = plots.GetBkgEstimate();
end_msg.ice_ring_score = plots.GetIceRingScoreArray();
end_msg.ice_ring_score_mean = plots.GetIceRingScore();
end_msg.az_int_result["dataset"] = plots.GetAzIntProfile();
end_msg.indexing_rate = result.indexing_rate;
if (full && !cancelled_) {
if (const auto rot = indexer->FinalizeRotationIndexing(); rot.has_value()) {
end_msg.rotation_lattice = rot->lattice;
// Write the refined geometry (not the nominal StartMessage values) into the _process.h5.
end_msg.refined_beam_center_x = rot->geom.GetBeamX_pxl();
end_msg.refined_beam_center_y = rot->geom.GetBeamY_pxl();
end_msg.refined_poni_rot1 = rot->geom.GetPoniRot1_rad();
end_msg.refined_poni_rot2 = rot->geom.GetPoniRot2_rad();
end_msg.refined_poni_rot3 = rot->geom.GetPoniRot3_rad();
if (rot->axis)
end_msg.refined_rotation_axis = rot->axis->GetAxis();
end_msg.rotation_lattice_type = LatticeMessage{
.centering = rot->search_result.centering,
.niggli_class = rot->search_result.niggli_class,
.crystal_system = rot->search_result.system
};
result.rotation_lattice_found = true;
}
result.consensus_cell = indexer->GetConsensusUnitCell();
end_msg.unit_cell = result.consensus_cell;
}
// Scaling and merging (full analysis only).
if (full && !cancelled_ && result.indexing_rate.has_value() && result.indexing_rate > 0
&& (config_.run_scaling || !config_.reference_data.empty())) {
// Scaling/merging is a long post-pass; report each sub-step as a phase so the GUI progress
// bar reflects what is happening instead of freezing on one label.
auto phase = [&](const std::string &p) {
if (observer) observer->OnPhase(p);
};
phase("Scaling and merging");
// Ice-ring handling (--detect-ice-rings): flag reflections sitting on a hexagonal-ice powder
// ring. Their integrated intensity is contaminated by the strong, variable ice background, so
// they drag the per-image scale fit. Flag them here so scaling (the per-image G fit and the
// fulls refit) skips them, while the combine, the merge and the statistics keep them - dropping
// them outright guts low/mid-resolution completeness on crystals that merge them fine.
// The eleven hexagonal-ice bands cover 16-26% of the unique reflections at typical resolutions
// whether or not the crystal has any ice, so handling ice unconditionally taxes clean data for
// nothing. Gate the flagging on the run's own measured ice strength (1 = no ice; the strongest
// ring over the smooth radial background, measured on the peak-excluded profile so Bragg peaks
// cannot pose as ice). The score is logged either way, so a run that gates off still reports why.
// Ice is looked for on TWO channels, because it arrives in two forms. Fine-grained ice is a
// smooth powder ring and shows up in the radial profile (ice_ring_score). Ice in large
// crystallites arrives as discrete reflections instead, leaves the profile flat, and is only
// visible as a pile-up of found SPOTS on the ring positions - measured against the same q width
// of ice-free flanks beside them, so the fall-off of spot density with resolution cancels. On
// the rotation battery the two channels barely overlap: the crystals with smooth ice read
// 2.1-2.4 / ~1.0, the ones whose spots are mostly ice read ~1.1 / 3.8-14.6, and a clean crystal
// reads 1.04 on both. Either channel alone is blind to half the cases.
const auto ice_score = plots.GetIceRingScore();
const auto ice_spot_ratio = plots.GetIceRingSpotRatio();
const float ice_min_score = experiment_.GetScalingSettings().GetIceMinScore();
const float ice_min_spot_ratio = experiment_.GetScalingSettings().GetIceMinSpotRatio();
const bool ice_present = (!ice_score.has_value() || *ice_score >= ice_min_score)
|| (ice_min_spot_ratio > 0.0f && ice_spot_ratio.has_value()
&& *ice_spot_ratio >= ice_min_spot_ratio);
if (experiment_.IsDetectIceRings() && !ice_present) {
logger.Info("Ice-ring handling: measured ice score {:.2f} < {:.2f} and spot ratio {:.2f} < "
"{:.2f}, no ice detected - ice-ring handling skipped entirely",
*ice_score, ice_min_score, ice_spot_ratio.value_or(NAN), ice_min_spot_ratio);
} else if (experiment_.IsDetectIceRings()) {
const float ice_width = config_.spot_finding.ice_ring_width_Q_recipA;
size_t total = 0, flagged = 0;
for (auto &outcome : indexer->GetIntegrationOutcome()) {
for (auto &r : outcome.reflections) {
++total;
r.on_ice_ring = IsOnIceRing(r.d, ice_width);
if (r.on_ice_ring)
++flagged;
}
}
logger.Info("Ice-ring handling: ice score {:.2f} (gate {:.2f}), spot ratio {:.2f} (gate "
"{:.2f}); flagged {} of {} reflections on ice rings (half-width {:.3f} A^-1); "
"excluded from scaling, kept for merging",
ice_score.value_or(NAN), ice_min_score, ice_spot_ratio.value_or(NAN),
ice_min_spot_ratio, flagged, total, ice_width);
}
// Scale the images and merge. Factored so it can run twice: first in P1 to give the
// space-group search a merged dataset, then again in the determined space group so the
// scaling sees symmetry equivalents and the final statistics are in the right symmetry.
// (With a user-fixed space group it simply runs once, already in that symmetry.)
struct ScaleMergeResult {
std::vector<MergedReflection> merged;
MergeStatistics statistics;
};
// The reference path computes each image's G once (per-image scaling against the
// reference); the scaling loop below is skipped, so G is stable across the two passes.
// Smoothing it more than once would compound the correction, so do it only on the first
// pass. The no-reference path recomputes G from scratch each pass and re-smooths correctly.
// Rotation self-scaling + 3D combine + merge is done by the dedicated RotationScaleMerge (a single
// allocate-once engine that recomputes partiality from the fitted mosaicity, combines the per-frame
// partials, scales the fulls and merges; the whole hot path runs on the GPU when one is present).
// Ingested once here (after ice flagging) and reused across both space-group passes; a fixed
// (forced) mosaicity is handled by the recompute. The decay + absorption correction surfaces run as
// post-scale-fulls stages (on by default, see ScalingSettings::CorrectionSurfaces); external-
// reference scaling, the stills -B (per-image B-factor) and wedge refinement are unsupported here.
// Two-pass second pass: if the de-novo indexer DEMOTED to a triclinic (primitive) cell while the reused
// merge space group is higher-symmetry (a huge oblique cell whose constrained centred refine was
// ill-posed, so the pseudo-symmetry guard kept the primitive), the reflections are in the primitive
// frame but the merge expects the group's CONVENTIONAL setting - reindex them into it (the same
// primitive->conventional change of basis the space-group search applies to a demoted result), else the
// merge folds the wrong equivalents. No-op when the indexer already produced the group's setting.
if (prepass_merge_sg_ && end_msg.rotation_lattice.has_value() && end_msg.rotation_lattice_type.has_value()
&& end_msg.rotation_lattice_type->crystal_system == gemmi::CrystalSystem::Triclinic) {
const auto *msg_sg = gemmi::find_spacegroup_by_number(static_cast<int>(*prepass_merge_sg_));
if (msg_sg && msg_sg->crystal_system() != gemmi::CrystalSystem::Triclinic) {
const auto cand = LatticeSearch(*end_msg.rotation_lattice);
// Exact integer reindex prim->conv: P[i][j] = conv_real[i] . rot_reciprocal[j] (as the
// centred-lattice test), robust to the refined cell not being exactly LatticeSearch's Niggli cell.
const Coord cv[3] = {cand.conventional.Vec0(), cand.conventional.Vec1(), cand.conventional.Vec2()};
const Coord rs[3] = {end_msg.rotation_lattice->Astar(), end_msg.rotation_lattice->Bstar(),
end_msg.rotation_lattice->Cstar()};
gemmi::Mat33 reindex;
double reindex_res = 0.0;
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j) {
const double v = cv[i] * rs[j];
reindex.a[i][j] = std::round(v);
reindex_res = std::max(reindex_res, std::fabs(v - reindex.a[i][j]));
}
// Only when LatticeSearch recovers the reused group's own metric (system + centring): otherwise
// the reused group and the indexed metric disagree and reindexing would be wrong - leave it.
if (cand.system == msg_sg->crystal_system() && cand.centering == msg_sg->centring_type()
&& reindex_res < 0.1) {
for (auto &io : indexer->GetIntegrationOutcome()) {
io.latt = io.latt.Multiply(reindex);
for (auto &r : io.reflections) {
const double h = r.h, k = r.k, l = r.l;
r.h = static_cast<int32_t>(std::lround(reindex.a[0][0]*h + reindex.a[0][1]*k + reindex.a[0][2]*l));
r.k = static_cast<int32_t>(std::lround(reindex.a[1][0]*h + reindex.a[1][1]*k + reindex.a[1][2]*l));
r.l = static_cast<int32_t>(std::lround(reindex.a[2][0]*h + reindex.a[2][1]*k + reindex.a[2][2]*l));
}
}
result.consensus_cell = cand.conventional.GetUnitCell();
end_msg.unit_cell = result.consensus_cell;
end_msg.rotation_lattice = end_msg.rotation_lattice->Multiply(reindex);
end_msg.rotation_lattice_type = LatticeMessage{ .centering = cand.centering,
.niggli_class = end_msg.rotation_lattice_type->niggli_class, .crystal_system = cand.system };
logger.Info("Two-pass: reindexed the de-novo primitive cell into the space-group {} "
"conventional setting for the merge", static_cast<int>(*prepass_merge_sg_));
}
}
}
// Two-pass second pass only: the indexer above ran DE NOVO (Run() cleared the space group so its
// pseudo-symmetry safeguards find the true cell); now reinstate pass-1's space group for the merge,
// exactly as a user -S would, so the symmetry is fixed without a re-search. No-op on pass 1 / single pass.
//
// Only when the lattice pass 2 actually found HAS that group's centring. Re-indexing de novo can
// land in a different setting from pass 1 - most often the primitive sub-cell of a centred
// lattice, which is what the reindex above exists to undo, and which it declines to do when the
// metric does not match. Stamping a centred group onto a primitive cell is not a small error:
// the centring absence rule then removes half the reflections that genuinely exist, so the merge
// reports more unique reflections than its own cell can hold (measured: a C2 group on a cell of
// exactly half the C-centred volume, reported as 117% complete with CC1/2 0.62 against the first
// pass's 92.6% and 0.98). When they disagree, let pass 2 search for itself rather than trust a
// group that belongs to a different lattice.
if (prepass_merge_sg_ && end_msg.rotation_lattice_type.has_value()) {
const auto *reuse_sg = gemmi::find_spacegroup_by_number(static_cast<int>(*prepass_merge_sg_));
if (reuse_sg && end_msg.rotation_lattice_type->centering != reuse_sg->centring_type()) {
logger.Warning("Two-pass: the second pass indexed a {}-centred lattice, but the first pass "
"determined {}, whose lattice is {}-centred - the group cannot describe this "
"lattice, so this pass will not be adopted",
end_msg.rotation_lattice_type->centering, reuse_sg->xhm(),
reuse_sg->centring_type());
result.lattice_conflicts_with_prepass_sg = true;
}
}
if (prepass_merge_sg_)
experiment_.SpaceGroupNumber(*prepass_merge_sg_);
const auto &rot_ss = experiment_.GetScalingSettings();
const bool is_rotation = experiment_.IsRotationIndexing(); // rotation indexing -> rotation scaling/merge
std::optional<RotationScaleMerge> rsm;
if (is_rotation) {
if (experiment_.GetRefineRotationWedgeInScaling()
|| rot_ss.GetRotationWedgeForScaling().has_value())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Rotation scaling/merging (RotationScaleMerge) does not support "
"wedge refinement");
// A reference MTZ is allowed for rotation: it fixes the space group / cell (on the CLI) and
// resolves the indexing ambiguity (below), but is NOT used to scale - the rotation merge stays
// self-consistent.
rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
static_cast<int>(config_.scaling_iter),
config_.nthreads, logger, config_.observation_dump_path);
rsm->Ingest();
}
auto scale_and_merge = [&](const std::string &label, bool for_search) -> ScaleMergeResult {
if (rsm) {
phase("Scale/combine/merge (" + label + ")");
auto r = rsm->Run(for_search);
result.error_model_isa = r.isa;
result.error_model_isa_asymptotic = r.isa_asymptotic;
result.error_model_a = r.error_model_a;
result.error_model_b = r.error_model_b;
return ScaleMergeResult{std::move(r.merged), std::move(r.statistics)};
}
// Stills (rotation goes through RotationScaleMerge above): self-scale each image against the
// running merge with ScaleOnTheFly (fixed partiality), then merge directly. This runs even
// with an external reference: the reference is used only to break the per-image indexing
// ambiguity and to report CCref / inherit the R-free set, NOT as a scale anchor (scaling each
// image against a foreign dataset injects cross-dataset systematics and is a worse reference
// than the data's own merge). One pass: the per-image scale G is the exact one-pass solution,
// so iterating would only rebuild the reference from the freshly-scaled (noisy) data and
// degrade weak stills (measured CC1/2 collapse at the default 3 iters).
{
phase("Scaling images (" + label + ")");
auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome());
indexer->ScaleAllImages(merge_result);
}
// Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal
// orientation tilt against the running merge and recompute each reflection's partiality + scale
// correction (no re-integration). Never for the P1 search pass.
if (!for_search && experiment_.GetScalingSettings().GetStillsPartialityRefine()) {
phase("Partiality post-refine (" + label + ")");
StillsPartialityRefine refiner(experiment_);
const double mean_tilt = refiner.Run(indexer->GetIntegrationOutcome(), config_.nthreads);
logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt);
// The refine fits its OWN per-image scale, with no floor, and adopts it unconditionally
// when the image has no prior CC - which is exactly the state the guard above leaves. So
// the guard was protecting nothing: measured on a jet dataset, of 367 images it left
// unscaled only 8 were still unscaled in the output, and 53 reached the merge at or below
// a fiftieth of the run median, the worst at a 4525th. Whatever produced the final scale,
// it is the final scale that has to be credible.
ScaleOnTheFly::RejectCollapsedScales(indexer->GetIntegrationOutcome());
}
const std::vector<IntegrationOutcome> &merge_input = indexer->GetIntegrationOutcome();
phase("Merging");
MergeOnTheFly merge_engine(experiment_);
// For the de-novo P1 pass (for_search) drop the ice-ring reflections from the merged
// intensities and the error model, so the space-group search sees clean data; the final
// in-symmetry merge keeps them so completeness is not lost.
merge_engine.ExcludeIceRings(for_search);
if (result.consensus_cell.has_value())
merge_engine.ReferenceCell(*result.consensus_cell);
// Merge-consistency filter: on floods with many spurious lattices (e.g. XFEL large cells)
// most indexed crystals do not correlate with the true structure; --min-image-cc drops the
// crystals whose per-image CC to the reference is below the limit, so the merge keeps only
// the coherent (real) ones. Set before the error model so the model, the merge and the
// reported statistics are all fitted over the same set of images. The P1 search pass keeps
// every image: its job is to find the symmetry, not to produce final intensities.
merge_engine.FilterByImageCC(!for_search
&& experiment_.GetScalingSettings().GetMinCCForImage() > 0.0);
merge_engine.RefineErrorModel(merge_input);
if (merge_engine.ErrorModelActive())
{
// XDS convention, so the pair reads directly against a CORRECT.LP; the stills ISa is
// already the whole-range 1/b, there being no asymptote tier on this path.
const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB());
logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} chi2={:.2f}",
em.a, em.b, em.isa, merge_engine.ErrorModelChi2());
}
for (size_t i = 0; i < merge_input.size(); ++i)
merge_engine.AddImage(merge_input[i], static_cast<int64_t>(i));
ScaleMergeResult out;
out.merged = merge_engine.ExportReflections();
// Automatic high-resolution cutoff (post-merge): trim the written reflections + reported
// shells to the CC1/2 fall-off. The merge, scaling and error model above ran over the full
// range, and the _process.h5 is written from the per-image outcomes, so no data is lost. A
// manual --scaling-high-resolution wins; the P1 search merge (for_search) is never cut so the
// space-group search still sees the full range. (Rotation is cut inside RotationScaleMerge.)
const auto &cut_ss = experiment_.GetScalingSettings();
const std::optional<double> effective_d_min = ApplyResolutionCutoff(
out.merged, cut_ss.GetHighResolutionLimit_A(), cut_ss.GetResolutionCutoff(),
cut_ss.GetResolutionCCTarget(), for_search, logger);
phase("Computing statistics");
out.statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data, effective_d_min);
{ // stills: no asymptote tier, and its ISa already IS the whole-range 1/b
const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB());
result.error_model_isa = em.isa;
result.error_model_a = em.a;
result.error_model_b = em.b;
}
logger.Info("Merge complete ({} unique reflections, {})", out.merged.size(), label);
return out;
};
// First pass: P1 when searching, or directly the user-fixed space group.
const auto initial_sg = experiment_.GetGemmiSpaceGroup();
auto sm = scale_and_merge(initial_sg ? initial_sg->short_name() : "P1", !initial_sg.has_value());
// Rotation two-pass geometry pre-pass: this first scale/merge has now fitted a frame-order-smoothed
// mosaicity and written it back onto the per-frame outcomes. Capture it for the second pass's Bragg
// prediction (so prediction uses the smoothed value, not a fresh per-image estimate) and post-refine +
// apply the detector geometry from this integration. The pre-pass then CONTINUES through the space-
// group search (below) so the second pass can reuse the same space group and lattice - re-searching at
// the slightly changed geometry could flip a borderline determination and mix the two passes.
if (geometry_prepass && rsm && experiment_.IsRotationIndexing()) {
const auto &outcomes = indexer->GetIntegrationOutcome();
prepass_mosaicity_.assign(outcomes.size(), NAN);
for (size_t o = 0; o < outcomes.size(); ++o)
if (outcomes[o].mosaicity_deg) prepass_mosaicity_[o] = *outcomes[o].mosaicity_deg;
if (result.consensus_cell.has_value()) {
if (const auto rot = indexer->FinalizeRotationIndexing(); rot && rot->axis) {
prepass_result_ = *rot; // full result (lattice+geom+axis), for the supercell fallback
PostRefineSettings prs;
prs.crystal_system = rot->search_result.system;
prs.num_threads = static_cast<int>(config_.nthreads);
prs.refine_geometry = true;
const auto pr = PostRefineRotationGeometry(outcomes, *rot->axis,
experiment_.GetDiffractionGeometry(), rot->lattice, prs, logger);
result.post_refine = pr;
if (pr.ok && pr.detector_refined) {
// Store, don't apply: this pass continues to write its own (header-geometry) output;
// Run() applies the refined geometry before the second pass.
prepass_detector_geometry_ = std::array<float, 3>{
static_cast<float>(pr.beam_x_after_px),
static_cast<float>(pr.beam_y_after_px),
static_cast<float>(pr.distance_after_mm)};
logger.Info("Two-pass: post-refined detector geometry (distance {:.3f} mm, "
"beam {:.2f},{:.2f} px) will drive the second integration pass",
pr.distance_after_mm, pr.beam_x_after_px, pr.beam_y_after_px);
} else {
logger.Info("Two-pass: geometry post-refine committed no detector change "
"- the second pass reproduces the first");
}
// DECISION POINT for the goniometer rotation scale. It does not ride on pr.ok - the
// scale is its own cross-validated fit and the crystals that have a stage fault are
// exactly the ones whose cell and detector steps do NOT pass, because the angle error
// is what their residual is made of. A calibration fault is rare (36 of 37 rotation
// datasets sit at 1.0000) and applying a 1 % angle correction to a healthy dataset
// would silently damage it, so the asymmetry is deliberate: committed only when the
// fit is both cross-validated and outside the tolerance. A manual --rotation-scale is
// already on the goniometer and is left alone.
if (pr.rotation_scale_suspect && !config_.rotation_scale.has_value())
prepass_rotation_scale_ = static_cast<float>(pr.rotation_scale);
}
}
}
// For a de-novo search, gate the P1 scaling+merge (the pass that feeds the space-group search) at
// <I/sigma> >= 1 over thin resolution shells: noise-dominated high-res shells otherwise mislead
// the determination, where a single artifact (a zinger the box sum integrates whole, or a
// weak-reflection profile-fit runaway) has an astronomical resolution-normalised E and wrecks the
// intensity-correlation / second-moment statistics. The final in-symmetry merge keeps the full
// range (only this P1 search pass is cut); a manual --scaling-high-resolution, if coarser, wins.
double d_min_search = 0.0;
if (!initial_sg.has_value()) {
std::vector<std::pair<float, float>> rs; // (d, I/sigma) over the full P1 merge
rs.reserve(sm.merged.size());
for (const auto &m : sm.merged)
if (std::isfinite(m.I) && std::isfinite(m.sigma) && m.sigma > 0.0f &&
std::isfinite(m.d) && m.d > 0.0f)
rs.emplace_back(m.d, m.I / m.sigma);
if (rs.size() >= 400) {
std::sort(rs.begin(), rs.end(),
[](const auto &a, const auto &b) { return a.first > b.first; }); // low -> high res
const int bins = std::clamp(static_cast<int>(rs.size() / 100), 1, 40);
const size_t per = (rs.size() + bins - 1) / static_cast<size_t>(bins);
for (size_t b = 0; b * per < rs.size(); ++b) {
const size_t lo = b * per, hi = std::min(rs.size(), lo + per);
double sum = 0.0;
for (size_t j = lo; j < hi; ++j) sum += rs[j].second;
if (sum / static_cast<double>(hi - lo) < 1.0) {
// Cut the noise-dominated high-res shells. When even the lowest-res shell fails,
// keep that shell alone rather than abandoning the cut. The bound is absolute
// while the merged I/sigma it tests is capped by the merge's own ISa, so a merge
// whose ISa has fallen below 1 cannot satisfy it in ANY shell - and that is
// exactly the merge the search must not be handed whole, since the noise which
// drove ISa down is what the cut exists to remove. Abandoning it feeds back:
// more noise into the search, a weaker operator statistic, and a lost symmetry.
if (b > 0)
d_min_search = rs[lo].first;
else if (per >= 400)
d_min_search = rs[per - 1].first;
break;
}
}
}
const double d_manual = experiment_.GetScalingSettings().GetHighResolutionLimit_A().value_or(0.0);
if (rsm && d_min_search > d_manual) {
logger.Info("Space-group search: re-scaling the P1 pass at {:.2f} A (<I/sigma> >= 1)", d_min_search);
rsm->SetDMinLimit(d_min_search);
sm = scale_and_merge("P1, space-group search", true);
rsm->SetDMinLimit(d_manual > 0.0 ? std::optional<double>(d_manual) : std::nullopt);
}
}
std::ostringstream stats_text;
// True when the point group below was chosen BY THE SEARCH (rather than given by the user) and is
// above triclinic - i.e. a promotion was made, which is what makes a later "the Laue class is
// holohedral so there is no twin law" conclusion circular.
// Pass 2 reuses pass-1's group instead of searching, so the promotion it needs to report was
// made there; a first / single pass has no carried-over group and decides this below.
bool promoted_point_group = prepass_merge_sg_.has_value() && prepass_promoted_point_group_;
if (!experiment_.GetGemmiSpaceGroup().has_value()) {
SearchSpaceGroupOptions sg_opts;
sg_opts.merge_friedel = experiment_.GetScalingSettings().GetMergeFriedel();
sg_opts.d_min_limit_A = std::max<double>(
d_min_search, experiment_.GetScalingSettings().GetHighResolutionLimit_A().value_or(0.0));
// Constrain the search to subgroups of the lattice (metric) symmetry found by rotation
// indexing. Centering is not constrained here - it is determined from the absences.
if (end_msg.rotation_lattice_type.has_value())
sg_opts.lattice_system = end_msg.rotation_lattice_type->crystal_system;
auto sg_search = SearchSpaceGroup(sm.merged, sg_opts);
// Second opinion from a merge that keeps only well-measured observations (see
// RotationScaleMerge::search_min_zeta). Where the two disagree the all-observation merge
// decides. The filter earns its place because a reflection that crosses the Ewald sphere
// near-tangentially is measured worst, and a handful of bad measurements can make a real
// operator look like a twin law; but it discards 40-80% of the observations, which starves
// the operator correlations. Both failure directions are real - the filtered arm has been
// seen to lose an operator the full merge confirms, and, at geometries away from the ones
// this comparison was first checked on, to confirm one the full merge refuses. Only one of
// the two has every observation behind it, so the filtered merge reports and the
// all-observation merge decides. Both searches keep their vetoes, so a merohedral twin
// stays refused in each (verified on a real twin at the three frame ranges where it
// over-promotes: refused in both arms every time).
if (rsm && rsm->GetSearchMinZeta() > 0.0 && !sg_search.point_group_hm.empty()) {
const double zeta = rsm->GetSearchMinZeta();
rsm->SetSearchMinZeta(0.0);
auto sm_all = scale_and_merge("P1, all observations", true);
rsm->SetSearchMinZeta(zeta);
const auto alt = SearchSpaceGroup(sm_all.merged, sg_opts);
// The order of the point group each arm confirmed. Read from the search, not from the
// space group it went on to pick: Stage B leaves best_space_group unset when no
// candidate has enough absences to be eligible, which is common on the filtered arm
// precisely because the filter removes the weak reflections the absences live in. An
// arm that confirmed 422 would then have scored 0 and lost to a P2 arm - a demotion,
// reported as "taking the higher symmetry".
const auto order_of = [](const SearchSpaceGroupResult &r) { return r.point_group_order; };
logger.Info("Space-group search: all-observation merge -> {} (order {}), "
"Lorentz-filtered -> {} (order {})",
alt.point_group_hm.empty() ? "?" : alt.point_group_hm, order_of(alt),
sg_search.point_group_hm.empty() ? "?" : sg_search.point_group_hm,
order_of(sg_search));
// The screw/centering determination in particular comes from the merge with all the
// observations, because systematic absences are decided by the WEAK reflections and the
// filter throws most of them away. Preferring the filtered arm on a tie cost four
// crystals their screw axes (P2(1) read as P2) for exactly that reason, with the point
// group and every intensity statistic identical.
// Checked while both arms are still intact - below, the all-observation arm is adopted
// by value and the two become the same result.
// Which point group, by its symmorphic representative rather than by name - gemmi calls
// both P321 and P312 "32", so two arms that disagree about WHICH 2-folds are real read
// as agreeing. 0 when no point group was chosen.
const auto pg_of = [](const SearchSpaceGroupResult &r) {
return r.point_group_representative.has_value() ? r.point_group_representative->number : 0;
};
if (order_of(alt) > 0 && order_of(alt) == order_of(sg_search)
&& pg_of(alt) != pg_of(sg_search)) {
// Same order, different symmetry: the two merges disagree about WHICH operators are
// real, and neither is higher, so there is nothing to prefer. Say so instead of
// picking silently - the two imply different molecular-replacement searches, and
// trying both is cheap next to processing the data again.
const std::string a = alt.best_space_group.has_value()
? alt.best_space_group->short_name() : alt.point_group_hm;
const std::string b = sg_search.best_space_group.has_value()
? sg_search.best_space_group->short_name() : sg_search.point_group_hm;
const std::string msg = fmt::format(
"Space group is AMBIGUOUS between {} and {} (point groups {} and {}, the same "
"order): the merge of all observations and the merge of only the well-measured "
"ones each support one of them, and neither is the higher symmetry, so the data "
"do not decide. Processing continues in {} (the all-observation choice) - try "
"BOTH in molecular replacement, or re-run with -S {} to force the other one.",
a, b, alt.point_group_hm, sg_search.point_group_hm, a, b);
logger.Warning("{}", msg);
stats_text << " !! " << msg << "\n\n";
result.warnings.push_back(msg);
}
if (order_of(alt) > order_of(sg_search))
logger.Info("Space-group search: all-observation merge supports {} where the "
"Lorentz-filtered one supports {} - taking the higher symmetry",
alt.point_group_hm, sg_search.point_group_hm);
else if (order_of(alt) < order_of(sg_search))
logger.Info("Space-group search: Lorentz-filtered merge supports {} where the "
"all-observation one supports {} - keeping the all-observation "
"answer, which has every observation behind it",
sg_search.point_group_hm, alt.point_group_hm);
sg_search = alt;
// The all-observation merge is the one carried forward: it is what the
// pre-promotion twinning numbers below are measured on, and the filtered merge is the
// second opinion reported above.
sm = std::move(sm_all);
}
// Twinning evidence measured BEFORE any promotion, on the P1/subgroup merge the search was
// given. Reported alongside the post-adoption analysis: once a higher Laue class has been
// adopted, its own twinning test can only say "no twin law exists within this class", which
// is circular - the promotion is what a twin would have caused. These numbers are not.
const auto pre_promotion_twinning = AnalyzeTwinning(sm.merged, nullptr);
if (pre_promotion_twinning.twinning_suspected)
logger.Warning("Twinning indicators BEFORE the space-group decision (subgroup merge): "
"<|L|> = {:.3f}, second moment = {:.3f} (untwinned 0.500 / 2.00) - if the "
"search promoted the point group, treat the promotion with suspicion",
pre_promotion_twinning.mean_abs_l, pre_promotion_twinning.second_moment);
promoted_point_group = !sg_search.point_group_hm.empty() && sg_search.point_group_hm != "1";
if (!sg_search.refused_point_group_hm.empty())
logger.Warning("Higher symmetry {} was confirmed by the operator correlations but refused: "
"{}. Processing in the lower symmetry (the recoverable direction).",
sg_search.refused_point_group_hm, sg_search.refused_reason);
// Miller-index reindex under a change of basis: (h,k,l)_conv = reindex * (h,k,l)_prim.
const auto reindex_hkl = [](auto &r, const gemmi::Mat33 &m) {
const double h = r.h, k = r.k, l = r.l;
r.h = static_cast<int32_t>(std::lround(m.a[0][0] * h + m.a[0][1] * k + m.a[0][2] * l));
r.k = static_cast<int32_t>(std::lround(m.a[1][0] * h + m.a[1][1] * k + m.a[1][2] * l));
r.l = static_cast<int32_t>(std::lround(m.a[2][0] * h + m.a[2][1] * k + m.a[2][2] * l));
};
// Intensity-based centred-lattice test (XDS CORRECT / POINTLESS style). When the indexer
// demoted to a primitive/triclinic cell but the metric is pseudo-symmetric, the point-group
// search above tested no operators (its candidate set is bounded by the indexer's lattice).
// A genuinely centred lattice and a merely pseudo-symmetric triclinic one are indistinguishable
// from spot POSITIONS - the 2-fold is an INTENSITY property - so decide it here, after P1
// integration: reindex the P1-merged data to the metric candidate's conventional setting and
// re-run the search. The I(h) vs I(Rh) correlation confirms a real centred lattice (high CC ->
// commit the centring) and rejects a pseudo-symmetric triclinic (low CC -> stays P1, no reindex).
std::optional<gemmi::Mat33> commit_reindex;
std::optional<UnitCell> commit_cell;
std::optional<CrystalLattice> commit_lattice;
std::optional<LatticeMessage> commit_lattice_type;
// Only a demoted (triclinic) indexing result reaches here with a primitive cell - that is both
// the case the point-group search could not test (a triclinic holohedry offers no operators)
// and the case where LatticeSearch's Niggli reduction is valid (it assumes a primitive input).
const bool searched_trivial = !sg_search.best_space_group || sg_search.best_space_group->number == 1;
const bool demoted_triclinic = end_msg.rotation_lattice_type.has_value()
&& end_msg.rotation_lattice_type->crystal_system == gemmi::CrystalSystem::Triclinic;
if (searched_trivial && demoted_triclinic && end_msg.rotation_lattice.has_value()) {
// The metric only generates the candidate; the intensity CC gate below makes the decision
// (a false pseudo-symmetry is rejected by the correlation).
const auto cand = LatticeSearch(*end_msg.rotation_lattice);
// Exact integer reindex from the indexed setting straight to the candidate conventional
// setting: P[i][j] = conv_real[i] . rot_reciprocal[j]. Derived from the two lattices, so it
// holds even if the refined indexed cell is not exactly LatticeSearch's Niggli cell (whose
// primitive_reduced->conventional reindex would otherwise be applied in the wrong setting).
const Coord cv[3] = {cand.conventional.Vec0(), cand.conventional.Vec1(), cand.conventional.Vec2()};
const Coord rs[3] = {end_msg.rotation_lattice->Astar(), end_msg.rotation_lattice->Bstar(),
end_msg.rotation_lattice->Cstar()};
gemmi::Mat33 reindex;
double reindex_res = 0.0;
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j) {
const double v = cv[i] * rs[j];
reindex.a[i][j] = std::round(v);
reindex_res = std::max(reindex_res, std::fabs(v - reindex.a[i][j]));
}
if (cand.system != gemmi::CrystalSystem::Triclinic && reindex_res < 0.1) {
std::vector<MergedReflection> merged_c = sm.merged;
for (auto &m : merged_c)
reindex_hkl(m, reindex);
SearchSpaceGroupOptions o2 = sg_opts;
o2.lattice_system = cand.system;
const auto s2 = SearchSpaceGroup(merged_c, o2);
if (s2.best_space_group.has_value() && s2.best_space_group->number > 1) {
// The centring cannot be confirmed from absences here (integrated in the primitive
// cell, so the centring-absent reflections do not exist) - it is metric-determined.
// Among the point-group-equivalent candidates the intensities cannot separate
// (e.g. P2 / P21 / C2), pick the one matching the metric candidate's centring.
gemmi::SpaceGroup chosen = *s2.best_space_group;
if (chosen.centring_type() != cand.centering)
for (const auto &alt : s2.alternatives)
if (alt.centring_type() == cand.centering) { chosen = alt; break; }
// Commit only if the adopted centring matches the metric candidate whose conventional
// cell we reindex into - otherwise the space group and the written cell disagree.
if (chosen.centring_type() == cand.centering) {
sg_search = s2;
// Swap the metric-matching candidate into the answer and leave the one it
// displaced among the alternatives, so the reported "or ..." list stays the
// set of groups the data cannot separate, each named once.
if (chosen.number != s2.best_space_group->number) {
sg_search.alternatives.push_back(*s2.best_space_group);
std::erase_if(sg_search.alternatives, [&](const gemmi::SpaceGroup &a) {
return a.number == chosen.number;
});
}
sg_search.best_space_group = chosen;
commit_reindex = reindex;
commit_cell = cand.conventional.GetUnitCell();
commit_lattice = end_msg.rotation_lattice->Multiply(reindex);
commit_lattice_type = LatticeMessage{ .centering = cand.centering,
.niggli_class = end_msg.rotation_lattice_type->niggli_class,
.crystal_system = cand.system };
}
}
}
}
// Adopt the determined space group and re-scale/merge in it, so scaling uses symmetry
// equivalents and the statistics come out in the right symmetry. P1 stands when nothing
// is determined.
if (sg_search.best_space_group.has_value()) {
const auto &sg = *sg_search.best_space_group;
logger.Info("Adopting space group {} (number {})", sg.short_name(), sg.number);
// A reindex is committed only when the intensity re-test confirmed a higher symmetry in a
// centred setting: bring the integrated reflections + cell into that setting and re-ingest
// the rotation merge, so the final scaling/merging folds the equivalents correctly.
if (commit_reindex) {
for (auto &io : indexer->GetIntegrationOutcome()) {
io.latt = io.latt.Multiply(*commit_reindex);
for (auto &r : io.reflections)
reindex_hkl(r, *commit_reindex);
}
result.consensus_cell = commit_cell;
end_msg.unit_cell = commit_cell;
// Keep the lattice metadata (hence the master-file UB matrix and indexed-lattice
// vectors) in the same conventional setting as the cell, reflections and space group.
end_msg.rotation_lattice = commit_lattice;
end_msg.rotation_lattice_type = commit_lattice_type;
if (rsm) {
rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
static_cast<int>(config_.scaling_iter),
config_.nthreads, logger, config_.observation_dump_path);
rsm->Ingest();
}
}
experiment_.SpaceGroupNumber(sg.number);
end_msg.space_group_number = sg.number;
result.space_group_number = sg.number;
phase("Re-scaling in space group " + sg.short_name());
sm = scale_and_merge(sg.short_name(), false);
}
result.space_group_search = sg_search;
} else {
// A space group was fixed by the user; surface it so the viewer/CLI can still show it.
result.space_group_number = experiment_.GetSpaceGroupNumber();
}
// Rotation two-pass geometry pre-pass: the determined space group is now fixed on experiment_, so the
// second pass REUSES it instead of re-searching at the slightly changed geometry (which could flip a
// borderline determination and mix the two passes). The lattice is deliberately NOT forced - the second
// pass re-indexes at the refined geometry so the CELL comes out self-consistent with it (forcing the
// pre-pass lattice would pin the nominal cell and undo the geometry refinement). Unlike a throwaway
// pre-pass, this one now CONTINUES to the write below so its own (header-geometry) result is saved
// as the first pass.
// Reference-based indexing-ambiguity resolution (rotation). When a reference MTZ is supplied and
// the crystal has an indexing ambiguity (merohedral: lattice symmetry higher than the Laue group),
// pick the reindexing that best correlates with the reference intensities and re-merge in it. The
// twin-law reindex is metric-preserving, so only the hkl labels change (the cell is unchanged).
// Stills resolve the ambiguity per image with ReindexAmbiguityResolver, not here.
if (rsm && !config_.reference_data.empty() && result.consensus_cell
&& experiment_.GetSpaceGroupNumber().has_value()) {
const int sg_num = static_cast<int>(*experiment_.GetSpaceGroupNumber());
const auto choice = ChooseReindex(
sm.merged, *result.consensus_cell, sg_num,
[&](const std::vector<MergedReflection> &m) {
return ReferenceIntensityCC(m, config_.reference_data, sg_num);
});
if (!choice.is_identity) {
logger.Info("Reference: resolved indexing ambiguity by reindexing to match the reference "
"(CC {:.3f}, vs {:.3f} unchanged)", choice.score, choice.identity_score);
for (auto &io : indexer->GetIntegrationOutcome())
for (auto &r : io.reflections) {
const gemmi::Op::Miller h = choice.op.apply_to_hkl({{r.h, r.k, r.l}});
r.h = h[0]; r.k = h[1]; r.l = h[2];
}
rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
static_cast<int>(config_.scaling_iter),
config_.nthreads, logger, config_.observation_dump_path);
rsm->Ingest();
phase("Re-merging in the reference indexing");
sm = scale_and_merge(experiment_.GetGemmiSpaceGroup()->short_name(), false);
}
}
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;
result.twinning = AnalyzeTwinning(sm.merged, twin_sg);
// Mark the conclusion as non-authoritative when the Laue class was reached by a promotion the
// search itself made, so the text cannot claim "no twin law exists" on its own say-so.
result.twinning.laue_class_was_chosen_by_promotion = promoted_point_group;
stats_text << TwinningAnalysisToText(result.twinning) << "\n";
// Symmetry axis vs spindle. Reported always on rotation data (it is a property of how the
// crystal was mounted, which the user can change), warned about when the two nearly coincide.
if (experiment_.IsRotationIndexing() && twin_sg && end_msg.rotation_lattice.has_value()) {
if (const auto gonio = experiment_.GetGoniometer()) {
const auto closest = ClosestSymmetryAxisToSpindle(*twin_sg, *end_msg.rotation_lattice,
gonio->GetAxis());
if (closest.has_value()) {
// Practical bound, not a derived one: the blind cone's half-angle is the maximum
// Bragg angle (~15 deg for 2 A data at 1 A), and measured on this battery a 13.6 deg
// case shows the loss while a 16.2 deg one is 99.7% complete.
constexpr double WARN_DEG = 15.0;
const auto [angle, order] = *closest;
if (angle < WARN_DEG) {
const std::string msg = fmt::format(
"The crystal's {}-fold axis is only {:.1f} deg from the spindle. A rotation sweep "
"never records the reflections whose reciprocal vector lies within the Bragg angle "
"of the spindle, and symmetry normally supplies them from an equivalent elsewhere; "
"it cannot here, because that axis maps the blind region onto itself. Those "
"reflections stay missing however long the sweep runs. The loss is confined to "
"that cone rather than spread over the data, so overall completeness may still "
"look reasonable - check it near the spindle direction. A second sweep on a "
"different axis, or re-mounting, recovers them.",
order, angle);
logger.Warning("{}", msg);
stats_text << " !! " << msg << "\n\n";
result.warnings.push_back(msg);
} else {
stats_text << "Closest symmetry axis to the spindle: " << order << "-fold at "
<< std::fixed << std::setprecision(1) << angle << " deg\n\n";
}
}
}
}
// Indexing-ambiguity (alternative-indexing) advisory. When the lattice metric symmetry exceeds
// the Laue symmetry the crystal can be validly indexed in several hands related by twin-law
// operators. For an obvious merohedral case (P3/P4/P6...) users expect this; but a PSEUDO-merohedral
// metric (e.g. a C2 crystal whose beta makes it pseudo-F-orthorhombic) is easy to miss, so surface it. The
// reindex operator is applied directly to (h,k,l), so its triplet reads as an h,k,l transform.
if (result.consensus_cell && twin_sg_number) {
const auto twin_ops = ReindexAmbiguityOperators(*result.consensus_cell,
static_cast<int>(*twin_sg_number), 2.0);
if (!twin_ops.empty()) {
std::string laws;
for (const auto &op : twin_ops) {
std::string t = op.triplet();
std::replace(t.begin(), t.end(), 'x', 'h');
std::replace(t.begin(), t.end(), 'y', 'k');
std::replace(t.begin(), t.end(), 'z', 'l');
laws += (laws.empty() ? "" : " ; ") + t;
}
if (config_.reference_data.empty()) {
const std::string msg = fmt::format(
"Indexing ambiguity: this cell / space group admits alternative indexing "
"(reindex operator(s): {}). {} are indexed in one hand at random, and rugnux can only "
"break this against an external reference. WITHOUT one the merge mixes the hands and "
"CC1/2 is degraded - supply a reference (-z reference.mtz or --model) to resolve it.",
laws, experiment_.IsRotationIndexing() ? "Lattices" : "Serial-stills crystals");
logger.Warning("{}", msg);
stats_text << " !! " << msg << "\n\n";
result.warnings.push_back(msg);
} else {
logger.Info("Indexing ambiguity present (reindex operator(s): {}); resolved against the "
"supplied reference.", laws);
}
}
}
// Dataset-wide Wilson B-factor estimate (like XDS's WILSON LINE B). Diagnostic only - it is not
// fed back into scaling; it just lands in the printed statistics, the mmCIF, and the log.
{
const GlobalWilsonB wilson = CalcGlobalWilsonB(sm.merged);
sm.statistics.wilson_b = wilson.b;
sm.statistics.wilson_b_correlation = wilson.correlation;
if (std::isfinite(wilson.b) && wilson.b > 0.0)
logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells)",
wilson.b, wilson.correlation, wilson.n_shells);
}
stats_text << sm.statistics;
result.merge_statistics_text = stats_text.str();
result.has_merge_statistics = true;
result.merge_statistics = sm.statistics;
// Per-image form of the sweep-quality ranges, for the _process.h5: one code per image, 0 where
// the image is in no flagged range, plus the vocabulary the codes index. Only filled when the
// diagnostic ran, so absent datasets mean "not looked for" rather than "all clean".
if (sm.statistics.sweep_quality.measured) {
end_msg.sweep_quality.assign(end_msg.max_image_number, 0);
for (const auto &r : sm.statistics.sweep_quality.ranges)
for (int64_t i = std::max<int64_t>(0, r.first_image);
i <= r.last_image && i < static_cast<int64_t>(end_msg.sweep_quality.size()); ++i)
end_msg.sweep_quality[i] = static_cast<uint8_t>(r.reason) + 1;
for (int r = 0; r <= static_cast<int>(SweepQualityReason::RadiationDamage); ++r)
end_msg.sweep_quality_reasons.emplace_back(
SweepQualityReasonCode(static_cast<SweepQualityReason>(r)));
}
{
// Stride rather than take the head: the merged list is ordered by hkl, so the first N
// reflections are one corner of reciprocal space and would not show the intensity range.
constexpr size_t MAX_POINTS = 4000;
const size_t n = sm.merged.size();
const size_t step = std::max<size_t>(1, n / MAX_POINTS);
result.merged_i_sigma.reserve(std::min(n, MAX_POINTS + 1));
for (size_t i = 0; i < n; i += step)
if (std::isfinite(sm.merged[i].I) && std::isfinite(sm.merged[i].sigma)
&& sm.merged[i].sigma > 0.0f)
result.merged_i_sigma.emplace_back(sm.merged[i].I, sm.merged[i].sigma);
}
result.has_reference = !config_.reference_data.empty();
// Inherit the campaign's shared R-free test set from the reference MTZ (overriding the
// per-hkl hash that the merge assigned), so every dataset flags the same free reflections.
if (config_.reference_has_free_flags && !config_.reference_data.empty() && !sm.merged.empty()) {
const auto sg = experiment_.GetSpaceGroupNumber().value_or(1);
const size_t matched = ApplyReferenceFreeFlags(sm.merged, static_cast<int32_t>(sg),
config_.reference_data);
logger.Info("R-free flags: inherited the reference test set ({} of {} merged reflections matched)",
matched, sm.merged.size());
}
// Radiation-damage report (rotation): the per-image scale CC-to-merge and mosaicity across the
// sweep (the per-image scaling the rotation merge already fits, binned by frame = dose). A
// per-image CC that falls, and/or a mosaicity that rises, with frame number is the classic
// radiation-damage signature - a data-quality-vs-dose read complementary to the fitted decay
// correction. The full per-image table is written to <prefix>_scaling.txt for detail.
if (experiment_.IsRotationIndexing()) {
const auto &outs = indexer->GetIntegrationOutcome();
const int nf = static_cast<int>(outs.size());
constexpr int nb = 10;
std::array<double, nb> cc_sum{}, mos_sum{};
std::array<int, nb> cc_n{}, mos_n{};
for (int f = 0; f < nf; ++f) {
const int b = nf > 1 ? std::min(nb - 1, f * nb / nf) : 0;
if (outs[f].image_scale_cc && std::isfinite(*outs[f].image_scale_cc)) {
cc_sum[b] += *outs[f].image_scale_cc; cc_n[b]++;
}
if (outs[f].mosaicity_deg && std::isfinite(*outs[f].mosaicity_deg)) {
mos_sum[b] += *outs[f].mosaicity_deg; mos_n[b]++;
}
}
int tot = 0;
for (int n : cc_n) tot += n;
if (tot >= 20) {
std::ostringstream os;
os << fmt::format("Radiation-damage report (per-image scale over {} frames):\n", nf);
os << " dose <CC to merge> <mosaicity deg>\n";
for (int b = 0; b < nb; ++b) {
if (cc_n[b] == 0 && mos_n[b] == 0) continue;
os << fmt::format(" {:3d}-{:3d}% {:>7} {:>7}\n", b * 10, (b + 1) * 10,
cc_n[b] ? fmt::format("{:.3f}", cc_sum[b] / cc_n[b]) : std::string("-"),
mos_n[b] ? fmt::format("{:.3f}", mos_sum[b] / mos_n[b]) : std::string("-"));
}
if (cc_n[0] > 0 && cc_n[nb - 1] > 0) {
const double cc0 = cc_sum[0] / cc_n[0], cc9 = cc_sum[nb - 1] / cc_n[nb - 1];
os << fmt::format(" => per-image CC to merge {:.3f} (first 10%) -> {:.3f} (last 10%){}",
cc0, cc9,
cc9 < cc0 - 0.05 ? " (falling: radiation damage / crystal decay)" : "");
}
// Relative B-factor across the sweep (measured before any decay correction, against the
// low-dose start): the resolution-dependent complement to the CC/mosaicity read above.
// Radiation damage fades the high-resolution intensity, so only a POSITIVE change is dose.
// A curve no straight line describes gets no first->last number at all, and a batch whose
// data could not measure one prints "-" rather than a value.
if (!sm.statistics.radiation_damage_b_batch.empty()) {
if (std::isfinite(sm.statistics.radiation_damage_delta_b))
os << fmt::format("\n => relative B-factor change over run = {:+.2f} A^2 (first->last){}",
sm.statistics.radiation_damage_delta_b,
sm.statistics.radiation_damage_delta_b > 5.0
? " (significant dose-dependent scaling: radiation damage)" : "");
else
os << "\n => relative B-factor: the per-batch curve below is not a trend, so no "
"first->last number describes it - dose does not come back, so whatever moved "
"here was not dose; see the sweep-quality report";
os << fmt::format("\n per-batch relative-B (A^2, {:.0f} deg/batch):",
sm.statistics.radiation_damage_batch_deg);
for (float bb : sm.statistics.radiation_damage_b_batch)
os << (std::isfinite(bb) ? fmt::format(" {:.1f}", bb) : std::string(" -"));
}
logger.Info("{}", os.str());
result.radiation_damage_text = os.str();
}
// Sweep-quality report: the stretches of the sweep over which the crystal delivered much
// less than the rest of the run, and what each one looks like. Nothing is excluded because
// of it - the frames still carry signal, and this is a message for the beamline, not a
// filter. Frame numbers are processed-image ordinals, inclusive, as in <prefix>_image.dat.
const auto &sq = sm.statistics.sweep_quality;
if (sq.measured) {
std::ostringstream os;
os << fmt::format("Sweep quality over {:.0f} deg (per-image scale with the incident flux, "
"which varied {:.2f}x, already divided out):\n",
sq.sweep_deg, sq.flux_peak_to_trough);
if (sq.ranges.empty()) {
os << " no stretch of the sweep is materially worse than the run";
} else {
os << " frames rotation diagnosis severity scale CC indexed\n";
for (const auto &r : sq.ranges)
os << fmt::format(" {:<17s} {:6.1f} deg {:<16s} {:5.2f} {:5.2f} {:5.2f} {:4.0f}%\n",
fmt::format("{}-{}", r.first_image, r.last_image), r.rotation_deg,
SweepQualityReasonText(r.reason), r.severity, r.mean_relative_scale,
r.mean_relative_cc, 100.0 * r.indexed_fraction);
os << " => severity is the fraction of the run's typical diffracting power missing "
"over the range";
if (sq.modulation_peak_to_trough >= 1.05f)
os << fmt::format("\n => once-per-revolution modulation of the per-image scale: "
"{:.1f}x peak to trough", sq.modulation_peak_to_trough);
}
logger.Info("{}", os.str());
}
}
if (result.consensus_cell && write_files && config_.write_merged) {
phase("Writing reflections");
const ErrorModelReport em_report{
result.error_model_isa > 0 ? fmt::format("{:.2f}", result.error_model_isa) : "?",
result.error_model_isa_asymptotic > 0 ? fmt::format("{:.2f}", result.error_model_isa_asymptotic)
: std::string(),
result.error_model_a > 0 ? fmt::format("{:.3f}", result.error_model_a) : std::string(),
result.error_model_b > 0 ? fmt::format("{:.4e}", result.error_model_b) : std::string()};
WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics,
em_report, result.twinning, config_.output_prefix);
// Per-image scaling table (G, B-factor, mosaicity, wedge, CC) for inspection / XDS
// comparison. The offline self-scaling result is otherwise not exposed (process.h5's
// per-image arrays are only filled on the online per-image path). Sourced from the
// partials, which carry the first-pass per-image scale.
ScalingResult(indexer->GetIntegrationOutcome()).SaveToFile(config_.output_prefix);
}
if (result.consensus_cell && write_files && !config_.model_path.empty()) {
phase("Validating against model");
const auto data_sg = experiment_.GetSpaceGroupNumber();
// With a reference MTZ the merohedral indexing was already resolved against it (rotation
// merge / stills scaling), so trust that; only probe indexing by R-free when model-only.
ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path,
config_.output_prefix, logger,
data_sg ? std::optional<int>(static_cast<int>(*data_sg)) : std::nullopt,
/*probe_indexing_ambiguity=*/config_.reference_data.empty());
}
}
// When merging is on, the _process.h5 is skipped because the merged reflections are the wanted
// output and the per-image file is large. If nothing indexed there are no merged reflections, so
// the run would finish successfully having written no file at all. Write one now: the per-image
// messages went past while there was no writer, but the end message still carries the per-image
// scalars (spot counts, background, resolution estimate) and the dataset metadata - which is what
// anyone asking "why did nothing index" needs, and it is small precisely because nothing indexed.
if (!writer && full && write_files && !cancelled_ && result.indexing_rate.value_or(0.0f) <= 0.0f) {
logger.Warning("No image indexed, so there are no merged reflections to write - writing "
"{}_process.h5 with the per-image analysis instead", config_.output_prefix);
writer = std::make_unique<FileWriter>(start_message, /*check_overwrite_at_start=*/true,
/*trusted_path=*/true);
}
if (writer) {
writer->WriteHDF5(end_msg);
writer->Finalize();
result.written_master_path = config_.output_prefix + "_process.h5";
}
if (observer)
observer->OnPhase(cancelled_ ? "Cancelled" : "Done");
// Total wall time for the whole run, including scaling/merging/writing above (not just the
// per-image pass), so the reported rate reflects everything rugnux did.
result.processing_time_s = std::chrono::duration<double>(std::chrono::steady_clock::now() - start_time).count();
if (result.processing_time_s > 0.0) {
result.frame_rate_hz = static_cast<double>(result.images_processed) / result.processing_time_s;
result.throughput_MBs = static_cast<double>(total_uncompressed_bytes) / (result.processing_time_s * 1e6);
}
logger.Info("{} {} images in {:.2f} s ({:.2f} Hz)", cancelled_ ? "Cancelled after" : "Processed",
result.images_processed, result.processing_time_s, result.frame_rate_hz);
// The geometry this pass integrated at, for the report. Taken from experiment_ rather than from the
// caller's copy: the second pass of the rotation two-pass runs at the post-refined geometry.
result.used_beam_x_pxl = experiment_.GetBeamX_pxl();
result.used_beam_y_pxl = experiment_.GetBeamY_pxl();
result.used_distance_mm = experiment_.GetDetectorDistance_mm();
return result;
}
Reference in New Issue View Git Blame Copy Permalink