leonarski_fandClaude Opus 5.5
2b476336ae
Every lattice decision counts spots or frames, and the sub-lattice's strong reflections win every count, so a crystal whose cell is doubled by a weak superstructure class (9min, 6z9g) is adopted at the half cell. This asks the question in intensities instead. On 60 frames spread over the sweep, after each frame's own integration, the 2a x 2b x 2c supercell of its primitive lattice is predicted to 3 A with the frame's own refined orientation and geometry and integrated on the same engine; the reflections are summed per parity class in two shells (20-5, 5-3 A), with a fit of intensity against partiality (I = a + b p) that separates what rocks like a Bragg reflection from what sits at the node whatever the rocking. Only the tested frames are predicted and nothing is retained, so memory is bounded (the previous prototype predicted the whole run through the merge and ran out of GPU memory). The probe's integrations are kept out of the engine's own counts, which the two-pass stencil guard reads. Results are bit-identical with and without it. REPORT ONLY - it decides nothing, because on the battery it does not yet separate a weak real class from what sits at the half-integer nodes of crystals whose cell is right. Real classes: 6z9g class 101 at 24 % of the lattice's intensity (29 % rocking), 9min 100 at 19 % (4 % rocking - its real class does not rock like the lattice either). On correct cells the largest classes reach 9-12 % raw (7n2s, 9i0a, 7os3) and 3-4 % rocking (7dkp, 7os3), and 7mzt reads 40 % / 19 % on a class the deposition does not have. The log line is the population a decision has to be calibrated on. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
897 lines
45 KiB
C++
897 lines
45 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),
|
|
indexer_(indexer),
|
|
rotation_indexer_counter(x) {
|
|
if (indexer && x.IsRotationIndexing())
|
|
rotation_indexer = std::make_unique<RotationIndexer>(x, *indexer, real_time);
|
|
// 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;
|
|
std::optional<float> spindle_blind_fraction;
|
|
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();
|
|
// Read per call, for the same reason: a beam-centre hypothesis moves the geometry after this
|
|
// object is built, and a spot mapped into reciprocal space through the centre the run started at
|
|
// is not the spot the trial centre sees. Held at construction, this scored every hypothesis of the
|
|
// beam-centre ladder against the file's geometry - which is why a centre scored 42/60 inside the
|
|
// ladder and 5/60 when the same value was given on the command line.
|
|
const DiffractionGeometry geom = experiment.GetDiffractionGeometry();
|
|
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;
|
|
// Kept from whichever run could measure it - the largest seed that answered - rather than from
|
|
// the run whose lattice won: the severity describes the frame's lattice rows, not the cell that
|
|
// closed on them, so a frame that indexed nothing still has one. When the escalation never
|
|
// leaves a seed below the score's spot floor it stays absent here, and the severity-only pass
|
|
// after the loop supplies it instead.
|
|
if (res.spindle_blind_fraction)
|
|
spindle_blind_fraction = res.spindle_blind_fraction;
|
|
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;
|
|
|
|
// The severity's spot floor (SPINDLE_MIN_SPOTS) was calibrated on a frame's whole spot list,
|
|
// but the escalation stops at the leanest seed that indexes - 30 spots on precisely the clean
|
|
// frames a grid scan produces - so riding on the indexing seed left the score absent exactly
|
|
// where automation most needs it, and absence maps to "engage" (see SpindleBlindFraction.h).
|
|
// Decouple the two: when no escalation pass could answer, spend one severity-only row pass over
|
|
// the full spot list (no reduction, no refinement); where that path does not exist - an indexer
|
|
// without a row search - or refuses, read the severity off the rows of the winning lattice,
|
|
// which any indexer produces.
|
|
if (!spindle_blind_fraction && experiment.GetGoniometer().has_value()) {
|
|
std::vector<Coord> recip;
|
|
recip.reserve(std::min<size_t>(msg.spots.size(), FFT_MAX_SPOTS));
|
|
for (const auto &i : msg.spots) {
|
|
if (index_ice_rings || !i.ice_ring) {
|
|
recip.push_back(i.ReciprocalCoord(geom));
|
|
if (recip.size() >= FFT_MAX_SPOTS)
|
|
break;
|
|
}
|
|
}
|
|
if (recip.size() >= SPINDLE_MIN_SPOTS) {
|
|
const auto algorithm = experiment.GetIndexingAlgorithm();
|
|
if (algorithm == IndexingAlgorithmEnum::FFT || algorithm == IndexingAlgorithmEnum::FFTW) {
|
|
const auto res = indexer_->Run(experiment, recip, /*severity_only=*/true);
|
|
spindle_blind_fraction = res.spindle_blind_fraction;
|
|
}
|
|
if (!spindle_blind_fraction && !indexer_result.lattice.empty()) {
|
|
const float theta_max_deg = SpindleThetaMax_deg(experiment.GetWavelength_A(),
|
|
experiment.GetDetectorMaxResolution_A());
|
|
if (const auto severity = SpindleBlindFractionFromLattice(
|
|
indexer_result.lattice.front(), experiment.GetGoniometer()->GetAxis(),
|
|
theta_max_deg))
|
|
spindle_blind_fraction = severity->score;
|
|
}
|
|
}
|
|
}
|
|
|
|
msg.spindle_blind_fraction = spindle_blind_fraction;
|
|
|
|
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_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
|
|
};
|
|
|
|
// This refinement holds one frame, so it deliberately does not back-rotate - the frame's angle is
|
|
// a gauge the orientation block absorbs, and no `axis` is passed above. The spots on that frame
|
|
// still diffracted at different angles across the exposure, though, and that spread does not go
|
|
// away with the gauge: hand the rocking geometry over so the acceptance gate can profile it out.
|
|
if (const auto &gon = experiment.GetGoniometer(); gon.has_value() && gon->IsScanning()) {
|
|
data.rocking_spindle = gon->GetAxis();
|
|
data.rocking_wedge_deg = gon->GetWedge_deg();
|
|
}
|
|
|
|
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_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,
|
|
const BraggIntegrateFn &probe_integrate) {
|
|
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;
|
|
// The width the PREDICTION window opens to, which is a different question from the width the
|
|
// partiality divides by. Prediction only decides membership, and a reflection just inside a
|
|
// generous window arrives with a partiality near zero and is weighted accordingly, so erring wide
|
|
// there is close to free (measured: a 28% wider window moves the merged statistics by under half a
|
|
// percent). The divisor is not free - it multiplies every partial - so it must stay on the width
|
|
// this frame actually measured.
|
|
std::optional<float> mos_predict;
|
|
|
|
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_predict = std::max(mos_measured.value_or(0.0f), prediction_mosaicity_override_[msg.number]);
|
|
}
|
|
float mos_deg = mos_predict.value_or(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 - unless this pass measures geometry and
|
|
// discards its intensities (predict_all_centring_nodes_), where the absences only halve the
|
|
// events the geometry is fitted from.
|
|
.centering = experiment.GetGemmiSpaceGroup().has_value() && !predict_all_centring_nodes_
|
|
? 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);
|
|
prediction_centring_ = settings_prediction.centering;
|
|
|
|
// 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();
|
|
const std::vector<Reflection> integrated = integrate(prediction.GetReflections(), nrefl, msg.number);
|
|
msg.integrated_reflections = integrated.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();
|
|
|
|
CalcISigmaAndWilsonBFactor(msg, integrated);
|
|
|
|
// A retained outcome is kept in the pass's arena (see ReflectionArena), not beside the worker's
|
|
// other allocations.
|
|
i_outcome.reflections = ReflectionVector(integrated.begin(), integrated.end(),
|
|
ArenaAllocator<Reflection>(retain_outcomes_ ? &reflection_arena
|
|
: nullptr));
|
|
|
|
ScaleImage(msg, i_outcome);
|
|
|
|
// Copy reflections to outgoing message
|
|
if (keep_reflections_in_message_)
|
|
msg.reflections.assign(i_outcome.reflections.begin(), i_outcome.reflections.end());
|
|
|
|
// 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);
|
|
}
|
|
|
|
if (probe_integrate && std::binary_search(supercell_probe_frames_.begin(), supercell_probe_frames_.end(),
|
|
msg.number))
|
|
ProbeSupercellFrame(msg, prediction, probe_integrate, outcome, settings_prediction);
|
|
}
|
|
|
|
SupercellProbeFit FitSupercellProbe(const SupercellProbeClass &c) {
|
|
SupercellProbeFit f;
|
|
if (c.n < 3)
|
|
return f;
|
|
const double n = static_cast<double>(c.n);
|
|
const double mp = c.sum_p / n, mi = c.sum_i / n;
|
|
const double spp = c.sum_pp - n * mp * mp, spi = c.sum_pi - n * mp * mi, sii = c.sum_ii - n * mi * mi;
|
|
f.mean_i = mi;
|
|
f.mean_p = mp;
|
|
f.mean_i_over_sigma = c.sum_i_over_sigma / n;
|
|
if (spp > 0.0) {
|
|
f.b = spi / spp;
|
|
f.b_se = std::sqrt(std::max(0.0, sii - f.b * spi) / (n - 2.0) / spp);
|
|
}
|
|
f.a = mi - f.b * mp;
|
|
return f;
|
|
}
|
|
|
|
void IndexAndRefine::ProbeSupercell(std::vector<int64_t> image_numbers) {
|
|
std::sort(image_numbers.begin(), image_numbers.end());
|
|
supercell_probe_frames_ = std::move(image_numbers);
|
|
supercell_probe_ = {};
|
|
}
|
|
|
|
SupercellProbe IndexAndRefine::GetSupercellProbe() {
|
|
const std::unique_lock ul(supercell_probe_mutex_);
|
|
return supercell_probe_;
|
|
}
|
|
|
|
void IndexAndRefine::ProbeSupercellFrame(const DataMessage &msg, BraggPrediction &prediction,
|
|
const BraggIntegrateFn &integrate, const IndexingOutcome &outcome,
|
|
const BraggPredictionSettings &settings) {
|
|
// The band the classes are read over. A superstructure class is weakest at low resolution (the two
|
|
// halves of the cell differ in detail, not in position), so the band reaches to 3 A; 20 A keeps the
|
|
// strongest, most overlap-prone low-order reflections of the lattice out of the denominator. The
|
|
// prediction stops at 3 A as well, which is what keeps eight times the nodes of a frame affordable.
|
|
constexpr float BAND_LOW_RES_A = 20.0f;
|
|
constexpr float BAND_HIGH_RES_A = 3.0f;
|
|
|
|
const CrystalLattice prim = outcome.lattice_candidate->ToPrimitive(outcome.symmetry.centering);
|
|
const CrystalLattice super(prim.Vec0() * 2.0f, prim.Vec1() * 2.0f, prim.Vec2() * 2.0f);
|
|
BraggPredictionSettings s = settings;
|
|
s.centering = 'P';
|
|
s.high_res_A = std::max(settings.high_res_A, BAND_HIGH_RES_A);
|
|
ApplyPredictionRange(s, experiment, super);
|
|
const auto nrefl = prediction.Calc(outcome.experiment, super, s);
|
|
const std::vector<Reflection> integrated = integrate(prediction.GetReflections(), nrefl, msg.number);
|
|
|
|
SupercellProbe frame{};
|
|
for (const auto &r : integrated) {
|
|
if (!(r.d > BAND_HIGH_RES_A && r.d <= BAND_LOW_RES_A) || !(r.sigma > 0.0f) || !std::isfinite(r.I)
|
|
|| !std::isfinite(r.partiality))
|
|
continue;
|
|
auto &c = frame[4 * (r.h & 1) + 2 * (r.k & 1) + (r.l & 1)][r.d > 5.0f ? 0 : 1];
|
|
c.n++;
|
|
c.sum_i += r.I;
|
|
c.sum_i_over_sigma += r.I / r.sigma;
|
|
c.sum_p += r.partiality;
|
|
c.sum_pp += r.partiality * r.partiality;
|
|
c.sum_pi += r.partiality * r.I;
|
|
c.sum_ii += static_cast<double>(r.I) * r.I;
|
|
}
|
|
const std::unique_lock ul(supercell_probe_mutex_);
|
|
for (int i = 0; i < 8; i++)
|
|
for (int j = 0; j < 2; j++)
|
|
supercell_probe_[i][j] += frame[i][j];
|
|
}
|
|
|
|
std::optional<IndexAndRefine::IndexingOutcome>
|
|
IndexAndRefine::DetermineRefineAnalyze(DataMessage &msg, const SpotFindingSettings &spot_finding_settings,
|
|
int64_t *spots_on_lattice) {
|
|
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,
|
|
spots_on_lattice, integrate_every_frame_))
|
|
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,
|
|
const BraggIntegrateFn &probe_integrate) {
|
|
auto outcome = DetermineRefineAnalyze(msg, spot_finding_settings);
|
|
if (outcome && spot_finding_settings.quick_integration)
|
|
QuickPredictAndIntegrate(msg, spot_finding_settings, prediction, integrate, *outcome, probe_integrate);
|
|
}
|
|
|
|
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);
|
|
}
|
|
|
|
int64_t IndexAndRefine::CountSpotsOnLattice(DataMessage &msg, const SpotFindingSettings &spot_finding_settings) {
|
|
int64_t spots_on_lattice = 0;
|
|
DetermineRefineAnalyze(msg, spot_finding_settings, &spots_on_lattice);
|
|
return spots_on_lattice;
|
|
}
|
|
|
|
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 {};
|
|
}
|
|
|
|
std::optional<RotationIndexer::ConstrainedRefit>
|
|
IndexAndRefine::RefineRotationLatticeConstrained(const CrystalLattice &latt,
|
|
gemmi::CrystalSystem system) const {
|
|
if (!rotation_indexer)
|
|
return {};
|
|
return rotation_indexer->RefineConstrained(latt, system);
|
|
}
|
|
|
|
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);
|
|
}
|