diff --git a/image_analysis/indexing/EigenRefine.h b/image_analysis/indexing/EigenRefine.h index a4511841..4d549be2 100644 --- a/image_analysis/indexing/EigenRefine.h +++ b/image_analysis/indexing/EigenRefine.h @@ -16,6 +16,8 @@ struct RefineParameters { std::optional 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 Refine(const std::vector &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 Refine(const std::vector &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; diff --git a/image_analysis/indexing/FFBIDXIndexer.cpp b/image_analysis/indexing/FFBIDXIndexer.cpp index 47adf527..d09d58c0 100644 --- a/image_analysis/indexing/FFBIDXIndexer.cpp +++ b/image_analysis/indexing/FFBIDXIndexer.cpp @@ -40,12 +40,14 @@ std::vector FFBIDXIndexer::Run(const std::vector &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); diff --git a/image_analysis/indexing/FFTIndexer.cpp b/image_analysis/indexing/FFTIndexer.cpp index 2e5852b6..f3432490 100644 --- a/image_analysis/indexing/FFTIndexer.cpp +++ b/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(M_PI) / settings.GetFFT_HighResolution_A(); @@ -100,8 +102,17 @@ std::vector FFTIndexer::ReduceResults(const std::vector & 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 FFTIndexer::Run(const std::vector &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 }; diff --git a/image_analysis/indexing/FFTIndexer.h b/image_analysis/indexing/FFTIndexer.h index 5cd5a190..0f38ea73 100644 --- a/image_analysis/indexing/FFTIndexer.h +++ b/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;