diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp
index a21b1ca1e..488a8f1fb 100644
--- a/common/DiffractionExperiment.cpp
+++ b/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;
diff --git a/frontend/src/components/IndexingSettings.tsx b/frontend/src/components/IndexingSettings.tsx
index 93dbbb849..adae9df31 100644
--- a/frontend/src/components/IndexingSettings.tsx
+++ b/frontend/src/components/IndexingSettings.tsx
@@ -104,7 +104,7 @@ function IndexingSettings({s: serverS, status}: MyProps) {
value={s.algorithm}
onChange={handleAlgorithmChange}
>
-
+
diff --git a/image_analysis/indexing/CMakeLists.txt b/image_analysis/indexing/CMakeLists.txt
index 5f7451c55..212810ce4 100644
--- a/image_analysis/indexing/CMakeLists.txt
+++ b/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)
diff --git a/image_analysis/indexing/FFBIDXIndexerCPU.cpp b/image_analysis/indexing/FFBIDXIndexerCPU.cpp
new file mode 100644
index 000000000..244c149a8
--- /dev/null
+++ b/image_analysis/indexing/FFBIDXIndexerCPU.cpp
@@ -0,0 +1,431 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute
+// 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
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#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, 3> dir; // Fibonacci points on the half sphere z >= 0
+ std::vector 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 dl2n;
+ for (int n = 0; n < 32; n++)
+ dl2n[n] = static_cast(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(i) / static_cast(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(std::cos(alpha)));
+ t.sin_alpha.push_back(static_cast(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(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(product);
+ exponent += static_cast(bits >> 23) - 127;
+ product = std::bit_cast((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(nspots)) - DELTA - static_cast(good + 1);
+ return std::isnan(s) ? std::numeric_limits::infinity() : s;
+ }
+
+ // Indices of the k lowest scores, best first
+ std::vector Best(const std::vector &score, int k) {
+ std::vector 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 FindCandidates(const Tables &t, const std::vector &spots, float length) {
+ const int nspots = static_cast(spots.size());
+ std::vector score(NUM_HALFSPHERE_POINTS);
+ std::array product;
+ std::array 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(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 value{}, rest{};
+
+ void Add(size_t spot, float x, float y) {
+ KahanAdd(value[spot % 32], rest[spot % 32], x, y);
+ }
+
+ float Reduce() {
+ std::array 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 &A, const std::array &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 &spots, const Eigen::Vector3f &z0) {
+ std::array STS;
+ std::array 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 A;
+ std::array 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::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 &spots, float *score) {
+ const int nspots = static_cast(spots.size());
+ const float *cos_a = t.cos_alpha.data();
+ const float *sin_a = t.sin_alpha.data();
+ std::array product;
+ std::array 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(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, size_t nspots) {
+ const Tables &t = GetTables();
+
+ std::vector 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 group{};
+ std::vector 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(it - group_length.begin());
+ if (it == group_length.end())
+ group_length.push_back(length);
+ }
+
+ std::vector> 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 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 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 &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 FFBIDXIndexerCPU::RunInternal(const std::vector &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);
+}
diff --git a/image_analysis/indexing/FFBIDXIndexerCPU.h b/image_analysis/indexing/FFBIDXIndexerCPU.h
new file mode 100644
index 000000000..2136bbf15
--- /dev/null
+++ b/image_analysis/indexing/FFBIDXIndexerCPU.h
@@ -0,0 +1,30 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include
+#include
+
+#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 cells;
+ Eigen::VectorX 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 &spots, size_t nspots);
+
+class FFBIDXIndexerCPU : public Indexer {
+ Eigen::Matrix3f input_cell;
+
+ void SetupUnitCell(const std::optional &cell) override;
+public:
+ std::vector RunInternal(const std::vector &coord, size_t nspots) override;
+};
diff --git a/image_analysis/indexing/IndexerFactory.cpp b/image_analysis/indexing/IndexerFactory.cpp
index cf8adecb2..52be026d3 100644
--- a/image_analysis/indexing/IndexerFactory.cpp
+++ b/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 CreateIndexer(const DiffractionExperiment& experiment)
switch (experiment.GetIndexingAlgorithm()) {
case IndexingAlgorithmEnum::FFBIDX:
#ifdef JFJOCH_USE_CUDA
- indexer = std::make_unique();
- indexer->Setup(experiment);
+ if (get_gpu_count() > 0)
+ indexer = std::make_unique();
#endif
+ if (!indexer)
+ indexer = std::make_unique();
+ indexer->Setup(experiment);
break;
case IndexingAlgorithmEnum::FFTW:
#ifdef JFJOCH_USE_FFTW
diff --git a/image_analysis/indexing/IndexerThreadPool.cpp b/image_analysis/indexing/IndexerThreadPool.cpp
index bd8e2696d..ea4150583 100644
--- a/image_analysis/indexing/IndexerThreadPool.cpp
+++ b/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 MakeIndexer(IndexingAlgorithmEnum algorithm, con
return std::make_unique();
}
#endif
+ if (algorithm == IndexingAlgorithmEnum::FFBIDX)
+ return std::make_unique();
#ifdef JFJOCH_USE_FFTW
if (algorithm == IndexingAlgorithmEnum::FFTW)
return std::make_unique(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
diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp
index a1401382a..e7aa66f59 100644
--- a/rugnux/Rugnux.cpp
+++ b/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(indexing_settings, IndexerConstruction::OnFirstUse,
diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h
index 20d350226..9b539dfc2 100644
--- a/rugnux/Rugnux.h
+++ b/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
diff --git a/rugnux/RugnuxStills.cpp b/rugnux/RugnuxStills.cpp
index 26d4f3d08..29872e86b 100644
--- a/rugnux/RugnuxStills.cpp
+++ b/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
diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp
index c6c8468be..9a9f45e01 100644
--- a/rugnux/rugnux_cli.cpp
+++ b/rugnux/rugnux_cli.cpp
@@ -138,7 +138,7 @@ void print_usage() {
std::cout << " -N, --threads 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 Start image number (default: 0)" << std::endl;
std::cout << " -e, --end-image End image number, exclusive: -s 0 -e 100 processes images 0-99 (default: all)" << std::endl;
std::cout << " -t, --stride Image stride (default: 1)" << std::endl;
@@ -197,7 +197,9 @@ void print_usage() {
std::cout << " --prepass-fraction 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 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 Space group number (96) or symbol (P43212) - for indexing and scaling" << std::endl;
std::cout << " -C, --unit-cell | Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl;
std::cout << " --fft-min-unit-cell 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);
}
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 70744345f..cecfe9c5a 100644
--- a/tests/CMakeLists.txt
+++ b/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)
diff --git a/tests/FFBIDXIndexerCPUTest.cpp b/tests/FFBIDXIndexerCPUTest.cpp
new file mode 100644
index 000000000..3c91c1f1e
--- /dev/null
+++ b/tests/FFBIDXIndexerCPUTest.cpp
@@ -0,0 +1,264 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+#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 SyntheticStill(const UnitCell &uc, uint32_t seed, size_t noise_spots) {
+ std::mt19937 rng(seed);
+ auto draw = [&rng]() {
+ return 2.0f * static_cast(rng()) / static_cast(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 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 &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(n) / static_cast(spots.size());
+ }
+
+ bool SameLengths(const UnitCell &a, const UnitCell &b, float rel_tol) {
+ std::array 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 &TestCells() {
+ static const std::vector 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(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(std::chrono::steady_clock::now() - t0).count() / calls;
+
+ const int threads = 16;
+ t0 = std::chrono::steady_clock::now();
+ std::vector 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(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 &spots) {
+ fast_feedback::config_runtime crt{};
+ fast_feedback::config_persistent cpers{
+ .max_output_cells = 32,
+ .max_input_cells = 1,
+ .max_spots = static_cast(spots.size()),
+ .num_candidate_vectors = 32,
+ .redundant_computations = true,
+ };
+ fast_feedback::refine::indexer 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> ReadSimpleData(const std::string &path) {
+ std::ifstream f(path);
+ REQUIRE(f.good());
+ std::vector 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 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 &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
diff --git a/viewer/widgets/JFJochViewerSettingsDock.cpp b/viewer/widgets/JFJochViewerSettingsDock.cpp
index b953a1d5d..0b66f75ca 100644
--- a/viewer/widgets/JFJochViewerSettingsDock.cpp
+++ b/viewer/widgets/JFJochViewerSettingsDock.cpp
@@ -371,7 +371,7 @@ QWidget *JFJochViewerSettingsDock::BuildMXPage() {
algo_->setSizeAdjustPolicy(QComboBox::AdjustToMinimumContentsLengthWithIcon);
algo_->setMinimumContentsLength(12);
algo_->addItem("Auto", static_cast(IndexingAlgorithmEnum::Auto));
- algo_->addItem("FFBIDX (GPU, known cell)", static_cast(IndexingAlgorithmEnum::FFBIDX));
+ algo_->addItem("FFBIDX (known cell)", static_cast(IndexingAlgorithmEnum::FFBIDX));
algo_->addItem("FFT (GPU, de-novo)", static_cast(IndexingAlgorithmEnum::FFT));
algo_->addItem("FFTW (CPU, de-novo)", static_cast(IndexingAlgorithmEnum::FFTW));
algo_->addItem("None", static_cast(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)")
|