FFBIDXIndexerCPU: the fast-feedback indexer on the CPU, used where there is no GPU

A port of the fast-feedback indexer's CUDA kernels to plain C++/Eigen: the same half-sphere and
rotation samples, the same scores and the same 32 output cells for Refine(). The loops are
reorganised so the inner one runs over samples and vectorises (GCC/Clang, no intrinsics), and the
cell search separates the rotation about the candidate vector out of the per-spot terms. The
vector refinement reproduces the library's warp sums - partial sums its shuffle reduction leaves in
the other lanes carry into the next round - and its float arithmetic (FMA contraction, the
Fibonacci points in float), because the candidates searched depend on it: on the library's 27 test
spot lists and 4 synthetic stills the CPU and GPU best cells agree to ~1e-6 relative, the top-32
sets in 30 of 31 cases completely.

FFBIDX now resolves to the GPU library where a GPU is visible and to the CPU port otherwise (no
GPU, --no-gpu, or a build without CUDA); Auto takes FFBIDX for stills whenever a cell is known,
also without a GPU, instead of FFTW. With a GPU nothing changes (a stills set and a rotation set
reproduce the baseline exactly). Without a GPU the stills indexer pool gets one indexer per worker
instead of eight, since the CPU indexer holds no card memory.

Single core, synthetic still, loaded 5950X: ~8 ms per call at 30 spots, ~15 ms at 80, ~35 ms at
250; 16 threads ~1300/700/330 calls/s.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit 9a8eb8669d1be6c3ec16db6786280affd6830ee5)

common/DiffractionExperiment.cpp

