FFTIndexer: Fixes to handle better multiple lattices

common/Coord.cpp

@@ -152,12 +152,6 @@ void Coord::swap(Coord &other) noexcept {
     std::swap(z, other.z);
 }
 
-// Then outside the class but in the same namespace:
-inline void swap(Coord& a, Coord& b) noexcept {
-    a.swap(b);
-}
-
-
 RotMatrix::RotMatrix() {
     for (int i = 0; i < 3; i++) {
         for (int j = 0; j < 3; j++)

common/Coord.h

@@ -59,4 +59,8 @@ public:
     std::vector<float> arr() const;
 };
 
+inline void swap(Coord& a, Coord& b) noexcept {
+    a.swap(b);
+}
+
 #endif //INDEX_COORD_H

image_analysis/indexing/FFTIndexer.cpp

@@ -89,7 +89,6 @@ std::vector<CrystalLattice> FFTIndexer::ReduceResults(const std::vector<Coord> &
             for (int k = 0; k < 3; k++) {
                 if (i + j + k + 2 >= results.size())
                     break;
-
                 Coord A = results[i];
                 Coord B = results[(i + j + 1)];
                 Coord C = results[(i + j + 1) + k + 1];
@@ -133,7 +132,7 @@ std::vector<Coord> FFTIndexer::FilterFFTResults() const {
          it != fft_result_map.rend() && count < max_vectors;
          ++it, ++count) {
         fft_result_filtered.emplace_back(it->second);
-         }
+    }
 
     std::vector<Coord> ret;
 
@@ -180,6 +179,7 @@ std::vector<Coord> FFTIndexer::FilterFFTResults() const {
                     // If it's less than 25%, the shorter peak is likely noise/aliasing,
                     // and the longer vector is the true primitive cell.
                     if (magnitude_ratio > MIN_FUNDAMENTAL_PEAK_RATIO) {
+                        dir_i = dir_j;
                         len_i = fft_result_filtered[j].length;
                         best_idx = j;
                     }
@@ -189,6 +189,12 @@ std::vector<Coord> FFTIndexer::FilterFFTResults() const {
         Coord best_dir = direction_vectors.at(fft_result_filtered[best_idx].direction);
         ret.push_back(best_dir * fft_result_filtered[best_idx].length);
     }
+
+    // Sort filtered vectors by magnitude
+    std::sort(ret.begin(), ret.end(), [](const Coord &A, const Coord &B) {
+       return A.Length() < B.Length();
+    });
+
     return ret;
 }
 
@@ -204,8 +210,8 @@ std::vector<CrystalLattice> FFTIndexer::RunInternal(const std::vector<Coord> &co
     assert(coord.size() <= FFT_MAX_SPOTS);
 
     ExecuteFFT(coord, nspots);
-    auto f = FilterFFTResults();
-    auto r = ReduceResults(f);
+    const auto f = FilterFFTResults();
+    const auto r = ReduceResults(f);
 
     Eigen::MatrixX3<float> oCell(r.size() * 3u, 3u);
     Eigen::VectorX<float> scores(r.size());
@@ -238,13 +244,5 @@ std::vector<CrystalLattice> FFTIndexer::RunInternal(const std::vector<Coord> &co
         .indexing_tolerance = indexing_tolerance
     };
 
-    auto ref_latt = Refine(coord, nspots, oCell, scores, parameters);
-    if (ref_latt.size() >= 1) {
-        auto uc = ref_latt.at(0).GetUnitCell();
-        if (uc.alpha < min_angle_deg || uc.alpha > max_angle_deg
-            || uc.beta < min_angle_deg || uc.beta > max_angle_deg
-            || uc.gamma < min_angle_deg || uc.gamma > max_angle_deg)
-            return {};
-    }
-    return ref_latt;
+    return Refine(coord, nspots, oCell, scores, parameters);
 }

image_analysis/indexing/PostIndexingRefinement.cpp

@@ -3,6 +3,8 @@
 
 #include "PostIndexingRefinement.h"
 
+#include <iostream>
+
 namespace {
     struct config_ifssr final {
         float threshold_contraction = .8; // contract error threshold by this value in every iteration
@@ -170,17 +172,18 @@ std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
         } else {
             if (row_norms.minCoeff() < p.min_length_A || row_norms.maxCoeff() > p.max_length_A)
                 continue;
-
-            float alpha = std::acos(cell_rows.row(1).normalized().dot(cell_rows.row(2).normalized())) * 180.0f / M_PI;
-            float beta = std::acos(cell_rows.row(0).normalized().dot(cell_rows.row(2).normalized())) * 180.0f / M_PI;
-            float gamma = std::acos(cell_rows.row(0).normalized().dot(cell_rows.row(1).normalized())) * 180.0f / M_PI;
-
-            if (alpha < p.min_angle_deg || alpha > p.max_angle_deg ||
-                beta < p.min_angle_deg || beta > p.max_angle_deg ||
-                gamma < p.min_angle_deg || gamma > p.max_angle_deg)
-                continue;
         }
 
+        // Filter for wrong angles
+        float alpha = std::acos(cell_rows.row(1).normalized().dot(cell_rows.row(2).normalized())) * 180.0f / M_PI;
+        float beta = std::acos(cell_rows.row(0).normalized().dot(cell_rows.row(2).normalized())) * 180.0f / M_PI;
+        float gamma = std::acos(cell_rows.row(0).normalized().dot(cell_rows.row(1).normalized())) * 180.0f / M_PI;
+
+        if (alpha < p.min_angle_deg || alpha > p.max_angle_deg ||
+            beta < p.min_angle_deg || beta > p.max_angle_deg ||
+            gamma < p.min_angle_deg || gamma > p.max_angle_deg)
+            continue;
+
         int64_t indexed_spot_count = 0;
         auto indexed_mask = ComputeIndexedMask(spots.topRows(nspots), cell_cols, p.indexing_tolerance, indexed_spot_count);
 
@@ -223,25 +226,21 @@ std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
     std::vector<RefinedCandidate> accepted;
 
     for (const auto &candidate: candidates) {
-        bool too_similar = false;
+        int64_t overlap = 0;
 
+        // Check all already selected lattices and see how many spots are already indexed for the candidate
+        // If the overlap is more than 40% of indexed spots - we assume the lattice doesn't bring anything new
         for (const auto &selected: accepted) {
-            int64_t overlap = 0;
             for (size_t i = 0; i < candidate.indexed_mask.size(); ++i) {
                 if (candidate.indexed_mask[i] && selected.indexed_mask[i])
                     overlap++;
             }
-
-            const int64_t max_set_size = std::max(candidate.indexed_spot_count, selected.indexed_spot_count);
-            if (overlap > static_cast<int64_t>(REFINE_CANDIDATE_OVERLAP_RATIO_THRESHOLD
-                                               * static_cast<float>(max_set_size))) {
-                too_similar = true;
-                break;
-            }
         }
 
-        if (!too_similar)
+        if (overlap < static_cast<int64_t>(REFINE_CANDIDATE_OVERLAP_RATIO_THRESHOLD
+                                           * static_cast<float>(candidate.indexed_spot_count))) {
             accepted.emplace_back(candidate);
+        }
     }
 
     ret.reserve(accepted.size());

image_analysis/indexing/PostIndexingRefinement.h

@@ -14,7 +14,7 @@
 
 constexpr float REFINE_CANDIDATE_SPOT_COUNT_RATIO_THRESHOLD = 0.9f;
 constexpr float REFINE_CANDIDATE_VOLUME_RATIO_THRESHOLD = 1.05f;
-constexpr float REFINE_CANDIDATE_OVERLAP_RATIO_THRESHOLD = 0.2f;
+constexpr float REFINE_CANDIDATE_OVERLAP_RATIO_THRESHOLD = 0.4f;
 constexpr float REFINE_MIN_VOLUME_EPSILON = 1e-12f;
 constexpr float REFINE_MIN_REFERENCE_LENGTH_EPSILON = 1e-6f;