diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 9c910064d..b814fbcf0 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -9,6 +9,7 @@ #include "IndexAndRefine.h" #include "bragg_integration/CalcISigma.h" +#include "bragg_integration/SystematicAbsence.h" #include "geom_refinement/XtalOptimizer.h" #include "indexing/AnalyzeIndexing.h" #include "indexing/FFTIndexer.h" @@ -213,6 +214,11 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data // 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; + // A spot counts only where the lattice centring allows its index. FFBIDX returns the given + // conventional cell, where the centring holds; a de-novo indexer returns a reduced primitive cell, + // where there is none to apply. + const auto seed_algorithm = experiment.GetIndexingAlgorithm(); + const char seed_centering = (seed_algorithm == IndexingAlgorithmEnum::FFBIDX) ? experiment.GetCentering() : 'P'; constexpr float SEED_STOP_FRACTION = 0.9f; // seed explained this well -> stop escalating IndexerResult indexer_result; std::optional spindle_blind_fraction; @@ -265,7 +271,12 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data // 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; + if (dh * dh + dk * dk + dl * dl < idx_tol_sq + && (seed_centering == 'P' || !systematic_absence(static_cast(std::rint(hf)), + static_cast(std::rint(kf)), + static_cast(std::rint(lf)), + seed_centering))) + ++n; } const float frac = recip.empty() ? 0.0f : static_cast(n) / recip.size(); if (frac > best_frac) { best_frac = frac; indexer_result = std::move(res); } diff --git a/image_analysis/indexing/FFBIDXIndexer.cpp b/image_analysis/indexing/FFBIDXIndexer.cpp index fb62a7040..ad0122663 100644 --- a/image_analysis/indexing/FFBIDXIndexer.cpp +++ b/image_analysis/indexing/FFBIDXIndexer.cpp @@ -11,9 +11,10 @@ void FFBIDXIndexer::SetupUnitCell(const std::optional &cell) { if (!cell.has_value()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "FFBIDX requires unit cell"); - reference_unit_cell = cell; - - CrystalLattice l1(cell.value()); + // Given the primitive cell of a centred lattice (IndexingCell); the lattices found are put back + // into the conventional cell in RunInternal. + const CrystalLattice l1 = IndexingCell(cell.value()); + reference_unit_cell = l1.GetUnitCell(); Eigen::Matrix3f m; indexer.iCellM() << l1.Vec0().x, l1.Vec0().y, l1.Vec0().z, @@ -53,5 +54,6 @@ std::vector FFBIDXIndexer::RunInternal(const std::vector .indexing_tolerance = indexing_tolerance }; - return Refine(coord, nspots, indexer.oCellM(), indexer.oScoreV(), parameters); + auto lattices = Refine(coord, nspots, indexer.oCellM(), indexer.oScoreV(), parameters); + return ToConventional(lattices); } diff --git a/image_analysis/indexing/FFBIDXIndexerCPU.cpp b/image_analysis/indexing/FFBIDXIndexerCPU.cpp index 244c149a8..31bd83e3a 100644 --- a/image_analysis/indexing/FFBIDXIndexerCPU.cpp +++ b/image_analysis/indexing/FFBIDXIndexerCPU.cpp @@ -399,9 +399,10 @@ void FFBIDXIndexerCPU::SetupUnitCell(const std::optional &cell) { if (!cell.has_value()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "FFBIDX requires unit cell"); - reference_unit_cell = cell; - - CrystalLattice l(cell.value()); + // Given the primitive cell of a centred lattice (IndexingCell); the lattices found are put back + // into the conventional cell in RunInternal. + const CrystalLattice l = IndexingCell(cell.value()); + reference_unit_cell = l.GetUnitCell(); input_cell << 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; @@ -427,5 +428,6 @@ std::vector FFBIDXIndexerCPU::RunInternal(const std::vectorGetAxis(); spindle_theta_max_deg = SpindleThetaMax_deg(wavelength_A, experiment.GetDetectorMaxResolution_A()); + centering = experiment.GetCentering(); // before SetupUnitCell, which reads it SetupUnitCell(experiment.GetUnitCell()); } @@ -41,3 +42,36 @@ IndexerResult Indexer::Run(const std::vector &coord, bool severity_only) } 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 Indexer::ToConventional(const std::vector &lattices) const { + std::vector ret; + ret.reserve(lattices.size()); + for (const auto &l: lattices) + ret.push_back(FromRows(to_conventional * Rows(l))); + return ret; +} diff --git a/image_analysis/indexing/Indexer.h b/image_analysis/indexing/Indexer.h index 3a3d2de88..903a38737 100644 --- a/image_analysis/indexing/Indexer.h +++ b/image_analysis/indexing/Indexer.h @@ -5,6 +5,7 @@ #include #include +#include #include "../../common/CrystalLattice.h" @@ -48,6 +49,19 @@ protected: DiffractionGeometry geom; std::optional reference_unit_cell; + // Lattice centring of the given space group ('P' without one). + char centering = 'P'; + // Rows of the given conventional cell = to_conventional * rows of the cell from IndexingCell. + Eigen::Matrix3f to_conventional = Eigen::Matrix3f::Identity(); + + // The cell a known-cell indexer is given for a centred lattice: the Niggli-reduced primitive cell, so + // that every node of it is a reflection the crystal can have (the conventional cell of a centred + // lattice has 2x - A, B, C, I - 3x - R - or 4x - F - as many nodes, on which a wrong orientation fits + // as well as the true one), and so that its angles stay near 90 deg (the plain primitive cell of an + // elongated C cell has one far outside the window Refine accepts). The given cell itself for 'P'. + CrystalLattice IndexingCell(const UnitCell &cell); + // The lattices found on IndexingCell put back into the given conventional cell. + std::vector ToConventional(const std::vector &lattices) const; // The rotation axis a subsequent sweep would use. Present on a grid scan too, where the axis is // defined but stationary - that is exactly the case the spindle severity is for. std::optional spindle_axis; diff --git a/tests/FFBIDXIndexerCPUTest.cpp b/tests/FFBIDXIndexerCPUTest.cpp index 3c91c1f1e..0f0071716 100644 --- a/tests/FFBIDXIndexerCPUTest.cpp +++ b/tests/FFBIDXIndexerCPUTest.cpp @@ -14,6 +14,7 @@ #include "../common/CUDAWrapper.h" #include "../common/DiffractionExperiment.h" #include "../image_analysis/indexing/FFBIDXIndexerCPU.h" +#include "../image_analysis/bragg_integration/SystematicAbsence.h" #ifdef JFJOCH_USE_CUDA #include "../image_analysis/indexing/FFBIDXIndexer.h" @@ -24,7 +25,7 @@ namespace { // within 0.001 1/A of the Ewald sphere of a 1 A beam along z, and then random spots in the same // resolution range. The seed is fixed, and the engine scaled by hand, so the frame is the same // with every standard library. - std::vector SyntheticStill(const UnitCell &uc, uint32_t seed, size_t noise_spots) { + std::vector SyntheticStill(const UnitCell &uc, uint32_t seed, size_t noise_spots, char centering = 'P') { std::mt19937 rng(seed); auto draw = [&rng]() { return 2.0f * static_cast(rng()) / static_cast(std::mt19937::max()) - 1.0f; @@ -43,6 +44,8 @@ namespace { for (int h = -40; h <= 40; h++) for (int k = -40; k <= 40; k++) for (int m = -40; m <= 40; m++) { + if (centering != 'P' && systematic_absence(h, k, m, centering)) + continue; const Coord r = h * astar + k * bstar + m * cstar; if (r.Length() > 0.5f || r.Length() < 0.02f) continue; @@ -120,6 +123,50 @@ TEST_CASE("Portable_FFBIDXIndexerCPU_SyntheticStill", "[portable][Indexing]") { } } +TEST_CASE("Portable_FFBIDXIndexerCPU_CentredLattices", "[portable][Indexing]") { + // Centred lattices given as their conventional cells. For R3 in the hexagonal setting the obverse + // orientation and the reverse one (60 deg about c) share every node of the primitive hexagonal + // lattice; only the centring tells which of the two the spots are on, and the reverse one indexes + // only the l = 3n layers. The lattice must come back in the conventional cell, with its spots on + // nodes the centring allows. + struct Case { UnitCell uc; int space_group; char centering; }; + const std::vector cases = { + {{60, 60, 140, 90, 90, 120}, 146, 'R'}, + {{70, 40, 50, 90, 105, 90}, 5, 'C'}, + // a C cell with b >> a: its plain primitive cell has a 156 deg angle, the reduced one does not + {{35, 170, 120, 90, 90, 90}, 21, 'C'}, + {{60, 70, 80, 90, 90, 90}, 23, 'I'}, + {{100, 100, 100, 90, 90, 90}, 196, 'F'} + }; + for (const auto &c: cases) { + for (uint32_t seed = 0; seed < 10; seed++) { + const auto spots = SyntheticStill(c.uc, 4321 + seed, 0, c.centering); + INFO("centring " << c.centering << ", seed " << seed << ", " << spots.size() << " spots"); + REQUIRE(spots.size() > 30); + + DiffractionExperiment experiment; + experiment.SetUnitCell(c.uc).SpaceGroupNumber(c.space_group); + REQUIRE(experiment.GetCentering() == c.centering); + FFBIDXIndexerCPU indexer; + indexer.Setup(experiment); + const auto result = indexer.Run(spots); + REQUIRE(!result.lattice.empty()); + const CrystalLattice &l = result.lattice[0]; + CHECK(SameLengths(l.GetUnitCell(), c.uc, 0.01f)); + size_t n = 0; + for (const auto &s: spots) { + const float h = s * l.Vec0(), k = s * l.Vec1(), m = s * l.Vec2(); + const float dh = h - std::round(h), dk = k - std::round(k), dm = m - std::round(m); + if (dh * dh + dk * dk + dm * dm < 0.01f + && !systematic_absence(static_cast(std::round(h)), static_cast(std::round(k)), + static_cast(std::round(m)), c.centering)) + n++; + } + CHECK(static_cast(n) / static_cast(spots.size()) > 0.95f); + } + } +} + // Time per call of the CPU indexer, alone and with Refine(), on one thread and on 16 at once TEST_CASE("FFBIDXIndexerCPU_Benchmark", "[.][ffbidx_bench]") { const UnitCell uc(40, 45, 80, 90, 82.5, 90); diff --git a/tests/IndexingUnitTest.cpp b/tests/IndexingUnitTest.cpp index b00f37c2c..46d3ff7cc 100644 --- a/tests/IndexingUnitTest.cpp +++ b/tests/IndexingUnitTest.cpp @@ -9,6 +9,7 @@ #include "../image_analysis/indexing/PostIndexingRefinement.h" #include "../image_analysis/IndexAndRefine.h" #include "../image_analysis/bragg_prediction/BraggPrediction.h" +#include "../image_analysis/bragg_integration/SystematicAbsence.h" #include "../common/Logger.h" inline double round_err(double x) { @@ -1000,3 +1001,155 @@ TEST_CASE("PostIndexingRefinement_CoplanarCandidateIsRejected", "[Indexing]") { for (const auto &l: lattices) CHECK_NOTHROW(l.GetUBMatrix()); } + +namespace { + // Spots on a centred lattice: every reciprocal node of the conventional cell up to 0.3 1/A that the + // centring allows. + std::vector CentredReciprocalCloud(const CrystalLattice &l, char centering) { + std::vector spots; + for (int h = -25; h <= 25; h++) + for (int k = -25; k <= 25; k++) + for (int m = -50; m <= 50; m++) { + if (systematic_absence(h, k, m, centering)) + continue; + const Coord r = h * l.Astar() + k * l.Bstar() + m * l.Cstar(); + if (r.Length() < 0.3f) + spots.push_back(r); + } + return spots; + } + + // Spots within 0.1 of a node of the lattice + size_t IndexedCount(const std::vector &spots, const CrystalLattice &l) { + size_t n = 0; + for (const auto &s: spots) { + const float h = s * l.Vec0(), k = s * l.Vec1(), m = s * l.Vec2(); + const float dh = h - std::round(h), dk = k - std::round(k), dm = m - std::round(m); + if (dh * dh + dk * dk + dm * dm < 0.01f) + n++; + } + return n; + } + + // Refine on two candidates in the PRIMITIVE cell, as FFBIDX hands them over: a wrong orientation, + // passed first, and the true one. + std::vector RefineTrueAndWrong(const std::vector &spots, + const CrystalLattice &truth, const CrystalLattice &wrong) { + Eigen::MatrixX3 cells(6, 3); + int row = 0; + for (const auto *l: {&wrong, &truth}) + for (const Coord &v: {l->Vec0(), l->Vec1(), l->Vec2()}) + cells.row(row++) << v.x, v.y, v.z; + Eigen::VectorX scores(2); + scores << 0.2f, 0.2f; + + RefineParameters p{ + .viable_cell_min_spots = 6, + .dist_tolerance_vs_reference = 0.05f, + .reference_unit_cell = truth.GetUnitCell(), + .min_length_A = 1.0f, + .max_length_A = 1000.0f, + .min_angle_deg = 30.0f, + .max_angle_deg = 150.0f, + .indexing_tolerance = 0.1f + }; + return Refine(spots, spots.size(), cells, scores, p); + } + + bool SameBasis(const CrystalLattice &a, const CrystalLattice &b) { + return (a.Vec0() - b.Vec0()).Length() < 0.5f && (a.Vec1() - b.Vec1()).Length() < 0.5f + && (a.Vec2() - b.Vec2()).Length() < 0.5f; + } +} + +TEST_CASE("PostIndexingRefinement_RhombohedralObverseNotReverse", "[Indexing]") { + // An R lattice in the hexagonal setting: the reverse setting (the obverse one turned 180 deg about c) + // has every node of the same primitive hexagonal lattice, so in the conventional cell both index all + // of the spots. In the primitive (rhombohedral) cell the reverse one keeps only the l = 3n layers. + const UnitCell uc(50.0, 50.0, 120.0, 90, 90, 120); + const CrystalLattice obverse(uc); + const CrystalLattice reverse(obverse.Vec0() * -1.0f, obverse.Vec1() * -1.0f, obverse.Vec2()); + const auto spots = CentredReciprocalCloud(obverse, 'R'); + REQUIRE(spots.size() > 1000); + + CHECK(IndexedCount(spots, reverse) == spots.size()); + const CrystalLattice obverse_p = obverse.ToPrimitive('R'), reverse_p = reverse.ToPrimitive('R'); + CHECK(IndexedCount(spots, obverse_p) == spots.size()); + CHECK(static_cast(IndexedCount(spots, reverse_p)) / static_cast(spots.size()) + == Catch::Approx(1.0 / 3.0).margin(0.02)); + + const auto out = RefineTrueAndWrong(spots, obverse_p, reverse_p); + REQUIRE(!out.empty()); + CHECK(SameBasis(out[0], obverse_p)); + CHECK(SameBasis(out[0].FromPrimitive('R'), obverse)); +} + +TEST_CASE("PostIndexingRefinement_CentredFalseOrientationScoresLower", "[Indexing]") { + // A C-centred cell with b = c: turning it 90 deg about a keeps the metric and every node of the + // primitive lattice of the conventional cell, but not the centring - in the conventional cell that + // orientation indexes every spot, in the primitive cell only half of them. + const UnitCell uc(70.0, 50.0, 50.0, 90, 90, 90); + const CrystalLattice truth(uc); + const CrystalLattice wrong(truth.Vec0(), truth.Vec2(), truth.Vec1() * -1.0f); + const auto spots = CentredReciprocalCloud(truth, 'C'); + REQUIRE(spots.size() > 1000); + + CHECK(IndexedCount(spots, wrong) == spots.size()); + const CrystalLattice truth_p = truth.ToPrimitive('C'), wrong_p = wrong.ToPrimitive('C'); + CHECK(IndexedCount(spots, truth_p) == spots.size()); + CHECK(static_cast(IndexedCount(spots, wrong_p)) / static_cast(spots.size()) + == Catch::Approx(0.5).margin(0.02)); + + const auto out = RefineTrueAndWrong(spots, truth_p, wrong_p); + REQUIRE(!out.empty()); + CHECK(SameBasis(out[0], truth_p)); +} + +TEST_CASE("IndexAndRefine_StillRhombohedralFFBIDX", "[Indexing]") { + // The per-frame entry point every stills flow goes through (the broker's FPGA and SIMPLON receivers, + // live rugnux stills, offline rugnux): an R3 crystal in the hexagonal setting, given as its space group + // and conventional cell, indexed by FFBIDX on predicted spots. The frame must come back on the obverse + // lattice - every indexed spot obeying -h+k+l = 3n - not on the reverse one, which shares the l = 3n + // layers only. + const UnitCell uc(60.0, 60.0, 140.0, 90, 90, 120); + DiffractionExperiment experiment(DetJF4M()); + experiment.DetectorDistance_mm(120).BeamX_pxl(1090).BeamY_pxl(1136).IncidentEnergy_keV(12.4); + experiment.SetUnitCell(uc).SpaceGroupNumber(146); + experiment.IndexingAlgorithm(IndexingAlgorithmEnum::FFBIDX); + REQUIRE(experiment.GetCentering() == 'R'); + + IndexerThreadPool pool(experiment.GetIndexingSettings()); + IndexAndRefine index_and_refine(experiment, &pool); + BraggPrediction prediction; + const BraggPredictionSettings pred_settings{.high_res_A = 2.0f, .ewald_dist_cutoff = 0.001f, .centering = 'R'}; + SpotFindingSettings spot_settings; + + const std::vector> orientations = { + {10, 18, 27}, {66, -14, 101}, {-35, 40, 5}, {120, 75, -60}, {5, -80, 33}, {-150, 20, 170}}; + for (const auto &o: orientations) { + const CrystalLattice truth = RotateLattice(CrystalLattice(uc), MakeRotation(o[0], o[1], o[2])); + const int n = prediction.Calc(experiment, truth, pred_settings); + DataMessage msg{}; + for (int i = 0; i < n; i++) { + SpotToSave s{}; + s.x = prediction.GetReflections().at(i).predicted_x; + s.y = prediction.GetReflections().at(i).predicted_y; + s.intensity = 100.0f; + msg.spots.push_back(s); + } + INFO("orientation " << o[0] << " " << o[1] << " " << o[2] << ", " << msg.spots.size() << " spots"); + REQUIRE(msg.spots.size() > 30); + + REQUIRE(index_and_refine.IndexFrameOnly(msg, spot_settings)); + size_t indexed = 0, obeying = 0; + for (const auto &s: msg.spots) { + if (!s.indexed) + continue; + indexed++; + if (!systematic_absence(static_cast(s.h), static_cast(s.k), static_cast(s.l), 'R')) + obeying++; + } + CHECK(indexed > 0.75 * msg.spots.size()); + CHECK(obeying == indexed); + } +}