diff --git a/image_analysis/IndexerWrapper.cpp b/image_analysis/IndexerWrapper.cpp
index e37a4d71..ec11922e 100644
--- a/image_analysis/IndexerWrapper.cpp
+++ b/image_analysis/IndexerWrapper.cpp
@@ -9,11 +9,27 @@ void IndexerWrapper::Setup(const UnitCell &cell) {
 #endif
 }
 
+// Select spots that belong to indexing solution
+// - cell       cell in real space
+// - spots      spots in reciprocal space
+// - threshold  radius around approximated miller indices
+template <typename Mat3, typename MatX3, typename float_type=typename Mat3::Scalar>
+inline auto select_indexed_spots(const Eigen::MatrixBase<Mat3>& cell,
+                                 const Eigen::MatrixBase<MatX3>& spots,
+                                 float_type threshold=.02f)
+{
+    using M3x = Eigen::MatrixX3<float_type>;
+    M3x resid = spots * cell.transpose();
+    const M3x miller = round(resid.array());
+    resid -= miller;
+    return resid.rowwise().norm().array() < threshold;
+}
+
 std::vector<IndexingResult> IndexerWrapper::Run(const std::vector<Coord> &coord) {
 #ifdef JFJOCH_USE_CUDA
     std::vector<IndexingResult> ret;
 
-    if (coord.size() <  MIN_SPOTS_TO_INDEX)
+    if (coord.size() <= min_spots)
         return ret;
 
     size_t nspots = std::min<size_t>(MAX_SPOTS_TO_INDEX, coord.size());
@@ -30,10 +46,21 @@ std::vector<IndexingResult> IndexerWrapper::Run(const std::vector<Coord> &coord)
     // Get best cell
     auto id = fast_feedback::refine::best_cell(indexer.oScoreV());
 
-    // Check if is viable
-    if (fast_feedback::refine::is_viable_cell(indexer.oCell(id), indexer.Spots(), 0.05f, 9u)) {
+    // Get indexed spots
+    auto arr = select_indexed_spots(indexer.oCell(id), indexer.Spots(), threshold);
+
+    auto indexed_spot_count = arr.count();
+    // Check if result is viable
+    if (indexed_spot_count > min_spots) {
         IndexingResult result;
         result.l = CrystalLattice(indexer.oCell(id));
+
+        result.indexed_spots.resize(coord.size());
+        for (int i = 0; i < arr.size(); i++)
+            result.indexed_spots[i] = arr[i];
+
+        result.indexed_spots_count = indexed_spot_count;
+
         ret.emplace_back(result);
     }
 
diff --git a/image_analysis/IndexerWrapper.h b/image_analysis/IndexerWrapper.h
index c9890a38..e2afb100 100644
--- a/image_analysis/IndexerWrapper.h
+++ b/image_analysis/IndexerWrapper.h
@@ -18,6 +18,8 @@
 
 struct IndexingResult {
     CrystalLattice l;
+    std::vector<bool> indexed_spots;
+    uint64_t indexed_spots_count;
 };
 
 class IndexerWrapper {
@@ -33,6 +35,8 @@ class IndexerWrapper {
     fast_feedback::refine::config_ifss<float> conf_ifss{};
     fast_feedback::refine::indexer_ifss<float> indexer{cpers, crt, conf_ifss};
 #endif
+    constexpr const static uint64_t min_spots = 9;
+    constexpr const static float threshold = 0.05f;
 public:
     void Setup(const UnitCell &cell);
     std::vector<IndexingResult> Run(const std::vector<Coord> &coord);
diff --git a/tests/IndexingUnitTest.cpp b/tests/IndexingUnitTest.cpp
index a89b11af..01243b71 100644
--- a/tests/IndexingUnitTest.cpp
+++ b/tests/IndexingUnitTest.cpp
@@ -94,6 +94,7 @@ TEST_CASE("FastFeedbackIndexer","[Indexing]") {
         //REQUIRE(c.a == Approx(uc.a));
         //REQUIRE(c.b == Approx(uc.b));
         //REQUIRE(c.c == Approx(uc.c));
+        REQUIRE(ret[0].indexed_spots_count == recip.size());
 
         double err[3] = {0.0, 0.0, 0.0};
         for (const auto &iter: recip) {