diff --git a/image_analysis/indexing/FFTIndexer.cpp b/image_analysis/indexing/FFTIndexer.cpp index 273ff977..87b33ff8 100644 --- a/image_analysis/indexing/FFTIndexer.cpp +++ b/image_analysis/indexing/FFTIndexer.cpp @@ -135,6 +135,7 @@ std::vector FFTIndexer::ReduceResults(const std::vector & // most triples of one lattice Niggli-reduce to the same cell (keeps the refine set small). const size_t n = std::min(results.size(), 64); const size_t n_short = std::min(n, 12); + std::vector candidate_cells; // parallel to `candidates`, see the dedup scan below for (size_t i = 0; i < n_short; i++) { for (size_t j = i + 1; j < n_short; j++) { for (size_t k = j + 1; k < n; k++) { @@ -151,10 +152,15 @@ std::vector FFTIndexer::ReduceResults(const std::vector & continue; bool duplicate = false; - for (const auto &c : candidates) - if (c.GetUnitCell().is_close(uc, 0.02f, 1.0f)) { duplicate = true; break; } - if (!duplicate) + for (const auto &c : candidate_cells) + if (c.is_close(uc, 0.02f, 1.0f)) { duplicate = true; break; } + if (!duplicate) { + // Keep each accepted candidate's cell rather than re-deriving it on the next + // triple: GetUnitCell costs three acos, and this scan runs over every candidate + // accepted so far, for every triple. + candidate_cells.push_back(uc); candidates.emplace_back(std::move(reduced)); + } } } } diff --git a/image_analysis/indexing/PostIndexingRefinement.cpp b/image_analysis/indexing/PostIndexingRefinement.cpp index 8725c75f..e65f6486 100644 --- a/image_analysis/indexing/PostIndexingRefinement.cpp +++ b/image_analysis/indexing/PostIndexingRefinement.cpp @@ -149,6 +149,10 @@ std::vector Refine(const std::vector &in_spots, std::vector candidates; + // Angle bounds as cosines, once, for the per-candidate test below. + const float cos_min_angle = std::cos(p.min_angle_deg * PI / 180.0f); + const float cos_max_angle = std::cos(p.max_angle_deg * PI / 180.0f); + for (int i = 0; i < scores.size(); i++) { Eigen::Matrix3f cell_rows = oCell.block(3u * i, 0u, 3u, 3u); Eigen::Matrix3f cell_cols = cell_rows.transpose(); @@ -208,14 +212,18 @@ std::vector Refine(const std::vector &in_spots, continue; } - // Filter for wrong angles - float alpha = std::acos(cell_rows.row(1).normalized().dot(cell_rows.row(2).normalized())) * 180.0f / PI; - float beta = std::acos(cell_rows.row(0).normalized().dot(cell_rows.row(2).normalized())) * 180.0f / PI; - float gamma = std::acos(cell_rows.row(0).normalized().dot(cell_rows.row(1).normalized())) * 180.0f / PI; + // Filter for wrong angles. Compared as COSINES, not angles: acos is strictly decreasing on + // [-1, 1], so "angle outside [min_angle, max_angle]" is exactly "cosine outside + // [cos(max_angle), cos(min_angle)]" with the ends swapped - and the three acos calls the + // comparison needed disappear. They were not cheap: this runs per candidate cell per image, + // and on a serial-stills run acos was 41% of the whole process. + const float cos_alpha = cell_rows.row(1).normalized().dot(cell_rows.row(2).normalized()); + const float cos_beta = cell_rows.row(0).normalized().dot(cell_rows.row(2).normalized()); + const float cos_gamma = cell_rows.row(0).normalized().dot(cell_rows.row(1).normalized()); - 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) + if (cos_alpha > cos_min_angle || cos_alpha < cos_max_angle || + cos_beta > cos_min_angle || cos_beta < cos_max_angle || + cos_gamma > cos_min_angle || cos_gamma < cos_max_angle) continue; int64_t indexed_spot_count = 0;