From a27f5dc482d05bc5361108354e49ff339681bee9 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Thu, 12 Mar 2026 14:16:19 +0100 Subject: [PATCH] PostIndexingRefinement: Update name --- image_analysis/indexing/CMakeLists.txt | 5 +- image_analysis/indexing/EigenRefine.h | 284 ------------------ image_analysis/indexing/FFBIDXIndexer.cpp | 2 +- image_analysis/indexing/FFTIndexer.cpp | 2 +- image_analysis/indexing/FFTIndexerCPU.cpp | 1 - .../indexing/PostIndexingRefinement.cpp | 258 ++++++++++++++++ .../indexing/PostIndexingRefinement.h | 36 +++ 7 files changed, 299 insertions(+), 289 deletions(-) delete mode 100644 image_analysis/indexing/EigenRefine.h create mode 100644 image_analysis/indexing/PostIndexingRefinement.cpp create mode 100644 image_analysis/indexing/PostIndexingRefinement.h diff --git a/image_analysis/indexing/CMakeLists.txt b/image_analysis/indexing/CMakeLists.txt index cba6ed81..fbbb2df7 100644 --- a/image_analysis/indexing/CMakeLists.txt +++ b/image_analysis/indexing/CMakeLists.txt @@ -9,10 +9,11 @@ ADD_LIBRARY(JFJochIndexing STATIC AnalyzeIndexing.h FitProfileRadius.cpp FitProfileRadius.h - EigenRefine.h + PostIndexingRefinement.h FFTResult.h FFTIndexer.cpp - FFTIndexer.h) + FFTIndexer.h + PostIndexingRefinement.cpp) TARGET_LINK_LIBRARIES(JFJochIndexing JFJochCommon) IF (JFJOCH_CUDA_AVAILABLE) diff --git a/image_analysis/indexing/EigenRefine.h b/image_analysis/indexing/EigenRefine.h deleted file mode 100644 index a72f59cf..00000000 --- a/image_analysis/indexing/EigenRefine.h +++ /dev/null @@ -1,284 +0,0 @@ -// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#pragma once - -#include -#include -#include -#include -#include -#include - -#include "../common/CrystalLattice.h" - -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_MIN_VOLUME_EPSILON = 1e-12f; -constexpr float REFINE_MIN_REFERENCE_LENGTH_EPSILON = 1e-6f; - -struct RefineParameters { - int64_t viable_cell_min_spots; - float dist_tolerance_vs_reference; - std::optional reference_unit_cell; - float min_length_A; - float max_length_A; - float min_angle_deg; - float max_angle_deg; - float indexing_tolerance; -}; - - -struct config_ifssr final { - float threshold_contraction=.8; // contract error threshold by this value in every iteration - float max_distance=.00075; // max distance to reciprocal spots for inliers - unsigned min_spots=8; // minimum number of spots to fit against - unsigned max_iter=32; // max number of iterations -}; - -static std::pair score_parts (float score) noexcept -{ - float nsp = -std::floor(score); - float s = score + nsp; - return std::make_pair(nsp-1, s); -} - -struct RefinedCandidate { - Eigen::Matrix3f cell; - float score; - float volume; - int64_t indexed_spot_count; - std::vector indexed_mask; -}; - -static inline Eigen::MatrixX3 CalculateResiduals( - const Eigen::Ref> &spots, - const Eigen::Matrix3f &cell) { - Eigen::MatrixX3 miller = (spots * cell).array().round().matrix(); - Eigen::MatrixX3 resid = miller * cell.inverse(); - resid -= spots; - return resid; -} - -static inline std::vector ComputeIndexedMask( - const Eigen::Ref> &spots, - const Eigen::Matrix3f &cell, - float indexing_tolerance, - int64_t &indexed_spot_count) { - const float indexing_tolerance_sq = indexing_tolerance * indexing_tolerance; - const Eigen::MatrixX3 resid = CalculateResiduals(spots, cell); - - std::vector mask(spots.rows(), 0); - indexed_spot_count = 0; - - for (int i = 0; i < spots.rows(); ++i) { - if (resid.row(i).squaredNorm() < indexing_tolerance_sq) { - mask[i] = 1; - indexed_spot_count++; - } - } - - return mask; -} - -template -static void RefineCandidateCells(const Eigen::Ref> &spots, - Eigen::DenseBase &cells, - Eigen::DenseBase &scores, - const config_ifssr &cifssr, - unsigned block = 0, unsigned nblocks = 1) { - using namespace Eigen; - using Mx3 = MatrixX3; - using M3 = Matrix3; - - const unsigned nspots = spots.rows(); - const unsigned ncells = scores.rows(); - VectorX below{nspots}; - MatrixX3 sel{nspots, 3u}; - Mx3 resid{nspots, 3u}; - Mx3 miller{nspots, 3u}; - M3 cell; - - const unsigned blocksize = (ncells + nblocks - 1u) / nblocks; - const unsigned startcell = block * blocksize; - const unsigned endcell = std::min(startcell + blocksize, ncells); - - for (unsigned j = startcell; j < endcell; j++) { - if (nspots < cifssr.min_spots) { - scores(j) = float{1.}; - continue; - } - - cell = cells.block(3u * j, 0u, 3u, 3u).transpose(); // cell: col vectors - const float scale = cell.colwise().norm().minCoeff(); - float threshold = score_parts(scores[j]).second / scale; - - for (unsigned niter = 1; niter < cifssr.max_iter && threshold > cifssr.max_distance; niter++) { - miller = (spots * cell).array().round().matrix(); - resid = miller * cell.inverse(); - resid -= spots; - - below = (resid.rowwise().norm().array() < threshold); - if (below.count() < cifssr.min_spots) - break; - - threshold *= cifssr.threshold_contraction; - sel.colwise() = below; - HouseholderQR qr{sel.select(spots, .0f)}; - cell = qr.solve(sel.select(miller, .0f)); - } - - resid = CalculateResiduals(spots, cell); - - ArrayX dist = resid.rowwise().norm(); - auto nth = std::begin(dist) + (cifssr.min_spots - 1); - std::nth_element(std::begin(dist), nth, std::end(dist)); - scores(j) = *nth; - - cells.block(3u * j, 0u, 3u, 3u) = cell.transpose(); - } -} - -inline std::vector Refine(const std::vector &in_spots, - size_t nspots, - Eigen::MatrixX3 &oCell, - Eigen::VectorX &scores, - RefineParameters &p) { - std::vector ret; - - Eigen::MatrixX3 spots(in_spots.size(), 3u); - - for (int i = 0; i < in_spots.size(); i++) { - spots(i, 0u) = in_spots[i].x; - spots(i, 1u) = in_spots[i].y; - spots(i, 2u) = in_spots[i].z; - } - - config_ifssr cifssr{ - .min_spots = static_cast(p.viable_cell_min_spots) - }; - - RefineCandidateCells(spots.topRows(nspots), oCell, scores, cifssr); - - std::vector candidates; - - for (int i = 0; i < scores.size(); i++) { - auto cell_block = oCell.block(3u * i, 0u, 3u, 3u); - Eigen::Matrix3f cell = cell_block.transpose(); - Eigen::Vector3f row_norms = cell_block.rowwise().norm(); - - if (p.reference_unit_cell) { - std::array obs = {row_norms(0), row_norms(1), row_norms(2)}; - std::array ref = { - static_cast(p.reference_unit_cell->a), - static_cast(p.reference_unit_cell->b), - static_cast(p.reference_unit_cell->c) - }; - std::sort(obs.begin(), obs.end()); - std::sort(ref.begin(), ref.end()); - - bool lengths_ok = true; - for (int k = 0; k < 3; ++k) { - const float denom = std::max(ref[k], REFINE_MIN_REFERENCE_LENGTH_EPSILON); - const float rel_dev = std::abs(obs[k] - ref[k]) / denom; - if (rel_dev > p.dist_tolerance_vs_reference) { - lengths_ok = false; - break; - } - } - if (!lengths_ok) - continue; - } else { - if (row_norms.minCoeff() < p.min_length_A || row_norms.maxCoeff() > p.max_length_A) - continue; - - float alpha = std::acos(cell_block.row(1).normalized().dot(cell_block.row(2).normalized())) * 180.0f / M_PI; - float beta = std::acos(cell_block.row(0).normalized().dot(cell_block.row(2).normalized())) * 180.0f / M_PI; - float gamma = std::acos(cell_block.row(0).normalized().dot(cell_block.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, p.indexing_tolerance, indexed_spot_count); - - if (indexed_spot_count < p.viable_cell_min_spots) - continue; - - candidates.emplace_back(RefinedCandidate{ - .cell = cell, - .score = scores(i), - .volume = std::abs(cell.determinant()), - .indexed_spot_count = indexed_spot_count, - .indexed_mask = std::move(indexed_mask) - }); - } - - std::sort(candidates.begin(), candidates.end(), - [](const RefinedCandidate &a, const RefinedCandidate &b) { - const auto max_spots = std::max(a.indexed_spot_count, b.indexed_spot_count); - const auto min_spots = std::min(a.indexed_spot_count, b.indexed_spot_count); - const bool spot_counts_close = (max_spots > 0) - && (static_cast(min_spots) / static_cast(max_spots) - >= REFINE_CANDIDATE_SPOT_COUNT_RATIO_THRESHOLD); - - if (!spot_counts_close) - return a.indexed_spot_count > b.indexed_spot_count; - - const float max_volume = std::max(a.volume, b.volume); - const float min_volume = std::max(std::min(a.volume, b.volume), REFINE_MIN_VOLUME_EPSILON); - const bool volume_differs = (max_volume / min_volume) > REFINE_CANDIDATE_VOLUME_RATIO_THRESHOLD; - - if (volume_differs) - return a.volume < b.volume; - - if (a.score != b.score) - return a.score < b.score; - - return a.indexed_spot_count > b.indexed_spot_count; - }); - - std::vector accepted; - - for (const auto &candidate : candidates) { - bool too_similar = false; - - 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) - accepted.emplace_back(candidate); - } - - ret.reserve(accepted.size()); - - for (auto &candidate : accepted) { - auto cell = candidate.cell; - if (cell.determinant() < .0f) - cell = -cell; - - ret.emplace_back( - Coord(cell(0,0), cell(0,1), cell(0,2)), - Coord(cell(1,0), cell(1,1), cell(1,2)), - Coord(cell(2,0), cell(2,1), cell(2,2)) - ); - } - - return ret; -} diff --git a/image_analysis/indexing/FFBIDXIndexer.cpp b/image_analysis/indexing/FFBIDXIndexer.cpp index d3a6ed97..731efa36 100644 --- a/image_analysis/indexing/FFBIDXIndexer.cpp +++ b/image_analysis/indexing/FFBIDXIndexer.cpp @@ -2,7 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include "FFBIDXIndexer.h" -#include "EigenRefine.h" +#include "PostIndexingRefinement.h" void FFBIDXIndexer::SetupUnitCell(const std::optional &cell) { if (!cell.has_value()) diff --git a/image_analysis/indexing/FFTIndexer.cpp b/image_analysis/indexing/FFTIndexer.cpp index 925f5c2a..d4c8d848 100644 --- a/image_analysis/indexing/FFTIndexer.cpp +++ b/image_analysis/indexing/FFTIndexer.cpp @@ -3,7 +3,7 @@ #include "FFTIndexer.h" #include -#include "EigenRefine.h" +#include "PostIndexingRefinement.h" FFTIndexer::FFTIndexer(const IndexingSettings &settings) : max_length_A(settings.GetFFT_MaxUnitCell_A()), diff --git a/image_analysis/indexing/FFTIndexerCPU.cpp b/image_analysis/indexing/FFTIndexerCPU.cpp index 2b909f29..abe6cc91 100644 --- a/image_analysis/indexing/FFTIndexerCPU.cpp +++ b/image_analysis/indexing/FFTIndexerCPU.cpp @@ -2,7 +2,6 @@ // SPDX-License-Identifier: GPL-3.0-only #include "FFTIndexerCPU.h" -#include "EigenRefine.h" #include #include diff --git a/image_analysis/indexing/PostIndexingRefinement.cpp b/image_analysis/indexing/PostIndexingRefinement.cpp new file mode 100644 index 00000000..51260f54 --- /dev/null +++ b/image_analysis/indexing/PostIndexingRefinement.cpp @@ -0,0 +1,258 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "PostIndexingRefinement.h" + +namespace { + struct config_ifssr final { + float threshold_contraction = .8; // contract error threshold by this value in every iteration + float max_distance = .00075; // max distance to reciprocal spots for inliers + unsigned min_spots = 8; // minimum number of spots to fit against + unsigned max_iter = 32; // max number of iterations + }; + + static std::pair score_parts(float score) noexcept { + float nsp = -std::floor(score); + float s = score + nsp; + return std::make_pair(nsp - 1, s); + } + + struct RefinedCandidate { + Eigen::Matrix3f cell; + float score; + float volume; + int64_t indexed_spot_count; + std::vector indexed_mask; + }; + + static inline Eigen::MatrixX3 CalculateResiduals( + const Eigen::Ref> &spots, + const Eigen::Matrix3f &cell) { + Eigen::MatrixX3 miller = (spots * cell).array().round().matrix(); + Eigen::MatrixX3 resid = miller * cell.inverse(); + resid -= spots; + return resid; + } + + static inline std::vector ComputeIndexedMask( + const Eigen::Ref> &spots, + const Eigen::Matrix3f &cell, + float indexing_tolerance, + int64_t &indexed_spot_count) { + const float indexing_tolerance_sq = indexing_tolerance * indexing_tolerance; + const Eigen::MatrixX3 resid = CalculateResiduals(spots, cell); + + std::vector mask(spots.rows(), 0); + indexed_spot_count = 0; + + for (int i = 0; i < spots.rows(); ++i) { + if (resid.row(i).squaredNorm() < indexing_tolerance_sq) { + mask[i] = 1; + indexed_spot_count++; + } + } + + return mask; + } + + template + static void RefineCandidateCells(const Eigen::Ref> &spots, + Eigen::DenseBase &cells, + Eigen::DenseBase &scores, + const config_ifssr &cifssr, + unsigned block = 0, unsigned nblocks = 1) { + using namespace Eigen; + using Mx3 = MatrixX3; + using M3 = Matrix3; + + const unsigned nspots = spots.rows(); + const unsigned ncells = scores.rows(); + VectorX below{nspots}; + MatrixX3 sel{nspots, 3u}; + Mx3 resid{nspots, 3u}; + Mx3 miller{nspots, 3u}; + M3 cell; + + const unsigned blocksize = (ncells + nblocks - 1u) / nblocks; + const unsigned startcell = block * blocksize; + const unsigned endcell = std::min(startcell + blocksize, ncells); + + for (unsigned j = startcell; j < endcell; j++) { + if (nspots < cifssr.min_spots) { + scores(j) = float{1.}; + continue; + } + + cell = cells.block(3u * j, 0u, 3u, 3u).transpose(); // cell: col vectors + const float scale = cell.colwise().norm().minCoeff(); + float threshold = score_parts(scores[j]).second / scale; + + for (unsigned niter = 1; niter < cifssr.max_iter && threshold > cifssr.max_distance; niter++) { + miller = (spots * cell).array().round().matrix(); + resid = miller * cell.inverse(); + resid -= spots; + + below = (resid.rowwise().norm().array() < threshold); + if (below.count() < cifssr.min_spots) + break; + + threshold *= cifssr.threshold_contraction; + sel.colwise() = below; + HouseholderQR qr{sel.select(spots, .0f)}; + cell = qr.solve(sel.select(miller, .0f)); + } + + resid = CalculateResiduals(spots, cell); + + ArrayX dist = resid.rowwise().norm(); + auto nth = std::begin(dist) + (cifssr.min_spots - 1); + std::nth_element(std::begin(dist), nth, std::end(dist)); + scores(j) = *nth; + + cells.block(3u * j, 0u, 3u, 3u) = cell.transpose(); + } + } +} + +std::vector Refine(const std::vector &in_spots, + size_t nspots, + Eigen::MatrixX3 &oCell, + Eigen::VectorX &scores, + RefineParameters &p) { + std::vector ret; + + Eigen::MatrixX3 spots(in_spots.size(), 3u); + + for (int i = 0; i < in_spots.size(); i++) { + spots(i, 0u) = in_spots[i].x; + spots(i, 1u) = in_spots[i].y; + spots(i, 2u) = in_spots[i].z; + } + + config_ifssr cifssr{ + .min_spots = static_cast(p.viable_cell_min_spots) + }; + + RefineCandidateCells(spots.topRows(nspots), oCell, scores, cifssr); + + std::vector candidates; + + for (int i = 0; i < scores.size(); i++) { + auto cell_block = oCell.block(3u * i, 0u, 3u, 3u); + Eigen::Matrix3f cell = cell_block.transpose(); + Eigen::Vector3f row_norms = cell_block.rowwise().norm(); + + if (p.reference_unit_cell) { + std::array obs = {row_norms(0), row_norms(1), row_norms(2)}; + std::array ref = { + static_cast(p.reference_unit_cell->a), + static_cast(p.reference_unit_cell->b), + static_cast(p.reference_unit_cell->c) + }; + std::sort(obs.begin(), obs.end()); + std::sort(ref.begin(), ref.end()); + + bool lengths_ok = true; + for (int k = 0; k < 3; ++k) { + const float denom = std::max(ref[k], REFINE_MIN_REFERENCE_LENGTH_EPSILON); + const float rel_dev = std::abs(obs[k] - ref[k]) / denom; + if (rel_dev > p.dist_tolerance_vs_reference) { + lengths_ok = false; + break; + } + } + if (!lengths_ok) + continue; + } else { + if (row_norms.minCoeff() < p.min_length_A || row_norms.maxCoeff() > p.max_length_A) + continue; + + float alpha = std::acos(cell_block.row(1).normalized().dot(cell_block.row(2).normalized())) * 180.0f / M_PI; + float beta = std::acos(cell_block.row(0).normalized().dot(cell_block.row(2).normalized())) * 180.0f / M_PI; + float gamma = std::acos(cell_block.row(0).normalized().dot(cell_block.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, p.indexing_tolerance, indexed_spot_count); + + if (indexed_spot_count < p.viable_cell_min_spots) + continue; + + candidates.emplace_back(RefinedCandidate{ + .cell = cell, + .score = scores(i), + .volume = std::abs(cell.determinant()), + .indexed_spot_count = indexed_spot_count, + .indexed_mask = std::move(indexed_mask) + }); + } + + std::sort(candidates.begin(), candidates.end(), + [](const RefinedCandidate &a, const RefinedCandidate &b) { + const auto max_spots = std::max(a.indexed_spot_count, b.indexed_spot_count); + const auto min_spots = std::min(a.indexed_spot_count, b.indexed_spot_count); + const bool spot_counts_close = (max_spots > 0) + && (static_cast(min_spots) / static_cast(max_spots) + >= REFINE_CANDIDATE_SPOT_COUNT_RATIO_THRESHOLD); + + if (!spot_counts_close) + return a.indexed_spot_count > b.indexed_spot_count; + + const float max_volume = std::max(a.volume, b.volume); + const float min_volume = std::max(std::min(a.volume, b.volume), REFINE_MIN_VOLUME_EPSILON); + const bool volume_differs = (max_volume / min_volume) > REFINE_CANDIDATE_VOLUME_RATIO_THRESHOLD; + + if (volume_differs) + return a.volume < b.volume; + + if (a.score != b.score) + return a.score < b.score; + + return a.indexed_spot_count > b.indexed_spot_count; + }); + + std::vector accepted; + + for (const auto &candidate: candidates) { + bool too_similar = false; + + 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) + accepted.emplace_back(candidate); + } + + ret.reserve(accepted.size()); + + for (auto &candidate: accepted) { + auto cell = candidate.cell; + if (cell.determinant() < .0f) + cell = -cell; + + ret.emplace_back( + Coord(cell(0, 0), cell(0, 1), cell(0, 2)), + Coord(cell(1, 0), cell(1, 1), cell(1, 2)), + Coord(cell(2, 0), cell(2, 1), cell(2, 2)) + ); + } + + return ret; +} diff --git a/image_analysis/indexing/PostIndexingRefinement.h b/image_analysis/indexing/PostIndexingRefinement.h new file mode 100644 index 00000000..e69c22f6 --- /dev/null +++ b/image_analysis/indexing/PostIndexingRefinement.h @@ -0,0 +1,36 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "../common/CrystalLattice.h" + +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_MIN_VOLUME_EPSILON = 1e-12f; +constexpr float REFINE_MIN_REFERENCE_LENGTH_EPSILON = 1e-6f; + +struct RefineParameters { + int64_t viable_cell_min_spots; + float dist_tolerance_vs_reference; + std::optional reference_unit_cell; + float min_length_A; + float max_length_A; + float min_angle_deg; + float max_angle_deg; + float indexing_tolerance; +}; + +std::vector Refine(const std::vector &in_spots, + size_t nspots, + Eigen::MatrixX3 &oCell, + Eigen::VectorX &scores, + RefineParameters &p); \ No newline at end of file