FFTIndexer: Limit cell angles to 30-150 range to avoid colinear vectors

image_analysis/indexing/EigenRefine.h

@@ -16,6 +16,8 @@ struct RefineParameters {
     std::optional<UnitCell> reference_unit_cell;
     float min_length_A;
     float max_length_A;
+    float min_angle_deg;
+    float max_angle_deg;
     float indexing_tolerance;
 };
 
@@ -112,7 +114,10 @@ inline std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
     int64_t id = -1;
 
     for (int i = 0; i < scores.size(); i++) {
-        Eigen::Vector3f row_norms = oCell.block(3u * i, 0u, 3u, 3u).rowwise().norm();
+        // Get cell vectors
+        auto cell = oCell.block(3u * i, 0u, 3u, 3u);
+
+        Eigen::Vector3f row_norms = cell.rowwise().norm();
 
         // Check for distance vs. reference unit cell
         if (p.reference_unit_cell) {
@@ -138,13 +143,22 @@ inline std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
             }
             if (!lengths_ok) continue;
         } else {
-            // Don't include lattices with weird lengths (only if no reference unit cell)
+            // Check lengths (A, B, C)
             if (row_norms.minCoeff() < p.min_length_A || row_norms.maxCoeff() > p.max_length_A)
                 continue;
+
+            // Calculate angles (alpha, beta, gamma) in degrees
+            float alpha = std::acos(cell.row(1).normalized().dot(cell.row(2).normalized())) * 180.0f / M_PI;
+            float beta = std::acos(cell.row(0).normalized().dot(cell.row(2).normalized())) * 180.0f / M_PI;
+            float gamma = std::acos(cell.row(0).normalized().dot(cell.row(1).normalized())) * 180.0f / M_PI;
+
+            // Check if angles are within allowed range
+            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;
         }
 
-        // Count spots that match indexing tolerance
-        auto cell = oCell.block(3u * i, 0u, 3u, 3u);
         M3x resid = spots.topRows(nspots) * cell.transpose();
         const M3x miller = round(resid.array());
         resid -= miller;

image_analysis/indexing/FFBIDXIndexer.cpp

@@ -40,12 +40,14 @@ std::vector<CrystalLattice> FFBIDXIndexer::Run(const std::vector<Coord> &coord,
     indexer.index(1, nspots);
 
     RefineParameters parameters{
-            .viable_cell_min_spots = viable_cell_min_spots,
-            .dist_tolerance_vs_reference = dist_tolerance_vs_reference,
-            .reference_unit_cell = reference_unit_cell,
-            .min_length_A = 1, // doesn't matter
-            .max_length_A = 1000, // doesn't matter
-            .indexing_tolerance = indexing_tolerance
+        .viable_cell_min_spots = viable_cell_min_spots,
+        .dist_tolerance_vs_reference = dist_tolerance_vs_reference,
+        .reference_unit_cell = reference_unit_cell,
+        .min_length_A = 1, // doesn't matter
+        .max_length_A = 1000, // doesn't matter
+        .min_angle_deg = 30,
+        .max_angle_deg = 150,
+        .indexing_tolerance = indexing_tolerance
     };
 
     return Refine(coord, nspots, indexer.oCellM(), indexer.oScoreV(), parameters);

image_analysis/indexing/FFTIndexer.cpp

@@ -6,10 +6,12 @@
 #include "EigenRefine.h"
 
 FFTIndexer::FFTIndexer(const IndexingSettings &settings)
-      : max_length_A(settings.GetFFT_MaxUnitCell_A()),
-        min_length_A(settings.GetFFT_MinUnitCell_A()),
-        nDirections(settings.GetFFT_NumVectors()),
-result_fft(nDirections) {
+    : max_length_A(settings.GetFFT_MaxUnitCell_A()),
+      min_length_A(settings.GetFFT_MinUnitCell_A()),
+      min_angle_deg(30),
+      max_angle_deg(150),
+      nDirections(settings.GetFFT_NumVectors()),
+      result_fft(nDirections) {
 
     float maxQ = 2.0f * static_cast<float>(M_PI) / settings.GetFFT_HighResolution_A();
 
@@ -100,8 +102,17 @@ std::vector<CrystalLattice> FFTIndexer::ReduceResults(const std::vector<Coord> &
                 C = C - std::round(C * A / (A * A)) * A;
                 C = C - std::round(C * B / (B * B)) * B;
 
+                float alpha = angle_deg(B, C);
+                float beta = angle_deg(A, C);
+                float gamma = angle_deg(A, B);
+
                 // Check if values are OK after lattice reduction
-                if (A.Length() >= min_length_A && B.Length() >= min_length_A && C.Length() >= min_length_A)
+                if (A.Length() >= min_length_A
+                    && B.Length() >= min_length_A
+                    && C.Length() >= min_length_A
+                    && alpha >= min_angle_deg && alpha <= max_angle_deg
+                    && beta >= min_angle_deg && beta <= max_angle_deg
+                    && gamma >= min_angle_deg && gamma <= max_angle_deg)
                     candidates.emplace_back(A, B, C);
             }
         }
@@ -188,6 +199,8 @@ std::vector<CrystalLattice> FFTIndexer::Run(const std::vector<Coord> &coord, siz
         .reference_unit_cell = reference_unit_cell,
         .min_length_A = min_length_A,
         .max_length_A = max_length_A,
+        .min_angle_deg = min_angle_deg,
+        .max_angle_deg = max_angle_deg,
         .indexing_tolerance = indexing_tolerance
     };
 

image_analysis/indexing/FFTIndexer.h

@@ -18,6 +18,9 @@ class FFTIndexer : public Indexer {
 protected:
     const float min_length_A;
     const float max_length_A;
+    const float min_angle_deg;
+    const float max_angle_deg;
+
     const int nDirections;
     float histogram_spacing;
     int64_t histogram_size;