leonarski_fandClaude Opus 5.5
42719adbf1
Known-cell (FFBIDX) indexing was given the conventional cell of a centred lattice. Its own candidate scoring and the indexed-spot count of Refine() then worked on the primitive lattice of that cell, which is 2x (A, B, C, I), 3x (R) or 4x (F) denser than the crystal's, so a wrong orientation that fits only those extra nodes indexed as many spots as the true one: - R in the hexagonal setting: the obverse orientation and the reverse one (60 deg about c) share every node of the primitive hexagonal lattice. About half of the frames came back reverse; only the l = 3n layers agreed between them, and the merge showed a false translational pseudo-symmetry. - C/I/F: false orientations that break the centring were accepted. The choice is made inside the indexer's top-candidate scoring, so a centring condition on the count in Refine() alone does not always settle it: on synthetic R3 stills the obverse orientation was sometimes not among the candidates at all. FFBIDX (GPU and CPU) is now given the Niggli-reduced primitive cell of the given lattice (Indexer::IndexingCell, with the centring of the given space group), and every lattice it finds is put back into the given conventional cell by the integer matrix between the two (Indexer::ToConventional), so everything downstream sees the conventional basis as before. Every node of the primitive cell is a possible reflection, so no centring count is needed in Refine(). The reduction matters: the plain primitive cell of an elongated C cell has an angle far outside the window Refine() accepts, and such a set then indexed nothing. The seed-escalation fraction in IndexAndRefine counts only spots the centring allows on the FFBIDX path. A primitive space group, a cell without a space group, and the de-novo (FFT) indexers are unchanged. All known-cell stills flows reach this through IndexAndRefine and the indexer thread pool: the broker's FPGA and SIMPLON receivers, live rugnux stills and offline rugnux. Rotation sweeps of R3/R32 crystals processed as stills with the known cell: reverse frames 47-50% before, 0% after; indexed spots obeying -h+k+l = 3n 0.67-0.68 before, >= 0.997 after; merged <I>(l=3n)/<I>(l!=3n) 1.2-11.9 before, 0.95-1.01 after. Still battery: primitive control unchanged; on an I222 set fewer frames index (67% -> 58%) because the extra ones were on a wrong lattice (indexed spots obeying the centring 0.54 -> 0.99), with R_split 0.088 -> 0.078, CC1/2 0.962 -> 0.981, R-free 0.162 -> 0.156; I23 indexes 3% more frames; C2221 unchanged. Tests: synthetic R3 obverse vs reverse and C-centred false orientation in Refine() on primitive candidates; FFBIDX CPU end-to-end on R, C (including an elongated one), I and F cells; IndexAndRefine on predicted R3 stills with FFBIDX. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UqXYdJKPUFLoeAchjzv3x4
78 lines
3.1 KiB
C++
78 lines
3.1 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "Indexer.h"
|
|
|
|
void Indexer::Setup(const DiffractionExperiment& experiment) {
|
|
geom = experiment.GetDiffractionGeometry();
|
|
indexing_tolerance = experiment.GetIndexingSettings().GetTolerance();
|
|
dist_tolerance_vs_reference = experiment.GetIndexingSettings().GetUnitCellDistTolerance();
|
|
step_down_to_reference = experiment.GetIndexingSettings().GetRotationIndexing();
|
|
search_short_rows = experiment.IsRotationIndexing();
|
|
viable_cell_min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots();
|
|
index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings();
|
|
wavelength_A = experiment.GetWavelength_A();
|
|
spindle_axis = {};
|
|
if (const auto goniometer = experiment.GetGoniometer())
|
|
spindle_axis = goniometer->GetAxis();
|
|
spindle_theta_max_deg = SpindleThetaMax_deg(wavelength_A, experiment.GetDetectorMaxResolution_A());
|
|
centering = experiment.GetCentering(); // before SetupUnitCell, which reads it
|
|
SetupUnitCell(experiment.GetUnitCell());
|
|
}
|
|
|
|
IndexerResult Indexer::Run(const std::vector<Coord> &coord, bool severity_only) {
|
|
IndexerResult ret;
|
|
auto start = std::chrono::steady_clock::now();
|
|
std::optional<SpindleSeverity> severity;
|
|
if (severity_only) {
|
|
severity = RunSeverityOnly(coord);
|
|
} else {
|
|
ret.lattice = RunInternal(coord, coord.size());
|
|
ret.executed = true;
|
|
severity = GetSpindleSeverity();
|
|
}
|
|
|
|
auto end = std::chrono::steady_clock::now();
|
|
std::chrono::duration<float> duration = end - start;
|
|
ret.indexing_time_s = duration.count();
|
|
if (severity) {
|
|
ret.spindle_blind_fraction = severity->score;
|
|
ret.spindle_row_length_A = severity->row_length_A;
|
|
ret.spindle_miss_angle_deg = severity->miss_angle_deg;
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
namespace {
|
|
Eigen::Matrix3f Rows(const CrystalLattice &l) {
|
|
Eigen::Matrix3f m;
|
|
m << l.Vec0().x, l.Vec0().y, l.Vec0().z, l.Vec1().x, l.Vec1().y, l.Vec1().z,
|
|
l.Vec2().x, l.Vec2().y, l.Vec2().z;
|
|
return m;
|
|
}
|
|
|
|
CrystalLattice FromRows(const Eigen::Matrix3f &m) {
|
|
return {Coord(m(0, 0), m(0, 1), m(0, 2)), Coord(m(1, 0), m(1, 1), m(1, 2)), Coord(m(2, 0), m(2, 1), m(2, 2))};
|
|
}
|
|
}
|
|
|
|
CrystalLattice Indexer::IndexingCell(const UnitCell &cell) {
|
|
const CrystalLattice conventional(cell);
|
|
to_conventional = Eigen::Matrix3f::Identity();
|
|
if (centering == 'P')
|
|
return conventional;
|
|
Eigen::Matrix3f reduced = Rows(conventional.ToPrimitive(centering).NiggliReduce());
|
|
if (reduced.determinant() < 0.0f)
|
|
reduced = -reduced; // right-handed, as Refine returns its cells
|
|
to_conventional = (Rows(conventional) * reduced.inverse()).array().round().matrix();
|
|
return FromRows(reduced);
|
|
}
|
|
|
|
std::vector<CrystalLattice> Indexer::ToConventional(const std::vector<CrystalLattice> &lattices) const {
|
|
std::vector<CrystalLattice> ret;
|
|
ret.reserve(lattices.size());
|
|
for (const auto &l: lattices)
|
|
ret.push_back(FromRows(to_conventional * Rows(l)));
|
|
return ret;
|
|
}
|