IndexAndRefine: Stabilize profile radius and mosaicity, by using moving average

image_analysis/CMakeLists.txt

@@ -28,7 +28,9 @@ ADD_LIBRARY(JFJochImageAnalysis STATIC
         dark_mask_analysis/DarkMaskAnalysis.cpp
         dark_mask_analysis/DarkMaskAnalysis.h
         RotationIndexer.cpp
-        RotationIndexer.h)
+        RotationIndexer.h
+        RotationParameters.cpp
+        RotationParameters.h)
 
 FIND_PACKAGE(Eigen3 3.4 REQUIRED NO_MODULE)  # provides Eigen3::Eigen
 

image_analysis/IndexAndRefine.cpp

@@ -20,8 +20,6 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool
         rotation_indexer = std::make_unique<RotationIndexer>(x, *indexer);
 }
 
-#include <iostream>
-
 IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(DataMessage &msg) {
     IndexingOutcome outcome(experiment);
 
@@ -131,13 +129,32 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
     if (!outcome.lattice_candidate)
         return;
 
-    float ewald_dist_cutoff = 0.001f;
+    CrystalLattice latt = outcome.lattice_candidate.value();
 
+    if (rotation_indexer) {
+        // Rotate lattice_candidate to the current image angle
+        auto gon = experiment.GetGoniometer();
+        if (gon)
+            latt = outcome.lattice_candidate->Multiply(gon->GetTransformation(-msg.number));
+        // Use moving average for mosaicity and profile_radius (also add beam center later)
+        if (msg.mosaicity_deg)
+            msg.mosaicity_deg = rotation_parameters.Mosaicity(msg.mosaicity_deg.value());
+        if (msg.profile_radius) {
+            msg.profile_radius = rotation_parameters.ProfileRadius(msg.profile_radius.value());
+        }
+    }
+
+    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;
 
+
+
+
+
     BraggPredictionSettings settings_prediction{
         .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(),
         .ewald_dist_cutoff = ewald_dist_cutoff,
@@ -146,15 +163,6 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
         .image_number = static_cast<int>(msg.number)
     };
 
-    CrystalLattice latt = outcome.lattice_candidate.value();
-
-    if (rotation_indexer) {
-        // Rotate lattice_candidate to the current image angle
-        auto gon = experiment.GetGoniometer();
-        if (gon)
-            latt = outcome.lattice_candidate->Multiply(gon->GetTransformation(-msg.number));
-    }
-
     auto nrefl = prediction.Calc(outcome.experiment, latt, settings_prediction);
     auto refl_ret = BraggIntegrate2D(outcome.experiment, image, prediction.GetReflections(), nrefl, msg.number);
 

image_analysis/IndexAndRefine.h

@@ -13,6 +13,7 @@
 #include "indexing/IndexerThreadPool.h"
 #include "lattice_search/LatticeSearch.h"
 #include "RotationIndexer.h"
+#include "RotationParameters.h"
 
 class IndexAndRefine {
     const bool index_ice_rings;
@@ -26,6 +27,8 @@ class IndexAndRefine {
     IndexerThreadPool *indexer_;
     std::unique_ptr<RotationIndexer> rotation_indexer;
 
+    RotationParameters rotation_parameters;
+
     struct IndexingOutcome {
         std::optional<CrystalLattice> lattice_candidate;
         DiffractionExperiment experiment;

image_analysis/RotationParameters.cpp

@@ -0,0 +1,27 @@
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "RotationParameters.h"
+
+RotationParameters::RotationParameters()
+    : profile_radius(100), beam_x(100), beam_y(100), mosaicity(100) {}
+
+float RotationParameters::ProfileRadius(float input) {
+    profile_radius.Add(input);
+    return profile_radius.Read().value_or(input);
+}
+
+float RotationParameters::Mosaicity(float input) {
+    mosaicity.Add(input);
+    return mosaicity.Read().value_or(input);
+}
+
+float RotationParameters::BeamX(float input) {
+    beam_x.Add(input);
+    return beam_x.Read().value_or(input);
+}
+
+float RotationParameters::BeamY(float input) {
+    beam_y.Add(input);
+    return beam_y.Read().value_or(input);
+}

image_analysis/RotationParameters.h

@@ -0,0 +1,20 @@
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#ifndef JFJOCH_ROTATIONPARAMETERS_H
+#define JFJOCH_ROTATIONPARAMETERS_H
+
+#include "../common/MovingAverage.h"
+
+class RotationParameters {
+    MovingAverage profile_radius, beam_x, beam_y, mosaicity;
+public:
+    RotationParameters();
+    float ProfileRadius(float input);
+    float Mosaicity(float input);
+    float BeamX(float input);
+    float BeamY(float input);
+};
+
+
+#endif //JFJOCH_ROTATIONPARAMETERS_H