RotationIndexer: Work in progress on rotation data indexer

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)
 

image_analysis/RotationIndexer.cpp

@@ -0,0 +1,87 @@
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// 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<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;
+}
+
+std::optional<CrystalLattice> 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 (!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 {};
+}

image_analysis/RotationIndexer.h

@@ -0,0 +1,59 @@
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#ifndef JFJOCH_SPOTROTATIONSTORE_H
+#define JFJOCH_SPOTROTATIONSTORE_H
+
+#include <vector>
+#include <mutex>
+
+#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<size_t> image_spot0;
+    std::vector<size_t> image_nspots;
+
+    std::vector<SpotToSave> v_;
+    std::vector<Coord> coords_;
+    std::optional<GoniometerAxis> 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<CrystalLattice> indexed_lattice;
+
+public:
+    RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer);
+    void SetLattice(const CrystalLattice &lattice);
+    std::optional<CrystalLattice> ProcessImage(int64_t image, const std::vector<SpotToSave>& spots);
+};
+
+
+#endif //JFJOCH_SPOTROTATIONSTORE_H