@@ -1714,11 +1714,10 @@ IndexingAlgorithmEnum DiffractionExperiment::GetIndexingAlgorithm() const {
                 return IndexingAlgorithmEnum::None;
             return IndexingAlgorithmEnum::FFBIDX;
         case IndexingAlgorithmEnum::Auto:
-            if (get_gpu_count() == 0)
-                return IndexingAlgorithmEnum::FFTW;
-            if (!cell)
-                return IndexingAlgorithmEnum::FFT;
-            return IndexingAlgorithmEnum::FFBIDX;
+            // FFBIDX runs on the GPU where there is one and on the CPU otherwise
+            if (cell)
+                return IndexingAlgorithmEnum::FFBIDX;
+            return get_gpu_count() == 0 ? IndexingAlgorithmEnum::FFTW : IndexingAlgorithmEnum::FFT;
         case IndexingAlgorithmEnum::FFT:
             // No silent switch to FFTW, too slow for live stills: a missing GPU turns indexing off
             return get_gpu_count() == 0 ? IndexingAlgorithmEnum::None : IndexingAlgorithmEnum::FFT;

frontend/src/components/IndexingSettings.tsx

@@ -104,7 +104,7 @@ function IndexingSettings({s: serverS, status}: MyProps) {
                         value={s.algorithm}
                         onChange={handleAlgorithmChange}
                     >
-                        <MenuItem value={indexing_algorithm.FFBIDX}>FFBIDX (known unit cell) - GPU</MenuItem>
+                        <MenuItem value={indexing_algorithm.FFBIDX}>FFBIDX (known unit cell) - GPU, or CPU without one</MenuItem>
                         <MenuItem value={indexing_algorithm.FFT}>FFT (guess unit cell) - GPU with CuFFT library</MenuItem>
                         <MenuItem value={indexing_algorithm.FFTW}>FFT (guess unit cell) - CPU with FFTW library</MenuItem>
                         <MenuItem value={indexing_algorithm.AUTO}>Auto (FFBIDX if unit cell provided, FFT otherwise)</MenuItem>

image_analysis/indexing/CMakeLists.txt

@@ -69,6 +69,13 @@ ELSE()
     TARGET_LINK_LIBRARIES(JFJochIndexing Eigen3::Eigen)
 ENDIF()
 
+# The CPU fast-feedback indexer is always built. Without errno and FP traps to preserve, GCC and Clang
+# vectorise its scoring loops (sqrt, rint, min/max); neither option changes a computed value.
+TARGET_SOURCES(JFJochIndexing PRIVATE FFBIDXIndexerCPU.cpp FFBIDXIndexerCPU.h)
+IF (NOT MSVC)
+    SET_SOURCE_FILES_PROPERTIES(FFBIDXIndexerCPU.cpp PROPERTIES COMPILE_OPTIONS "-fno-math-errno;-fno-trapping-math")
+ENDIF()
+
 # FFTW (fftw3f) is always available via FetchContent -> the CPU FFT indexer is always built.
 TARGET_SOURCES(JFJochIndexing PRIVATE FFTIndexerCPU.cpp FFTIndexerCPU.h)
 TARGET_LINK_LIBRARIES(JFJochIndexing fftw3f)

image_analysis/indexing/FFBIDXIndexerCPU.cpp

@@ -0,0 +1,431 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+// Ported from the CUDA kernels of the fast-feedback indexer,
+// https://github.com/paulscherrerinstitute/fast-feedback-indexer (indexer/src/indexer_gpu.cu)
+// (c) 2022 Paul Scherrer Institute, BSD-3-Clause
+//
+// The samples and the scores are the library's; the loops are reorganised for a CPU. Each block of
+// samples is scored against one spot at a time, so the inner loop runs over samples and vectorises,
+// and in the cell search the rotation about the candidate vector is separated out of the spot terms.
+
+#include "FFBIDXIndexerCPU.h"
+
+#include <algorithm>
+#include <array>
+#include <bit>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <numbers>
+#include <numeric>
+
+#include "PostIndexingRefinement.h"
+
+namespace {
+    // The GPU library's runtime defaults (config_runtime in ffbidx/indexer.h) and rugnux's persistent
+    // configuration of it (FFBIDXIndexer.h). Keep them in step: the CPU vs GPU test compares the two.
+    constexpr int NUM_HALFSPHERE_POINTS = 32 * 1024;
+    constexpr int NUM_ANGLE_POINTS = 1024;  // the library's own choice for 32k half-sphere points
+    constexpr int NUM_CANDIDATE_VECTORS = 32;
+    constexpr int NUM_OUTPUT_CELLS = 32;
+    constexpr int MIN_SPOTS = 6;            // vector refinement
+    constexpr float TRIML = 0.001f;
+    constexpr float TRIMH = 0.3f;
+    constexpr float DELTA = 0.1f;
+    constexpr float DIST1 = 0.2f;
+    constexpr float DIST3 = 0.15f;
+    constexpr float LENGTH_THRESHOLD = 1e-9f;
+
+    // Half-sphere directions are scored in blocks of this many, so a block's accumulators stay in L1
+    constexpr int BLOCK = 1024;
+    static_assert(NUM_HALFSPHERE_POINTS % BLOCK == 0);
+
+    // Every factor of a score's product is in [TRIML + DELTA, TRIMH + DELTA]; after 8 of them the
+    // product is still a normal float
+    constexpr int RENORMALISE_EVERY = 8;
+
+    struct Tables {
+        std::array<std::vector<float>, 3> dir;      // Fibonacci points on the half sphere z >= 0
+        std::vector<float> cos_alpha, sin_alpha;    // rotation angles about a cell vector
+    };
+
+    // (a, rest) = a + x * y + rest, a Kahan sum, with the product fused into the first addition as
+    // nvcc compiles the library's
+    inline void KahanAdd(float &a, float &rest, float x, float y = 1.0f) {
+        const float s = std::fma(x, y, rest);
+        const float t = a;
+        a = t + s;
+        rest = s - (a - t);
+    }
+
+    Tables MakeTables() {
+        constexpr double pi2 = 2.0 * std::numbers::pi;
+        const double dl = 3.0 - std::sqrt(5.0);
+
+        Tables t;
+        for (auto &v: t.dir)
+            v.resize(NUM_HALFSPHERE_POINTS);
+        std::array<float, 32> dl2n;
+        for (int n = 0; n < 32; n++)
+            dl2n[n] = static_cast<float>(std::fmod(dl * std::ldexp(1.0, n), pi2));
+        for (int i = 0; i < NUM_HALFSPHERE_POINTS; i++) {
+            const float z = 1.0f - static_cast<float>(i) / static_cast<float>(NUM_HALFSPHERE_POINTS - 1);
+            const float r = std::sqrt(std::fma(-z, z, 1.0f));
+            float l = 0, rest = 0;
+            for (int n = 0; n < 32; n++)
+                if (i & (1 << n))
+                    KahanAdd(l, rest, dl2n[n]);
+            t.dir[0][i] = std::cos(l) * r;
+            t.dir[1][i] = std::sin(l) * r;
+            t.dir[2][i] = z;
+        }
+        for (int q = 0; q < NUM_ANGLE_POINTS; q++) {
+            const double alpha = q * pi2 / NUM_ANGLE_POINTS;
+            t.cos_alpha.push_back(static_cast<float>(std::cos(alpha)));
+            t.sin_alpha.push_back(static_cast<float>(std::sin(alpha)));
+        }
+        return t;
+    }
+
+    const Tables &GetTables() {
+        static const Tables t = MakeTables();
+        return t;
+    }
+
+    inline float Dist2Int(float x) {
+        return std::abs(x - std::rint(x));
+    }
+
+    // One spot's share of a sample's score: a factor of the product, and whether it counts as good
+    inline void AddSpot(float d, float good_limit, float &product, int &good) {
+        good += static_cast<int>(d < good_limit);
+        product *= std::clamp(d, TRIML, TRIMH) + DELTA;
+    }
+
+    // Moves the power of two of the product into the exponent, leaving the product in [1, 2). The
+    // library takes a log2 every 8 spots instead; this is the same sum, and it vectorises.
+    inline void Renormalise(float &product, int &exponent) {
+        const auto bits = std::bit_cast<uint32_t>(product);
+        exponent += static_cast<int>(bits >> 23) - 127;
+        product = std::bit_cast<float>((bits & 0x007fffffu) | 0x3f800000u);
+    }
+
+    // exp2(mean log2(trim(d) + delta)) - delta - (good + 1): lower is better. A NaN, from a degenerate
+    // candidate, goes last instead of breaking the sort.
+    inline float Score(float log2_sum, int good, int nspots) {
+        const float s = std::exp2(log2_sum / static_cast<float>(nspots)) - DELTA - static_cast<float>(good + 1);
+        return std::isnan(s) ? std::numeric_limits<float>::infinity() : s;
+    }
+
+    // Indices of the k lowest scores, best first
+    std::vector<int> Best(const std::vector<float> &score, int k) {
+        std::vector<int> idx(score.size());
+        std::iota(idx.begin(), idx.end(), 0);
+        auto lower = [&score](int a, int b) { return score[a] < score[b] || (score[a] == score[b] && a < b); };
+        std::nth_element(idx.begin(), idx.begin() + k, idx.end(), lower);
+        idx.resize(k);
+        std::sort(idx.begin(), idx.end(), lower);
+        return idx;
+    }
+
+    // Stage 1: every half-sphere direction as the direction of a cell vector of this length, scored on
+    // how close the spots' projections on it come to integers
+    std::vector<int> FindCandidates(const Tables &t, const std::vector<Eigen::Vector3f> &spots, float length) {
+        const int nspots = static_cast<int>(spots.size());
+        std::vector<float> score(NUM_HALFSPHERE_POINTS);
+        std::array<float, BLOCK> product;
+        std::array<int, BLOCK> exponent, good;
+
+        for (int b = 0; b < NUM_HALFSPHERE_POINTS; b += BLOCK) {
+            const float *vx = t.dir[0].data() + b;
+            const float *vy = t.dir[1].data() + b;
+            const float *vz = t.dir[2].data() + b;
+            product.fill(1.0f);
+            exponent.fill(0);
+            good.fill(0);
+            for (int i = 0; i < nspots; i++) {
+                const Eigen::Vector3f s = length * spots[i];
+                for (int k = 0; k < BLOCK; k++)
+                    AddSpot(Dist2Int(s.x() * vx[k] + s.y() * vy[k] + s.z() * vz[k]), DIST1, product[k], good[k]);
+                if ((i + 1) % RENORMALISE_EVERY == 0)
+                    for (int k = 0; k < BLOCK; k++)
+                        Renormalise(product[k], exponent[k]);
+            }
+            for (int k = 0; k < BLOCK; k++)
+                score[b + k] = Score(static_cast<float>(exponent[k]) + std::log2(product[k]), good[k], nspots);
+        }
+        return Best(score, NUM_CANDIDATE_VECTORS);
+    }
+
+    // The library's dot product: rounding matters here, as a spot on the edge of a tolerance decides
+    // which spots the candidate is refined against
+    inline float Dot(const Eigen::Vector3f &a, const Eigen::Vector3f &b) {
+        return std::fma(a[0], b[0], std::fma(a[1], b[1], a[2] * b[2]));
+    }
+
+    // A sum as the library forms it on a warp of 32 lanes: spot i is added on lane i % 32, and the
+    // lanes are reduced by shuffles that leave partial sums behind in every lane but the first. Those
+    // stay in the lanes and enter the next round's sum. Reproduced as it is, because the refined
+    // candidates - and so the cells searched - depend on it.
+    struct WarpSum {
+        std::array<float, 32> value{}, rest{};
+
+        void Add(size_t spot, float x, float y) {
+            KahanAdd(value[spot % 32], rest[spot % 32], x, y);
+        }
+
+        float Reduce() {
+            std::array<float, 32> r{};
+            for (int d = 1; d < 32; d *= 2) {
+                const auto old_value = value;
+                const auto old_r = r;
+                for (int l = 0; l < 32; l++) {
+                    const int from = (l % (2 * d) < d) ? l + d : l;
+                    r[l] = old_r[from];
+                    KahanAdd(value[l], r[l], old_value[from]);
+                }
+            }
+            return value[0];
+        }
+    };
+
+    // x = A^-1 v for symmetric A given by its upper triangle (00 01 02 11 12 22), by cofactors, as
+    // the library solves it (products fused as nvcc compiles them)
+    Eigen::Vector3f SymSolve(const std::array<float, 6> &A, const std::array<float, 3> &v) {
+        auto cofactor = [](float a, float b, float c, float d) { return std::fma(a, b, -(c * d)); };
+        const float co00 = cofactor(A[3], A[5], A[4], A[4]);
+        const float co01 = cofactor(A[1], A[5], A[2], A[4]);
+        const float co02 = cofactor(A[1], A[4], A[2], A[3]);
+        const float co11 = cofactor(A[0], A[5], A[2], A[2]);
+        const float co12 = cofactor(A[0], A[4], A[1], A[2]);
+        const float co22 = cofactor(A[0], A[3], A[1], A[1]);
+        // a0 * b0 + a1 * b1 + a2 * b2
+        auto sum3 = [](float a0, float b0, float a1, float b1, float a2, float b2) {
+            float t = a0 * b0;
+            float rest = 0;
+            KahanAdd(t, rest, a1, b1);
+            KahanAdd(t, rest, a2, b2);
+            return t;
+        };
+        const float r_det = 1.0f / sum3(A[0], co00, -A[1], co01, A[2], co02);
+        return {sum3(co00, v[0], -co01, v[1], co02, v[2]) * r_det,
+                sum3(-co01, v[0], co11, v[1], -co12, v[2]) * r_det,
+                sum3(co02, v[0], -co12, v[1], co22, v[2]) * r_det};
+    }
+
+    // Stage 2: a candidate refined by least squares against the spots it nearly indexes, in four rounds
+    // of tightening tolerance; the Miller indices stay those of the unrefined vector
+    Eigen::Vector3f RefineCandidate(const std::vector<Eigen::Vector3f> &spots, const Eigen::Vector3f &z0) {
+        std::array<WarpSum, 6> STS;
+        std::array<WarpSum, 3> STm;
+        Eigen::Vector3f z = z0;
+        for (const float beta: {0.1f, 0.05f, 0.025f, 0.01f}) {
+            int good = 0;
+            for (size_t i = 0; i < spots.size(); i++) {
+                const Eigen::Vector3f &s = spots[i];
+                const float m = std::rint(Dot(s, z0));
+                if (std::abs(Dot(s, z) - m) >= beta)
+                    continue;
+                good++;
+                int e = 0;
+                for (int j = 0; j < 3; j++) {
+                    STm[j].Add(i, s[j], m);
+                    for (int k = j; k < 3; k++)
+                        STS[e++].Add(i, s[j], s[k]);
+                }
+            }
+            if (good < MIN_SPOTS)
+                break;
+            std::array<float, 6> A;
+            std::array<float, 3> v;
+            for (int e = 0; e < 3; e++)
+                v[e] = STm[e].Reduce();
+            for (int e = 0; e < 6; e++)
+                A[e] = STS[e].Reduce();
+            z = SymSolve(A, v);
+        }
+        return z;
+    }
+
+    // The input cell turned so its vector cell_vec lies along c, in the frame (x, y, z = c/|c|) the
+    // rotation about c is sampled in: a(alpha) = a_z z + a_xy (cos alpha x + sin alpha y), and b the
+    // same at alpha + delta.
+    struct CellFrame {
+        int cell_vec;
+        Eigen::Vector3f x, y, z;
+        float c_length;
+        float a_z, a_xy, b_z, b_xy;
+        float cos_delta, sin_delta;
+
+        Eigen::Vector3f A(float cos_a, float sin_a) const {
+            return a_z * z + a_xy * (cos_a * x + sin_a * y);
+        }
+
+        Eigen::Vector3f B(float cos_a, float sin_a) const {
+            const float c = cos_a * cos_delta - sin_a * sin_delta;
+            const float s = sin_a * cos_delta + cos_a * sin_delta;
+            return b_z * z + b_xy * (c * x + s * y);
+        }
+    };
+
+    CellFrame MakeFrame(const Eigen::Matrix3f &cell, int cell_vec, const Eigen::Vector3f &c) {
+        CellFrame f;
+        f.cell_vec = cell_vec;
+        f.c_length = c.norm();
+        f.z = c / f.c_length;
+
+        // Mirroring on the bisector of the input vector and z takes one onto the other
+        const Eigen::Vector3f v = cell.row(cell_vec);
+        Eigen::Vector3f a = cell.row((cell_vec + 1) % 3);
+        Eigen::Vector3f b = cell.row((cell_vec + 2) % 3);
+        Eigen::Vector3f t = v + v.norm() * f.z;
+        if (t.norm() >= std::numeric_limits<float>::epsilon())
+            t.normalize();
+        else
+            t = f.z.cross(a).normalized();
+        a = 2.0f * a.dot(t) * t - a;
+        b = 2.0f * b.dot(t) * t - b;
+
+        f.a_z = a.dot(f.z);
+        f.x = (a - f.a_z * f.z).normalized();
+        f.a_xy = a.dot(f.x);
+        f.y = f.z.cross(f.x);
+
+        f.b_z = b.dot(f.z);
+        const Eigen::Vector3f b_dir = (b - f.b_z * f.z).normalized();
+        f.b_xy = b.dot(b_dir);
+        f.cos_delta = b_dir.dot(f.x);
+        f.sin_delta = b_dir.dot(f.y);
+        return f;
+    }
+
+    // Stage 3: the cell of one frame at every rotation angle, scored on how close each spot comes to
+    // a lattice point. The spot's projections on z, x and y are taken once; c does not turn at all.
+    void ScoreAngles(const Tables &t, const CellFrame &f, const std::vector<Eigen::Vector3f> &spots, float *score) {
+        const int nspots = static_cast<int>(spots.size());
+        const float *cos_a = t.cos_alpha.data();
+        const float *sin_a = t.sin_alpha.data();
+        std::array<float, NUM_ANGLE_POINTS> product;
+        std::array<int, NUM_ANGLE_POINTS> exponent, good;
+        product.fill(1.0f);
+        exponent.fill(0);
+        good.fill(0);
+
+        for (int i = 0; i < nspots; i++) {
+            const float zs = f.z.dot(spots[i]);
+            const float xs = f.x.dot(spots[i]);
+            const float ys = f.y.dot(spots[i]);
+            const float dc = Dist2Int(f.c_length * zs);
+            const float a0 = f.a_z * zs;
+            const float a1 = f.a_xy * xs;
+            const float a2 = f.a_xy * ys;
+            const float b0 = f.b_z * zs;
+            const float b1 = f.b_xy * (f.cos_delta * xs + f.sin_delta * ys);
+            const float b2 = f.b_xy * (f.cos_delta * ys - f.sin_delta * xs);
+            for (int q = 0; q < NUM_ANGLE_POINTS; q++) {
+                const float da = Dist2Int(a0 + a1 * cos_a[q] + a2 * sin_a[q]);
+                const float db = Dist2Int(b0 + b1 * cos_a[q] + b2 * sin_a[q]);
+                AddSpot(std::sqrt(dc * dc + da * da + db * db), DIST3, product[q], good[q]);
+            }
+            if ((i + 1) % RENORMALISE_EVERY == 0)
+                for (int q = 0; q < NUM_ANGLE_POINTS; q++)
+                    Renormalise(product[q], exponent[q]);
+        }
+        // The library starts this sum at 1, not 0; kept, so the scores agree with the GPU
+        for (int q = 0; q < NUM_ANGLE_POINTS; q++)
+            score[q] = Score(1.0f + static_cast<float>(exponent[q]) + std::log2(product[q]), good[q], nspots);
+    }
+}
+
+// Following the fast-feedback indexer (H.-C. Stadler, PSI) and TORO, Gasparotto et al. (2024)
+// J. Appl. Cryst. 57, 931-944
+FFBIDXCells FFBIDXFindCells(const Eigen::Matrix3f &cell, const std::vector<Coord> &coord, size_t nspots) {
+    const Tables &t = GetTables();
+
+    std::vector<Eigen::Vector3f> spots(nspots);
+    for (size_t i = 0; i < nspots; i++)
+        spots[i] = Eigen::Vector3f(coord[i].x, coord[i].y, coord[i].z);
+
+    // Cell vectors of the same length share one set of candidates
+    std::array<int, 3> group{};
+    std::vector<float> group_length;
+    for (int v = 0; v < 3; v++) {
+        const float length = cell.row(v).norm();
+        auto it = std::find_if(group_length.begin(), group_length.end(),
+                               [length](float l) { return std::abs(l - length) < LENGTH_THRESHOLD; });
+        group[v] = static_cast<int>(it - group_length.begin());
+        if (it == group_length.end())
+            group_length.push_back(length);
+    }
+
+    std::vector<std::vector<Eigen::Vector3f>> candidates(group_length.size());
+    for (size_t g = 0; g < group_length.size(); g++) {
+        for (const int i: FindCandidates(t, spots, group_length[g])) {
+            const Eigen::Vector3f dir(t.dir[0][i], t.dir[1][i], t.dir[2][i]);
+            candidates[g].push_back(RefineCandidate(spots, group_length[g] * dir));
+        }
+    }
+
+    // Every cell vector along every candidate of its length, both ways round
+    std::vector<CellFrame> frames;
+    for (int v = 0; v < 3; v++) {
+        for (const auto &c: candidates[group[v]]) {
+            frames.push_back(MakeFrame(cell, v, c));
+            frames.push_back(MakeFrame(cell, v, -c));
+        }
+    }
+
+    std::vector<float> score(frames.size() * NUM_ANGLE_POINTS);
+    for (size_t i = 0; i < frames.size(); i++)
+        ScoreAngles(t, frames[i], spots, score.data() + i * NUM_ANGLE_POINTS);
+
+    FFBIDXCells ret;
+    ret.cells.resize(3 * NUM_OUTPUT_CELLS, 3);
+    ret.scores.resize(NUM_OUTPUT_CELLS);
+    const auto best = Best(score, NUM_OUTPUT_CELLS);
+    for (int k = 0; k < NUM_OUTPUT_CELLS; k++) {
+        const CellFrame &f = frames[best[k] / NUM_ANGLE_POINTS];
+        const int q = best[k] % NUM_ANGLE_POINTS;
+        ret.scores(k) = score[best[k]];
+        ret.cells.row(3 * k + f.cell_vec) = f.c_length * f.z;
+        ret.cells.row(3 * k + (f.cell_vec + 1) % 3) = f.A(t.cos_alpha[q], t.sin_alpha[q]);
+        ret.cells.row(3 * k + (f.cell_vec + 2) % 3) = f.B(t.cos_alpha[q], t.sin_alpha[q]);
+    }
+    return ret;
+}
+
+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());
+    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;
+}
+
+std::vector<CrystalLattice> FFBIDXIndexerCPU::RunInternal(const std::vector<Coord> &coord, size_t nspots) {
+    if (nspots > coord.size())
+        nspots = coord.size();
+
+    if (nspots < viable_cell_min_spots)
+        return {};
+
+    auto found = FFBIDXFindCells(input_cell, coord, nspots);
+
+    RefineParameters parameters{
+        .viable_cell_min_spots = viable_cell_min_spots,
+        .dist_tolerance_vs_reference = dist_tolerance_vs_reference,
+        .reference_unit_cell = reference_unit_cell,
+        .min_length_A = 1, // doesn't matter
+        .max_length_A = 1000, // doesn't matter
+        .min_angle_deg = 30,
+        .max_angle_deg = 150,
+        .indexing_tolerance = indexing_tolerance
+    };
+
+    return Refine(coord, nspots, found.cells, found.scores, parameters);
+}

