From 118fed1d8392145e880920ffaf2b2e81af7ede4a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 29 Nov 2025 21:02:14 +0100 Subject: [PATCH] IndexAndRefine: Refactor integration --- image_analysis/CMakeLists.txt | 6 +- image_analysis/IndexAndRefine.cpp | 277 ++++++++++++++++++ .../{RotationIndexer.h => IndexAndRefine.h} | 34 +-- image_analysis/MXAnalysisAfterFPGA.cpp | 16 +- image_analysis/MXAnalysisAfterFPGA.h | 10 +- image_analysis/MXAnalysisWithoutFPGA.cpp | 13 +- image_analysis/MXAnalysisWithoutFPGA.h | 6 +- image_analysis/RotationIndexer.cpp | 163 ----------- image_analysis/SpotAnalyze.cpp | 128 -------- image_analysis/SpotAnalyze.h | 23 -- receiver/JFJochReceiver.cpp | 16 +- receiver/JFJochReceiver.h | 3 +- receiver/JFJochReceiverFPGA.cpp | 3 +- receiver/JFJochReceiverLite.cpp | 2 +- tools/jfjoch_indexing_test.cpp | 31 +- 15 files changed, 344 insertions(+), 387 deletions(-) create mode 100644 image_analysis/IndexAndRefine.cpp rename image_analysis/{RotationIndexer.h => IndexAndRefine.h} (78%) delete mode 100644 image_analysis/RotationIndexer.cpp delete mode 100644 image_analysis/SpotAnalyze.cpp delete mode 100644 image_analysis/SpotAnalyze.h diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 5680f254..6cdda11c 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -3,10 +3,8 @@ ADD_LIBRARY(JFJochImageAnalysis STATIC MXAnalysisWithoutFPGA.h MXAnalysisAfterFPGA.h MXAnalysisAfterFPGA.cpp - SpotAnalyze.cpp - SpotAnalyze.h - RotationIndexer.cpp - RotationIndexer.h + IndexAndRefine.cpp + IndexAndRefine.h dark_mask_analysis/DarkMaskAnalysis.cpp dark_mask_analysis/DarkMaskAnalysis.h) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp new file mode 100644 index 00000000..4bfc522b --- /dev/null +++ b/image_analysis/IndexAndRefine.cpp @@ -0,0 +1,277 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "IndexAndRefine.h" + +#include "bragg_integration/BraggIntegrate2D.h" +#include "bragg_integration/CalcISigma.h" +#include "geom_refinement/XtalOptimizer.h" +#include "indexing/AnalyzeIndexing.h" +#include "indexing/FFTIndexer.h" +#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()), + 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; + } + } + } +} + +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; + } + } +} + +void IndexAndRefine::ProcessImage(DataMessage &msg, + const SpotFindingSettings &spot_finding_settings, + const CompressedImage &image, + BraggPrediction &prediction) { + if (!indexer_) + return; + + msg.indexing_result = false; + std::optional lattice_candidate; + DiffractionExperiment experiment_copy(experiment); + + LatticeMessage symmetry{ + .centering = 'P', + .niggli_class = 0, + .crystal_system = gemmi::CrystalSystem::Triclinic + }; + + bool beam_center_updated = false; + + if (!axis_) { + // 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_)); + } + + auto indexer_result = indexer_->Run(experiment, recip).get(); + msg.indexing_time_s = indexer_result.indexing_time_s; + + if (!indexer_result.lattice.empty()) { + auto latt = indexer_result.lattice[0]; + + auto sg = experiment.GetGemmiSpaceGroup(); + + // If space group provided => always enforce symmetry in refinement + // If space group not provided => guess symmetry + if (sg) { + // If space group provided but no unit cell fixed, it is better to keep triclinic for now + if (experiment.GetUnitCell()) { + symmetry = LatticeMessage{ + .centering = sg->centring_type(), + .niggli_class = 0, + .crystal_system = sg->crystal_system() + }; + } + } else { + auto sym_result = LatticeSearch(latt); + symmetry = LatticeMessage{ + .centering = sym_result.centering, + .niggli_class = sym_result.niggli_class, + .crystal_system = sym_result.system + }; + + latt = sym_result.conventional; + } + + lattice_candidate = latt; + } + } else { + std::unique_lock ul(m); + const auto rot = axis_->GetTransformation(msg.number); + + 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_)); + } + } + + accumulated_images++; + last_accumulated_image = msg.number; + + 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; + } + } + + 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 (XtalOptimizer(data, msg.spots)) { + if (AnalyzeIndexing(msg, experiment, data.latt)) { + msg.lattice_type = symmetry; + + float ewald_dist_cutoff = 0.001f; + + if (msg.profile_radius) + ewald_dist_cutoff = msg.profile_radius.value() * 2.0f; + if (experiment.GetBraggIntegrationSettings().GetFixedProfileRadius_recipA()) + 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); + + constexpr size_t kMaxReflections = 10000; + if (res.size() > kMaxReflections) { + msg.reflections.assign(res.begin(), res.begin() + kMaxReflections); + } else + msg.reflections = res; + + CalcISigma(msg); + CalcWilsonBFactor(msg); + } + } + } + } +} + +std::optional IndexAndRefine::Finalize(bool retry) { + if (axis_) { + std::unique_lock ul(m); + + if (!indexed_lattice || retry) + TryIndexRotationData(); + return RotationIndexerResult{ + .lattice = indexed_lattice.value(), + .search_result = search_result_, + .geom = geom_ + }; + } + return {}; +} diff --git a/image_analysis/RotationIndexer.h b/image_analysis/IndexAndRefine.h similarity index 78% rename from image_analysis/RotationIndexer.h rename to image_analysis/IndexAndRefine.h index 44339ec6..cc1be5c6 100644 --- a/image_analysis/RotationIndexer.h +++ b/image_analysis/IndexAndRefine.h @@ -9,6 +9,7 @@ #include "../common/DiffractionSpot.h" #include "../common/DiffractionExperiment.h" +#include "bragg_integration/BraggPrediction.h" #include "indexing/IndexerThreadPool.h" #include "lattice_search/LatticeSearch.h" @@ -24,42 +25,39 @@ struct RotationIndexerResult { DiffractionGeometry geom; }; -class RotationIndexer { - mutable std::mutex m; - const DiffractionExperiment& experiment; - +class IndexAndRefine { 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 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_; - 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; + std::optional indexed_lattice; + + std::optional axis_; + + IndexerThreadPool *indexer_; int64_t image_stride; int64_t first_image_to_try_indexing; - std::optional indexed_lattice; - - - void TryIndex(); + void TryIndexRotationData(); public: - RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer); + IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer); void SetLattice(const CrystalLattice &lattice); - std::optional ProcessImage(int64_t image, const std::vector& spots); + void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, const CompressedImage &image, BraggPrediction &prediction); std::optional Finalize(bool retry); }; -#endif //JFJOCH_SPOTROTATIONSTORE_H \ No newline at end of file +#endif //JFJOCH_SPOTROTATIONSTORE_H diff --git a/image_analysis/MXAnalysisAfterFPGA.cpp b/image_analysis/MXAnalysisAfterFPGA.cpp index f6599166..a09a00d0 100644 --- a/image_analysis/MXAnalysisAfterFPGA.cpp +++ b/image_analysis/MXAnalysisAfterFPGA.cpp @@ -4,7 +4,7 @@ #include "MXAnalysisAfterFPGA.h" #include "spot_finding/DetModuleSpotFinder_cpu.h" #include "../common/CUDAWrapper.h" -#include "SpotAnalyze.h" +#include "spot_finding/SpotUtils.h" #include "bragg_integration/BraggPredictionFactory.h" double stddev(const std::vector &v) { @@ -25,19 +25,14 @@ double stddev(const std::vector &v) { } -MXAnalysisAfterFPGA::MXAnalysisAfterFPGA(const DiffractionExperiment &in_experiment) -: experiment(in_experiment) { +MXAnalysisAfterFPGA::MXAnalysisAfterFPGA(const DiffractionExperiment &in_experiment, IndexAndRefine &indexer) +: experiment(in_experiment), indexer(indexer) { if (experiment.IsSpotFindingEnabled()) find_spots = true; prediction = CreateBraggPrediction(); } -MXAnalysisAfterFPGA &MXAnalysisAfterFPGA::SetIndexer(IndexerThreadPool *input) { - indexer = input; - return *this; -} - void MXAnalysisAfterFPGA::ReadFromFPGA(const DeviceOutput *output, const SpotFindingSettings &settings, size_t module_number) { if (!find_spots || !settings.enable) return; @@ -93,7 +88,10 @@ void MXAnalysisAfterFPGA::Process(DataMessage &message, const SpotFindingSetting if (!find_spots) return; - SpotAnalyze(experiment, spot_finding_settings, spots, message.image, *prediction, indexer, message); + SpotAnalyze(experiment, spot_finding_settings, spots, message); + + if (spot_finding_settings.indexing) + indexer.ProcessImage(message, spot_finding_settings, message.image, *prediction); spots.clear(); } diff --git a/image_analysis/MXAnalysisAfterFPGA.h b/image_analysis/MXAnalysisAfterFPGA.h index b4df6935..33b521df 100644 --- a/image_analysis/MXAnalysisAfterFPGA.h +++ b/image_analysis/MXAnalysisAfterFPGA.h @@ -8,11 +8,12 @@ #include "bragg_integration/BraggPrediction.h" #include "indexing/IndexerThreadPool.h" #include "spot_finding/StrongPixelSet.h" +#include "IndexAndRefine.h" class MXAnalysisAfterFPGA { mutable std::mutex read_from_cpu_mutex; const DiffractionExperiment &experiment; - IndexerThreadPool *indexer = nullptr; + IndexAndRefine &indexer; std::unique_ptr prediction; bool find_spots = false; @@ -25,7 +26,7 @@ class MXAnalysisAfterFPGA { std::vector arr_strong_pixel; public: - explicit MXAnalysisAfterFPGA(const DiffractionExperiment& experiment); + MXAnalysisAfterFPGA(const DiffractionExperiment& experiment, IndexAndRefine &indexer); void ReadFromFPGA(const DeviceOutput* output, const SpotFindingSettings& settings, @@ -35,10 +36,7 @@ public: const SpotFindingSettings &settings, size_t module_number); - void Process(DataMessage &message, - const SpotFindingSettings& settings); - - MXAnalysisAfterFPGA& SetIndexer(IndexerThreadPool *input); + void Process(DataMessage &message, const SpotFindingSettings& settings); }; diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 567345eb..5dc2fa55 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -6,14 +6,14 @@ #include "spot_finding/StrongPixelSet.h" #include "../compression/JFJochDecompress.h" -#include "SpotAnalyze.h" +#include "spot_finding/SpotUtils.h" #include "spot_finding/ImageSpotFinderFactory.h" #include "bragg_integration/BraggPredictionFactory.h" MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_experiment, const AzimuthalIntegration &in_integration, const PixelMask &in_mask, - IndexerThreadPool *in_indexer) + IndexAndRefine &in_indexer) : experiment(in_experiment), integration(in_integration), roi_map(experiment.ExportROIMap()), @@ -160,9 +160,12 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, const std::vector spots = spotFinder->Run(settings, mask_resolution); - SpotAnalyze(experiment, settings, spots, - CompressedImage(updated_image, experiment.GetXPixelsNum(), experiment.GetYPixelsNum()), - *prediction, indexer, output); + SpotAnalyze(experiment, settings, spots, output); + + if (settings.indexing) + indexer.ProcessImage(output, settings, + CompressedImage(updated_image, experiment.GetXPixelsNum(), experiment.GetYPixelsNum()), + *prediction); } profile.Add(azim_sum, azim_count); diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index 2389c5e8..205618b8 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -14,6 +14,7 @@ #include "bragg_integration/BraggPrediction.h" #include "spot_finding/ImageSpotFinder.h" #include "indexing/IndexerThreadPool.h" +#include "IndexAndRefine.h" class MXAnalysisWithoutFPGA { std::mutex m; @@ -30,7 +31,7 @@ class MXAnalysisWithoutFPGA { std::vector mask_1bit; std::unique_ptr spotFinder; - IndexerThreadPool *indexer; + IndexAndRefine &indexer; std::unique_ptr prediction; std::vector &updated_image; @@ -49,8 +50,7 @@ class MXAnalysisWithoutFPGA { void Analyze(DataMessage &output, const uint8_t *image, T err_pixel_val, T sat_pixel_val, AzimuthalIntegrationProfile &profile, const SpotFindingSettings &settings); public: MXAnalysisWithoutFPGA(const DiffractionExperiment &experiment, const AzimuthalIntegration &integration, - const PixelMask &mask, - IndexerThreadPool *indexer = nullptr); + const PixelMask &mask, IndexAndRefine &indexer); void Analyze(DataMessage &output, std::vector &buffer, AzimuthalIntegrationProfile &profile, const SpotFindingSettings &spot_finding_settings); }; diff --git a/image_analysis/RotationIndexer.cpp b/image_analysis/RotationIndexer.cpp deleted file mode 100644 index b654ade7..00000000 --- a/image_analysis/RotationIndexer.cpp +++ /dev/null @@ -1,163 +0,0 @@ -// 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/SpotAnalyze.cpp b/image_analysis/SpotAnalyze.cpp deleted file mode 100644 index ea71fb6e..00000000 --- a/image_analysis/SpotAnalyze.cpp +++ /dev/null @@ -1,128 +0,0 @@ -// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#include "SpotAnalyze.h" - -#include "geom_refinement/XtalOptimizer.h" -#include "spot_finding/SpotUtils.h" -#include "spot_finding/StrongPixelSet.h" -#include "bragg_integration/BraggIntegrate2D.h" -#include "indexing/AnalyzeIndexing.h" -#include "bragg_integration/CalcISigma.h" -#include "lattice_search/LatticeSearch.h" - -void SpotAnalyze(const DiffractionExperiment &experiment, - const SpotFindingSettings &spot_finding_settings, - const std::vector &spots, - const CompressedImage &image, - BraggPrediction &prediction, - IndexerThreadPool *indexer, - DataMessage &output) { - auto geom = experiment.GetDiffractionGeometry(); - - SpotAnalyze(experiment, spot_finding_settings, spots, output); - - if ((indexer != nullptr) && spot_finding_settings.indexing) { - std::vector recip; - recip.reserve(output.spots.size()); - bool index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(); - for (const auto &i: output.spots) { - if (index_ice_rings || !i.ice_ring) - recip.push_back(i.ReciprocalCoord(geom)); - } - - auto indexer_result = indexer->Run(experiment, recip).get(); - output.indexing_time_s = indexer_result.indexing_time_s; - - if (indexer_result.lattice.empty()) - output.indexing_result = false; - else { - auto latt = indexer_result.lattice[0]; - - bool beam_center_updated = false; - - auto sg = experiment.GetGemmiSpaceGroup(); - LatticeMessage symmetry{ - .centering = 'P', - .niggli_class = 0, - .crystal_system = gemmi::CrystalSystem::Triclinic - }; - - // If space group provided => always enforce symmetry in refinement - // If space group not provided => guess symmetry - - if (sg) { - // If space group provided but no unit cell fixed, it is better to keep triclinic for now - if (experiment.GetUnitCell()) { - symmetry = LatticeMessage{ - .centering = sg->centring_type(), - .niggli_class = 0, - .crystal_system = sg->crystal_system() - }; - } - } else { - auto sym_result = LatticeSearch(latt); - symmetry = LatticeMessage{ - .centering = sym_result.centering, - .niggli_class = sym_result.niggli_class, - .crystal_system = sym_result.system, - .primitive = sym_result.primitive_reduced - }; - output.lattice_type = symmetry; - latt = sym_result.conventional; - } - - DiffractionExperiment experiment_copy(experiment); - XtalOptimizerData data{ - .geom = experiment_copy.GetDiffractionGeometry(), - .latt = latt, - .crystal_system = symmetry.crystal_system, - .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), - }; - - 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, output.spots)) { - experiment_copy.BeamX_pxl(data.geom.GetBeamX_pxl()).BeamY_pxl(data.geom.GetBeamY_pxl()); - beam_center_updated = true; - latt = data.latt; - } - break; - } - - if (AnalyzeIndexing(output, experiment_copy, latt)) { - float ewald_dist_cutoff = 0.001f; - - if (output.profile_radius) - ewald_dist_cutoff = output.profile_radius.value() * 2.0f; - if (experiment.GetBraggIntegrationSettings().GetFixedProfileRadius_recipA()) - ewald_dist_cutoff = experiment.GetBraggIntegrationSettings().GetFixedProfileRadius_recipA().value() - * 3.0f; - - if (beam_center_updated) { - output.beam_corr_x = data.beam_corr_x; - output.beam_corr_y = data.beam_corr_y; - } - - if (spot_finding_settings.quick_integration) { - auto res = BraggIntegrate2D(experiment_copy, image, latt, - prediction, ewald_dist_cutoff, output.number, symmetry.centering); - - constexpr size_t kMaxReflections = 10000; - if (res.size() > kMaxReflections) { - output.reflections.assign(res.begin(), res.begin() + kMaxReflections); - } else - output.reflections = res; - - CalcISigma(output); - CalcWilsonBFactor(output); - } - } - } - } -} diff --git a/image_analysis/SpotAnalyze.h b/image_analysis/SpotAnalyze.h deleted file mode 100644 index 1c7ecf81..00000000 --- a/image_analysis/SpotAnalyze.h +++ /dev/null @@ -1,23 +0,0 @@ -// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#ifndef JFJOCH_SPOTANALYZE_H -#define JFJOCH_SPOTANALYZE_H - -#include -#include "../common/DiffractionSpot.h" -#include "../common/DiffractionExperiment.h" -#include "../common/JFJochMessages.h" -#include "bragg_integration/BraggPrediction.h" -#include "spot_finding/SpotFindingSettings.h" -#include "indexing/IndexerThreadPool.h" - -void SpotAnalyze(const DiffractionExperiment &experiment, - const SpotFindingSettings &spot_finding_settings, - const std::vector &spots, - const CompressedImage &image, - BraggPrediction &prediction, - IndexerThreadPool *indexer, - DataMessage &output); - -#endif //JFJOCH_SPOTANALYZE_H diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp index 19623a80..1e0e7a43 100644 --- a/receiver/JFJochReceiver.cpp +++ b/receiver/JFJochReceiver.cpp @@ -19,22 +19,22 @@ JFJochReceiver::JFJochReceiver(const DiffractionExperiment &in_experiment, ZMQPreviewSocket *in_zmq_preview_socket, ZMQMetadataSocket *in_zmq_metadata_socket, IndexerThreadPool *indexing_thread_pool) - : experiment(in_experiment), + : logger(logger), + experiment(in_experiment), + spot_finding_settings(spot_finding_settings), image_buffer(in_image_buffer), image_pusher(in_image_pusher), preview_image(in_preview_image), - current_status(in_current_status), - plots(in_plots), - spot_finding_settings(spot_finding_settings), - logger(logger), - numa_policy(in_numa_policy), zmq_preview_socket(in_zmq_preview_socket), zmq_metadata_socket(in_zmq_metadata_socket), + current_status(in_current_status), + plots(in_plots), + scan_result(in_experiment), serialmx_filter(in_experiment), + numa_policy(in_numa_policy), pixel_mask(in_pixel_mask), az_int_mapping(experiment, pixel_mask), - scan_result(in_experiment), - indexer_thread_pool(indexing_thread_pool) { + indexer(experiment, indexing_thread_pool) { image_buffer.StartMeasurement(experiment); current_status.SetProgress(0); diff --git a/receiver/JFJochReceiver.h b/receiver/JFJochReceiver.h index b82ce908..2544f414 100644 --- a/receiver/JFJochReceiver.h +++ b/receiver/JFJochReceiver.h @@ -23,6 +23,7 @@ #include "../common/NUMAHWPolicy.h" #include "../common/ScanResultGenerator.h" #include "../image_analysis/indexing/IndexerThreadPool.h" +#include "../image_analysis/IndexAndRefine.h" class JFJochReceiver { protected: @@ -84,7 +85,7 @@ protected: std::optional max_delay; std::mutex max_delay_mutex; - IndexerThreadPool *indexer_thread_pool; + IndexAndRefine indexer; std::string writer_error; diff --git a/receiver/JFJochReceiverFPGA.cpp b/receiver/JFJochReceiverFPGA.cpp index 589cbda4..97053239 100644 --- a/receiver/JFJochReceiverFPGA.cpp +++ b/receiver/JFJochReceiverFPGA.cpp @@ -298,8 +298,7 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { try { numa_policy.Bind(threadid); - analyzer = std::make_unique(experiment); - analyzer->SetIndexer(indexer_thread_pool); + analyzer = std::make_unique(experiment, indexer); } catch (const JFJochException &e) { frame_transformation_ready.count_down(); logger.Error("Thread setup error {}", e.what()); diff --git a/receiver/JFJochReceiverLite.cpp b/receiver/JFJochReceiverLite.cpp index 0ebe36fd..f2e952b6 100644 --- a/receiver/JFJochReceiverLite.cpp +++ b/receiver/JFJochReceiverLite.cpp @@ -243,7 +243,7 @@ void JFJochReceiverLite::DataAnalysisThread(uint32_t id) { measurement_started.wait(); try { - analysis = std::make_unique(experiment, az_int_mapping, pixel_mask, indexer_thread_pool); + analysis = std::make_unique(experiment, az_int_mapping, pixel_mask, indexer); } catch (const JFJochException &e) { Cancel(e); return; diff --git a/tools/jfjoch_indexing_test.cpp b/tools/jfjoch_indexing_test.cpp index b8bb3909..53ee0000 100644 --- a/tools/jfjoch_indexing_test.cpp +++ b/tools/jfjoch_indexing_test.cpp @@ -4,7 +4,8 @@ #include "../reader/JFJochHDF5Reader.h" #include "../image_analysis/indexing/IndexerFactory.h" #include "../common/Logger.h" -#include "../image_analysis/RotationIndexer.h" +#include "../image_analysis/IndexAndRefine.h" +#include "../image_analysis/bragg_integration/BraggPredictionFactory.h" #include "../writer/FileWriter.h" #include "../image_analysis/indexing/AnalyzeIndexing.h" @@ -30,12 +31,16 @@ int main(int argc, char** argv) { // FileWriter fw(msg); IndexerThreadPool indexer(experiment.GetIndexingSettings()); + SpotFindingSettings settings{}; + settings.quick_integration = false; + settings.indexing = true; auto geom = experiment.GetDiffractionGeometry(); - RotationIndexer rot_index(experiment, indexer); - int cnt = 0; - int cnt_2 = 0; + IndexAndRefine rot_index(experiment, &indexer); + auto prediction = CreateBraggPrediction(); + + int cnt = 0; for (int i = 0; i < reader.GetNumberOfImages(); i++) { auto img = reader.LoadImage(i); @@ -44,19 +49,13 @@ int main(int argc, char** argv) { DataMessage msg = img->ImageData(); - auto img_res = rot_index.ProcessImage(i, msg.spots); - if (img_res ) { - auto img_latt = img_res->lattice; - DiffractionExperiment x(experiment); - x.DetectorDistance_mm(img_res->geom.GetDetectorDistance_mm()) - .BeamX_pxl(img_res->geom.GetBeamX_pxl()) - .BeamY_pxl(img_res->geom.GetBeamY_pxl()); - cnt++; - if (AnalyzeIndexing(msg, x, img_latt)) - cnt_2++; - } + rot_index.ProcessImage(msg, settings, msg.image, *prediction); + + if (msg.indexing_result.value_or(false)) + cnt++; + if (i % 200 == 0) - logger.Info("Image {} {} {}", i, cnt, cnt_2); + logger.Info("Image {} {}", i, cnt); } auto latt_1 = rot_index.Finalize(false);