FFTIndexer: Use Niggli Reduction

image_analysis/indexing/FFTIndexer.cpp

@@ -97,23 +97,20 @@ std::vector<CrystalLattice> FFTIndexer::ReduceResults(const std::vector<Coord> &
                 // sort vectors by length for reduction
                 Sort(A,B,C);
 
-                // Lattice reduction
-                B = B - std::round(B * A / (A * A)) * A;
-                C = C - std::round(C * A / (A * A)) * A;
-                C = C - std::round(C * B / (B * B)) * B;
+                CrystalLattice raw(A, B, C);
+                CrystalLattice reduced = raw.NiggliReduce(); // Reduce cell
+                const auto uc = reduced.GetUnitCell();
 
-                float alpha = angle_deg(B, C);
-                float beta = angle_deg(A, C);
-                float gamma = angle_deg(A, B);
+                if (uc.a < min_length_A || uc.b < min_length_A || uc.c < min_length_A)
+                    continue;
 
-                // Check if values are OK after lattice reduction
-                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);
+                float alpha = uc.alpha, beta = uc.beta, gamma = uc.gamma;
+                if (alpha < min_angle_deg || alpha > max_angle_deg ||
+                    beta  < min_angle_deg || beta  > max_angle_deg ||
+                    gamma < min_angle_deg || gamma > max_angle_deg)
+                    continue;
+
+                candidates.emplace_back(std::move(reduced));
             }
         }
     }
@@ -141,8 +138,7 @@ std::vector<Coord> FFTIndexer::FilterFFTResults() const {
     std::vector<Coord> ret;
 
     // Remove vectors less than 5 deg apart, as most likely these are colinear
-    constexpr float COS_5_DEG = 0.996194;
-
+    const float COS_5_DEG = std::cos(5.0f * M_PI / 180.0f);
 
     // Minimum relative amplitude to accept a shorter vector (fundamental)
     // over a longer one (harmonic).