image_analysis/indexing/FFBIDXIndexerCPU.h

@@ -0,0 +1,30 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <vector>
+#include <Eigen/Dense>
+
+#include "../../common/Coord.h"
+#include "Indexer.h"
+
+// The best cells of the fast-feedback indexer for one input cell, before Refine(): three rows per
+// cell (the real-space cell vectors), scores ascending - lower is better. The layout of the GPU
+// library's oCellM()/oScoreV(), so both feed Refine() the same way.
+struct FFBIDXCells {
+    Eigen::MatrixX3<float> cells;
+    Eigen::VectorX<float> scores;
+};
+
+// The fast-feedback indexer on the CPU: same samples, scores and output as the GPU library.
+// cell holds the real-space cell vectors as rows; the first nspots of spots are used.
+FFBIDXCells FFBIDXFindCells(const Eigen::Matrix3f &cell, const std::vector<Coord> &spots, size_t nspots);
+
+class FFBIDXIndexerCPU : public Indexer {
+    Eigen::Matrix3f input_cell;
+
+    void SetupUnitCell(const std::optional<UnitCell> &cell) override;
+public:
+    std::vector<CrystalLattice> RunInternal(const std::vector<Coord> &coord, size_t nspots) override;
+};

image_analysis/indexing/IndexerFactory.cpp

