diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 33ed4ed60..5680f254b 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -5,6 +5,8 @@ ADD_LIBRARY(JFJochImageAnalysis STATIC MXAnalysisAfterFPGA.cpp SpotAnalyze.cpp SpotAnalyze.h + RotationIndexer.cpp + RotationIndexer.h dark_mask_analysis/DarkMaskAnalysis.cpp dark_mask_analysis/DarkMaskAnalysis.h) diff --git a/image_analysis/RotationIndexer.cpp b/image_analysis/RotationIndexer.cpp new file mode 100644 index 000000000..ebe5a83e0 --- /dev/null +++ b/image_analysis/RotationIndexer.cpp @@ -0,0 +1,87 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "RotationIndexer.h" + +RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPool &indexer) + : experiment(x), + index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()), + image_spot0(x.GetImageNum(), 0), + image_nspots(x.GetImageNum(), 0), + axis_(x.GetGoniometer()), + geom_(x.GetDiffractionGeometry()), + 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; +} + +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 (!indexed_lattice && image >= last_accumulated_image + image_stride) { + v_.reserve(v_.size() + spots.size()); + coords_.reserve(coords_.size() + spots.size()); + + image_spot0[image] = v_.size() - 1; + int i = 0; + for (const auto &s: spots) { + if (index_ice_rings || !s.ice_ring) { + v_.emplace_back(s); + coords_.emplace_back(rot * s.ReciprocalCoord(geom_)); + i++; + } + } + image_nspots[image] = i; + + accumulated_images++; + last_accumulated_image = image; + + if (accumulated_images >= min_images_for_indexing && image >= first_image_to_try_indexing) { + auto indexer_result = indexer_.Run(experiment, coords_).get(); + if (!indexer_result.lattice.empty()) + indexed_lattice = indexer_result.lattice[0]; + + // Find lattice type + // Run refinement + } + } + + if (indexed_lattice) + return indexed_lattice->Multiply(rot); + return {}; +} diff --git a/image_analysis/RotationIndexer.h b/image_analysis/RotationIndexer.h new file mode 100644 index 000000000..c72f7d862 --- /dev/null +++ b/image_analysis/RotationIndexer.h @@ -0,0 +1,59 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#ifndef JFJOCH_SPOTROTATIONSTORE_H +#define JFJOCH_SPOTROTATIONSTORE_H + +#include +#include + +#include "../common/DiffractionSpot.h" +#include "../common/DiffractionExperiment.h" +#include "indexing/IndexerThreadPool.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 + +class RotationIndexer { + mutable std::mutex m; + const DiffractionExperiment& experiment; + + constexpr static float min_accum_angle_deg = 10.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 image_spot0; + std::vector image_nspots; + + std::vector v_; + std::vector coords_; + std::optional axis_; + const DiffractionGeometry geom_; + + 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; + +public: + RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer); + void SetLattice(const CrystalLattice &lattice); + std::optional ProcessImage(int64_t image, const std::vector& spots); +}; + + +#endif //JFJOCH_SPOTROTATIONSTORE_H \ No newline at end of file