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

RotationIndexer: Work in progress

Browse Source
This commit is contained in:
leonarski_f
2025-11-29 17:28:33 +01:00
parent 163518f724
commit 549bb9fb5d
4 changed files with 109 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/RotationIndexer.cpp
+78 -37
View File
@@ -1,8 +1,10 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "RotationIndexer.h"
#include <iostream>
#include "geom_refinement/XtalOptimizer.h"
#include "indexing/FFTIndexer.h"
#include "lattice_search/LatticeSearch.h"
RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPool &indexer)
@@ -38,6 +40,71 @@ void RotationIndexer::SetLattice(const CrystalLattice &lattice) {
indexed_lattice = lattice;
}
#include <iostream>
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
auto sym_result = LatticeSearch(indexer_result.lattice[0]);
indexed_lattice = indexer_result.lattice[0];
/*
// Run refinement
DiffractionExperiment experiment_copy(experiment);
XtalOptimizerData data{
.geom = experiment_copy.GetDiffractionGeometry(),
.latt = sym_result.conventional,
.crystal_system = sym_result.system,
.min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(),
.refine_beam_center = true,
.refine_distance_mm = false,
.axis = axis_
};
if (data.crystal_system == gemmi::CrystalSystem::Trigonal)
data.crystal_system = gemmi::CrystalSystem::Hexagonal;
if (XtalOptimizer(data, v_sel)) {
indexed_lattice = data.latt;
}*/
}
}
std::optional<CrystalLattice> RotationIndexer::ProcessImage(int64_t image, const std::vector<SpotToSave> &spots) {
std::unique_lock ul(m);
@@ -47,7 +114,7 @@ std::optional<CrystalLattice> RotationIndexer::ProcessImage(int64_t image, const
const auto rot = axis_->GetTransformation(image);
if (!indexed_lattice && image >= last_accumulated_image + image_stride) {
if (image >= last_accumulated_image + image_stride) {
v_.reserve(v_.size() + spots.size());
coords_.reserve(coords_.size() + spots.size());
@@ -61,44 +128,18 @@ std::optional<CrystalLattice> RotationIndexer::ProcessImage(int64_t image, const
accumulated_images++;
last_accumulated_image = image;
if (accumulated_images >= min_images_for_indexing && image >= first_image_to_try_indexing) {
// Index
auto indexer_result = indexer_.Run(experiment, coords_).get();
if (!indexer_result.lattice.empty()) {
// Find lattice type
auto sym_result = LatticeSearch(indexer_result.lattice[0]);
// Run refinement
DiffractionExperiment experiment_copy(experiment);
XtalOptimizerData data{
.geom = experiment_copy.GetDiffractionGeometry(),
.latt = sym_result.conventional,
.crystal_system = sym_result.system,
.min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(),
.refine_beam_center = true,
.refine_distance_mm = false,
.axis = axis_
};
if (data.crystal_system == gemmi::CrystalSystem::Trigonal)
data.crystal_system = gemmi::CrystalSystem::Hexagonal;
if (XtalOptimizer(data, v_)) {
indexed_lattice = data.latt;
std::cout << sym_result.niggli_class << std::endl;
auto uc = indexed_lattice->GetUnitCell();
std::cout << uc.a << " " << uc.b << " " << uc.c << " " << uc.alpha << " " << uc.beta << " " << uc.gamma
<< std::endl;
std::cout << indexed_lattice->Vec0().x << " " << indexed_lattice->Vec0().y << " " << indexed_lattice->Vec0().z << std::endl;
std::cout << indexed_lattice->Vec1().x << " " << indexed_lattice->Vec1().y << " " << indexed_lattice->Vec1().z << std::endl;
std::cout << indexed_lattice->Vec2().x << " " << indexed_lattice->Vec2().y << " " << indexed_lattice->Vec2().z << std::endl;
}
}
}
if (!indexed_lattice && accumulated_images >= min_images_for_indexing && image >= first_image_to_try_indexing)
TryIndex();
}
if (indexed_lattice)
return indexed_lattice->Multiply(rot);
return {};
}
std::optional<CrystalLattice> RotationIndexer::Finalize(bool retry) {
if (!indexed_lattice || retry)
TryIndex();
return indexed_lattice;
}
image_analysis/RotationIndexer.h
+3 -2
View File
@@ -11,7 +11,6 @@
#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
@@ -22,6 +21,7 @@ class RotationIndexer {
mutable std::mutex m;
const DiffractionExperiment& experiment;
constexpr static int64_t max_spots = 32768;
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;
@@ -42,11 +42,12 @@ class RotationIndexer {
int64_t first_image_to_try_indexing;
std::optional<CrystalLattice> indexed_lattice;
void TryIndex();
public:
RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer);
void SetLattice(const CrystalLattice &lattice);
std::optional<CrystalLattice> ProcessImage(int64_t image, const std::vector<SpotToSave>& spots);
std::optional<CrystalLattice> Finalize(bool retry);
};
tools/CMakeLists.txt
+1 -1
View File
@@ -36,7 +36,7 @@ TARGET_LINK_LIBRARIES(jfjoch_lite_perf_test JFJochReceiver JFJochWriter)
INSTALL(TARGETS jfjoch_lite_perf_test RUNTIME COMPONENT jfjoch)
ADD_EXECUTABLE(jfjoch_indexing_test jfjoch_indexing_test.cpp)
TARGET_LINK_LIBRARIES(jfjoch_indexing_test JFJochReader JFJochImageAnalysis)
TARGET_LINK_LIBRARIES(jfjoch_indexing_test JFJochReader JFJochImageAnalysis JFJochWriter)
INSTALL(TARGETS jfjoch_indexing_test RUNTIME COMPONENT jfjoch)
ADD_EXECUTABLE(jfjoch_extract_hkl jfjoch_extract_hkl.cpp)
tools/jfjoch_indexing_test.cpp
+27 -12
View File
@@ -5,6 +5,7 @@
#include "../image_analysis/indexing/IndexerFactory.h"
#include "../common/Logger.h"
#include "../image_analysis/RotationIndexer.h"
#include "../writer/FileWriter.h"
int main(int argc, char** argv) {
Logger logger("jfjoch_indexing_test");
@@ -19,7 +20,13 @@ int main(int argc, char** argv) {
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Error opening file");
DiffractionExperiment experiment(reader.GetDataset()->experiment);
experiment.IndexingAlgorithm(IndexingAlgorithmEnum::Auto);
experiment.SetUnitCell({});
experiment.IndexingAlgorithm(IndexingAlgorithmEnum::FFT);
experiment.FilePrefix("output");
StartMessage msg;
experiment.FillMessage(msg);
// FileWriter fw(msg);
IndexerThreadPool indexer(experiment.GetIndexingSettings());
@@ -34,17 +41,25 @@ int main(int argc, char** argv) {
DataMessage msg = img->ImageData();
auto output = rot_index.ProcessImage(i, msg.spots);
logger.Info("Result {} {}", msg.number, output.has_value());
if (output.has_value()) {
auto a = output->Vec0();
auto b = output->Vec1();
auto c = output->Vec2();
logger.Info("Lattice {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f}",
a[0],a[1],a[2],b[0],b[1],b[2],c[0],c[1], c[2]);
}
rot_index.ProcessImage(i, msg.spots);
if (i % 200 == 0)
logger.Info("Image {}", i);
}
auto latt_1 = rot_index.Finalize(false);
if (latt_1) {
auto uc = latt_1->GetUnitCell();
logger.Info("Unit cell {} {} {} {} {} {}", uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
logger.Info("Latt1 {} {} {}", latt_1->Vec0().x, latt_1->Vec0().y, latt_1->Vec0().z);
logger.Info("Latt1 {} {} {}", latt_1->Vec1().x, latt_1->Vec1().y, latt_1->Vec1().z);
logger.Info("Latt1 {} {} {}", latt_1->Vec2().x, latt_1->Vec2().y, latt_1->Vec2().z);
}
auto latt_2 = rot_index.Finalize(true);
if (latt_2) {
auto uc = latt_2->GetUnitCell();
logger.Info("Unit cell {} {} {} {} {} {}", uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
logger.Info("Latt1 {} {} {}", latt_2->Vec0().x, latt_2->Vec0().y, latt_2->Vec0().z);
logger.Info("Latt1 {} {} {}", latt_2->Vec1().x, latt_2->Vec1().y, latt_2->Vec1().z);
logger.Info("Latt1 {} {} {}", latt_2->Vec2().x, latt_2->Vec2().y, latt_2->Vec2().z);
}
}