mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests 1 Actions Packages Projects Releases 55 Wiki Activity

RotationIndexer: Remove logic controlling image accumulation

Browse Source
This commit is contained in:
leonarski_f
2026-05-29 16:27:22 +02:00
parent e3924a08dd
commit 7b20eedcbd
4 changed files with 34 additions and 189 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/IndexAndRefine.cpp
+8 -3
View File
@@ -26,10 +26,10 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool
}
IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotation(DataMessage &msg) {
const auto result = rotation_indexer->GetLattice();
IndexingOutcome outcome(experiment);
rotation_indexer->ProcessImage(msg.number, msg.spots);
auto result = rotation_indexer->GetLattice();
if (result) {
// For rotation indexing, indexing rate is calculated only for frames, where "global" rotation indexing solution was found
msg.indexing_result = false;
@@ -292,8 +292,13 @@ void IndexAndRefine::ProcessImage(DataMessage &msg,
}
std::optional<RotationIndexerResult> IndexAndRefine::Finalize() {
if (rotation_indexer)
if (rotation_indexer) {
if (const auto latt = rotation_indexer->GetLattice())
return latt;
rotation_indexer->RunIndexing();
return rotation_indexer->GetLattice();
}
return {};
}
image_analysis/RotationIndexer.cpp
+13 -57
View File
@@ -16,30 +16,11 @@ RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPo
geom_(x.GetDiffractionGeometry()),
updated_geom_(geom_),
indexer_(indexer) {
if (axis_) {
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 = std::max<int64_t>(0, x.GetImageNum() - 1);
next_image_to_try_indexing = first_image_to_try_indexing;
image_stride = 1;
} else {
first_image_to_try_indexing = std::max<int64_t>(min_images_for_indexing,
x.GetIndexingSettings().GetRotationIndexingMinAngularRange_deg() / angle_norm_deg);
next_image_to_try_indexing = first_image_to_try_indexing;
image_stride = std::ceil(x.GetIndexingSettings().GetRotationIndexingAngularStride_deg() / angle_norm_deg);
if (image_stride == 0)
image_stride = 1;
}
}
}
}
IndexerResult RotationIndexer::RunIndexing() const {
void RotationIndexer::RunIndexing() {
std::unique_lock ul(m);
std::vector<Coord> coords;
coords.reserve(max_spots_per_image * v_.size());
@@ -50,13 +31,9 @@ IndexerResult RotationIndexer::RunIndexing() const {
for (const auto &s: v_[i])
coords.emplace_back(rot * s.ReciprocalCoord(geom_));
}
const auto indexer_result = indexer_.Run(experiment, coords);
return indexer_.Run(experiment, coords);
}
void RotationIndexer::TryIndex() {
if (auto indexer_result = RunIndexing();
!indexer_result.lattice.empty() && indexer_result.lattice[0].CalcVolume() > 1.0) {
if (!indexer_result.lattice.empty() && indexer_result.lattice[0].CalcVolume() > 1.0) {
auto sg = experiment.GetGemmiSpaceGroup();
if (sg) {
search_result_ = LatticeSearchResult{
@@ -97,15 +74,6 @@ void RotationIndexer::TryIndex() {
axis_ = data.axis;
}
}
if (!indexed_lattice) {
if (indexing_range_multiplier < max_indexing_range_multiplier) {
indexing_range_multiplier++;
next_image_to_try_indexing = first_image_to_try_indexing * indexing_range_multiplier;
} else {
next_image_to_try_indexing = INT64_MAX;
}
}
}
void RotationIndexer::ProcessImage(int64_t image, const std::vector<SpotToSave> &spots) {
@@ -115,15 +83,12 @@ void RotationIndexer::ProcessImage(int64_t image, const std::vector<SpotToSave>
if (!axis_)
return;
if (accumulated_images >= max_images_for_indexing)
if (accumulated_spots >= max_spots)
return;
if (indexed_lattice)
return;
if (image < last_accumulated_image + image_stride)
return;
v_[image].reserve(spots.size());
for (const auto &s: spots) {
@@ -131,26 +96,20 @@ void RotationIndexer::ProcessImage(int64_t image, const std::vector<SpotToSave>
v_[image].emplace_back(s);
}
if (v_[image].size() > max_spots_per_image) {
std::ranges::nth_element(v_[image], v_[image].begin() + max_spots_per_image,
// truncate spots, so we don't get above max_spots (total) and max_spots_per_image (for this image)
size_t max_spots_limit = std::min(max_spots_per_image, max_spots - accumulated_spots);
if (v_[image].size() > max_spots_limit) {
std::ranges::nth_element(v_[image], v_[image].begin() + max_spots_limit,
[](const SpotToSave &a, const SpotToSave &b) {
return a.intensity > b.intensity;
}
);
v_[image].resize(max_spots_per_image);
v_[image].resize(max_spots_limit);
}
accumulated_images++;
last_accumulated_image = image;
const bool short_scan_last_image =
(experiment.GetImageNum() < min_images_for_indexing) &&
(image >= experiment.GetImageNum() - 1);
if ((accumulated_images >= min_images_for_indexing || short_scan_last_image) &&
image >= next_image_to_try_indexing)
TryIndex();
accumulated_spots += v_[image].size();
}
std::optional<RotationIndexerResult> RotationIndexer::GetLattice() const {
@@ -166,6 +125,3 @@ std::optional<RotationIndexerResult> RotationIndexer::GetLattice() const {
};
}
int64_t RotationIndexer::GetAccumulatedImages() const {
return accumulated_images;
}
image_analysis/RotationIndexer.h
+6 -26
View File
@@ -1,8 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#ifndef JFJOCH_ROTATIONINDEXER_H
#define JFJOCH_ROTATIONINDEXER_H
#pragma once
#include <vector>
#include <mutex>
@@ -12,12 +11,6 @@
#include "indexing/IndexerThreadPool.h"
#include "lattice_search/LatticeSearch.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
struct RotationIndexerResult {
CrystalLattice lattice;
LatticeSearchResult search_result;
@@ -29,13 +22,10 @@ class RotationIndexer {
mutable std::mutex m;
const DiffractionExperiment& experiment;
constexpr static int64_t max_spots_per_image = 200;
constexpr static int64_t max_images_for_indexing = 100;
constexpr static int64_t min_images_for_indexing = 10;
constexpr static int64_t max_indexing_range_multiplier = 3;
constexpr static size_t max_spots = 10000;
constexpr static size_t max_spots_per_image = 200;
const bool index_ice_rings;
float angle_norm_deg = 1.0f;
std::vector<std::vector<SpotToSave>> v_;
@@ -46,24 +36,14 @@ class RotationIndexer {
IndexerThreadPool &indexer_;
int64_t last_accumulated_image = -1;
int64_t accumulated_images = 0;
int64_t image_stride = 1;
int64_t first_image_to_try_indexing = INT64_MAX;
int64_t next_image_to_try_indexing = INT64_MAX;
int64_t indexing_range_multiplier = 1;
size_t accumulated_spots = 0;
std::optional<CrystalLattice> indexed_lattice;
IndexerResult RunIndexing() const;
void TryIndex();
public:
RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer);
void ProcessImage(int64_t image, const std::vector<SpotToSave>& spots);
void RunIndexing();
std::optional<RotationIndexerResult> GetLattice() const;
int64_t GetAccumulatedImages() const;
size_t GetAccumulatedImages() const;
};
#endif //JFJOCH_ROTATIONINDEXER_H
tests/RotationIndexerTest.cpp
+7 -103
View File
@@ -70,7 +70,11 @@ TEST_CASE("RotationIndexer") {
s.indexed = true;
spots.push_back(s);
}
indexer.ProcessImage(img, spots);
if (img == 30)
indexer.RunIndexing();
auto result = indexer.GetLattice();
if (result.has_value())
cnt++;
@@ -82,112 +86,12 @@ TEST_CASE("RotationIndexer") {
REQUIRE(ret.has_value());
auto uc = ret->lattice.GetUnitCell();
auto uc_ref = latt_base.GetUnitCell();
REQUIRE(std::fabs(uc.a - uc_ref.a) < 0.02 );
REQUIRE(std::fabs(uc.b - uc_ref.b) < 0.02 );
REQUIRE(std::fabs(uc.c - uc_ref.c) < 0.02 );
REQUIRE(std::fabs(uc.a - uc_ref.a) < 0.1);
REQUIRE(std::fabs(uc.b - uc_ref.b) < 0.1);
REQUIRE(std::fabs(uc.c - uc_ref.c) < 0.1);
REQUIRE(std::fabs(uc.alpha - uc_ref.alpha) < 0.1);
REQUIRE(std::fabs(uc.beta - uc_ref.beta) < 0.1);
REQUIRE(std::fabs(uc.gamma - uc_ref.gamma) < 0.1);
CHECK(ret->search_result.centering == 'P');
CHECK(ret->search_result.system == gemmi::CrystalSystem::Orthorhombic);
}
TEST_CASE("RotationIndexer short scan indexes only at the end") {
DiffractionExperiment exp_i;
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
.BeamX_pxl(1000)
.BeamY_pxl(1000)
.PoniRot1_rad(0.01)
.PoniRot2_rad(0.02)
.DetectorDistance_mm(200)
.ImagesPerTrigger(6);
IndexingSettings settings;
#ifdef JFJOCH_USE_CUDA
settings.Algorithm(IndexingAlgorithmEnum::FFT);
#elif JFJOCH_USE_FFTW
settings.Algorithm(IndexingAlgorithmEnum::FFTW);
#else
return;
#endif
settings.RotationIndexing(true).RotationIndexingAngularStride_deg(1.0).RotationIndexingMinAngularRange_deg(30.0);
exp_i.ImportIndexingSettings(settings);
CrystalLattice latt_base(40, 50, 80, 90, 90, 90);
latt_base = latt_base.Multiply(RotMatrix(2.0, Coord(sqrt(3)/3,sqrt(3)/3,sqrt(3)/3)));
GoniometerAxis axis("omega", 0.0f, 30.0f, Coord(1,0,0), std::nullopt);
exp_i.Goniometer(axis);
BraggPredictionSettings prediction_settings{
.high_res_A = 1.3,
.ewald_dist_cutoff = 0.002
};
IndexerThreadPool indexer_thread_pool(exp_i.GetIndexingSettings());
RotationIndexer indexer(exp_i, indexer_thread_pool);
BraggPrediction prediction;
for (int img = 0; img < 5; ++img) {
std::vector<SpotToSave> spots;
const float angle_deg = axis.GetAngle_deg(img) + axis.GetWedge_deg() / 2.0f;
const RotMatrix rot = axis.GetTransformationAngle(angle_deg);
const CrystalLattice latt_img = latt_base.Multiply(rot.transpose());
const auto n = prediction.Calc(exp_i, latt_img, prediction_settings);
for (int i = 0; i < n; ++i) {
const auto& r = prediction.GetReflections().at(i);
SpotToSave s{};
s.x = r.predicted_x;
s.y = r.predicted_y;
s.image = img;
s.intensity = 1.0f;
s.phi = angle_deg;
s.ice_ring = false;
s.indexed = true;
spots.push_back(s);
}
indexer.ProcessImage(img, spots);
auto result = indexer.GetLattice();
CHECK_FALSE(result.has_value());
}
std::vector<SpotToSave> last_spots;
const int last_img = 5;
const float angle_deg = axis.GetAngle_deg(last_img) + axis.GetWedge_deg() / 2.0f;
const RotMatrix rot = axis.GetTransformationAngle(angle_deg);
const CrystalLattice latt_img = latt_base.Multiply(rot.transpose());
const auto n = prediction.Calc(exp_i, latt_img, prediction_settings);
for (int i = 0; i < n; ++i) {
const auto& r = prediction.GetReflections().at(i);
SpotToSave s{};
s.x = r.predicted_x;
s.y = r.predicted_y;
s.image = last_img;
s.intensity = 1.0f;
s.phi = angle_deg;
s.ice_ring = false;
s.indexed = true;
last_spots.push_back(s);
}
indexer.ProcessImage(last_img, last_spots);
auto ret_last = indexer.GetLattice();
REQUIRE(ret_last.has_value());
auto ret = indexer.GetLattice();
REQUIRE(ret.has_value());
auto uc = ret->lattice.GetUnitCell();
auto uc_ref = latt_base.GetUnitCell();
REQUIRE(std::fabs(uc.a - uc_ref.a) < 0.02 );
REQUIRE(std::fabs(uc.b - uc_ref.b) < 0.02 );
REQUIRE(std::fabs(uc.c - uc_ref.c) < 0.02 );
REQUIRE(std::fabs(uc.alpha - uc_ref.alpha) < 0.1);
REQUIRE(std::fabs(uc.beta - uc_ref.beta) < 0.1);
REQUIRE(std::fabs(uc.gamma - uc_ref.gamma) < 0.1);
CHECK(ret->search_result.centering == 'P');
CHECK(ret->search_result.system == gemmi::CrystalSystem::Orthorhombic);
}