Merge branch 'ffbidx-centring' into rc182

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UqXYdJKPUFLoeAchjzv3x4

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<float> 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<int>(std::rint(hf)),
+                                                                     static_cast<int>(std::rint(kf)),
+                                                                     static_cast<int>(std::rint(lf)),
+                                                                     seed_centering)))
+                    ++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); }

image_analysis/indexing/FFBIDXIndexer.cpp

@@ -11,9 +11,10 @@ void FFBIDXIndexer::SetupUnitCell(const std::optional<UnitCell> &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<CrystalLattice> FFBIDXIndexer::RunInternal(const std::vector<Coord>
         .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);
 }

image_analysis/indexing/FFBIDXIndexerCPU.cpp

@@ -399,9 +399,10 @@ void FFBIDXIndexerCPU::SetupUnitCell(const std::optional<UnitCell> &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<CrystalLattice> FFBIDXIndexerCPU::RunInternal(const std::vector<Coor
         .indexing_tolerance = indexing_tolerance
     };
 
-    return Refine(coord, nspots, found.cells, found.scores, parameters);
+    auto lattices = Refine(coord, nspots, found.cells, found.scores, parameters);
+    return ToConventional(lattices);
 }

image_analysis/indexing/Indexer.cpp

@@ -16,6 +16,7 @@ void Indexer::Setup(const DiffractionExperiment& experiment) {
     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());
 }
 
@@ -41,3 +42,36 @@ IndexerResult Indexer::Run(const std::vector<Coord> &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<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;
+}

image_analysis/indexing/Indexer.h

@@ -5,6 +5,7 @@
 
 #include <exception>
 #include <optional>
+#include <Eigen/Dense>
 
 #include "../../common/CrystalLattice.h"
 
@@ -48,6 +49,19 @@ protected:
 
     DiffractionGeometry geom;
     std::optional<UnitCell> 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<CrystalLattice> ToConventional(const std::vector<CrystalLattice> &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<Coord> spindle_axis;

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<Coord> SyntheticStill(const UnitCell &uc, uint32_t seed, size_t noise_spots) {
+    std::vector<Coord> 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<float>(rng()) / static_cast<float>(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<Case> 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<int>(std::round(h)), static_cast<int>(std::round(k)),
+                                           static_cast<int>(std::round(m)), c.centering))
+                    n++;
+            }
+            CHECK(static_cast<float>(n) / static_cast<float>(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);

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<Coord> CentredReciprocalCloud(const CrystalLattice &l, char centering) {
+        std::vector<Coord> 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<Coord> &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<CrystalLattice> RefineTrueAndWrong(const std::vector<Coord> &spots,
+                                                   const CrystalLattice &truth, const CrystalLattice &wrong) {
+        Eigen::MatrixX3<float> 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<float> 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<float>(IndexedCount(spots, reverse_p)) / static_cast<float>(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<float>(IndexedCount(spots, wrong_p)) / static_cast<float>(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<std::array<float, 3>> 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<int>(s.h), static_cast<int>(s.k), static_cast<int>(s.l), 'R'))
+                obeying++;
+        }
+        CHECK(indexed > 0.75 * msg.spots.size());
+        CHECK(obeying == indexed);
+    }
+}