v1.0.0-rc.122 (#29)

This is an UNSTABLE release.

* jfjoch_broker: Add thresholding to prefer shorter vectors after FFT
* jfjoch_broker: Add experimental mosaicity estimation for rotation experiments (this is work in progress)
* jfjoch_viewer: Display file opening errors
* jfjoch_viewer: When loading files over DBus add retry/back-off till the file is available

Reviewed-on: #29
Co-authored-by: Filip Leonarski <filip.leonarski@psi.ch>
Co-committed-by: Filip Leonarski <filip.leonarski@psi.ch>

image_analysis/IndexAndRefine.cpp

@@ -129,7 +129,7 @@ void IndexAndRefine::ProcessImage(DataMessage &msg,
             msg.beam_corr_y = data.beam_corr_y;
         }
 
-        if (AnalyzeIndexing(msg, experiment_copy, *lattice_candidate)) {
+        if (AnalyzeIndexing(msg, experiment_copy, *lattice_candidate, experiment_copy.GetGoniometer())) {
             msg.lattice_type = symmetry;
 
             float ewald_dist_cutoff = 0.001f;

image_analysis/bragg_integration/BraggPrediction.cpp

@@ -2,10 +2,7 @@
 // SPDX-License-Identifier: GPL-3.0-only
 
 #include "BraggPrediction.h"
-
-inline bool odd(const int val) {
-    return (val & 1) != 0;
-}
+#include "SystematicAbsence.h"
 
 BraggPrediction::BraggPrediction(int max_reflections)
     : max_reflections(max_reflections), reflections(max_reflections) {}
@@ -55,52 +52,7 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
             const float AhBk_z = Ah_z + Bstar.z * k;
 
             for (int l = -settings.max_hkl; l < settings.max_hkl; l++) {
-                bool absent = false;
-                // (000) is not too interesting
-                if (h == 0 && k == 0 && l == 0)
-                    absent = true;
-
-                switch (settings.centering) {
-                    case 'I':
-                        // Body-centered: h + k + l must be even
-                        if (odd(h + k + l))
-                            absent = true;
-                        break;
-                    case 'A':
-                        // A-centered: k + l must be even
-                        if (odd(k + l))
-                            absent = true;
-                        break;
-                    case 'B':
-                        // B-centered: h + l must be even
-                        if (odd(h + l))
-                            absent = true;
-                        break;
-                    case 'C':
-                        // C-centered: h + k must be even
-                        if (odd(h + k))
-                            absent = true;
-                        break;
-                    case 'F': { // Face-centered: h, k, l all even or all odd
-                        if (odd(h+k) || odd(h+l) || odd(k+l))
-                            absent = true;
-                        break;
-                    }
-                    case 'R': {
-                        // Rhombohedral in hexagonal setting (hR, a_h=b_h, γ=120°):
-                        // Reflection condition: -h + k + l = 3n  (equivalently h - k + l = 3n)
-                        int mod = (-h + k + l) % 3;
-                        if (mod < 0) mod += 3;
-                        if (mod != 0)
-                            absent = true;
-                        break;
-                    }
-                    default:
-                        // P or unspecified: no systematic absences
-                        break;
-                }
-
-                if (absent)
+                if (systematic_absence(h, k, l, settings.centering))
                     continue;
 
                 if (i >= max_reflections)

image_analysis/bragg_integration/BraggPredictionRotation.cpp

@@ -0,0 +1,209 @@
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "BraggPredictionRotation.h"
+
+#include <cmath>
+#include <algorithm>
+
+#include "../../common/DiffractionGeometry.h"
+#include "../../common/GoniometerAxis.h"
+#include "../../common/JFJochException.h"
+#include "SystematicAbsence.h"
+
+namespace {
+    inline float deg_to_rad(float deg) {
+        return deg * (static_cast<float>(M_PI) / 180.0f);
+    }
+    inline float rad_to_deg(float rad) {
+        return rad * (180.0f / static_cast<float>(M_PI));
+    }
+
+    // Solve A cos(phi) + B sin(phi) + D = 0, return solutions in [phi0, phi1]
+    // Returns 0..2 solutions.
+    inline int solve_trig(float A, float B, float D,
+                          float phi0, float phi1,
+                          float out_phi[2]) {
+        const float R = std::sqrt(A*A + B*B);
+        if (!(R > 0.0f))
+            return 0;
+
+        const float rhs = -D / R;
+        if (rhs < -1.0f || rhs > 1.0f)
+            return 0;
+
+        const float phi_ref = std::atan2(B, A);
+        const float delta = std::acos(rhs);
+
+        float s1 = phi_ref + delta;
+        float s2 = phi_ref - delta;
+
+        // Shift candidates by +/- 2*pi so they land near the interval.
+        const float two_pi = 2.0f * static_cast<float>(M_PI);
+        auto shift_near = [&](float x) {
+            while (x < phi0 - two_pi) x += two_pi;
+            while (x > phi1 + two_pi) x -= two_pi;
+            return x;
+        };
+
+        s1 = shift_near(s1);
+        s2 = shift_near(s2);
+
+        int n = 0;
+        if (s1 >= phi0 && s1 <= phi1) out_phi[n++] = s1;
+        if (s2 >= phi0 && s2 <= phi1) {
+            if (n == 0 || std::fabs(s2 - out_phi[0]) > 1e-6f) out_phi[n++] = s2;
+        }
+        return n;
+    }
+
+    // Convert angle to fractional image_number using goniometer start/increment
+    inline float phi_deg_to_image_number(float phi_deg, float start_deg, float increment_deg) {
+        return (phi_deg - start_deg) / increment_deg;
+    }
+} // namespace
+
+std::vector<Reflection> BraggPredictionRotation::Calc(const DiffractionExperiment& experiment,
+                                                        const CrystalLattice& lattice,
+                                                        const BraggPredictionRotationSettings& settings) const {
+    std::vector<Reflection> out;
+    out.reserve(200000); // tune
+
+    const auto gon_opt = experiment.GetGoniometer();
+    if (!gon_opt.has_value())
+        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+                              "BraggPredictionRotationCPU requires a goniometer axis");
+
+    const GoniometerAxis& gon = *gon_opt;
+    const auto geom = experiment.GetDiffractionGeometry();
+
+    // Determine prediction interval in spindle angle.
+    // We treat each image as [angle(image), angle(image)+wedge)
+    const float start_deg = gon.GetStart_deg();
+    const float inc_deg   = gon.GetIncrement_deg();
+    const float wedge_deg = gon.GetWedge_deg();
+
+    float phi0_deg = gon.GetAngle_deg(static_cast<float>(settings.image_first));
+    float phi1_deg = gon.GetAngle_deg(static_cast<float>(settings.image_last)) + wedge_deg;
+
+    if (settings.pad_one_wedge) {
+        phi0_deg -= wedge_deg;
+        phi1_deg += wedge_deg;
+    }
+
+    const float phi0 = deg_to_rad(phi0_deg);
+    const float phi1 = deg_to_rad(phi1_deg);
+
+    // Resolution cutoff in reciprocal space
+    const float one_over_dmax = 1.0f / settings.high_res_A;
+    const float one_over_dmax_sq = one_over_dmax * one_over_dmax;
+
+    // S0 has length 1/lambda (your convention)
+    const Coord S0 = geom.GetScatteringVector();
+    const float k2 = (S0 * S0); // (1/lambda)^2
+
+    // Rotation axis in lab frame (already normalized in ctor, but keep safe)
+    const Coord w = gon.GetAxis().Normalize();
+
+    const Coord Astar = lattice.Astar();
+    const Coord Bstar = lattice.Bstar();
+    const Coord Cstar = lattice.Cstar();
+
+    const float det_w = static_cast<float>(experiment.GetXPixelsNum());
+    const float det_h = static_cast<float>(experiment.GetYPixelsNum());
+
+    for (int h = -settings.max_hkl; h <= settings.max_hkl; ++h) {
+        const Coord Ah = Astar * static_cast<float>(h);
+
+        for (int k = -settings.max_hkl; k <= settings.max_hkl; ++k) {
+            const Coord AhBk = Ah + Bstar * static_cast<float>(k);
+
+            for (int l = -settings.max_hkl; l <= settings.max_hkl; ++l) {
+                if (systematic_absence(h, k, l, settings.centering))
+                    continue;
+
+                const Coord g0 = AhBk + Cstar * static_cast<float>(l);
+                const float g02 = g0 * g0;
+
+                if (!(g02 > 0.0f))
+                    continue;
+                if (g02 > one_over_dmax_sq)
+                    continue;
+
+                // g0 = g_par + g_perp wrt spindle axis
+                const float g_par_s = g0 * w;
+                const Coord g_par = w * g_par_s;
+                const Coord g_perp = g0 - g_par;
+
+                const float g_perp2 = g_perp * g_perp;
+                if (g_perp2 < 1e-12f) {
+                    // Rotation does not move this reciprocal vector: skip for now.
+                    // (Could be handled by checking whether it satisfies Ewald at all.)
+                    continue;
+                }
+
+                // Equation: |S0 + g(phi)|^2 = |S0|^2
+                // g(phi) = g_par + cos(phi) g_perp + sin(phi) (w x g_perp)
+                const Coord p = S0 + g_par;
+                const Coord w_x_gperp = w % g_perp;
+
+                const float A = 2.0f * (p * g_perp);
+                const float B = 2.0f * (p * w_x_gperp);
+                const float D = (p * p) + g_perp2 - k2;
+
+                float sols[2]{};
+                const int nsol = solve_trig(A, B, D, phi0, phi1, sols);
+                if (nsol == 0)
+                    continue;
+
+                for (int si = 0; si < nsol; ++si) {
+                    const float phi = sols[si];
+
+                    // Rotate g0 by phi about w
+                    const RotMatrix R(phi, w);
+                    const Coord g = R * g0;
+                    const Coord S = S0 + g;
+
+                    // Project to detector using canonical geometry code
+                    const auto [x, y] = geom.RecipToDector(g);
+                    if (!std::isfinite(x) || !std::isfinite(y))
+                        continue;
+                    if (x < 0.0f || x >= det_w || y < 0.0f || y >= det_h)
+                        continue;
+
+                    // Convert phi to fractional image number
+                    const float phi_deg = rad_to_deg(phi);
+
+                    Reflection r{};
+                    r.h = h;
+                    r.k = k;
+                    r.l = l;
+                    r.angle_deg = phi_deg;
+                    r.image_number = phi_deg_to_image_number(phi_deg, start_deg, inc_deg);
+                    r.predicted_x = x;
+                    r.predicted_y = y;
+                    r.d = 1.0f / std::sqrt(g02);
+
+                    // diagnostic: should be ~0
+                    r.dist_ewald = S.Length() - std::sqrt(k2);
+
+                    r.I = 0.0f;
+                    r.bkg = 0.0f;
+                    r.sigma = 0.0f;
+
+                    out.push_back(r);
+                }
+            }
+        }
+    }
+
+    std::sort(out.begin(), out.end(), [](const Reflection &a, const Reflection &b) {
+        if (a.angle_deg != b.angle_deg)
+            return a.angle_deg < b.angle_deg;
+        if (a.h != b.h) return a.h < b.h;
+        if (a.k != b.k) return a.k < b.k;
+        return a.l < b.l;
+    });
+
+    return out;
+}

image_analysis/bragg_integration/BraggPredictionRotation.h

@@ -0,0 +1,38 @@
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#ifndef JFJOCH_BRAGGPREDICTIONROTATION_H
+#define JFJOCH_BRAGGPREDICTIONROTATION_H
+
+#include <vector>
+#include <cstdint>
+
+#include "../../common/CrystalLattice.h"
+#include "../../common/DiffractionExperiment.h"
+#include "../../common/Reflection.h"
+
+struct BraggPredictionRotationSettings {
+    float high_res_A = 1.5f;
+    int max_hkl = 100;
+    char centering = 'P';
+
+    // Predict only in this image interval (inclusive)
+    int64_t image_first = 0;
+    int64_t image_last = 0;
+
+    // If true, expands [phi_start, phi_end] by one wedge on each side
+    // (useful if you later want to model partials / edge effects).
+    bool pad_one_wedge = false;
+};
+
+class BraggPredictionRotation {
+public:
+    BraggPredictionRotation() = default;
+
+    std::vector<Reflection> Calc(const DiffractionExperiment& experiment,
+                                 const CrystalLattice& lattice,
+                                 const BraggPredictionRotationSettings& settings) const;
+};
+
+
+#endif //JFJOCH_BRAGGPREDICTIONROTATION_H

image_analysis/bragg_integration/CMakeLists.txt

@@ -8,6 +8,9 @@ ADD_LIBRARY(JFJochBraggIntegration STATIC
         CalcISigma.h
         BraggPredictionFactory.cpp
         BraggPredictionFactory.h
+        SystematicAbsence.h
+        BraggPredictionRotation.cpp
+        BraggPredictionRotation.h
 )
 
 TARGET_LINK_LIBRARIES(JFJochBraggIntegration JFJochCommon)

image_analysis/bragg_integration/SystematicAbsence.h

@@ -0,0 +1,26 @@
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#ifndef JFJOCH_SYSTEMATICABSENCE_H
+#define JFJOCH_SYSTEMATICABSENCE_H
+
+static inline bool odd_int(int v) { return (v & 1) != 0; }
+
+static inline bool systematic_absence(int h, int k, int l, char centering) {
+    if (h == 0 && k == 0 && l == 0) return true;
+    switch (centering) {
+        case 'I': return odd_int(h + k + l);
+        case 'A': return odd_int(k + l);
+        case 'B': return odd_int(h + l);
+        case 'C': return odd_int(h + k);
+        case 'F': return (odd_int(h + k) || odd_int(h + l) || odd_int(k + l));
+        case 'R': {
+            int mod = (-h + k + l) % 3;
+            if (mod < 0) mod += 3;
+            return mod != 0;
+        }
+        default:  return false; // P
+    }
+}
+
+#endif //JFJOCH_SYSTEMATICABSENCE_H

image_analysis/indexing/AnalyzeIndexing.cpp

@@ -7,23 +7,126 @@
 
 #include "FitProfileRadius.h"
 
-inline bool ok(float x) {
-    if (!std::isfinite(x))
-        return false;
-    if (x < 0.0)
-        return false;
-    return true;
-};
+namespace {
+    inline bool ok(float x) {
+        if (!std::isfinite(x))
+            return false;
+        if (x < 0.0)
+            return false;
+        return true;
+    }
 
-bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experiment, const CrystalLattice &latt) {
+    inline float deg_to_rad(float deg) {
+        return deg * (static_cast<float>(M_PI) / 180.0f);
+    }
 
+    inline float rad_to_deg(float rad) {
+        return rad * (180.0f / static_cast<float>(M_PI));
+    }
+
+    // Wrap to [-180, 180] (useful for residuals)
+    inline float wrap_deg_pm180(float deg) {
+        while (deg > 180.0f) deg -= 360.0f;
+        while (deg < -180.0f) deg += 360.0f;
+        return deg;
+    }
+
+    // Solve A cos(phi) + B sin(phi) + D = 0, return solutions in [phi0, phi1] (radians)
+    inline int solve_trig(float A, float B, float D,
+                          float phi0, float phi1,
+                          float out_phi[2]) {
+        const float R = std::sqrt(A * A + B * B);
+        if (!(R > 0.0f))
+            return 0;
+
+        const float rhs = -D / R;
+        if (rhs < -1.0f || rhs > 1.0f)
+            return 0;
+
+        const float phi_ref = std::atan2(B, A);
+        const float delta = std::acos(rhs);
+
+        float s1 = phi_ref + delta;
+        float s2 = phi_ref - delta;
+
+        const float two_pi = 2.0f * static_cast<float>(M_PI);
+        auto shift_near = [&](float x) {
+            while (x < phi0 - two_pi) x += two_pi;
+            while (x > phi1 + two_pi) x -= two_pi;
+            return x;
+        };
+
+        s1 = shift_near(s1);
+        s2 = shift_near(s2);
+
+        int n = 0;
+        if (s1 >= phi0 && s1 <= phi1) out_phi[n++] = s1;
+        if (s2 >= phi0 && s2 <= phi1) {
+            if (n == 0 || std::fabs(s2 - out_phi[0]) > 1e-6f) out_phi[n++] = s2;
+        }
+        return n;
+    }
+
+    // Find predicted phi (deg) for given g0 around phi_obs (deg) within +/- half_window_deg.
+    // Returns nullopt if no solution in the local window.
+    inline std::optional<float> predict_phi_deg_local(const Coord &g0,
+                                                      const Coord &S0,
+                                                      const Coord &w_unit,
+                                                      float phi_obs_deg,
+                                                      float half_window_deg) {
+        const float phi0 = deg_to_rad(phi_obs_deg - half_window_deg);
+        const float phi1 = deg_to_rad(phi_obs_deg + half_window_deg);
+
+        // Decompose g0 into parallel/perp to w
+        const float g_par_s = g0 * w_unit;
+        const Coord g_par = w_unit * g_par_s;
+        const Coord g_perp = g0 - g_par;
+
+        const float g_perp2 = g_perp * g_perp;
+        if (g_perp2 < 1e-12f)
+            return std::nullopt;
+
+        const float k2 = (S0 * S0); // |S0|^2 = (1/lambda)^2
+
+        // Equation: |S0 + g(phi)|^2 = |S0|^2
+        const Coord p = S0 + g_par;
+        const Coord w_x_gperp = w_unit % g_perp;
+
+        const float A = 2.0f * (p * g_perp);
+        const float B = 2.0f * (p * w_x_gperp);
+        const float D = (p * p) + g_perp2 - k2;
+
+        float sols[2]{};
+        const int nsol = solve_trig(A, B, D, phi0, phi1, sols);
+        if (nsol == 0)
+            return std::nullopt;
+
+        // Pick the solution closest to phi_obs
+        const float phi_obs = deg_to_rad(phi_obs_deg);
+        float best_phi = sols[0];
+        float best_err = std::fabs(sols[0] - phi_obs);
+
+        if (nsol == 2) {
+            const float err2 = std::fabs(sols[1] - phi_obs);
+            if (err2 < best_err) {
+                best_err = err2;
+                best_phi = sols[1];
+            }
+        }
+
+        return rad_to_deg(best_phi);
+    }
+} // namespace
+
+bool AnalyzeIndexing(DataMessage &message,
+                     const DiffractionExperiment &experiment,
+                     const CrystalLattice &latt,
+                     const std::optional<GoniometerAxis> &rotation_axis) {
     size_t nspots = message.spots.size();
 
     uint64_t indexed_spot_count = 0;
     std::vector<uint8_t> indexed_spots(nspots);
 
-    std::vector<float> distance_ewald_sphere(nspots);
-
     // Check spots
     const Coord a = latt.Vec0();
     const Coord b = latt.Vec1();
@@ -36,7 +139,8 @@ bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experime
     const auto geom = experiment.GetDiffractionGeometry();
     const auto indexing_tolerance = experiment.GetIndexingSettings().GetTolerance();
     const auto viable_cell_min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots();
-     // identify indexed spots
+
+    // identify indexed spots
     for (int i = 0; i < message.spots.size(); i++) {
         auto recip = message.spots[i].ReciprocalCoord(geom);
 
@@ -79,8 +183,54 @@ bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experime
         message.spot_count_indexed = indexed_spot_count;
         message.indexing_lattice = latt;
         message.indexing_unit_cell = latt.GetUnitCell();
+
+        message.mosaicity_deg = std::nullopt;
+        if (rotation_axis.has_value()) {
+            const auto gon_opt = experiment.GetGoniometer();
+            if (gon_opt.has_value()) {
+                const auto &gon = *gon_opt;
+
+                const Coord w = rotation_axis->GetAxis().Normalize();
+                const Coord S0 = geom.GetScatteringVector();
+                const float wedge_deg = gon.GetWedge_deg();
+                const float start_deg = gon.GetStart_deg();
+                const float inc_deg = gon.GetIncrement_deg();
+
+                double sum_sq = 0.0;
+                int count = 0;
+
+                for (const auto &s: message.spots) {
+                    if (!s.indexed)
+                        continue;
+
+                    // Observed angle: use frame center
+                    const float image_center = static_cast<float>(s.image) + 0.5f;
+                    const float phi_obs_deg = start_deg + inc_deg * image_center;
+
+                    // g0 at phi=0 assumption
+                    const Coord g0 = astar * static_cast<float>(s.h)
+                                     + bstar * static_cast<float>(s.k)
+                                     + cstar * static_cast<float>(s.l);
+
+                    // Local solve window: +/- 1 wedge (easy/robust first try)
+                    const auto phi_pred_deg_opt = predict_phi_deg_local(g0, S0, w, phi_obs_deg, wedge_deg);
+                    if (!phi_pred_deg_opt.has_value())
+                        continue;
+
+                    float dphi = wrap_deg_pm180(phi_obs_deg - phi_pred_deg_opt.value());
+                    sum_sq += static_cast<double>(dphi) * static_cast<double>(dphi);
+                    count++;
+                }
+
+                if (count > 0) {
+                    message.mosaicity_deg = static_cast<float>(std::sqrt(sum_sq / static_cast<double>(count)));
+                }
+            }
+        }
+
         return true;
     }
+
     message.indexing_result = false;
     return false;
 }

image_analysis/indexing/AnalyzeIndexing.h

@@ -10,7 +10,10 @@
 
 constexpr static float min_percentage_spots = 0.20f;
 
-bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experiment, const CrystalLattice &latt);
+bool AnalyzeIndexing(DataMessage &message,
+                     const DiffractionExperiment &experiment,
+                     const CrystalLattice &latt,
+                     const std::optional<GoniometerAxis> &rotation_axis);
 
 
 #endif //JFJOCH_ANALYZEINDEXING_H

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;
 }