IndexAndRefine: Refactor integration

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)
 

image_analysis/IndexAndRefine.cpp

@@ -0,0 +1,277 @@
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// 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<int64_t>(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<SpotToSave> v_sel;
+    std::vector<Coord> coords_sel;
+
+    if (coords_.size() > max_spots) {
+        // Indices into v_ / coords_
+        std::vector<std::size_t> 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<CrystalLattice> 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<Coord> 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<RotationIndexerResult> 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 {};
+}

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<SpotToSave> v_;
     std::vector<Coord> coords_;
-    std::optional<GoniometerAxis> 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<CrystalLattice> indexed_lattice;
+
+    std::optional<GoniometerAxis> axis_;
+
+    IndexerThreadPool *indexer_;
 
     int64_t image_stride;
     int64_t first_image_to_try_indexing;
 
-    std::optional<CrystalLattice> 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<RotationIndexerResult> ProcessImage(int64_t image, const std::vector<SpotToSave>& spots);
+    void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, const CompressedImage &image, BraggPrediction &prediction);
     std::optional<RotationIndexerResult> Finalize(bool retry);
 };
 
 
-#endif //JFJOCH_SPOTROTATIONSTORE_H
+#endif //JFJOCH_SPOTROTATIONSTORE_H

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<float> &v) {
@@ -25,19 +25,14 @@ double stddev(const std::vector<float> &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();
 }

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<BraggPrediction> prediction;
 
     bool find_spots = false;
@@ -25,7 +26,7 @@ class MXAnalysisAfterFPGA {
     std::vector<uint32_t> 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);
 };
 
 

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<DiffractionSpot> 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);

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<bool> mask_1bit;
 
     std::unique_ptr<ImageSpotFinder> spotFinder;
-    IndexerThreadPool *indexer;
+    IndexAndRefine &indexer;
     std::unique_ptr<BraggPrediction> prediction;
     std::vector<int32_t> &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<uint8_t> &buffer, AzimuthalIntegrationProfile &profile, const SpotFindingSettings &spot_finding_settings);
 };
 

image_analysis/RotationIndexer.cpp

@@ -1,163 +0,0 @@
-// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
-// 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<int64_t>(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<SpotToSave> v_sel;
-    std::vector<Coord> coords_sel;
-
-    if (coords_.size() > max_spots) {
-
-        // Indices into v_ / coords_
-        std::vector<std::size_t> 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<RotationIndexerResult> RotationIndexer::ProcessImage(int64_t image, const std::vector<SpotToSave> &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<RotationIndexerResult> RotationIndexer::Finalize(bool retry) {
-    if (!indexed_lattice || retry)
-        TryIndex();
-    return RotationIndexerResult{
-        .lattice = indexed_lattice.value(),
-        .search_result = search_result_,
-        .geom = geom_
-    };
-}

image_analysis/SpotAnalyze.cpp

@@ -1,128 +0,0 @@
-// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
-// 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<DiffractionSpot> &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<Coord> 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);
-                }
-            }
-        }
-    }
-}

image_analysis/SpotAnalyze.h

@@ -1,23 +0,0 @@
-// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
-// SPDX-License-Identifier: GPL-3.0-only
-
-#ifndef JFJOCH_SPOTANALYZE_H
-#define JFJOCH_SPOTANALYZE_H
-
-#include <vector>
-#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<DiffractionSpot> &spots,
-                 const CompressedImage &image,
-                 BraggPrediction &prediction,
-                 IndexerThreadPool *indexer,
-                 DataMessage &output);
-
-#endif //JFJOCH_SPOTANALYZE_H

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);

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<uint64_t> max_delay;
     std::mutex max_delay_mutex;
 
-    IndexerThreadPool *indexer_thread_pool;
+    IndexAndRefine indexer;
 
     std::string writer_error;
 

receiver/JFJochReceiverFPGA.cpp

@@ -298,8 +298,7 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) {
 
     try {
         numa_policy.Bind(threadid);
-        analyzer = std::make_unique<MXAnalysisAfterFPGA>(experiment);
-        analyzer->SetIndexer(indexer_thread_pool);
+        analyzer = std::make_unique<MXAnalysisAfterFPGA>(experiment, indexer);
     } catch (const JFJochException &e) {
         frame_transformation_ready.count_down();
         logger.Error("Thread setup error {}", e.what());

receiver/JFJochReceiverLite.cpp

@@ -243,7 +243,7 @@ void JFJochReceiverLite::DataAnalysisThread(uint32_t id) {
     measurement_started.wait();
 
     try {
-        analysis = std::make_unique<MXAnalysisWithoutFPGA>(experiment, az_int_mapping, pixel_mask, indexer_thread_pool);
+        analysis = std::make_unique<MXAnalysisWithoutFPGA>(experiment, az_int_mapping, pixel_mask, indexer);
     } catch (const JFJochException &e) {
         Cancel(e);
         return;

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);