@@ -2,6 +2,8 @@
 // SPDX-License-Identifier: GPL-3.0-only
 
 #include "IndexerFactory.h"
+#include "FFBIDXIndexerCPU.h"
+#include "../../common/CUDAWrapper.h"
 
 #ifdef JFJOCH_USE_CUDA
 #include "FFBIDXIndexer.h"
@@ -17,9 +19,12 @@ std::unique_ptr<Indexer> CreateIndexer(const DiffractionExperiment& experiment)
     switch (experiment.GetIndexingAlgorithm()) {
         case IndexingAlgorithmEnum::FFBIDX:
 #ifdef JFJOCH_USE_CUDA
-            indexer = std::make_unique<FFBIDXIndexer>();
-            indexer->Setup(experiment);
+            if (get_gpu_count() > 0)
+                indexer = std::make_unique<FFBIDXIndexer>();
 #endif
+            if (!indexer)
+                indexer = std::make_unique<FFBIDXIndexerCPU>();
+            indexer->Setup(experiment);
             break;
         case IndexingAlgorithmEnum::FFTW:
 #ifdef JFJOCH_USE_FFTW

image_analysis/indexing/IndexerThreadPool.cpp

@@ -11,6 +11,7 @@
 
 #include "../common/CUDAWrapper.h"
 #include "../common/Logger.h"
+#include "FFBIDXIndexerCPU.h"
 
 #ifdef JFJOCH_USE_CUDA
 #include "FFBIDXIndexer.h"
@@ -173,6 +174,8 @@ static std::unique_ptr<Indexer> MakeIndexer(IndexingAlgorithmEnum algorithm, con
             return std::make_unique<FFBIDXIndexer>();
     }
 #endif
+    if (algorithm == IndexingAlgorithmEnum::FFBIDX)
+        return std::make_unique<FFBIDXIndexerCPU>();
 #ifdef JFJOCH_USE_FFTW
     if (algorithm == IndexingAlgorithmEnum::FFTW)
         return std::make_unique<FFTIndexerCPU>(settings);
@@ -446,7 +449,7 @@ IndexerResult IndexerThreadPool::Run(const DiffractionExperiment &experiment, co
     const bool servable =
         (algorithm == IndexingAlgorithmEnum::FFT    && have_gpu &&
              (requested == IndexingAlgorithmEnum::Auto || requested == IndexingAlgorithmEnum::FFT)) ||
-        (algorithm == IndexingAlgorithmEnum::FFBIDX && have_gpu &&
+        (algorithm == IndexingAlgorithmEnum::FFBIDX &&
              (requested == IndexingAlgorithmEnum::Auto || requested == IndexingAlgorithmEnum::FFBIDX)) ||
         (algorithm == IndexingAlgorithmEnum::FFTW   && fftw_built &&
              ((requested == IndexingAlgorithmEnum::Auto && !have_gpu) || requested == IndexingAlgorithmEnum::FFTW));
@@ -457,8 +460,6 @@ IndexerResult IndexerThreadPool::Run(const DiffractionExperiment &experiment, co
                   "always uses the GPU FFT indexer here; select FFT or Auto, or run FFTW on a CPU-only node.";
         else if (algorithm == IndexingAlgorithmEnum::FFT && !have_gpu)
             msg = "FFT is the GPU indexer but no GPU is available. Select FFTW or Auto for CPU indexing.";
-        else if (algorithm == IndexingAlgorithmEnum::FFBIDX && !have_gpu)
-            msg = "FFBIDX is a GPU indexer but no GPU is available. Select FFTW or Auto for CPU indexing.";
         else if (algorithm == IndexingAlgorithmEnum::FFTW)
             msg = "FFTW (CPU) indexing was requested but this build has no FFTW indexer.";
         else

rugnux/Rugnux.cpp

@@ -593,6 +593,12 @@ int Rugnux::FirstPassIndexerThreads() const {
     return gpu_work_beside_ && gpu_takes_turns() ? 1 : 2;
 }
 
+int Rugnux::StillsIndexerThreads() const {
+    if (experiment_.GetIndexingAlgorithm() == IndexingAlgorithmEnum::FFBIDX && get_gpu_count() == 0)
+        return std::clamp(config_.nthreads, 1, 64);
+    return std::clamp(config_.nthreads, 1, 8);
+}
+
 int Rugnux::FirstPassSpotWorkers() const {
     const int gpus = get_gpu_count();
     if (gpus == 0)
@@ -1054,7 +1060,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
             // Measured on a 16 Mpx serial set, de-novo FFT: 83 frames/s with 4, 125 with 8, 133 with 16 -
             // by eight the card is what limits it. Each GPU FFT indexer holds ~0.5 GB of the card
             // (sixteen of them took the run to 12.7 GB), so eight rather than one per worker.
-            indexing_settings.IndexingThreads(std::clamp(config_.nthreads, 1, 8));
+            indexing_settings.IndexingThreads(StillsIndexerThreads());
         }
         indexing_bound_A = indexing_settings.GetFFT_MaxUnitCell_A();
         indexer_pool = std::make_unique<IndexerThreadPool>(indexing_settings, IndexerConstruction::OnFirstUse,

rugnux/Rugnux.h

@@ -1489,6 +1489,9 @@ class Rugnux {
     // An FFT indexer holds 0.3-0.6 GB of the card, and which one indexes a frame does not change what it
     // finds.
     [[nodiscard]] int FirstPassIndexerThreads() const;
+    // The indexers of a stills pool: eight where each holds card memory; one per worker for the CPU
+    // fast-feedback indexer, which holds next to nothing.
+    [[nodiscard]] int StillsIndexerThreads() const;
     // The first pass's spot-finding workers, each with GPU engines of its own where there is a card.
     [[nodiscard]] int FirstPassSpotWorkers() const;
     // What a probe pass that only indexes holds of a card at its most: its first pass's workers and FFT

rugnux/RugnuxStills.cpp

@@ -841,7 +841,7 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) {
     // As many indexer threads as the image pass has (RunPipeline): the hypotheses are indexed side by side, and
     // with the default four the workers waited on the pool.
     IndexingSettings pool_settings = experiment_.GetIndexingSettings();
-    pool_settings.IndexingThreads(std::clamp(config_.nthreads, 1, 8));
+    pool_settings.IndexingThreads(StillsIndexerThreads());
     IndexerThreadPool pool(pool_settings, IndexerConstruction::OnFirstUse, gpu_work_beside_);
 
     // A hypothesis is indexed once: the searches below meet the same one more than once - the decided centre

rugnux/rugnux_cli.cpp

@@ -138,7 +138,7 @@ void print_usage() {
     std::cout << "   -N, --threads <num>              Number of threads (default: every CPU the run may use - its affinity and cgroup quota - of which the "
                  "per-image loop uses at most 16 per GPU; an explicit value is used as given)" << std::endl;
     std::cout << "   --no-gpu                         Run on the CPU only, as on a machine without a GPU: CUDA is not "
-                 "touched and every GPU step takes its CPU path (-X Auto and FFT index with FFTW)" << std::endl;
+                 "touched and every GPU step takes its CPU path (FFT indexes with FFTW, FFBIDX with its CPU implementation)" << std::endl;
     std::cout << "   -s, --start-image <num>          Start image number (default: 0)" << std::endl;
     std::cout << "   -e, --end-image <num>            End image number, exclusive: -s 0 -e 100 processes images 0-99 (default: all)" << std::endl;
     std::cout << "   -t, --stride <num>               Image stride (default: 1)" << std::endl;
@@ -197,7 +197,9 @@ void print_usage() {
     std::cout << "   --prepass-fraction <f>           Integrate the first fraction f of the sweep, but at least 120 deg of it, in the two-pass rotation geometry pre-pass (0.05-1, default 0.5; 1 = the whole sweep); pass 2 reads every frame. The pre-scan and the first-pass indexing sample the same part of the sweep" << std::endl;
     std::cout << "   -X, --indexing-algorithm <txt>   Indexing algorithm (Auto|FFT|cuFFT|FFTW|FFBIDX|None, default Auto). FFT is "
                  "the de-novo FFT search, on the GPU (cuFFT) where there is one and --no-gpu is not given, otherwise on the "
-                 "CPU (FFTW); cuFFT and FFTW name one implementation each. cuFFT and FFBIDX (known cell) need a GPU" << std::endl;
+                 "CPU (FFTW); cuFFT and FFTW name one implementation each, and cuFFT needs a GPU. FFBIDX is the known-cell "
+                 "fast-feedback indexer for stills, likewise on the GPU where there is one and otherwise on the CPU; "
+                 "Auto picks it for stills whenever a cell is known" << std::endl;
     std::cout << "   -S, --space-group <num|symbol>   Space group number (96) or symbol (P43212) - for indexing and scaling" << std::endl;
     std::cout << "   -C, --unit-cell <cell>           Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl;
     std::cout << "   --fft-min-unit-cell <num>        Shortest cell axis the FFT search accepts, in A (default: 10). A candidate with a shorter axis is discarded, but a de-novo rotation run also tries a second first-pass hypothesis with the floor at 5 A, so a small-molecule cell is indexed without this flag; -C lowers it on its own to cover the cell given" << std::endl;
@@ -1789,10 +1791,8 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
         }
     }
 
-    if ((indexing_algorithm == IndexingAlgorithmEnum::FFT || indexing_algorithm == IndexingAlgorithmEnum::FFBIDX)
-        && get_gpu_count() == 0) {
-        logger.Error("-X {} runs on the GPU only, and {}; -X FFT or Auto index on the CPU",
-                     indexing_algorithm == IndexingAlgorithmEnum::FFT ? "cuFFT" : "FFBIDX",
+    if (indexing_algorithm == IndexingAlgorithmEnum::FFT && get_gpu_count() == 0) {
+        logger.Error("-X cuFFT runs on the GPU only, and {}; -X FFT or Auto index on the CPU",
                      no_gpu ? "--no-gpu is given" : "no GPU is visible");
         exit(EXIT_FAILURE);
     }

tests/CMakeLists.txt

@@ -35,6 +35,7 @@ ADD_EXECUTABLE(jfjoch_portable_test EXCLUDE_FROM_ALL
         LoadReferenceMtzFreeFlagsTest.cpp
         LoadReferenceSfCifTest.cpp
         RotationIndexerTest.cpp
+        FFBIDXIndexerCPUTest.cpp
 )
 target_link_libraries(jfjoch_portable_test Catch2WithMain JFJochRugnux JFJochReader JFJochImageAnalysis
                       JFJochWriter JFJochCommon)
@@ -44,6 +45,13 @@ IF (JFJOCH_PORTABLE_ONLY AND JFJOCH_CUDA_AVAILABLE AND TARGET CUDA::cufft_static
     SET_TARGET_PROPERTIES(jfjoch_portable_test PROPERTIES CUDA_RESOLVE_DEVICE_SYMBOLS ON)
 ENDIF()
 
+# The CPU vs GPU test of the fast-feedback indexer reads the library's own test spot lists
+IF (JFJOCH_CUDA_AVAILABLE)
+    FetchContent_GetProperties(fast-indexer SOURCE_DIR FFBIDX_SOURCE_DIR)
+    SET(FFBIDX_SIMPLE_DATA_DIR "${FFBIDX_SOURCE_DIR}/data/simple/files")
+    TARGET_COMPILE_DEFINITIONS(jfjoch_portable_test PRIVATE FFBIDX_SIMPLE_DATA_DIR="${FFBIDX_SIMPLE_DATA_DIR}")
+ENDIF()
+
 IF (JFJOCH_PORTABLE_ONLY)
     RETURN()
 ENDIF()
@@ -67,6 +75,7 @@ ADD_EXECUTABLE(jfjoch_test
         JFJochReceiverIntegrationTest.cpp
         AcquisitionCountersTest.cpp
         IndexingUnitTest.cpp
+        FFBIDXIndexerCPUTest.cpp
         ThreadSafeFIFOSetTest.cpp
         ReceiveRingBufferTest.cpp
         FPGANetworkTest.cpp
@@ -215,6 +224,9 @@ target_link_libraries(jfjoch_test Catch2WithMain JFJochBroker JFJochReceiver JFJ
                       JFJochRugnuxOutput JFJochRugnux JFJochModelWriter JFJochImageAnalysis JFJochCommon JFJochHLSSimulation JFJochPreview
                       jfjoch_xds_plugin ZLIB::ZLIB httplib::httplib)
 target_include_directories(jfjoch_test PRIVATE .)
+IF (JFJOCH_CUDA_AVAILABLE)
+    TARGET_COMPILE_DEFINITIONS(jfjoch_test PRIVATE FFBIDX_SIMPLE_DATA_DIR="${FFBIDX_SIMPLE_DATA_DIR}")
+ENDIF()
 
 find_package(Threads REQUIRED)
 

tests/FFBIDXIndexerCPUTest.cpp

@@ -0,0 +1,264 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+#include <chrono>
+#include <cmath>
+#include <fstream>
+#include <random>
+#include <sstream>
+#include <thread>
+
+#include <Eigen/Geometry>
+
+#include "../common/CUDAWrapper.h"
+#include "../common/DiffractionExperiment.h"
+#include "../image_analysis/indexing/FFBIDXIndexerCPU.h"
+
+#ifdef JFJOCH_USE_CUDA
+#include "../image_analysis/indexing/FFBIDXIndexer.h"
+#endif
+
+namespace {
+    // A still of a test cell in a fixed random orientation: the reciprocal lattice points to 2 A that lie
+    // 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::mt19937 rng(seed);
+        auto draw = [&rng]() {
+            return 2.0f * static_cast<float>(rng()) / static_cast<float>(std::mt19937::max()) - 1.0f;
+        };
+        const Eigen::Matrix3f rot = Eigen::Quaternionf(draw(), draw(), draw(), draw()).normalized().toRotationMatrix();
+
+        const CrystalLattice l(uc);
+        auto rotate = [&rot](const Coord &v) {
+            const Eigen::Vector3f x = rot * Eigen::Vector3f(v.x, v.y, v.z);
+            return Coord(x.x(), x.y(), x.z());
+        };
+        const Coord astar = rotate(l.Astar()), bstar = rotate(l.Bstar()), cstar = rotate(l.Cstar());
+
+        const Coord k0(0, 0, 1.0f);
+        std::vector<Coord> spots;
+        for (int h = -40; h <= 40; h++)
+            for (int k = -40; k <= 40; k++)
+                for (int m = -40; m <= 40; m++) {
+                    const Coord r = h * astar + k * bstar + m * cstar;
+                    if (r.Length() > 0.5f || r.Length() < 0.02f)
+                        continue;
+                    if (std::abs((r + k0).Length() - 1.0f) < 0.001f)
+                        spots.push_back(r);
+                }
+        while (noise_spots > 0) {
+            const Coord r(0.5f * draw(), 0.5f * draw(), 0.5f * draw());
+            if (r.Length() > 0.5f || r.Length() < 0.02f)
+                continue;
+            spots.push_back(r);
+            noise_spots--;
+        }
+        return spots;
+    }
+
+    // Fraction of spots within 0.1 of a lattice point (in Miller index units)
+    float IndexedFraction(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 static_cast<float>(n) / static_cast<float>(spots.size());
+    }
+
+    bool SameLengths(const UnitCell &a, const UnitCell &b, float rel_tol) {
+        std::array<double, 3> x = {a.a, a.b, a.c}, y = {b.a, b.b, b.c};
+        std::sort(x.begin(), x.end());
+        std::sort(y.begin(), y.end());
+        for (int i = 0; i < 3; i++)
+            if (std::abs(x[i] - y[i]) > rel_tol * y[i])
+                return false;
+        return true;
+    }
+
+    // The cell vectors as rows
+    Eigen::Matrix3f CellMatrix(const UnitCell &uc) {
+        const CrystalLattice l(uc);
+        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;
+    }
+
+    const std::vector<UnitCell> &TestCells() {
+        static const std::vector<UnitCell> cells = {
+            {30, 40, 50, 90, 90, 90},
+            {80, 80, 90, 90, 90, 120},
+            {40, 45, 80, 90, 82.5, 90},
+            {60, 60, 60, 90, 90, 90}
+        };
+        return cells;
+    }
+}
+
+TEST_CASE("Portable_FFBIDXIndexerCPU_SyntheticStill", "[portable][Indexing]") {
+    for (size_t i = 0; i < TestCells().size(); i++) {
+        const UnitCell &uc = TestCells()[i];
+        const auto lattice_spots = SyntheticStill(uc, 1234 + i, 0);
+        auto spots = SyntheticStill(uc, 1234 + i, lattice_spots.size() / 5);
+        INFO("cell " << i << ", " << lattice_spots.size() << " lattice spots");
+        REQUIRE(lattice_spots.size() > 30);
+
+        DiffractionExperiment experiment;
+        experiment.SetUnitCell(uc);
+        FFBIDXIndexerCPU indexer;
+        indexer.Setup(experiment);
+        const auto result = indexer.Run(spots);
+        REQUIRE(!result.lattice.empty());
+        CHECK(SameLengths(result.lattice[0].GetUnitCell(), uc, 0.01f));
+        CHECK(IndexedFraction(lattice_spots, result.lattice[0]) > 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);
+    DiffractionExperiment experiment;
+    experiment.SetUnitCell(uc);
+
+    for (const size_t n: {30, 80, 250}) {
+        auto spots = SyntheticStill(uc, 99, 0);
+        const auto lattice_spots = spots.size();
+        spots = SyntheticStill(uc, 99, n > lattice_spots ? n - lattice_spots : 0);
+        spots.resize(n);
+
+        const int calls = 20;
+        auto t0 = std::chrono::steady_clock::now();
+        for (int i = 0; i < calls; i++)
+            FFBIDXFindCells(CellMatrix(uc), spots, n);
+        const double find_ms = std::chrono::duration<double, std::milli>(std::chrono::steady_clock::now() - t0).count() / calls;
+
+        FFBIDXIndexerCPU indexer;
+        indexer.Setup(experiment);
+        t0 = std::chrono::steady_clock::now();
+        for (int i = 0; i < calls; i++)
+            indexer.Run(spots);
+        const double run_ms = std::chrono::duration<double, std::milli>(std::chrono::steady_clock::now() - t0).count() / calls;
+
+        const int threads = 16;
+        t0 = std::chrono::steady_clock::now();
+        std::vector<std::thread> workers;
+        for (int t = 0; t < threads; t++)
+            workers.emplace_back([&] {
+                FFBIDXIndexerCPU w;
+                w.Setup(experiment);
+                for (int i = 0; i < calls; i++)
+                    w.Run(spots);
+            });
+        for (auto &w: workers)
+            w.join();
+        const double wall_s = std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count();
+
+        WARN(n << " spots: FindCells " << find_ms << " ms, with Refine " << run_ms << " ms; "
+             << threads << " threads " << threads * calls / wall_s << " calls/s");
+    }
+}
+
+#ifdef JFJOCH_USE_CUDA
+
+namespace {
+    // The GPU library's raw output, before Refine(), with the configuration FFBIDXIndexer uses
+    FFBIDXCells RunGPU(const Eigen::Matrix3f &cell, const std::vector<Coord> &spots) {
+        fast_feedback::config_runtime<float> crt{};
+        fast_feedback::config_persistent<float> cpers{
+            .max_output_cells = 32,
+            .max_input_cells = 1,
+            .max_spots = static_cast<unsigned>(spots.size()),
+            .num_candidate_vectors = 32,
+            .redundant_computations = true,
+        };
+        fast_feedback::refine::indexer<float> indexer{cpers, crt};
+        indexer.iCellM() = cell;
+        for (size_t i = 0; i < spots.size(); i++) {
+            indexer.spotX(i) = spots[i].x;
+            indexer.spotY(i) = spots[i].y;
+            indexer.spotZ(i) = spots[i].z;
+        }
+        indexer.index(1, spots.size());
+        return {indexer.oCellM(), indexer.oScoreV()};
+    }
+
+    // Cell (rows) and spots of one of the library's test files: the first line is the cell, then a
+    // spot per line; '#' starts a comment
+    std::pair<Eigen::Matrix3f, std::vector<Coord>> ReadSimpleData(const std::string &path) {
+        std::ifstream f(path);
+        REQUIRE(f.good());
+        std::vector<float> v;
+        std::string line;
+        while (std::getline(f, line)) {
+            if (line.empty() || line[0] == '#')
+                continue;
+            std::istringstream s(line);
+            float x;
+            while (s >> x)
+                v.push_back(x);
+        }
+        REQUIRE(v.size() >= 12);
+        Eigen::Matrix3f cell;
+        cell << v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8];
+        std::vector<Coord> spots;
+        for (size_t i = 9; i + 2 < v.size(); i += 3)
+            spots.emplace_back(v[i], v[i + 1], v[i + 2]);
+        return {cell, spots};
+    }
+
+    float RelativeDifference(const Eigen::Matrix3f &a, const Eigen::Matrix3f &b) {
+        return (a - b).norm() / b.norm();
+    }
+
+    // CPU and GPU agree when the best scores agree and the CPU's best cell is one of the GPU's best
+    // cells - the samples are the same, so a matching cell agrees to rounding, and two cells within
+    // rounding of each other in score may come out in either order
+    void CompareWithGPU(const std::string &name, const Eigen::Matrix3f &cell, const std::vector<Coord> &spots) {
+        const auto gpu = RunGPU(cell, spots);
+        const auto cpu = FFBIDXFindCells(cell, spots, spots.size());
+
+        float best_match = INFINITY;
+        for (int k = 0; k < gpu.scores.size(); k++)
+            best_match = std::min(best_match, RelativeDifference(cpu.cells.block<3, 3>(0, 0), gpu.cells.block<3, 3>(3 * k, 0)));
+
+        int top_matched = 0;
+        for (int j = 0; j < cpu.scores.size(); j++)
+            for (int k = 0; k < gpu.scores.size(); k++)
+                if (RelativeDifference(cpu.cells.block<3, 3>(3 * j, 0), gpu.cells.block<3, 3>(3 * k, 0)) < 1e-3f) {
+                    top_matched++;
+                    break;
+                }
+
+        UNSCOPED_INFO(name << ": " << spots.size() << " spots, score CPU " << cpu.scores(0) << " GPU " << gpu.scores(0)
+                      << ", best cell to GPU " << best_match << ", " << top_matched << "/32 cells matched");
+        CHECK(std::abs(cpu.scores(0) - gpu.scores(0)) <= 1e-3f * std::abs(gpu.scores(0)));
+        CHECK(best_match < 1e-3f);
+        CHECK(top_matched >= 24);
+    }
+}
+
+TEST_CASE("FFBIDXIndexerCPU_MatchesGPU", "[Indexing][GPU]") {
+    if (get_gpu_count() == 0)
+        SKIP("no GPU");
+
+    for (int i = 0; i < 9; i++) {
+        for (const std::string kind: {"local", "peakfinder8", "radial"}) {
+            const std::string name = "image" + std::to_string(i) + "_" + kind + ".txt";
+            const auto [cell, spots] = ReadSimpleData(std::string(FFBIDX_SIMPLE_DATA_DIR) + "/" + name);
+            CompareWithGPU(name, cell, spots);
+        }
+    }
+
+    for (size_t i = 0; i < TestCells().size(); i++) {
+        const auto spots = SyntheticStill(TestCells()[i], 1234 + i, 10);
+        CompareWithGPU("synthetic " + std::to_string(i), CellMatrix(TestCells()[i]), spots);
+    }
+}
+
+#endif

viewer/widgets/JFJochViewerSettingsDock.cpp

@@ -371,7 +371,7 @@ QWidget *JFJochViewerSettingsDock::BuildMXPage() {
     algo_->setSizeAdjustPolicy(QComboBox::AdjustToMinimumContentsLengthWithIcon);
     algo_->setMinimumContentsLength(12);
     algo_->addItem("Auto", static_cast<int>(IndexingAlgorithmEnum::Auto));
-    algo_->addItem("FFBIDX (GPU, known cell)", static_cast<int>(IndexingAlgorithmEnum::FFBIDX));
+    algo_->addItem("FFBIDX (known cell)", static_cast<int>(IndexingAlgorithmEnum::FFBIDX));
     algo_->addItem("FFT (GPU, de-novo)", static_cast<int>(IndexingAlgorithmEnum::FFT));
     algo_->addItem("FFTW (CPU, de-novo)", static_cast<int>(IndexingAlgorithmEnum::FFTW));
     algo_->addItem("None", static_cast<int>(IndexingAlgorithmEnum::None));
@@ -765,14 +765,14 @@ void JFJochViewerSettingsDock::UpdateAlgorithmDescription() {
     QString text;
     switch (a) {
         case IndexingAlgorithmEnum::FFBIDX:
-            text = "GPU, needs a known cell — best for sparse serial stills"; break;
+            text = "needs a known cell, GPU or CPU — best for sparse serial stills"; break;
         case IndexingAlgorithmEnum::FFT:
             text = "GPU, de-novo — best for strong rotation data"; break;
         case IndexingAlgorithmEnum::FFTW:
             text = "CPU, de-novo — no GPU needed"; break;
         case IndexingAlgorithmEnum::Auto: {
-            const IndexingAlgorithmEnum r = !gpu ? IndexingAlgorithmEnum::FFTW
-                : (cell_known ? IndexingAlgorithmEnum::FFBIDX : IndexingAlgorithmEnum::FFT);
+            const IndexingAlgorithmEnum r = cell_known ? IndexingAlgorithmEnum::FFBIDX
+                : (gpu ? IndexingAlgorithmEnum::FFT : IndexingAlgorithmEnum::FFTW);
             const QString rn = r == IndexingAlgorithmEnum::FFBIDX ? "FFBIDX"
                 : r == IndexingAlgorithmEnum::FFT ? "FFT" : "FFTW";
             text = QString("resolves to %1 here (%2, cell %3)")