diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp
index 2bf090d1c..4a67037e1 100644
--- a/image_analysis/rotation_indexer/RotationIndexer.cpp
+++ b/image_analysis/rotation_indexer/RotationIndexer.cpp
@@ -99,17 +99,33 @@ namespace {
     // Miller indices within tol. Comparing a symmetry-constrained refinement against an
     // unconstrained (triclinic) one is a data-driven test for a false promotion: a wrong
     // higher-symmetry constraint snaps a pseudo cell onto ideal angles and misplaces most spots.
-    float IndexedFraction(const CrystalLattice &latt, const std::vector<Coord> &coords, float tol) {
+    //
+    // A centred lattice is counted with its centring: a spot whose index on the conventional cell is
+    // one the centring extinguishes is not on the lattice. Without that a primitive lattice labelled
+    // centred indexes exactly what the same cell labelled P does, at half the primitive volume - and
+    // the sub-cell rule of the candidate selection takes the smaller of two cells that index alike.
+    float IndexedFraction(const CrystalLattice &latt, const std::vector<Coord> &coords, float tol,
+                          char centering = 'P') {
         if (coords.empty())
             return 0.0f;
         const Coord a = latt.Vec0(), b = latt.Vec1(), c = latt.Vec2();
+        // Rows: the primitive cell's Miller indices as combinations of the conventional ones.
+        const gemmi::Mat33 to_prim = gemmi::rot_as_mat33(gemmi::centred_to_primitive(centering)).transpose();
         const float tol_sq = tol * tol;
         size_t indexed = 0;
         for (const Coord &s : coords) {
-            const float dh = a * s - std::round(a * s);   // Coord operator* = dot product = Miller index
-            const float dk = b * s - std::round(b * s);
-            const float dl = c * s - std::round(c * s);
-            if (dh * dh + dk * dk + dl * dl < tol_sq)
+            const float h = std::round(a * s), k = std::round(b * s), l = std::round(c * s);
+            const float dh = a * s - h;   // Coord operator* = dot product = Miller index
+            const float dk = b * s - k;
+            const float dl = c * s - l;
+            if (dh * dh + dk * dk + dl * dl >= tol_sq)
+                continue;
+            bool allowed = true;
+            for (int i = 0; i < 3 && centering != 'P'; i++) {
+                const double p = to_prim[i][0] * h + to_prim[i][1] * k + to_prim[i][2] * l;
+                allowed = allowed && std::abs(p - std::round(p)) < 0.01;
+            }
+            if (allowed)
                 ++indexed;
         }
         return static_cast<float>(indexed) / static_cast<float>(coords.size());
@@ -421,18 +437,20 @@ void RotationIndexer::RunIndexingLocked() {
         // Refine (constrained metric solve + score by the refined-geometry indexed fraction, the
         // reliable discriminator). Runs on its own thread per solve.
         struct Solved { bool ok = false; float frac = 0.0f; XtalOptimizerData data; };
-        auto solve = [&](XtalOptimizerData d) -> Solved {
+        auto solve = [&](XtalOptimizerData d, char centering) -> Solved {
             const bool ok = XtalOptimizer(d, v_, kCeresThreads);
-            const float frac = ok ? IndexedFraction(d.latt, accumulate(d.geom, d.axis), index_tol) : 0.0f;
+            const float frac = ok ? IndexedFraction(d.latt, accumulate(d.geom, d.axis), index_tol, centering)
+                                  : 0.0f;
             return {ok, frac, std::move(d)};
         };
         std::vector<std::future<Solved>> constrained_f(n_try), tri_f(n_try);
         for (size_t ci = 0; ci < n_try; ci++) {
             if (!work[ci].viable)
                 continue;
-            constrained_f[ci] = std::async(std::launch::async, [&, ci] { return solve(work[ci].constrained); });
+            constrained_f[ci] = std::async(std::launch::async,
+                                           [&, ci] { return solve(work[ci].constrained, work[ci].sr.centering); });
             if (work[ci].has_tri)
-                tri_f[ci] = std::async(std::launch::async, [&, ci] { return solve(work[ci].tri); });
+                tri_f[ci] = std::async(std::launch::async, [&, ci] { return solve(work[ci].tri, 'P'); });
         }
 
         // Assemble and select serially, in candidate order - identical to refining them one by one.