From cb983f5c742aae241950fca9cafebd2c81457c40 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 30 Nov 2025 13:20:26 +0100 Subject: [PATCH] IndexAndRefine: Another iteration --- image_analysis/CMakeLists.txt | 4 +- image_analysis/IndexAndRefine.cpp | 208 +++++++---------------------- image_analysis/IndexAndRefine.h | 37 +---- image_analysis/RotationIndexer.cpp | 163 ++++++++++++++++++++++ image_analysis/RotationIndexer.h | 64 +++++++++ 5 files changed, 280 insertions(+), 196 deletions(-) create mode 100644 image_analysis/RotationIndexer.cpp create mode 100644 image_analysis/RotationIndexer.h diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 6cdda11c..83e9bbda 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -6,7 +6,9 @@ ADD_LIBRARY(JFJochImageAnalysis STATIC IndexAndRefine.cpp IndexAndRefine.h dark_mask_analysis/DarkMaskAnalysis.cpp - dark_mask_analysis/DarkMaskAnalysis.h) + dark_mask_analysis/DarkMaskAnalysis.h + RotationIndexer.cpp + RotationIndexer.h) FIND_PACKAGE(Eigen3 3.4 REQUIRED NO_MODULE) # provides Eigen3::Eigen diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index c0867ec2..d52d90df 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -11,99 +11,17 @@ #include "lattice_search/LatticeSearch.h" IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer) - : min_accum_angle_deg(x.GetIndexingSettings().GetRotationIndexingMinAngularRange_deg()), - stride_angle_deg(x.GetIndexingSettings().GetRotationIndexingAngularStride_deg()), - index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()), + : index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()), experiment(x), geom_(x.GetDiffractionGeometry()), - updated_geom_(geom_), indexer_(indexer) { - auto axis = x.GetGoniometer(); - - if (indexer && axis && x.GetIndexingSettings().GetRotationIndexing()) { - float angle_norm_deg = std::fabs(axis->GetIncrement_deg()); - if (angle_norm_deg > 1e-6) { - axis_ = axis; - if (x.GetImageNum() > min_images_for_indexing) { - first_image_to_try_indexing = std::max(min_images_for_indexing, - min_accum_angle_deg / angle_norm_deg); - image_stride = std::ceil(stride_angle_deg / angle_norm_deg); - if (image_stride == 0) - image_stride = 1; - } - } - } + if (indexer && x.GetGoniometer() && x.GetIndexingSettings().GetRotationIndexing()) + rotation_indexer = std::make_unique(x, *indexer); } void IndexAndRefine::SetLattice(const CrystalLattice &lattice) { - if (axis_) { - std::unique_lock ul(m); - indexed_lattice = lattice; - } -} - -void IndexAndRefine::TryIndexRotationData() { - if (!indexer_) - return; - - // Index - std::vector v_sel; - std::vector coords_sel; - - if (coords_.size() > max_spots) { - // Indices into v_ / coords_ - std::vector idx(coords_.size()); - std::iota(idx.begin(), idx.end(), std::size_t{0}); - - // Sort indices by descending intensity - std::partial_sort( - idx.begin(), - idx.begin() + max_spots, - idx.end(), - [this](std::size_t a, std::size_t b) { - // Replace `.intensity` with the actual SpotToSave intensity member - return v_[a].intensity > v_[b].intensity; - } - ); - - v_sel.reserve(max_spots); - coords_sel.reserve(max_spots); - - for (std::size_t i = 0; i < max_spots; ++i) { - const std::size_t k = idx[i]; - v_sel.emplace_back(v_[k]); - coords_sel.emplace_back(coords_[k]); - } - } else { - v_sel = v_; - coords_sel = coords_; - } - - auto indexer_result = indexer_->Run(experiment, coords_sel).get(); - if (!indexer_result.lattice.empty()) { - // Find lattice type - search_result_ = LatticeSearch(indexer_result.lattice[0]); - // Run refinement - - DiffractionExperiment experiment_copy(experiment); - XtalOptimizerData data{ - .geom = experiment_copy.GetDiffractionGeometry(), - .latt = search_result_.conventional, - .crystal_system = search_result_.system, - .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), - .refine_beam_center = true, - .refine_distance_mm = true, - .axis = axis_ - }; - - if (data.crystal_system == gemmi::CrystalSystem::Trigonal) - data.crystal_system = gemmi::CrystalSystem::Hexagonal; - - if (XtalOptimizer(data, v_sel)) { - indexed_lattice = data.latt; - updated_geom_ = data.geom; - } - } + if (rotation_indexer) + rotation_indexer->SetLattice(lattice); } void IndexAndRefine::ProcessImage(DataMessage &msg, @@ -125,13 +43,24 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, bool beam_center_updated = false; - if (!axis_) { + if (rotation_indexer) { + auto result = rotation_indexer->ProcessImage(msg.number, msg.spots); + if (result) { + lattice_candidate = result->lattice; + experiment_copy.BeamX_pxl(result->geom.GetBeamX_pxl()) + .BeamY_pxl(result->geom.GetBeamY_pxl()) + .DetectorDistance_mm(result->geom.GetDetectorDistance_mm()); + symmetry.centering = result->search_result.centering; + symmetry.niggli_class = result->search_result.niggli_class; + symmetry.crystal_system = result->search_result.system; + } + } else { // Convert input spots to reciprocal space std::vector recip; recip.reserve(msg.spots.size()); for (const auto &i: msg.spots) { if (index_ice_rings || !i.ice_ring) - recip.push_back(i.ReciprocalCoord(updated_geom_)); + recip.push_back(i.ReciprocalCoord(geom_)); } auto indexer_result = indexer_->Run(experiment, recip).get(); @@ -163,72 +92,40 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, latt = sym_result.conventional; } + XtalOptimizerData data{ + .geom = geom_, + .latt = latt, + .crystal_system = symmetry.crystal_system, + .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), + .refine_beam_center = true, + .refine_distance_mm = false + }; - lattice_candidate = latt; - } - } else { - std::unique_lock ul(m); - const auto rot = axis_->GetTransformation(msg.number); + if (symmetry.crystal_system == gemmi::CrystalSystem::Trigonal) + data.crystal_system = gemmi::CrystalSystem::Hexagonal; - if (msg.number >= last_accumulated_image + image_stride) { - v_.reserve(v_.size() + msg.spots.size()); - coords_.reserve(coords_.size() + msg.spots.size()); - - for (const auto &s: msg.spots) { - if (index_ice_rings || !s.ice_ring) { - v_.emplace_back(s); - coords_.emplace_back(rot * s.ReciprocalCoord(geom_)); - } + switch (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm()) { + case GeomRefinementAlgorithmEnum::None: + break; + case GeomRefinementAlgorithmEnum::BeamCenter: + if (XtalOptimizer(data, msg.spots)) { + experiment_copy.BeamX_pxl(data.geom.GetBeamX_pxl()).BeamY_pxl(data.geom.GetBeamY_pxl()); + beam_center_updated = true; + } + break; } - accumulated_images++; - last_accumulated_image = msg.number; + lattice_candidate = data.latt; - if (!indexed_lattice - && accumulated_images >= min_images_for_indexing - && msg.number >= first_image_to_try_indexing) - TryIndexRotationData(); - } - - if (indexed_lattice) { - lattice_candidate = indexed_lattice->Multiply(rot.transpose()); - symmetry = LatticeMessage{ - .centering = search_result_.centering, - .niggli_class = search_result_.niggli_class, - .crystal_system = search_result_.system - }; - beam_center_updated = true; - experiment_copy.BeamX_pxl(updated_geom_.GetBeamX_pxl()) - .BeamY_pxl(updated_geom_.GetBeamY_pxl()) - .DetectorDistance_mm(updated_geom_.GetDetectorDistance_mm()); + if (beam_center_updated) { + msg.beam_corr_x = data.beam_corr_x; + msg.beam_corr_y = data.beam_corr_y; + } } } if (lattice_candidate) { - XtalOptimizerData data{ - .geom = updated_geom_, - .latt = lattice_candidate.value(), - .crystal_system = symmetry.crystal_system, - .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), - .refine_beam_center = true, - .refine_distance_mm = false - }; - - if (symmetry.crystal_system == gemmi::CrystalSystem::Trigonal) - data.crystal_system = gemmi::CrystalSystem::Hexagonal; - - switch (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm()) { - case GeomRefinementAlgorithmEnum::None: - break; - case GeomRefinementAlgorithmEnum::BeamCenter: - if (XtalOptimizer(data, msg.spots)) { - experiment_copy.BeamX_pxl(data.geom.GetBeamX_pxl()).BeamY_pxl(data.geom.GetBeamY_pxl()); - beam_center_updated = true; - lattice_candidate = data.latt; - } - break; - } - if (AnalyzeIndexing(msg, experiment, data.latt)) { + if (AnalyzeIndexing(msg, experiment, *lattice_candidate)) { msg.lattice_type = symmetry; float ewald_dist_cutoff = 0.001f; @@ -239,11 +136,6 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, ewald_dist_cutoff = experiment.GetBraggIntegrationSettings().GetFixedProfileRadius_recipA().value() * 3.0f; - if (beam_center_updated) { - msg.beam_corr_x = data.beam_corr_x; - msg.beam_corr_y = data.beam_corr_y; - } - if (spot_finding_settings.quick_integration) { auto res = BraggIntegrate2D(experiment_copy, image, *lattice_candidate, prediction, ewald_dist_cutoff, msg.number, symmetry.centering); @@ -262,17 +154,7 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, } std::optional IndexAndRefine::Finalize(bool retry) { - if (axis_) { - std::unique_lock ul(m); - - if (!indexed_lattice || retry) - TryIndexRotationData(); - if (indexed_lattice) - return RotationIndexerResult{ - .lattice = indexed_lattice.value(), - .search_result = search_result_, - .geom = geom_ - }; - } + if (rotation_indexer) + return rotation_indexer->Finalize(retry); return {}; } diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 005292d2..0f186ad6 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -1,8 +1,8 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only -#ifndef JFJOCH_SPOTROTATIONSTORE_H -#define JFJOCH_SPOTROTATIONSTORE_H +#ifndef JFJOCH_INDEXANDREFINE_H +#define JFJOCH_INDEXANDREFINE_H #include #include @@ -12,46 +12,19 @@ #include "bragg_integration/BraggPrediction.h" #include "indexing/IndexerThreadPool.h" #include "lattice_search/LatticeSearch.h" - -// RotationIndexer works as following: -// 1. First accumulates spot results from rotation images (only images within a certain stride are included) -// 2. When a minimum number of images is reached (at least 10 images AND at least 10 deg), indexing is attempted -// 3. If indexing is successful - lattice is provided that is used by subsequent images -// 4. If indexing is not-successful - accumulation procedure is continued - -struct RotationIndexerResult { - CrystalLattice lattice; - LatticeSearchResult search_result; - DiffractionGeometry geom; -}; +#include "RotationIndexer.h" class IndexAndRefine { - constexpr static int64_t max_spots = 60000; - constexpr static int64_t min_images_for_indexing = 10; - - const float min_accum_angle_deg = 20.0; - const float stride_angle_deg = 0.5; const bool index_ice_rings; const DiffractionExperiment& experiment; const DiffractionGeometry geom_; - mutable std::mutex m; - std::vector v_; - std::vector coords_; - DiffractionGeometry updated_geom_; - LatticeSearchResult search_result_; - int64_t last_accumulated_image = -1; - int64_t accumulated_images = 0; std::optional indexed_lattice; std::optional axis_; IndexerThreadPool *indexer_; - - int64_t image_stride = 1; - int64_t first_image_to_try_indexing = INT64_MAX; - - void TryIndexRotationData(); + std::unique_ptr rotation_indexer; public: IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer); void SetLattice(const CrystalLattice &lattice); @@ -60,4 +33,4 @@ public: }; -#endif //JFJOCH_SPOTROTATIONSTORE_H +#endif //JFJOCH_INDEXANDREFINE_H diff --git a/image_analysis/RotationIndexer.cpp b/image_analysis/RotationIndexer.cpp new file mode 100644 index 00000000..b654ade7 --- /dev/null +++ b/image_analysis/RotationIndexer.cpp @@ -0,0 +1,163 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "RotationIndexer.h" + +#include "geom_refinement/XtalOptimizer.h" +#include "indexing/FFTIndexer.h" +#include "lattice_search/LatticeSearch.h" + +RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPool &indexer) + : experiment(x), + index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()), + axis_(x.GetGoniometer()), + geom_(x.GetDiffractionGeometry()), + updated_geom_(geom_), + indexer_(indexer) { + if (axis_) { + float angle_norm_deg = std::fabs(axis_->GetIncrement_deg()); + if (angle_norm_deg < 1e-6) { + // Guard against rotation close to zero + axis_ = std::nullopt; + } else { + if (x.GetImageNum() < min_images_for_indexing) { + // For short measurements - only indexing at the end + first_image_to_try_indexing = INT64_MAX; + image_stride = 1; + } else { + first_image_to_try_indexing = std::max(min_images_for_indexing, + min_accum_angle_deg / angle_norm_deg); + image_stride = std::ceil(stride_angle_deg / angle_norm_deg); + if (image_stride == 0) + image_stride = 1; + } + } + } +} + +void RotationIndexer::SetLattice(const CrystalLattice &lattice) { + std::unique_lock ul(m); + indexed_lattice = lattice; +} + +void RotationIndexer::TryIndex() { + // Index + std::vector v_sel; + std::vector coords_sel; + + if (coords_.size() > max_spots) { + + // Indices into v_ / coords_ + std::vector idx(coords_.size()); + std::iota(idx.begin(), idx.end(), std::size_t{0}); + + // Sort indices by descending intensity + std::partial_sort( + idx.begin(), + idx.begin() + max_spots, + idx.end(), + [this](std::size_t a, std::size_t b) { + // Replace `.intensity` with the actual SpotToSave intensity member + return v_[a].intensity > v_[b].intensity; + } + ); + + v_sel.reserve(max_spots); + coords_sel.reserve(max_spots); + + for (std::size_t i = 0; i < max_spots; ++i) { + const std::size_t k = idx[i]; + v_sel.emplace_back(v_[k]); + coords_sel.emplace_back(coords_[k]); + } + } else { + v_sel = v_; + coords_sel = coords_; + } + + auto indexer_result = indexer_.Run(experiment, coords_sel).get(); + if (!indexer_result.lattice.empty()) { + // Find lattice type + search_result_ = LatticeSearch(indexer_result.lattice[0]); + // Run refinement + + DiffractionExperiment experiment_copy(experiment); + XtalOptimizerData data{ + .geom = experiment_copy.GetDiffractionGeometry(), + .latt = search_result_.conventional, + .crystal_system = search_result_.system, + .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), + .refine_beam_center = true, + .refine_distance_mm = true, + .axis = axis_ + }; + + if (data.crystal_system == gemmi::CrystalSystem::Trigonal) + data.crystal_system = gemmi::CrystalSystem::Hexagonal; + + if (XtalOptimizer(data, v_sel)) { + indexed_lattice = data.latt; + updated_geom_ = data.geom; + } + } +} + +std::optional RotationIndexer::ProcessImage(int64_t image, const std::vector &spots) { + std::unique_lock ul(m); + + // For non-rotation just ignore the whole procedure + if (!axis_) + return {}; + + const auto rot = axis_->GetTransformation(image); + + if (image >= last_accumulated_image + image_stride) { + v_.reserve(v_.size() + spots.size()); + coords_.reserve(coords_.size() + spots.size()); + + for (const auto &s: spots) { + if (index_ice_rings || !s.ice_ring) { + v_.emplace_back(s); + coords_.emplace_back(rot * s.ReciprocalCoord(geom_)); + } + } + + accumulated_images++; + last_accumulated_image = image; + + if (!indexed_lattice && accumulated_images >= min_images_for_indexing && image >= first_image_to_try_indexing) + TryIndex(); + } + + if (indexed_lattice) { + // Refine indexing lattice with spots found in the image + auto latt_image = indexed_lattice->Multiply(rot.transpose()); + XtalOptimizerData data{ + .geom = updated_geom_, + .latt = latt_image, + .crystal_system = search_result_.system, + .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), + .refine_beam_center = true, + .refine_distance_mm = false + }; + + if (XtalOptimizer(data, spots)) + return RotationIndexerResult{ + .lattice = data.latt, + .search_result = search_result_, + .geom = data.geom + }; + } + return {}; +} + + +std::optional RotationIndexer::Finalize(bool retry) { + if (!indexed_lattice || retry) + TryIndex(); + return RotationIndexerResult{ + .lattice = indexed_lattice.value(), + .search_result = search_result_, + .geom = geom_ + }; +} diff --git a/image_analysis/RotationIndexer.h b/image_analysis/RotationIndexer.h new file mode 100644 index 00000000..75bf3575 --- /dev/null +++ b/image_analysis/RotationIndexer.h @@ -0,0 +1,64 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#ifndef JFJOCH_ROTATIONINDEXER_H +#define JFJOCH_ROTATIONINDEXER_H + +#include +#include + +#include "../common/DiffractionSpot.h" +#include "../common/DiffractionExperiment.h" +#include "indexing/IndexerThreadPool.h" +#include "lattice_search/LatticeSearch.h" + +// RotationIndexer works as following: +// 1. First accumulates spot results from rotation images (only images within a certain stride are included) +// 2. When a minimum number of images is reached (at least 10 images AND at least 10 deg), indexing is attempted +// 3. If indexing is successful - lattice is provided that is used by subsequent images +// 4. If indexing is not-successful - accumulation procedure is continued + +struct RotationIndexerResult { + CrystalLattice lattice; + LatticeSearchResult search_result; + DiffractionGeometry geom; +}; + +class RotationIndexer { + mutable std::mutex m; + const DiffractionExperiment& experiment; + + constexpr static int64_t max_spots = 60000; + constexpr static float min_accum_angle_deg = 20.0; + constexpr static float stride_angle_deg = 0.5; + constexpr static int64_t min_images_for_indexing = 10; + + const bool index_ice_rings; + + std::vector v_; + std::vector coords_; + std::optional axis_; + const DiffractionGeometry geom_; + DiffractionGeometry updated_geom_; + LatticeSearchResult search_result_; + + IndexerThreadPool &indexer_; + + int64_t last_accumulated_image = -1; + int64_t accumulated_images = 0; + + int64_t image_stride; + int64_t first_image_to_try_indexing; + + std::optional indexed_lattice; + + + void TryIndex(); +public: + RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer); + void SetLattice(const CrystalLattice &lattice); + std::optional ProcessImage(int64_t image, const std::vector& spots); + std::optional Finalize(bool retry); +}; + +#endif //JFJOCH_ROTATIONINDEXER_H \ No newline at end of file