diff --git a/image_analysis/indexing/fast-feedback-indexer/indexer_gpu.cu b/image_analysis/indexing/fast-feedback-indexer/indexer_gpu.cu
index f21ed7c57..750fdfbc9 100644
--- a/image_analysis/indexing/fast-feedback-indexer/indexer_gpu.cu
+++ b/image_analysis/indexing/fast-feedback-indexer/indexer_gpu.cu
@@ -1228,10 +1228,18 @@ namespace {
         //                 A[a] = B[b0]; b0 += 1
         //         a += 1
         // ----------------------------
-        if (top_val[n_cand-1] <= cand[0].value)
+        // Candidates are ordered by value, equal values by sample index: a total order, so the global top
+        // list is the same whatever order the blocks merge in (a value tie would otherwise go to the block
+        // that happened to merge first). The block-local lists are in this order already: the block's
+        // radix sort is stable and its samples ascend with the thread index.
+        auto less = [](float_type va, unsigned sa, float_type vb, unsigned sb) {
+            return va < vb || (va == vb && sa < sb);
+        };
+
+        if (!less(cand[0].value, cand[0].sample, top_val[n_cand-1], top_sample[n_cand-1]))
             return;
 
-        if (cand[n_cand-1].value <= top_val[0]) {
+        if (!less(top_val[0], top_sample[0], cand[n_cand-1].value, cand[n_cand-1].sample)) {
             for (unsigned i=0; i<n_cand; ++i) {
                 top_val[i] = cand[i].value;
                 top_sample[i] = cand[i].sample;
@@ -1246,13 +1254,14 @@ namespace {
 
         while (top < n_cand) {
             auto val = top_val[top];
+            const unsigned val_sample = top_sample[top];
             if (cand1 < cand2) {
-                if (val < cand[cand1].value) {
+                if (less(val, val_sample, cand[cand1].value, cand[cand1].sample)) {
                     cand[cand1].value = val;
                     cand[cand1].sample = top_sample[top];
                     --cand1;
                 }
-                if (cand[cand0].value < cand[cand2].value) {
+                if (less(cand[cand0].value, cand[cand0].sample, cand[cand2].value, cand[cand2].sample)) {
                     top_val[top] = cand[cand0].value;
                     top_sample[top] = cand[cand0].sample;
                     ++cand0;
@@ -1261,8 +1270,8 @@ namespace {
                     top_sample[top] = cand[cand2].sample;
                     --cand2;
                 }
-            } else if (val > cand[cand0].value) {
-                if (val < cand[cand1].value) {
+            } else if (less(cand[cand0].value, cand[cand0].sample, val, val_sample)) {
+                if (less(val, val_sample, cand[cand1].value, cand[cand1].sample)) {
                     cand[cand1].value = val;
                     cand[cand1].sample = top_sample[top];
                     --cand1;
@@ -1289,10 +1298,41 @@ namespace {
                                                     cell_cand_t<float_type>* __restrict__ cand,
                                                     const unsigned n_cand) noexcept
     {
-        if (out.score[n_cand-1] <= cand[0].value)
+        // A total order, as in merge_top_vecs: equal values by cell vector (with its flip bit), candidate
+        // vector and rotation sample, so a tie does not go to whichever block merged first. Within a block
+        // only the rotation sample varies, and the block's stable radix sort leaves it ascending on ties.
+        auto key_less = [](unsigned ca, unsigned va, unsigned ra, unsigned cb, unsigned vb, unsigned rb) {
+            return ca != cb ? ca < cb : (va != vb ? va < vb : ra < rb);
+        };
+        auto out_key = [&](unsigned i, unsigned& c, unsigned& v, unsigned& r) {
+            c = *reinterpret_cast<unsigned*>(&out.z[3u * i]);
+            v = *reinterpret_cast<unsigned*>(&out.x[3u * i]);
+            r = *reinterpret_cast<unsigned*>(&out.y[3u * i]);
+        };
+        auto out_less_cand = [&](unsigned i, const cell_cand_t<float_type>& b) {
+            if (out.score[i] != b.value)
+                return out.score[i] < b.value;
+            unsigned c, v, r;
+            out_key(i, c, v, r);
+            return key_less(c, v, r, b.cell_vec, b.vsample, b.rsample);
+        };
+        auto cand_less_out = [&](const cell_cand_t<float_type>& a, unsigned i) {
+            if (a.value != out.score[i])
+                return a.value < out.score[i];
+            unsigned c, v, r;
+            out_key(i, c, v, r);
+            return key_less(a.cell_vec, a.vsample, a.rsample, c, v, r);
+        };
+        auto cand_less = [&](const cell_cand_t<float_type>& a, const cell_cand_t<float_type>& b) {
+            if (a.value != b.value)
+                return a.value < b.value;
+            return key_less(a.cell_vec, a.vsample, a.rsample, b.cell_vec, b.vsample, b.rsample);
+        };
+
+        if (!cand_less_out(cand[0], n_cand-1))
             return;
 
-        if (cand[n_cand-1].value <= out.score[0]) {
+        if (!out_less_cand(0, cand[n_cand-1])) {
             for (unsigned i=0u; i<n_cand; i++) {
                 out.score[i] = cand[i].value;
                 *reinterpret_cast<unsigned*>(&out.x[3u * i]) = cand[i].vsample;
@@ -1310,14 +1350,14 @@ namespace {
         while (top < n_cand) {
             auto val = out.score[top];
             if (cand1 < cand2) {
-                if (val < cand[cand1].value) {
+                if (out_less_cand(top, cand[cand1])) {
                     cand[cand1].value = val;
                     cand[cand1].vsample = *reinterpret_cast<unsigned*>(&out.x[3u * top]);
                     cand[cand1].rsample = *reinterpret_cast<unsigned*>(&out.y[3u * top]);
                     cand[cand1].cell_vec = *reinterpret_cast<unsigned*>(&out.z[3u * top]);
                     --cand1;
                 }
-                if (cand[cand0].value < cand[cand2].value) {
+                if (cand_less(cand[cand0], cand[cand2])) {
                     out.score[top] = cand[cand0].value;
                     *reinterpret_cast<unsigned*>(&out.x[3u * top]) = cand[cand0].vsample;
                     *reinterpret_cast<unsigned*>(&out.y[3u * top]) = cand[cand0].rsample;
@@ -1330,8 +1370,8 @@ namespace {
                     *reinterpret_cast<unsigned*>(&out.z[3u * top]) = cand[cand2].cell_vec;
                     --cand2;
                 }
-            } else if (val > cand[cand0].value) {
-                if (val < cand[cand1].value) {
+            } else if (cand_less_out(cand[cand0], top)) {
+                if (out_less_cand(top, cand[cand1])) {
                     cand[cand1].value = val;
                     cand[cand1].vsample = *reinterpret_cast<unsigned*>(&out.x[3u * top]);
                     cand[cand1].rsample = *reinterpret_cast<unsigned*>(&out.y[3u * top]);