diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 9c1ba7d3..21dda4a2 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -26,10 +26,10 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool } IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotation(DataMessage &msg) { + const auto result = rotation_indexer->GetLattice(); + IndexingOutcome outcome(experiment); - rotation_indexer->ProcessImage(msg.number, msg.spots); - auto result = rotation_indexer->GetLattice(); if (result) { // For rotation indexing, indexing rate is calculated only for frames, where "global" rotation indexing solution was found msg.indexing_result = false; @@ -292,8 +292,13 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, } std::optional IndexAndRefine::Finalize() { - if (rotation_indexer) + if (rotation_indexer) { + if (const auto latt = rotation_indexer->GetLattice()) + return latt; + + rotation_indexer->RunIndexing(); return rotation_indexer->GetLattice(); + } return {}; } diff --git a/image_analysis/RotationIndexer.cpp b/image_analysis/RotationIndexer.cpp index 9904b0b1..6abe30d2 100644 --- a/image_analysis/RotationIndexer.cpp +++ b/image_analysis/RotationIndexer.cpp @@ -16,30 +16,11 @@ RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPo geom_(x.GetDiffractionGeometry()), updated_geom_(geom_), indexer_(indexer) { - if (axis_) { - 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 = std::max(0, x.GetImageNum() - 1); - next_image_to_try_indexing = first_image_to_try_indexing; - image_stride = 1; - } else { - first_image_to_try_indexing = std::max(min_images_for_indexing, - x.GetIndexingSettings().GetRotationIndexingMinAngularRange_deg() / angle_norm_deg); - next_image_to_try_indexing = first_image_to_try_indexing; - image_stride = std::ceil(x.GetIndexingSettings().GetRotationIndexingAngularStride_deg() / angle_norm_deg); - if (image_stride == 0) - image_stride = 1; - } - } - } } -IndexerResult RotationIndexer::RunIndexing() const { +void RotationIndexer::RunIndexing() { + std::unique_lock ul(m); + std::vector coords; coords.reserve(max_spots_per_image * v_.size()); @@ -50,13 +31,9 @@ IndexerResult RotationIndexer::RunIndexing() const { for (const auto &s: v_[i]) coords.emplace_back(rot * s.ReciprocalCoord(geom_)); } + const auto indexer_result = indexer_.Run(experiment, coords); - return indexer_.Run(experiment, coords); -} - -void RotationIndexer::TryIndex() { - if (auto indexer_result = RunIndexing(); - !indexer_result.lattice.empty() && indexer_result.lattice[0].CalcVolume() > 1.0) { + if (!indexer_result.lattice.empty() && indexer_result.lattice[0].CalcVolume() > 1.0) { auto sg = experiment.GetGemmiSpaceGroup(); if (sg) { search_result_ = LatticeSearchResult{ @@ -97,15 +74,6 @@ void RotationIndexer::TryIndex() { axis_ = data.axis; } } - - if (!indexed_lattice) { - if (indexing_range_multiplier < max_indexing_range_multiplier) { - indexing_range_multiplier++; - next_image_to_try_indexing = first_image_to_try_indexing * indexing_range_multiplier; - } else { - next_image_to_try_indexing = INT64_MAX; - } - } } void RotationIndexer::ProcessImage(int64_t image, const std::vector &spots) { @@ -115,15 +83,12 @@ void RotationIndexer::ProcessImage(int64_t image, const std::vector if (!axis_) return; - if (accumulated_images >= max_images_for_indexing) + if (accumulated_spots >= max_spots) return; if (indexed_lattice) return; - if (image < last_accumulated_image + image_stride) - return; - v_[image].reserve(spots.size()); for (const auto &s: spots) { @@ -131,26 +96,20 @@ void RotationIndexer::ProcessImage(int64_t image, const std::vector v_[image].emplace_back(s); } - if (v_[image].size() > max_spots_per_image) { - std::ranges::nth_element(v_[image], v_[image].begin() + max_spots_per_image, + // truncate spots, so we don't get above max_spots (total) and max_spots_per_image (for this image) + size_t max_spots_limit = std::min(max_spots_per_image, max_spots - accumulated_spots); + + if (v_[image].size() > max_spots_limit) { + std::ranges::nth_element(v_[image], v_[image].begin() + max_spots_limit, [](const SpotToSave &a, const SpotToSave &b) { return a.intensity > b.intensity; } ); - v_[image].resize(max_spots_per_image); + v_[image].resize(max_spots_limit); } - accumulated_images++; - last_accumulated_image = image; - - const bool short_scan_last_image = - (experiment.GetImageNum() < min_images_for_indexing) && - (image >= experiment.GetImageNum() - 1); - - if ((accumulated_images >= min_images_for_indexing || short_scan_last_image) && - image >= next_image_to_try_indexing) - TryIndex(); + accumulated_spots += v_[image].size(); } std::optional RotationIndexer::GetLattice() const { @@ -166,6 +125,3 @@ std::optional RotationIndexer::GetLattice() const { }; } -int64_t RotationIndexer::GetAccumulatedImages() const { - return accumulated_images; -} diff --git a/image_analysis/RotationIndexer.h b/image_analysis/RotationIndexer.h index ccaca1b9..61dd947f 100644 --- a/image_analysis/RotationIndexer.h +++ b/image_analysis/RotationIndexer.h @@ -1,8 +1,7 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only -#ifndef JFJOCH_ROTATIONINDEXER_H -#define JFJOCH_ROTATIONINDEXER_H +#pragma once #include #include @@ -12,12 +11,6 @@ #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; @@ -29,13 +22,10 @@ class RotationIndexer { mutable std::mutex m; const DiffractionExperiment& experiment; - constexpr static int64_t max_spots_per_image = 200; - constexpr static int64_t max_images_for_indexing = 100; - constexpr static int64_t min_images_for_indexing = 10; - constexpr static int64_t max_indexing_range_multiplier = 3; + constexpr static size_t max_spots = 10000; + constexpr static size_t max_spots_per_image = 200; const bool index_ice_rings; - float angle_norm_deg = 1.0f; std::vector> v_; @@ -46,24 +36,14 @@ class RotationIndexer { IndexerThreadPool &indexer_; - int64_t last_accumulated_image = -1; - int64_t accumulated_images = 0; - - int64_t image_stride = 1; - int64_t first_image_to_try_indexing = INT64_MAX; - int64_t next_image_to_try_indexing = INT64_MAX; - int64_t indexing_range_multiplier = 1; + size_t accumulated_spots = 0; std::optional indexed_lattice; - IndexerResult RunIndexing() const; - - void TryIndex(); public: RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer); void ProcessImage(int64_t image, const std::vector& spots); + void RunIndexing(); std::optional GetLattice() const; - int64_t GetAccumulatedImages() const; + size_t GetAccumulatedImages() const; }; - -#endif //JFJOCH_ROTATIONINDEXER_H \ No newline at end of file diff --git a/tests/RotationIndexerTest.cpp b/tests/RotationIndexerTest.cpp index 029dbac7..c6634290 100644 --- a/tests/RotationIndexerTest.cpp +++ b/tests/RotationIndexerTest.cpp @@ -70,7 +70,11 @@ TEST_CASE("RotationIndexer") { s.indexed = true; spots.push_back(s); } + indexer.ProcessImage(img, spots); + if (img == 30) + indexer.RunIndexing(); + auto result = indexer.GetLattice(); if (result.has_value()) cnt++; @@ -82,112 +86,12 @@ TEST_CASE("RotationIndexer") { REQUIRE(ret.has_value()); auto uc = ret->lattice.GetUnitCell(); auto uc_ref = latt_base.GetUnitCell(); - REQUIRE(std::fabs(uc.a - uc_ref.a) < 0.02 ); - REQUIRE(std::fabs(uc.b - uc_ref.b) < 0.02 ); - REQUIRE(std::fabs(uc.c - uc_ref.c) < 0.02 ); + REQUIRE(std::fabs(uc.a - uc_ref.a) < 0.1); + REQUIRE(std::fabs(uc.b - uc_ref.b) < 0.1); + REQUIRE(std::fabs(uc.c - uc_ref.c) < 0.1); REQUIRE(std::fabs(uc.alpha - uc_ref.alpha) < 0.1); REQUIRE(std::fabs(uc.beta - uc_ref.beta) < 0.1); REQUIRE(std::fabs(uc.gamma - uc_ref.gamma) < 0.1); CHECK(ret->search_result.centering == 'P'); CHECK(ret->search_result.system == gemmi::CrystalSystem::Orthorhombic); } - -TEST_CASE("RotationIndexer short scan indexes only at the end") { - DiffractionExperiment exp_i; - exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV) - .BeamX_pxl(1000) - .BeamY_pxl(1000) - .PoniRot1_rad(0.01) - .PoniRot2_rad(0.02) - .DetectorDistance_mm(200) - .ImagesPerTrigger(6); - - IndexingSettings settings; -#ifdef JFJOCH_USE_CUDA - settings.Algorithm(IndexingAlgorithmEnum::FFT); -#elif JFJOCH_USE_FFTW - settings.Algorithm(IndexingAlgorithmEnum::FFTW); -#else - return; -#endif - - settings.RotationIndexing(true).RotationIndexingAngularStride_deg(1.0).RotationIndexingMinAngularRange_deg(30.0); - exp_i.ImportIndexingSettings(settings); - - CrystalLattice latt_base(40, 50, 80, 90, 90, 90); - latt_base = latt_base.Multiply(RotMatrix(2.0, Coord(sqrt(3)/3,sqrt(3)/3,sqrt(3)/3))); - - GoniometerAxis axis("omega", 0.0f, 30.0f, Coord(1,0,0), std::nullopt); - exp_i.Goniometer(axis); - - BraggPredictionSettings prediction_settings{ - .high_res_A = 1.3, - .ewald_dist_cutoff = 0.002 - }; - - IndexerThreadPool indexer_thread_pool(exp_i.GetIndexingSettings()); - RotationIndexer indexer(exp_i, indexer_thread_pool); - - BraggPrediction prediction; - - for (int img = 0; img < 5; ++img) { - std::vector spots; - const float angle_deg = axis.GetAngle_deg(img) + axis.GetWedge_deg() / 2.0f; - const RotMatrix rot = axis.GetTransformationAngle(angle_deg); - const CrystalLattice latt_img = latt_base.Multiply(rot.transpose()); - - const auto n = prediction.Calc(exp_i, latt_img, prediction_settings); - for (int i = 0; i < n; ++i) { - const auto& r = prediction.GetReflections().at(i); - SpotToSave s{}; - s.x = r.predicted_x; - s.y = r.predicted_y; - s.image = img; - s.intensity = 1.0f; - s.phi = angle_deg; - s.ice_ring = false; - s.indexed = true; - spots.push_back(s); - } - indexer.ProcessImage(img, spots); - auto result = indexer.GetLattice(); - CHECK_FALSE(result.has_value()); - } - - std::vector last_spots; - const int last_img = 5; - const float angle_deg = axis.GetAngle_deg(last_img) + axis.GetWedge_deg() / 2.0f; - const RotMatrix rot = axis.GetTransformationAngle(angle_deg); - const CrystalLattice latt_img = latt_base.Multiply(rot.transpose()); - - const auto n = prediction.Calc(exp_i, latt_img, prediction_settings); - for (int i = 0; i < n; ++i) { - const auto& r = prediction.GetReflections().at(i); - SpotToSave s{}; - s.x = r.predicted_x; - s.y = r.predicted_y; - s.image = last_img; - s.intensity = 1.0f; - s.phi = angle_deg; - s.ice_ring = false; - s.indexed = true; - last_spots.push_back(s); - } - - indexer.ProcessImage(last_img, last_spots); - auto ret_last = indexer.GetLattice(); - REQUIRE(ret_last.has_value()); - - auto ret = indexer.GetLattice(); - REQUIRE(ret.has_value()); - auto uc = ret->lattice.GetUnitCell(); - auto uc_ref = latt_base.GetUnitCell(); - REQUIRE(std::fabs(uc.a - uc_ref.a) < 0.02 ); - REQUIRE(std::fabs(uc.b - uc_ref.b) < 0.02 ); - REQUIRE(std::fabs(uc.c - uc_ref.c) < 0.02 ); - REQUIRE(std::fabs(uc.alpha - uc_ref.alpha) < 0.1); - REQUIRE(std::fabs(uc.beta - uc_ref.beta) < 0.1); - REQUIRE(std::fabs(uc.gamma - uc_ref.gamma) < 0.1); - CHECK(ret->search_result.centering == 'P'); - CHECK(ret->search_result.system == gemmi::CrystalSystem::Orthorhombic); -} \ No newline at end of file