From 47a066d631582243ed2b374ab6d8d9b60609fdc5 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Thu, 4 Jun 2026 18:44:36 +0200 Subject: [PATCH] FFTIndexer: Fixes to handle better multiple lattices --- common/Coord.cpp | 6 --- common/Coord.h | 4 ++ image_analysis/indexing/FFTIndexer.cpp | 24 ++++++------ .../indexing/PostIndexingRefinement.cpp | 37 +++++++++---------- .../indexing/PostIndexingRefinement.h | 2 +- 5 files changed, 34 insertions(+), 39 deletions(-) diff --git a/common/Coord.cpp b/common/Coord.cpp index cde960ea..54809077 100644 --- a/common/Coord.cpp +++ b/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++) diff --git a/common/Coord.h b/common/Coord.h index 96d0a320..33191d15 100644 --- a/common/Coord.h +++ b/common/Coord.h @@ -59,4 +59,8 @@ public: std::vector arr() const; }; +inline void swap(Coord& a, Coord& b) noexcept { + a.swap(b); +} + #endif //INDEX_COORD_H diff --git a/image_analysis/indexing/FFTIndexer.cpp b/image_analysis/indexing/FFTIndexer.cpp index 978ee639..a0f224d3 100644 --- a/image_analysis/indexing/FFTIndexer.cpp +++ b/image_analysis/indexing/FFTIndexer.cpp @@ -89,7 +89,6 @@ std::vector FFTIndexer::ReduceResults(const std::vector & 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 FFTIndexer::FilterFFTResults() const { it != fft_result_map.rend() && count < max_vectors; ++it, ++count) { fft_result_filtered.emplace_back(it->second); - } + } std::vector ret; @@ -180,6 +179,7 @@ std::vector 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 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 FFTIndexer::RunInternal(const std::vector &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 oCell(r.size() * 3u, 3u); Eigen::VectorX scores(r.size()); @@ -238,13 +244,5 @@ std::vector FFTIndexer::RunInternal(const std::vector &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); } \ No newline at end of file diff --git a/image_analysis/indexing/PostIndexingRefinement.cpp b/image_analysis/indexing/PostIndexingRefinement.cpp index 36b0b837..45960eca 100644 --- a/image_analysis/indexing/PostIndexingRefinement.cpp +++ b/image_analysis/indexing/PostIndexingRefinement.cpp @@ -3,6 +3,8 @@ #include "PostIndexingRefinement.h" +#include + namespace { struct config_ifssr final { float threshold_contraction = .8; // contract error threshold by this value in every iteration @@ -170,17 +172,18 @@ std::vector Refine(const std::vector &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 Refine(const std::vector &in_spots, std::vector 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(REFINE_CANDIDATE_OVERLAP_RATIO_THRESHOLD - * static_cast(max_set_size))) { - too_similar = true; - break; - } } - if (!too_similar) + if (overlap < static_cast(REFINE_CANDIDATE_OVERLAP_RATIO_THRESHOLD + * static_cast(candidate.indexed_spot_count))) { accepted.emplace_back(candidate); + } } ret.reserve(accepted.size()); diff --git a/image_analysis/indexing/PostIndexingRefinement.h b/image_analysis/indexing/PostIndexingRefinement.h index e69c22f6..df7d0329 100644 --- a/image_analysis/indexing/PostIndexingRefinement.h +++ b/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;