leonarski_f
a0604ff244
Prediction and partiality read one number, the per-image sigma_M, so a sigma_M that moves takes the integrated reflection population with it and there is no way to ask which of the two uses carries a downstream difference. --prediction-mosaicity fixes the width the prediction window opens to while the partiality keeps using the measured sigma_M, which separates them. Measured with it on a rotation crystal whose lattice search returns two different cells a few tens of microns of detector distance apart: widening the prediction window from the narrower branch's 0.203 deg to the wider branch's 0.272 deg adds 28% more partials and moves the merged statistics by less than half a percent (I/sigma 6.1 -> 6.2, R_meas 12.7 -> 12.6%, ISa 10.7 -> 10.9); narrowing the wide branch the other way removes 26% of its partials and recovers nothing. On a tetragonal reference crystal a 4.7x over-wide window costs 11%. The prediction window is not where a mosaicity difference turns into a merged-data difference - the captured-fraction gate and the partiality weighting downstream absorb a generous window. Diagnostic only; off by default, so nothing changes unless it is asked for.
739 lines
36 KiB
C++
739 lines
36 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
#include <limits>
|
|
#include <cstdlib>
|
|
#include "IndexAndRefine.h"
|
|
|
|
#include "bragg_integration/CalcISigma.h"
|
|
#include "geom_refinement/XtalOptimizer.h"
|
|
#include "indexing/AnalyzeIndexing.h"
|
|
#include "indexing/FFTIndexer.h"
|
|
#include "indexing/MultiLatticeSearch.h"
|
|
#include "lattice_search/LatticeSearch.h"
|
|
#include "scale_merge/ReindexAmbiguity.h"
|
|
#include "scale_merge/ScaleOnTheFly.h"
|
|
|
|
namespace {
|
|
// Iterations the offline geometry refinement is allowed, standing in for the 40 ms the online path
|
|
// spends. Ceres' own default is 50; the per-image problem is small and converges well inside that,
|
|
// so this bounds the pathological case rather than the normal one.
|
|
constexpr int OFFLINE_REFINE_ITERATIONS = 50;
|
|
|
|
// How far the predictor has to walk each index for THIS crystal. The predictor keeps only
|
|
// reflections with |q| <= 1/d_min, and h = a.q for the real-space axis a, so |h| <= a/d_min exactly
|
|
// (Cauchy-Schwarz, equality when a lies along q) - and independently |k| <= b/d_min, |l| <= c/d_min.
|
|
// Each index is therefore bounded by its OWN axis, which is why the limits are per-axis: a single
|
|
// cube would have to be sized for the longest axis and would walk the short ones far past anything
|
|
// the resolution cut can keep. One index of margin covers the rounding.
|
|
int MaxIndexForAxis(float axis_A, float d_min_A) {
|
|
return static_cast<int>(std::ceil(axis_A / d_min_A)) + 1;
|
|
}
|
|
|
|
// The resolution the prediction walks out to: an explicit setting, else as far as the detector
|
|
// reaches. The predictor drops any reflection that misses the detector anyway, so this bound cannot
|
|
// add reflections the geometry does not offer - it only decides how much of the lattice is examined,
|
|
// which is why the detector's own reach is the right default and a fixed number was not. A geometry
|
|
// with no scattering angle at all (no distance, no wavelength) reports 0, which is not a resolution
|
|
// to divide by; nothing can be predicted from it either way.
|
|
float IntegrationDMin_A(const DiffractionExperiment &experiment) {
|
|
if (const auto fixed = experiment.GetBraggIntegrationSettings().GetDMinLimit_A())
|
|
return *fixed;
|
|
const float detector_A = experiment.GetDetectorMaxResolution_A();
|
|
return detector_A > 0.0f ? detector_A : 1.0f;
|
|
}
|
|
|
|
// An explicit setting is enforced as given, on every index - it is one number, deliberately, because
|
|
// it exists to bound the work rather than to describe the crystal. Otherwise the cell decides. The
|
|
// online path always carries a value (the broker bootstraps one and the API can change it), so a
|
|
// live acquisition never has its per-frame cost decided by whichever crystal was mounted.
|
|
void ApplyPredictionRange(BraggPredictionSettings &settings, const DiffractionExperiment &experiment,
|
|
const CrystalLattice &latt) {
|
|
if (const auto fixed = experiment.GetBraggIntegrationSettings().GetMaxHKL()) {
|
|
settings.max_h = settings.max_k = settings.max_l = *fixed;
|
|
return;
|
|
}
|
|
const UnitCell cell = latt.GetUnitCell();
|
|
settings.max_h = MaxIndexForAxis(cell.a, settings.high_res_A);
|
|
settings.max_k = MaxIndexForAxis(cell.b, settings.high_res_A);
|
|
settings.max_l = MaxIndexForAxis(cell.c, settings.high_res_A);
|
|
}
|
|
}
|
|
|
|
IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer,
|
|
bool retain_outcomes, bool real_time)
|
|
: retain_outcomes_(retain_outcomes),
|
|
real_time(real_time),
|
|
experiment(x),
|
|
geom_(x.GetDiffractionGeometry()),
|
|
indexer_(indexer),
|
|
rotation_indexer_counter(x) {
|
|
if (indexer && x.IsRotationIndexing())
|
|
rotation_indexer = std::make_unique<RotationIndexer>(x, *indexer);
|
|
// Only retain the whole-run per-image reflections when a later scaling/merge pass will read them.
|
|
if (retain_outcomes_)
|
|
integration_outcome.resize(x.GetImageNum());
|
|
|
|
mosaicity.resize(x.GetImageNum(), NAN);
|
|
scale_cc.resize(x.GetImageNum(), 0);
|
|
unit_cells.resize(x.GetImageNum());
|
|
}
|
|
|
|
std::optional<float> IndexAndRefine::RotationAngle(int64_t image) const {
|
|
// Mid-exposure rotation angle for image index `image`, matching the angle used for prediction.
|
|
if (const auto g = experiment.GetGoniometer())
|
|
return g->GetAngle_deg(static_cast<float>(image)) + g->GetWedge_deg() / 2.0f;
|
|
return std::nullopt;
|
|
}
|
|
|
|
void IndexAndRefine::AddImageToRotationIndexer(DataMessage &msg) {
|
|
if (rotation_indexer)
|
|
rotation_indexer->ProcessImage(msg.number, msg.spots, RotationAngle(msg.number));
|
|
}
|
|
|
|
IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotation(DataMessage &msg) {
|
|
IndexingOutcome outcome(experiment);
|
|
|
|
if (!rotation_indexer)
|
|
return outcome;
|
|
|
|
auto result = rotation_indexer->GetLattice();
|
|
|
|
if (!result.has_value()) {
|
|
auto rot_cnt = rotation_indexer_counter.Process(msg.number);
|
|
if (rot_cnt.first)
|
|
rotation_indexer->ProcessImage(msg.number, msg.spots, RotationAngle(msg.number));
|
|
if (rot_cnt.second)
|
|
rotation_indexer->RunIndexing();
|
|
result = rotation_indexer->GetLattice();
|
|
}
|
|
|
|
if (result.has_value()) {
|
|
// For rotation indexing, indexing rate is calculated only for frames, where "global" rotation indexing solution was found
|
|
msg.indexing_result = false;
|
|
|
|
// get rotated lattice
|
|
auto gon = result->axis;
|
|
if (gon) {
|
|
const float angle_deg = gon->GetAngle_deg(msg.number) + gon->GetWedge_deg() / 2.0f;
|
|
const auto rot_to_image = gon->GetTransformationAngle(-angle_deg);
|
|
outcome.lattice_candidate = result->lattice.Multiply(rot_to_image);
|
|
|
|
outcome.extra_lattice_candidates.reserve(result->extra_lattices.size());
|
|
for (const auto &el : result->extra_lattices)
|
|
outcome.extra_lattice_candidates.push_back(el.Multiply(rot_to_image));
|
|
}
|
|
|
|
outcome.experiment.BeamX_pxl(result->geom.GetBeamX_pxl())
|
|
.BeamY_pxl(result->geom.GetBeamY_pxl())
|
|
.DetectorDistance_mm(result->geom.GetDetectorDistance_mm())
|
|
.PoniRot1_rad(result->geom.GetPoniRot1_rad())
|
|
.PoniRot2_rad(result->geom.GetPoniRot2_rad())
|
|
.Goniometer(result->axis);
|
|
outcome.symmetry.centering = result->search_result.centering;
|
|
outcome.symmetry.niggli_class = result->search_result.niggli_class;
|
|
outcome.symmetry.crystal_system = result->search_result.system;
|
|
}
|
|
return outcome;
|
|
}
|
|
|
|
IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(DataMessage &msg) {
|
|
auto indexing_start_time = std::chrono::steady_clock::now();
|
|
|
|
IndexingOutcome outcome(experiment);
|
|
|
|
// Seed the indexer with the strongest few spots first and escalate to more only if that fails.
|
|
// On flooded / noisy frames (XFEL, ice) a lean high-quality seed finds the lattice far more
|
|
// reliably than the full spot list, whose many spurious peaks derail the search; on clean frames
|
|
// the lean seed already works, so nothing is lost. The FULL spot list is still used for geometry
|
|
// refinement and integration downstream, so higher-resolution accuracy is preserved. Cost is
|
|
// ~1 indexer call on frames that index cleanly, up to 3 only on the hard ones. msg.spots is
|
|
// already ordered non-ice-first, strongest-first (FilterSpotsByCount), so the prefix IS the seed.
|
|
const float idx_tol = experiment.GetIndexingSettings().GetTolerance();
|
|
const float idx_tol_sq = idx_tol * idx_tol;
|
|
constexpr float SEED_STOP_FRACTION = 0.9f; // seed explained this well -> stop escalating
|
|
IndexerResult indexer_result;
|
|
bool any_executed = false;
|
|
float best_frac = -1.0f;
|
|
// Read per call rather than cached at construction: the run measures for itself whether the
|
|
// crystal has ice, and that verdict lands after this object is built.
|
|
const bool index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings();
|
|
for (size_t seed_cap : {size_t{30}, size_t{80}, std::numeric_limits<size_t>::max()}) {
|
|
std::vector<Coord> recip;
|
|
recip.reserve(std::min<size_t>(seed_cap, msg.spots.size()));
|
|
for (const auto &i: msg.spots) {
|
|
if (index_ice_rings || !i.ice_ring) {
|
|
recip.push_back(i.ReciprocalCoord(geom_));
|
|
if (recip.size() >= seed_cap)
|
|
break;
|
|
}
|
|
}
|
|
auto res = indexer_->Run(experiment, recip);
|
|
any_executed |= res.executed;
|
|
if (!res.lattice.empty()) {
|
|
// Keep the seed the lattice explains the largest FRACTION of: a lean clean seed a good
|
|
// lattice indexes almost fully beats a flooded seed it fits only in small part. This
|
|
// auto-selects the lean seed on noisy frames (XFEL) and the full seed where the extra spots
|
|
// are real signal (weak synchrotron) -- no per-dataset setting.
|
|
const Coord a = res.lattice[0].Vec0(), b = res.lattice[0].Vec1(), c = res.lattice[0].Vec2();
|
|
int n = 0;
|
|
for (const auto &q : recip) {
|
|
const float hf = q * a, kf = q * b, lf = q * c;
|
|
// std::rint, not std::round: rounding half away from zero has to be a libm call, half to
|
|
// even is inlined. Only the squared residual is used, and the rules can differ only at an
|
|
// exact .5, where either leaves |frac| = 0.5 - so the count is the same either way.
|
|
const float dh = hf - std::rint(hf), dk = kf - std::rint(kf), dl = lf - std::rint(lf);
|
|
if (dh * dh + dk * dk + dl * dl < idx_tol_sq) ++n;
|
|
}
|
|
const float frac = recip.empty() ? 0.0f : static_cast<float>(n) / recip.size();
|
|
if (frac > best_frac) { best_frac = frac; indexer_result = std::move(res); }
|
|
// A lattice that already explains nearly the whole seed is kept whatever a larger seed
|
|
// returns: the winner is the highest explained FRACTION, and adding weaker spots almost
|
|
// always lowers it. Stop here - this is what keeps clean frames at one indexer call.
|
|
if (frac >= SEED_STOP_FRACTION)
|
|
break;
|
|
}
|
|
if (recip.size() < seed_cap) // already fed every available spot; a larger cap won't add any
|
|
break;
|
|
}
|
|
|
|
if (any_executed)
|
|
msg.indexing_result = false;
|
|
|
|
if (!indexer_result.lattice.empty()) {
|
|
auto latt = indexer_result.lattice[0];
|
|
if (latt.CalcVolume() > 1.0) {
|
|
auto sg = experiment.GetGemmiSpaceGroup();
|
|
const auto algorithm = experiment.GetIndexingAlgorithm();
|
|
const bool de_novo = (algorithm == IndexingAlgorithmEnum::FFT
|
|
|| algorithm == IndexingAlgorithmEnum::FFTW);
|
|
|
|
// If space group and cell provided => enforce that symmetry in refinement.
|
|
// If not => detect the symmetry from the lattice.
|
|
if (sg && experiment.GetUnitCell()) {
|
|
outcome.symmetry = LatticeMessage{
|
|
.centering = sg->centring_type(),
|
|
.niggli_class = 0,
|
|
.crystal_system = sg->crystal_system()
|
|
};
|
|
// Place every frame's cell in ONE consistent setting for the whole dataset:
|
|
// mixed axis orders (e.g. [78,78,38] vs [38,78,78]) index the same reflection
|
|
// as different HKLs and cannot be merged. LatticeSearch gives the conventional
|
|
// setting when its detected symmetry agrees with the user's space group. On
|
|
// noisy frames it can instead pick an alternative Bravais setting (e.g. the
|
|
// sqrt2 C-centred description of a primitive tetragonal cell,
|
|
// [78,78,38]->[110,111,38]); there:
|
|
// - FFBIDX already returns the reference setting (c-last), consistent with the
|
|
// conventional frames, so its raw lattice is safe -> use it (FFBIDX neutral);
|
|
// - de-novo indexers (FFT/FFTW) return a Niggli-primitive cell with a DIFFERENT
|
|
// axis order (c-first) that would corrupt the merge -> reject the frame.
|
|
// niggli_class is left unassigned (0): it needs the primitive cell incl.
|
|
// centering, which LatticeSearch cannot recover from a (possibly centred, e.g.
|
|
// C2) user cell. A proper primitive-cell indexing path (CrystFEL-style) is deferred.
|
|
auto sym_result = LatticeSearch(latt);
|
|
if (sym_result.system == sg->crystal_system())
|
|
outcome.lattice_candidate = sym_result.conventional;
|
|
else if (!de_novo)
|
|
outcome.lattice_candidate = latt;
|
|
// else: de-novo + symmetry mismatch -> leave unset, frame is not indexed
|
|
} else {
|
|
auto sym_result = LatticeSearch(latt);
|
|
outcome.symmetry = LatticeMessage{
|
|
.centering = sym_result.centering,
|
|
.niggli_class = sym_result.niggli_class,
|
|
.crystal_system = sym_result.system
|
|
};
|
|
outcome.lattice_candidate = sym_result.conventional;
|
|
}
|
|
|
|
// Multi-lattice search for stills: store rotations that map the reference
|
|
// lattice to each accepted extra lattice. Candidates are materialized later
|
|
// in RefineGeometryIfNeeded so they're rooted in the refined (and, for
|
|
// monoclinic, reordered) main lattice.
|
|
if (outcome.lattice_candidate && indexer_result.lattice.size() > 1) {
|
|
auto ml_latt = MultiLatticeSearch(indexer_result.lattice);
|
|
for (auto &ml : ml_latt) {
|
|
if (outcome.extra_lattice_rotations.size() >= experiment.GetIndexingSettings().GetMaxExtraLattices())
|
|
break;
|
|
outcome.extra_lattice_rotations.push_back(ml.rotation_vector);
|
|
RotMatrix rot(ml.rotation_vector.Length(), ml.rotation_vector.Normalize());
|
|
outcome.extra_lattice_candidates.push_back(outcome.lattice_candidate->Multiply(rot));
|
|
}
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
auto indexing_end_time = std::chrono::steady_clock::now();
|
|
msg.indexing_time_s = std::chrono::duration<float>(indexing_end_time - indexing_start_time).count();
|
|
|
|
return outcome;
|
|
}
|
|
|
|
namespace {
|
|
// Count spots whose fractional Miller index falls within the indexing tolerance of an integer for a
|
|
// given lattice + geometry - the "how well does this model explain the spots" score used by -r multi.
|
|
int CountIndexedSpots(const DiffractionGeometry &geom, const CrystalLattice &latt,
|
|
const std::vector<SpotToSave> &spots, float tol_sq) {
|
|
const Coord a = latt.Vec0(), b = latt.Vec1(), c = latt.Vec2();
|
|
int n = 0;
|
|
for (const auto &s : spots) {
|
|
const Coord recip = s.ReciprocalCoord(geom);
|
|
const float hf = recip * a, kf = recip * b, lf = recip * c;
|
|
// std::rint rather than std::round - only the squared residual is used, so the tie rule cannot
|
|
// change the count (see DetermineLatticeAndSymmetry).
|
|
const float dh = hf - std::rint(hf), dk = kf - std::rint(kf), dl = lf - std::rint(lf);
|
|
if (dh * dh + dk * dk + dl * dl < tol_sq) ++n;
|
|
}
|
|
return n;
|
|
}
|
|
} // namespace
|
|
|
|
void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::IndexingOutcome &outcome) {
|
|
if (!outcome.lattice_candidate)
|
|
return;
|
|
|
|
auto start_time = std::chrono::steady_clock::now();
|
|
|
|
XtalOptimizerData data{
|
|
.geom = outcome.experiment.GetDiffractionGeometry(),
|
|
.latt = *outcome.lattice_candidate,
|
|
.crystal_system = outcome.symmetry.crystal_system,
|
|
.min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(),
|
|
// Match the [30,150] deg bound the indexers already use (FFBIDXIndexer, FFT settings):
|
|
// the struct default [60,120] clamps a monoclinic beta outside that window (e.g. a C2
|
|
// beta near 132 deg) to the boundary, corrupting the per-frame cell refinement.
|
|
.min_angle_deg = 30.0f,
|
|
.max_angle_deg = 150.0f,
|
|
.refine_beam_center = true,
|
|
.refine_distance_mm = false,
|
|
.refine_detector_angles = false,
|
|
.refine_unit_cell = !experiment.IsRotationIndexing(),
|
|
// The whole spot list is passed below, not the indexed subset, so on weak images most of what
|
|
// enters the fit at the loose first tolerance is arbitrarily indexed noise. Weight every spot by
|
|
// how strong it is for its resolution so those contribute without dragging the orientation.
|
|
.weight_spots_by_confidence = true,
|
|
// Online: 40 ms is the real budget per image, so the wall clock is the right bound even though
|
|
// it makes the answer depend on machine load. Offline: bound the same refinement by iterations
|
|
// instead, so reprocessing the same file twice gives the same lattice.
|
|
.max_time = 0.04,
|
|
.max_iterations = real_time ? 0 : OFFLINE_REFINE_ITERATIONS
|
|
};
|
|
|
|
|
|
if (outcome.symmetry.crystal_system == gemmi::CrystalSystem::Trigonal)
|
|
data.crystal_system = gemmi::CrystalSystem::Hexagonal;
|
|
|
|
switch (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm()) {
|
|
case GeomRefinementAlgorithmEnum::None:
|
|
break;
|
|
case GeomRefinementAlgorithmEnum::OrientationOnly:
|
|
XtalOptimizerRotationOnly(data, msg.spots, 0.2);
|
|
XtalOptimizerRotationOnly(data, msg.spots, 0.1);
|
|
XtalOptimizerRotationOnly(data, msg.spots, 0.05);
|
|
break;
|
|
case GeomRefinementAlgorithmEnum::BeamCenter:
|
|
if (XtalOptimizer(data, {msg.spots})) {
|
|
outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl())
|
|
.BeamY_pxl(data.geom.GetBeamY_pxl());
|
|
outcome.beam_center_updated = true;
|
|
}
|
|
break;
|
|
case GeomRefinementAlgorithmEnum::Flex: {
|
|
// Try all three refinements per image and keep whichever indexes the most spots. Beam+cell
|
|
// refinement helps some stills but diverges on sparse spot lists (few spots, long axes),
|
|
// where it pushes a good lattice out of tolerance; scoring by indexed-spot count lets each
|
|
// image fall back to orientation-only or no refinement when refinement would hurt. Ties
|
|
// prefer less refinement (strict >, not >=) to avoid overfitting.
|
|
const float tol = experiment.GetIndexingSettings().GetTolerance();
|
|
const float tol_sq = tol * tol;
|
|
|
|
XtalOptimizerData d_none = data;
|
|
XtalOptimizerData d_orient = data;
|
|
XtalOptimizerRotationOnly(d_orient, msg.spots, 0.2);
|
|
XtalOptimizerRotationOnly(d_orient, msg.spots, 0.1);
|
|
XtalOptimizerRotationOnly(d_orient, msg.spots, 0.05);
|
|
XtalOptimizerData d_beam = data;
|
|
const bool beam_ok = XtalOptimizer(d_beam, {msg.spots});
|
|
|
|
const int s_none = CountIndexedSpots(d_none.geom, d_none.latt, msg.spots, tol_sq);
|
|
const int s_orient = CountIndexedSpots(d_orient.geom, d_orient.latt, msg.spots, tol_sq);
|
|
const int s_beam = beam_ok ? CountIndexedSpots(d_beam.geom, d_beam.latt, msg.spots, tol_sq) : -1;
|
|
|
|
if (s_beam > s_none && s_beam > s_orient) {
|
|
data = d_beam;
|
|
outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl())
|
|
.BeamY_pxl(data.geom.GetBeamY_pxl());
|
|
outcome.beam_center_updated = true;
|
|
} else if (s_orient > s_none) {
|
|
data = d_orient;
|
|
} else {
|
|
data = d_none;
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
|
|
outcome.lattice_candidate = data.latt;
|
|
|
|
if (outcome.symmetry.crystal_system == gemmi::CrystalSystem::Monoclinic)
|
|
outcome.lattice_candidate->ReorderMonoclinic();
|
|
|
|
// Rebuild extra-lattice candidates from the refined (and possibly reordered) main
|
|
// lattice so they share its cell and obtuse-beta convention.
|
|
if (!outcome.extra_lattice_rotations.empty()) {
|
|
outcome.extra_lattice_candidates.clear();
|
|
outcome.extra_lattice_candidates.reserve(outcome.extra_lattice_rotations.size());
|
|
for (const auto &rv : outcome.extra_lattice_rotations) {
|
|
RotMatrix rot(rv.Length(), rv.Normalize());
|
|
outcome.extra_lattice_candidates.push_back(outcome.lattice_candidate->Multiply(rot));
|
|
}
|
|
}
|
|
|
|
// Quick orientation-only refinement of extra lattices (stills path).
|
|
// Cell, beam center, detector geometry are taken from the first lattice.
|
|
if (!experiment.IsRotationIndexing() && !outcome.extra_lattice_candidates.empty()) {
|
|
for (auto &el : outcome.extra_lattice_candidates) {
|
|
XtalOptimizerData data_extra{
|
|
.geom = data.geom,
|
|
.latt = el,
|
|
.crystal_system = data.crystal_system,
|
|
.min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(),
|
|
.refine_beam_center = false,
|
|
.refine_distance_mm = false,
|
|
.refine_detector_angles = false,
|
|
.refine_unit_cell = false,
|
|
.refine_rotation_axis = false,
|
|
.index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(),
|
|
.max_time = 0.02,
|
|
.max_iterations = real_time ? 0 : OFFLINE_REFINE_ITERATIONS / 2
|
|
};
|
|
XtalOptimizerRotationOnly(data_extra, msg.spots, 0.1);
|
|
el = data_extra.latt;
|
|
}
|
|
}
|
|
|
|
if (outcome.beam_center_updated) {
|
|
msg.beam_corr_x = data.beam_corr_x;
|
|
msg.beam_corr_y = data.beam_corr_y;
|
|
}
|
|
|
|
auto end_time = std::chrono::steady_clock::now();
|
|
msg.refinement_time_s = std::chrono::duration_cast<std::chrono::duration<double>>(end_time - start_time).count();
|
|
}
|
|
|
|
void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
|
|
const SpotFindingSettings &spot_finding_settings,
|
|
BraggPrediction &prediction,
|
|
const BraggIntegrateFn &integrate,
|
|
const IndexAndRefine::IndexingOutcome &outcome) {
|
|
if (!outcome.lattice_candidate)
|
|
return;
|
|
|
|
|
|
CrystalLattice latt = outcome.lattice_candidate.value();
|
|
|
|
// Prediction uses each frame's OWN mosaicity/profile_radius (image-local). We deliberately do NOT
|
|
// smooth them here with a running moving average: it averaged the last N *processed* frames, whose
|
|
// order under the parallel per-image loop is thread-arrival order, making the predicted rocking
|
|
// width - and hence which reflections are integrated - non-deterministic run-to-run. Prediction only
|
|
// decides membership (a reflection on the cutoff contributes ~nothing), so the per-frame value is
|
|
// fine here. The mosaicity smoothing that actually matters - keeping the partialities of one rocking
|
|
// event consistent so they tile the curve and sum toward 1 - is done deterministically in frame
|
|
// order before the 3D combine (RotationScaleMerge), where partiality is recomputed from it.
|
|
|
|
float ewald_dist_cutoff = 0.001f;
|
|
if (msg.profile_radius)
|
|
ewald_dist_cutoff = msg.profile_radius.value() * 2.0f;
|
|
|
|
if (experiment.GetBraggIntegrationSettings().GetFixedProfileRadius_recipA())
|
|
ewald_dist_cutoff = experiment.GetBraggIntegrationSettings().GetFixedProfileRadius_recipA().value() * 3.0f;
|
|
|
|
float wedge_deg = 0.0f;
|
|
// The rocking width, when this frame has one. A frame too sparse to fit its own (and, on rotation,
|
|
// one integrated purely from the sweep's lattice) predicts with the default below, but must NOT
|
|
// report that default onward: RotationScaleMerge averages the reported values in frame order to
|
|
// recompute every partiality, and a placeholder entered there is read as a measurement.
|
|
std::optional<float> mos_measured;
|
|
|
|
if (experiment.GetGoniometer().has_value()) {
|
|
// Full oscillation wedge of one frame; BraggPredictionRot halves it to the +/- half-wedge of the
|
|
// partiality erf pair (Kabsch). Passing the full increment gives a half-wedge of increment/2 -
|
|
// matching ScaleOnTheFly's RotationPartiality, so the predicted partiality is used directly there.
|
|
wedge_deg = experiment.GetGoniometer()->GetWedge_deg();
|
|
|
|
if (msg.mosaicity_deg) {
|
|
mos_measured = msg.mosaicity_deg.value();
|
|
mosaicity[msg.number] = *mos_measured;
|
|
}
|
|
// Second pass of the rotation two-pass: widen the prediction to the frame-order-smoothed mosaicity
|
|
// that RotationScaleMerge fitted in the first pass. Take the MAX with this frame's own estimate so
|
|
// prediction is never NARROWER than the first pass - a too-narrow smoothed value would otherwise drop
|
|
// reflections and collapse the multiplicity. (A wider value only helps prediction cover the spot.)
|
|
if (msg.number >= 0 && msg.number < static_cast<int64_t>(prediction_mosaicity_override_.size())
|
|
&& std::isfinite(prediction_mosaicity_override_[msg.number])
|
|
&& prediction_mosaicity_override_[msg.number] > 0.0f) {
|
|
mos_measured = std::max(mos_measured.value_or(0.0f), prediction_mosaicity_override_[msg.number]);
|
|
mosaicity[msg.number] = *mos_measured;
|
|
}
|
|
}
|
|
// The width the prediction window opens to. It is the measured sigma_M unless a fixed prediction
|
|
// mosaicity was asked for: the reported (and therefore integrated-against) sigma_M below is left
|
|
// alone either way, so the two uses of the number can be separated.
|
|
float mos_deg = mos_measured.value_or(0.1f);
|
|
if (const auto forced_pred = experiment.GetBraggIntegrationSettings().GetForcedPredictionMosaicity_deg())
|
|
mos_deg = *forced_pred;
|
|
|
|
IntegrationOutcome i_outcome{
|
|
.geom = outcome.experiment.GetDiffractionGeometry(),
|
|
.latt = latt,
|
|
.mosaicity_deg = mos_measured,
|
|
.image_scale_cc = msg.image_scale_cc,
|
|
};
|
|
|
|
BraggPredictionSettings settings_prediction{
|
|
.high_res_A = IntegrationDMin_A(experiment),
|
|
.ewald_dist_cutoff = ewald_dist_cutoff,
|
|
// Centering is a hypothesis to confirm, not assume: with no user-fixed space group, predict
|
|
// in P so the centering-absent reflections are integrated and the space-group search can
|
|
// confirm or disprove centering (and catch a missed superstructure). A user-fixed space
|
|
// group is trusted, so reject its absences here.
|
|
.centering = experiment.GetGemmiSpaceGroup().has_value() ? outcome.symmetry.centering : 'P',
|
|
.wedge_deg = std::fabs(wedge_deg),
|
|
.mosaicity_deg = std::fabs(mos_deg),
|
|
// FWHM -> sigma; 0 when monochromatic, leaving the prediction unchanged.
|
|
.bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f,
|
|
};
|
|
ApplyPredictionRange(settings_prediction, experiment, latt);
|
|
|
|
// Online is bounded by what the image-buffer slot can carry, offline by the prediction limit; both
|
|
// come from BraggPrediction so the cap, the prediction and the transport headroom cannot drift.
|
|
// The predictor applies it, so what is over the cap is never integrated.
|
|
prediction.output_limit = real_time ? BraggPrediction::kOnlineMaxReflections
|
|
: BraggPrediction::kPredictionOutput;
|
|
|
|
// Predict, then integrate with the selected integrator (box-sum or profile-fit).
|
|
auto pred_start_time = std::chrono::steady_clock::now();
|
|
auto nrefl = prediction.Calc(outcome.experiment, latt, settings_prediction);
|
|
auto pred_end_time = std::chrono::steady_clock::now();
|
|
msg.bragg_prediction_time_s = std::chrono::duration<float>(pred_end_time - pred_start_time).count();
|
|
|
|
// The engine picks box-sum vs profile-fit internally from the experiment's IntegratorMode; the
|
|
// caller's callback binds it to the right image (GPU-resident buffer, host buffer, or the assembled
|
|
// FPGA image read straight on the CPU).
|
|
auto integration_start_time = std::chrono::steady_clock::now();
|
|
i_outcome.reflections = integrate(prediction.GetReflections(), nrefl, msg.number);
|
|
msg.integrated_reflections = i_outcome.reflections.size();
|
|
auto integration_end_time = std::chrono::steady_clock::now();
|
|
msg.integration_time_s = std::chrono::duration<float>(integration_end_time - integration_start_time).count();
|
|
|
|
CalcISigma(msg, i_outcome.reflections);
|
|
CalcWilsonBFactor(msg, i_outcome.reflections);
|
|
|
|
ScaleImage(msg, i_outcome);
|
|
|
|
// Copy reflections to outgoing message
|
|
msg.reflections = i_outcome.reflections;
|
|
|
|
// Persist the per-image result for the whole-run scaling/merge pass, unless the caller opted out
|
|
// (viewer interactive use only needs the current image, returned above via msg).
|
|
if (retain_outcomes_) {
|
|
const std::unique_lock ul(reflections_mutex);
|
|
integration_outcome[msg.number] = std::move(i_outcome);
|
|
}
|
|
}
|
|
|
|
std::optional<IndexAndRefine::IndexingOutcome>
|
|
IndexAndRefine::DetermineRefineAnalyze(DataMessage &msg, const SpotFindingSettings &spot_finding_settings) {
|
|
if (!indexer_ || !spot_finding_settings.indexing)
|
|
return std::nullopt;
|
|
|
|
IndexingOutcome outcome(experiment);
|
|
|
|
if (rotation_indexer)
|
|
outcome = DetermineLatticeAndSymmetryRotation(msg);
|
|
else
|
|
outcome = DetermineLatticeAndSymmetry(msg);
|
|
|
|
if (!outcome.lattice_candidate)
|
|
return std::nullopt;
|
|
|
|
if (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm() != GeomRefinementAlgorithmEnum::None)
|
|
RefineGeometryIfNeeded(msg, outcome);
|
|
|
|
if (!outcome.lattice_candidate.has_value())
|
|
return std::nullopt;
|
|
|
|
// AnalyzeIndexing answers "is this frame worth integrating"; msg.indexing_result carries the
|
|
// stricter "does this frame index on its own", which on rotation is not the same question.
|
|
if (!AnalyzeIndexing(msg, outcome.experiment, *outcome.lattice_candidate, outcome.extra_lattice_candidates))
|
|
return std::nullopt;
|
|
|
|
{
|
|
std::unique_lock ul(reflections_mutex);
|
|
unit_cells[msg.number] = outcome.lattice_candidate->GetUnitCell();
|
|
}
|
|
msg.lattice_type = outcome.symmetry;
|
|
return outcome;
|
|
}
|
|
|
|
void IndexAndRefine::ProcessImage(DataMessage &msg,
|
|
const SpotFindingSettings &spot_finding_settings,
|
|
BraggPrediction &prediction,
|
|
const BraggIntegrateFn &integrate) {
|
|
auto outcome = DetermineRefineAnalyze(msg, spot_finding_settings);
|
|
if (outcome && spot_finding_settings.quick_integration)
|
|
QuickPredictAndIntegrate(msg, spot_finding_settings, prediction, integrate, *outcome);
|
|
}
|
|
|
|
bool IndexAndRefine::IndexFrameOnly(DataMessage &msg, const SpotFindingSettings &spot_finding_settings) {
|
|
// The per-frame verdict, not "was there anything to integrate": this is what the rotation first pass
|
|
// scores candidate lattices on, and on rotation a frame can be integrable without indexing on its
|
|
// own (AnalyzeIndexing), which would score a sparse frame for every candidate alike.
|
|
return DetermineRefineAnalyze(msg, spot_finding_settings).has_value()
|
|
&& msg.indexing_result.value_or(false);
|
|
}
|
|
|
|
std::optional<RotationIndexerResult> IndexAndRefine::FinalizeRotationIndexing() {
|
|
if (rotation_indexer) {
|
|
if (const auto latt = rotation_indexer->GetLattice())
|
|
return latt;
|
|
|
|
rotation_indexer->RunIndexing();
|
|
return rotation_indexer->GetLattice();
|
|
}
|
|
return {};
|
|
}
|
|
|
|
IndexAndRefine &IndexAndRefine::ReferenceIntensities(std::vector<MergedReflection> &reference) {
|
|
// An external reference is trusted to be in the correct hand, so use it to break the merohedral
|
|
// indexing ambiguity per image (serial stills index each crystal independently).
|
|
reindex_resolver = std::make_unique<ReindexAmbiguityResolver>(experiment, reference);
|
|
return *this;
|
|
}
|
|
|
|
void IndexAndRefine::ScaleImage(DataMessage &msg, IntegrationOutcome& outcome) {
|
|
if (!reindex_resolver)
|
|
return;
|
|
|
|
// The external reference fixes the cell/space group, breaks the indexing ambiguity and reports CCref,
|
|
// but is NEVER a scale anchor: scaling an image against a foreign dataset injects cross-dataset
|
|
// systematics and is a worse reference than the data's own merge, so scaling self-references at the
|
|
// post-measurement merge for both workflows. Rotation resolves the ambiguity globally and self-scales
|
|
// in RotationScaleMerge (ChooseReindex / ReferenceIntensityCC), so there is nothing to do per image.
|
|
// Stills resolve the merohedral ambiguity per image here (each crystal indexes in a random hand; pick
|
|
// the hand best-correlated with the reference, once and for good).
|
|
if (experiment.IsRotationIndexing())
|
|
return;
|
|
|
|
auto scaling_start_time = std::chrono::steady_clock::now();
|
|
reindex_resolver->Resolve(outcome.reflections);
|
|
auto scaling_end_time = std::chrono::steady_clock::now();
|
|
msg.image_scale_time_s = std::chrono::duration<float>(scaling_end_time - scaling_start_time).count();
|
|
}
|
|
|
|
ScalingResult IndexAndRefine::ScaleAllImages(const std::vector<MergedReflection> &reference, size_t nthreads) {
|
|
ScaleOnTheFly scaling(experiment, reference);
|
|
scaling.Scale(integration_outcome, nthreads);
|
|
scale_cc.resize(integration_outcome.size());
|
|
|
|
for (int i = 0; i < integration_outcome.size(); i++)
|
|
scale_cc.at(i) = integration_outcome[i].image_scale_cc.value_or(NAN);
|
|
|
|
return ScalingResult(integration_outcome);
|
|
}
|
|
|
|
const std::vector<float> &IndexAndRefine::GetImageCC() const {
|
|
return scale_cc;
|
|
}
|
|
|
|
const std::vector<std::optional<UnitCell> > & IndexAndRefine::GetUnitCells() const {
|
|
return unit_cells;
|
|
}
|
|
|
|
std::optional<UnitCell> IndexAndRefine::GetConsensusUnitCell() const {
|
|
const auto dist_tolerance = experiment.GetIndexingSettings().GetUnitCellDistTolerance();
|
|
const auto angle_tolerance = experiment.GetIndexingSettings().GetUnitCellAngleTolerance_deg();
|
|
|
|
if (rotation_indexer) {
|
|
auto result = rotation_indexer->GetLattice();
|
|
if (!result)
|
|
return {};
|
|
return result->lattice.GetUnitCell();
|
|
}
|
|
|
|
std::vector<UnitCell> cells;
|
|
{
|
|
std::unique_lock ul(reflections_mutex);
|
|
cells.reserve(unit_cells.size());
|
|
for (const auto &cell: unit_cells) {
|
|
if (cell && cell->is_finite())
|
|
cells.emplace_back(*cell);
|
|
}
|
|
}
|
|
|
|
if (cells.empty())
|
|
return {};
|
|
|
|
if (experiment.GetUnitCell()) {
|
|
std::vector<UnitCell> accepted;
|
|
accepted.reserve(cells.size());
|
|
|
|
for (const auto &cell: cells) {
|
|
if (cell.is_close(*experiment.GetUnitCell(), dist_tolerance, angle_tolerance))
|
|
accepted.emplace_back(cell);
|
|
}
|
|
|
|
return MeanUnitCell(accepted);
|
|
}
|
|
|
|
size_t best_count = 0;
|
|
UnitCell best_reference{};
|
|
|
|
for (const auto &ref: cells) {
|
|
size_t count = 0;
|
|
for (const auto &cell: cells) {
|
|
if (cell.is_close(ref, dist_tolerance, angle_tolerance))
|
|
++count;
|
|
}
|
|
|
|
if (count > best_count) {
|
|
best_count = count;
|
|
best_reference = ref;
|
|
}
|
|
}
|
|
|
|
if (best_count == 0)
|
|
return {};
|
|
|
|
std::vector<UnitCell> accepted;
|
|
accepted.reserve(best_count);
|
|
|
|
for (const auto &cell: cells) {
|
|
if (cell.is_close(best_reference, dist_tolerance, angle_tolerance))
|
|
accepted.emplace_back(cell);
|
|
}
|
|
|
|
return MeanUnitCell(accepted);
|
|
}
|
|
|
|
std::vector<IntegrationOutcome> &IndexAndRefine::GetIntegrationOutcome() {
|
|
return integration_outcome;
|
|
}
|
|
|
|
const std::vector<IntegrationOutcome> &IndexAndRefine::GetIntegrationOutcome() const {
|
|
return integration_outcome;
|
|
}
|
|
|
|
void IndexAndRefine::ForceRotationIndexerLattice(const CrystalLattice &lattice) {
|
|
if (rotation_indexer)
|
|
rotation_indexer->ForceLattice(lattice);
|
|
}
|
|
|
|
void IndexAndRefine::ForceRotationIndexerResult(const RotationIndexerResult &result) {
|
|
if (rotation_indexer)
|
|
rotation_indexer->ForceResult(result);
|
|
}
|