From 8056d48c2666b5015c25ebe8be56d66b38e3d2cd Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Tue, 16 Dec 2025 12:06:03 +0100
Subject: [PATCH] jfjoch_broker: Add thresholding to prefer shorter vectors
 after FFT

---
 docs/CHANGELOG.md                      |  1 +
 image_analysis/indexing/FFTIndexer.cpp | 42 ++++++++++++++++++++++++--
 2 files changed, 41 insertions(+), 2 deletions(-)

diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md
index 08de2078..9d2cfcec 100644
--- a/docs/CHANGELOG.md
+++ b/docs/CHANGELOG.md
@@ -3,6 +3,7 @@
 ### 1.0.0-rc.122
 This is an UNSTABLE release.
 
+* jfjoch_broker: Add thresholding to prefer shorter vectors after FFT
 * jfjoch_viewer: Display file opening errors
 * jfjoch_viewer: When loading files over DBus add retry/back-off till the file is available
 
diff --git a/image_analysis/indexing/FFTIndexer.cpp b/image_analysis/indexing/FFTIndexer.cpp
index 0465c3f5..925f5c2a 100644
--- a/image_analysis/indexing/FFTIndexer.cpp
+++ b/image_analysis/indexing/FFTIndexer.cpp
@@ -142,18 +142,56 @@ std::vector<Coord> FFTIndexer::FilterFFTResults() const {
 
     // Remove vectors less than 5 deg apart, as most likely these are colinear
     constexpr float COS_5_DEG = 0.996194;
+
+
+    // Minimum relative amplitude to accept a shorter vector (fundamental)
+    // over a longer one (harmonic).
+    // 0.25 means the fundamental frequency must have at least 25% of the
+    // intensity of the harmonic. If less, the odd-indexed spots are likely
+    // systematic absences or noise, and the longer vector is the true cell.
+    constexpr float MIN_FUNDAMENTAL_PEAK_RATIO = 0.25f;
+
     std::vector<bool> ignore(fft_result_filtered.size(), false);
 
     for (int i = 0; i < fft_result_filtered.size(); i++) {
         if (ignore[i])
             continue;
+
         Coord dir_i = direction_vectors.at(fft_result_filtered[i].direction);
+        float len_i = fft_result_filtered[i].length;
+        int best_idx = i; // Index of the vector we currently plan to keep
+
         for (int j = i + 1; j < fft_result_filtered.size(); j++) {
+            if (ignore[j]) continue;
+
             Coord dir_j = direction_vectors.at(fft_result_filtered[j].direction);
-            if (std::fabs(dir_i * dir_j) > COS_5_DEG)
+
+            // If vectors are colinear (angle < 5 deg)
+            if (std::fabs(dir_i * dir_j) > COS_5_DEG) {
                 ignore[j] = true;
+
+                // CHECK: Is the new candidate (j) shorter than current best (best_idx)?
+                // We prefer shorter vectors (fundamental periodicity) over longer ones (harmonics)
+                // BUT only if the shorter vector has "enough" amplitude to be real.
+
+                if (fft_result_filtered[j].length < len_i * 0.9f) {
+                    // Compare against 'i' (the strongest in the cluster) to define the noise floor.
+                    // Using 'best_idx' could allow stepping down into noise if we already swapped to a weak peak.
+                    float magnitude_ratio = fft_result_filtered[j].magnitude / fft_result_filtered[i].magnitude;
+
+                    // Heuristic: If the shorter vector has at least 25% of the amplitude of the
+                    // stronger (longer) vector, assume the shorter one is the true unit cell.
+                    // 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) {
+                        len_i = fft_result_filtered[j].length;
+                        best_idx = j;
+                    }
+                }
+            }
         }
-        ret.push_back(dir_i * fft_result_filtered[i].length);
+        Coord best_dir = direction_vectors.at(fft_result_filtered[best_idx].direction);
+        ret.push_back(best_dir * fft_result_filtered[best_idx].length);
     }
     return ret;
 }