diff --git a/image_analysis/RotationIndexer.cpp b/image_analysis/RotationIndexer.cpp index a54acf30..8ba6a04c 100644 --- a/image_analysis/RotationIndexer.cpp +++ b/image_analysis/RotationIndexer.cpp @@ -1,8 +1,10 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only + #include "RotationIndexer.h" -#include + #include "geom_refinement/XtalOptimizer.h" +#include "indexing/FFTIndexer.h" #include "lattice_search/LatticeSearch.h" RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPool &indexer) @@ -38,6 +40,71 @@ void RotationIndexer::SetLattice(const CrystalLattice &lattice) { indexed_lattice = lattice; } +#include + +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 + auto sym_result = LatticeSearch(indexer_result.lattice[0]); + indexed_lattice = indexer_result.lattice[0]; +/* + // Run refinement + + DiffractionExperiment experiment_copy(experiment); + XtalOptimizerData data{ + .geom = experiment_copy.GetDiffractionGeometry(), + .latt = sym_result.conventional, + .crystal_system = sym_result.system, + .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), + .refine_beam_center = true, + .refine_distance_mm = false, + .axis = axis_ + }; + + if (data.crystal_system == gemmi::CrystalSystem::Trigonal) + data.crystal_system = gemmi::CrystalSystem::Hexagonal; + + if (XtalOptimizer(data, v_sel)) { + indexed_lattice = data.latt; + }*/ + } +} + std::optional RotationIndexer::ProcessImage(int64_t image, const std::vector &spots) { std::unique_lock ul(m); @@ -47,7 +114,7 @@ std::optional RotationIndexer::ProcessImage(int64_t image, const const auto rot = axis_->GetTransformation(image); - if (!indexed_lattice && image >= last_accumulated_image + image_stride) { + if (image >= last_accumulated_image + image_stride) { v_.reserve(v_.size() + spots.size()); coords_.reserve(coords_.size() + spots.size()); @@ -61,44 +128,18 @@ std::optional RotationIndexer::ProcessImage(int64_t image, const accumulated_images++; last_accumulated_image = image; - if (accumulated_images >= min_images_for_indexing && image >= first_image_to_try_indexing) { - // Index - auto indexer_result = indexer_.Run(experiment, coords_).get(); - if (!indexer_result.lattice.empty()) { - // Find lattice type - auto sym_result = LatticeSearch(indexer_result.lattice[0]); - // Run refinement - - DiffractionExperiment experiment_copy(experiment); - XtalOptimizerData data{ - .geom = experiment_copy.GetDiffractionGeometry(), - .latt = sym_result.conventional, - .crystal_system = sym_result.system, - .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), - .refine_beam_center = true, - .refine_distance_mm = false, - .axis = axis_ - }; - - if (data.crystal_system == gemmi::CrystalSystem::Trigonal) - data.crystal_system = gemmi::CrystalSystem::Hexagonal; - - if (XtalOptimizer(data, v_)) { - indexed_lattice = data.latt; - std::cout << sym_result.niggli_class << std::endl; - auto uc = indexed_lattice->GetUnitCell(); - std::cout << uc.a << " " << uc.b << " " << uc.c << " " << uc.alpha << " " << uc.beta << " " << uc.gamma - << std::endl; - std::cout << indexed_lattice->Vec0().x << " " << indexed_lattice->Vec0().y << " " << indexed_lattice->Vec0().z << std::endl; - std::cout << indexed_lattice->Vec1().x << " " << indexed_lattice->Vec1().y << " " << indexed_lattice->Vec1().z << std::endl; - std::cout << indexed_lattice->Vec2().x << " " << indexed_lattice->Vec2().y << " " << indexed_lattice->Vec2().z << std::endl; - - } - } - } + if (!indexed_lattice && accumulated_images >= min_images_for_indexing && image >= first_image_to_try_indexing) + TryIndex(); } if (indexed_lattice) return indexed_lattice->Multiply(rot); return {}; } + + +std::optional RotationIndexer::Finalize(bool retry) { + if (!indexed_lattice || retry) + TryIndex(); + return indexed_lattice; +} diff --git a/image_analysis/RotationIndexer.h b/image_analysis/RotationIndexer.h index 5c6671be..9f8a8b62 100644 --- a/image_analysis/RotationIndexer.h +++ b/image_analysis/RotationIndexer.h @@ -11,7 +11,6 @@ #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 @@ -22,6 +21,7 @@ class RotationIndexer { mutable std::mutex m; const DiffractionExperiment& experiment; + constexpr static int64_t max_spots = 32768; 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; @@ -42,11 +42,12 @@ class RotationIndexer { 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); }; diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 815b7793..eeb04752 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -36,7 +36,7 @@ TARGET_LINK_LIBRARIES(jfjoch_lite_perf_test JFJochReceiver JFJochWriter) INSTALL(TARGETS jfjoch_lite_perf_test RUNTIME COMPONENT jfjoch) ADD_EXECUTABLE(jfjoch_indexing_test jfjoch_indexing_test.cpp) -TARGET_LINK_LIBRARIES(jfjoch_indexing_test JFJochReader JFJochImageAnalysis) +TARGET_LINK_LIBRARIES(jfjoch_indexing_test JFJochReader JFJochImageAnalysis JFJochWriter) INSTALL(TARGETS jfjoch_indexing_test RUNTIME COMPONENT jfjoch) ADD_EXECUTABLE(jfjoch_extract_hkl jfjoch_extract_hkl.cpp) diff --git a/tools/jfjoch_indexing_test.cpp b/tools/jfjoch_indexing_test.cpp index 00860282..27351559 100644 --- a/tools/jfjoch_indexing_test.cpp +++ b/tools/jfjoch_indexing_test.cpp @@ -5,6 +5,7 @@ #include "../image_analysis/indexing/IndexerFactory.h" #include "../common/Logger.h" #include "../image_analysis/RotationIndexer.h" +#include "../writer/FileWriter.h" int main(int argc, char** argv) { Logger logger("jfjoch_indexing_test"); @@ -19,7 +20,13 @@ int main(int argc, char** argv) { throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Error opening file"); DiffractionExperiment experiment(reader.GetDataset()->experiment); - experiment.IndexingAlgorithm(IndexingAlgorithmEnum::Auto); + experiment.SetUnitCell({}); + experiment.IndexingAlgorithm(IndexingAlgorithmEnum::FFT); + + experiment.FilePrefix("output"); + StartMessage msg; + experiment.FillMessage(msg); + // FileWriter fw(msg); IndexerThreadPool indexer(experiment.GetIndexingSettings()); @@ -34,17 +41,25 @@ int main(int argc, char** argv) { DataMessage msg = img->ImageData(); - auto output = rot_index.ProcessImage(i, msg.spots); - - logger.Info("Result {} {}", msg.number, output.has_value()); - if (output.has_value()) { - auto a = output->Vec0(); - auto b = output->Vec1(); - auto c = output->Vec2(); - - logger.Info("Lattice {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f}", - a[0],a[1],a[2],b[0],b[1],b[2],c[0],c[1], c[2]); - } + rot_index.ProcessImage(i, msg.spots); + if (i % 200 == 0) + logger.Info("Image {}", i); } + auto latt_1 = rot_index.Finalize(false); + if (latt_1) { + auto uc = latt_1->GetUnitCell(); + logger.Info("Unit cell {} {} {} {} {} {}", uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma); + logger.Info("Latt1 {} {} {}", latt_1->Vec0().x, latt_1->Vec0().y, latt_1->Vec0().z); + logger.Info("Latt1 {} {} {}", latt_1->Vec1().x, latt_1->Vec1().y, latt_1->Vec1().z); + logger.Info("Latt1 {} {} {}", latt_1->Vec2().x, latt_1->Vec2().y, latt_1->Vec2().z); + } + auto latt_2 = rot_index.Finalize(true); + if (latt_2) { + auto uc = latt_2->GetUnitCell(); + logger.Info("Unit cell {} {} {} {} {} {}", uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma); + logger.Info("Latt1 {} {} {}", latt_2->Vec0().x, latt_2->Vec0().y, latt_2->Vec0().z); + logger.Info("Latt1 {} {} {}", latt_2->Vec1().x, latt_2->Vec1().y, latt_2->Vec1().z); + logger.Info("Latt1 {} {} {}", latt_2->Vec2().x, latt_2->Vec2().y, latt_2->Vec2().z); + } } \ No newline at end of file