diff --git a/image_analysis/indexing/PostIndexingRefinement.cpp b/image_analysis/indexing/PostIndexingRefinement.cpp
index 59a51d709..20dd2823e 100644
--- a/image_analysis/indexing/PostIndexingRefinement.cpp
+++ b/image_analysis/indexing/PostIndexingRefinement.cpp
@@ -3,6 +3,7 @@
 
 #include "../../common/JFJochMath.h"
 #include "PostIndexingRefinement.h"
+#include "../../common/ParallelFor.h"
 
 #include <iostream>
 #include <limits>
@@ -335,7 +336,11 @@ std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
         reference_cells.push_back(reference.NiggliReduce().GetUnitCell());
     }
 
-   for (int i = 0; i < scores.size(); i++) {
+    // Each candidate is tested on its own - the step-down to the reference and the indexed mask are the
+    // costly parts - so they are spread over the same threads as the refinement above, each into its own
+    // slot, and kept in candidate order.
+    std::vector<std::optional<RefinedCandidate>> tested(scores.size());
+    ParallelFor(static_cast<int>(scores.size()), p.refine_threads, [&](int i) {
         Eigen::Matrix3f cell_rows = oCell.block(3u * i, 0u, 3u, 3u);
         Eigen::Matrix3f cell_cols = cell_rows.transpose();
         Eigen::Vector3f row_norms = cell_rows.rowwise().norm();
@@ -343,18 +348,18 @@ std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
         if (!reference_cells.empty()) {
             if (!MatchesReference(cell_rows, reference_cells, p.dist_tolerance_vs_reference)) {
                 if (!p.step_down_to_reference)
-                    continue;
+                    return;
                 const auto sub = StepDownToReference(cell_rows, reference_cells, p.dist_tolerance_vs_reference,
                                                      spots.topRows(nspots), p.indexing_tolerance);
                 if (!sub)
-                    continue;
+                    return;
                 cell_rows = *sub;
                 cell_cols = cell_rows.transpose();
                 row_norms = cell_rows.rowwise().norm();
             }
         } else {
             if (row_norms.minCoeff() < p.min_length_A || row_norms.maxCoeff() > p.max_length_A)
-                continue;
+                return;
         }
 
         // Filter for wrong angles. Compared as COSINES, not angles: acos is strictly decreasing on
@@ -369,28 +374,31 @@ std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
         if (cos_alpha > cos_min_angle || cos_alpha < cos_max_angle ||
             cos_beta > cos_min_angle || cos_beta < cos_max_angle ||
             cos_gamma > cos_min_angle || cos_gamma < cos_max_angle)
-            continue;
+            return;
 
         // Three coplanar rows pass both tests above - 120 deg apart in one plane is inside any angle
         // window - and the refinement can walk a candidate there. Such a cell has no reciprocal basis,
         // so it is not a lattice: the same volume test as CrystalLattice::VolumeFraction().
         if (std::abs(cell_rows.determinant()) < MIN_BASIS_VOLUME_FRACTION * row_norms.prod())
-            continue;
+            return;
 
         int64_t indexed_spot_count = 0;
         auto indexed_mask = ComputeIndexedMask(spots.topRows(nspots), cell_cols, p.indexing_tolerance, indexed_spot_count);
 
         if (indexed_spot_count < p.viable_cell_min_spots)
-            continue;
+            return;
 
-        candidates.emplace_back(RefinedCandidate{
+        tested[i] = RefinedCandidate{
             .cell = cell_rows,
             .score = scores(i),
             .volume = std::abs(cell_rows.determinant()),
             .indexed_spot_count = indexed_spot_count,
             .indexed_mask = std::move(indexed_mask)
-        });
-    }
+        };
+    });
+    for (auto &t : tested)
+        if (t)
+            candidates.push_back(std::move(*t));
 
     std::sort(candidates.begin(), candidates.end(),
               [](const RefinedCandidate &a, const RefinedCandidate &b) {